Dynamics of the interaction of a pair of thin evaporating droplets on compliant substrates

Abstract The dynamics of the interaction of a system of two thin volatile liquid droplets resting on a soft viscoelastic solid substrate are investigated theoretically. The developed model fully considers the effect of evaporative cooling and the generated Marangoni stresses due to the induced thermal gradients, while also accounting for the effect of the gas phase composition and the diffusion of vapour in the atmosphere of the droplets. Using the framework of lubrication theory, we derive evolution equations for both the droplet profile and the displacement of the elastic solid, which are solved in combination with Laplace's equation for the vapour concentration in the gas phase. A disjoining-pressure/precursor-film approach is used to describe contact-line motion. The evolution equations are solved numerically, using the finite-element method, and we present a thorough parametric analysis to investigate the physical properties and mechanisms that affect the dynamics of droplet interactions. The results show that the droplets interact through both the soft substrate and the gas phase. In the absence of thermocapillary phenomena, the combined effect of non-uniform evaporation due to the increased vapour concentration between the two droplets and elastocapillary phenomena determines whether the drop–drop interaction is attractive or repulsive. The Marangoni stresses suppress droplet attraction at the early stages of the drying process and lead to longer droplet lifetimes. For substrates with intermediate stiffness, the emergence of spontaneous symmetry breaking at late stages of evaporation is found. The rich dynamics of this complex system is explored by constructing a detailed map of the dynamic regimes.


Introduction
The dynamic wetting behavior of liquid droplets on soft solid substrates has received a great deal of attention lately, due to its relevance in diverse applications (Bico et al. 2018) ranging from biology, i.e. the inhibition of the dispersal of cancer cells (Douezan et al. 2012) or the control of medicine dispersal on tissues, to the control of the spreading of the deposited particles over a compliant substrate after the evaporation of ink-jetted microdroplets (Park & Moon 2006) and microfabrication of materials in technology (Kong et al. 2014;Bonaccurso et al. 2005;Pericet-Camara et al. 2007).
The evaporation of droplets on rigid substrates has been widely studied over the years, underlining various aspects of evaporation such as droplet lifetimes (Stauber et al. 2014(Stauber et al. , 2015)), the impact of capillary flow on the coffee stain effect (Deegan et al. 1997), or the effect of substrate properties (Erbil 2012).A key concept in droplet evaporation is the so-called shielding effect, where neighboring droplets interact with each other, since the presence of the vapour from adjacent droplets reduces the evaporation rate and increases the droplet lifetime in comparison to those of a single isolated droplet (Fabrikant 1985;Wray et al. 2020Wray et al. , 2021;;Masoud et al. 2021).On the contrary, the study of droplet evaporation on compliant solid substrates is insubstantial so far.
The unbalanced vertical forces acting on the contact line of the liquid result in a local deformation of the compliant substrate affecting droplet shape and, ultimately, the dynamics of the flow (Andreotti & Snoeijer 2020).This unique attribute of the compliant substrates is responsible for the creation of a macroscopic surface protrusion, most widely referred to as a wetting ridge (Shanahan 1988;Park et al. 2014;Chen et al. 2020).The structure of the wetting ridge is an immediate outcome of the balance between capillary and elastic forces, while a key role of the solid surface tension has been recently identified (Jerison et al. 2011;Style & Dufresne 2012;Marchand et al. 2012;Style et al. 2013).Besides slowing down the contact line, the wetting ridge can also constitute the reason for periodic stick-slip behavior, with a periodic depinning of the contact line (Kajiya et al. 2012;Lopes & Bonaccurso 2013;Yu et al. 2013;Kajiya et al. 2014;Karpitschka et al. 2016;van Gorcum et al. 2018;Mokbel et al. 2022).Starting from a rigid substrate, when we increase the substrate softness we observe an initial strong increase in the wetting ridge size while the droplet footprint is kept relatively constant, which is replaced by a strong depression of the substrate under the droplet in much softer substrates (Charitatos & Kumar 2020;Henkel et al. 2021).The arising solid angle is governed by a balance of surface tensions (Style & Dufresne 2012;Marchand et al. 2012;Style et al. 2013).The possible strain dependence of the solid surface tension (i.e., the Shuttleworth effect) (Andreotti & Snoeijer 2016;Xu et al. 2017;Pandey et al. 2020;Henkel et al. 2022) gives rise to another complexity.
Depending on a combination of capillarity and bulk elasticity, adjacent droplets on solid substrates can interact, leading either to droplet attraction and even coalescence on thick substrates, or to droplet repulsion on thinner substrates (Hernández-Sánchez et al. 2012;Karpitschka et al. 2016;Leong & Le 2020).Karpitschka et al. (2016) described the droplet interaction when deposited on soft solids as the "inverse Cheerios effect" with direct reference to the liquid-on-solid analog of the so-called "Cheerios effect" (Vella & Mahadevan 2005); the latter refers to solid particles attraction when floating on liquids, mediated by surface tension forces, and it has been named after the sticking of breakfast cereals either to the walls of a bowl or to each other.However, there are substantial differences between the "Cheerios effect" and the "inverted Cheerios effect" regarding the driving force and the mechanism which mediates the interaction.The shape of the liquid interface in the "Cheerios effect" is specified by the balance between surface tension and gravity, while the interaction is driven by a change in gravitational potential energy.On the contrary, in the "inverse Cheerios effect" there is no gravity involved and the shape of the solid interface is specified by elastocapillarity (Karpitschka et al. 2015).
Despite the innumerable studies concerning droplet coalescence (Eggers et al. 1999;Aarts et al. 2005), not many of them consider the role that a compliant substrate might play in the process.More recently, Henkel et al. (2021) investigated the two basic coarsening modes of two droplets on soft substrates without evaporation; the volume or mass transfer mode (also referred to as drop collapse mode, diffusion-controlled ripening or Ostwald ripening) and the translation mode (also referred to as coalescence, collision or migration mode).On the one hand, on the mass transfer mode, material is transferred from the small drop to the larger one until the smaller drop has completely vanished, while the centers of mass of each droplet remain constant.This mass transfer might occur either through the vapour phase in case of volatile droplets, or through an adsorption layer in case of non-evaporating partially wetting fluids.On the other hand, on the translation mode, there is droplet migration towards each other, until their contact lines touch, leading to coalescence (Henkel et al. 2021).
When compared to rigid substrates, early experimental studies (Lopes & Bonaccurso 2012;Pu & Severtson 2012;Lopes & Bonaccurso 2013;Yu et al. 2013;Chuang et al. 2014;Gerber et al. 2019) highlighted the faster evaporation of droplets observed on softer substrates, due to the longer pinning of the contact line throughout evaporation, caused by the formation of the wetting ridge.In addition, some of these experimental studies (Lopes & Bonaccurso 2013;Yu et al. 2013) revealed that while on the initial stages of evaporation, the droplet appears to remain pinned, while the contact angle is decreased.After the droplet depins, the opposite behavior is observed, i.e. the contact angle remains constant and the contact radius decreases.Then, at the late stages of droplet evaporation, the contact angle slightly increases, before ultimately decrease until the droplet completely evaporates.
As it has been established in studies for rigid substrates, the dynamics of the contact line plays a crucial role on droplet evaporation (Lopes & Bonaccurso 2012;Pu & Severtson 2012;Lopes & Bonaccurso 2013;Yu et al. 2013;Chuang et al. 2014;Gerber et al. 2019) and in order to model the contact line motion an approach, followed by several researchers, has been to introduce the effect of disjoining pressure, by assuming the presence of an adsorbed precursor film ahead of the droplet, which is stabilized by the action of intermolecular interactions between the two interfaces.The presence of the precursor film is evident in experimental studies with microscopic techniques (Kavehpour et al. 2003;Xu et al. 2004;Hoang & Kavehpour 2011) and constitutes the reason for the levelled transition from the contact line to the flat gas-solid interface, circumventing the singularity arising from the shear stress (Wang et al. 2021), due to the contradiction between the non-zero displacement on the contact line and the no-slip boundary condition on the liquid-solid interface on the same point.This approach has been widely implemented for the modelling of not only perfectly wetting fluids (Bonn et al. 2009), but also partially wetting fluids (Schwartz 1998;Schwartz & Eley 1998;Gomba & Homsy 2010).A similar approach has been also used by Charitatos & Kumar (2021) to examine droplet evaporation on partially wetted soft viscoelastic substrates.
To account for the droplet evaporation, two approximations have been mostly used (Wilson & D'Ambrosio 2023).In the so-called one-sided model, the attention is solely drawn to the liquid phase, since vapour viscosity, density and thermal conductivity are considered negligible.In this approach, evaporation is limited by the rate on which the molecules of the solvent are headed from the liquid towards the gas phase (Burelbach et al. 1988).On that principle, the work of Moosman & Homsy (1980), Ajaev & Homsy (2001) and later Ajaev (2005) considered an adsorbed thin film ahead of the evaporating liquid, with non-zero film thickness, which is in thermodynamic equilibrium with both the solid and the gas phases.The work of Ajaev has been the grounding principle for many researchers studying qualitatively droplet spreading and evaporation of more complex systems, such as droplets with nanoparticles (Matar et al. 2007), the deposition of particles while in the presence of surfactants (Karapetsas et al. 2016), as well as the evaporation of droplets which consist of ethanol-water or other binary mixtures (Williams et al. 2021), or the vapour adsorption of hygroscopic aqueous solution droplets (Wang et al. 2021).Concerning droplet evaporation on compliant substrates, the only theoretical work so far is the work of Charitatos & Kumar (2021), where a one-sided model is developed to study the evaporation of a single droplet.
However, when the evaporation is diffusion-limited and, hence, the vapour phase can not be considered irrelevant, quantitative results can only be achieved by employing a two-sided approximation (Sultan et al. 2005;Schofield et al. 2020;Wray et al. 2020).Typically, one-sided models consider that the evaporation flux is only a function of pressure and temperature differences in the liquid-gas interface and the transport equations in the liquid phase are the only prerequisite for the system modelling.On the contrary, the diffusion-limited model includes the simultaneous solution of a diffusion equation concerning the vapor gas phase concentration.The second approach entails a higher computational cost but provides a more accurate description of phase change phenomena (Deegan et al. 1997;Hu & Larson 2005;Masoud & Felske 2009;Cazabat & Guéna 2010;Mikishev & Nepomnyashchy 2013;Larson 2014).A common assumption to reduce the computational cost, is to consider that droplet retains a spherical-cap shape.This assumption, though, is not always safe, since the droplet shape might be significantly distorted by forces, such as gravitational forces (e.g.evaporation on inclined substrates), the effect of Marangoni or elastic stresses, etc. Hartmann et al. (2023) recently developed a long-wave model of a sessile shallow droplet of evaporating partially wetting fluid on a rigid substrate, which, similarly to the earlier work of Sultan et al. (2005), captured the transition between the diffusion-limited and the one-sided model.
Several authors have also investigated the effect of Marangoni stresses on droplet spreading and evaporation.Marangoni flows can be induced by thermal gradients, a variation in the concentration of surfactants, or by the presence of binary mixtures (Dunn et al. 2009a;Williams et al. 2021;Wang et al. 2021).The coffee-ring effect may be suppressed by the action of Marangoni stresses (Karapetsas et al. 2016;Seo et al. 2017), since they counter the outward liquid flow facilitated by contact-line pinning.Additionally, the coalescence of merging droplets with different surface tensions (such as water and ethanol) has been shown to be strongly delayed (Karpitschka & Riegler 2010;Chen et al. 2021).Concerning evaporating droplets, Talbot et al. (2012), Schofield et al. (2018) and Dunn et al. (2009b) showed that the thermal effects can significantly extend droplet lifetime.
There have been different approaches in the literature to describe soft substrate wetting, ranging from the development of long-wave models (Kumar & Matar 2004;Matar et al. 2005;Gielok et al. 2017;Gomez & Velay-Lizancos 2020) to full-scale computational studies (Bueno et al. 2017(Bueno et al. , 2018)).Kumar & Matar (2004) and Matar et al. (2005) first developed a long-wave approach to study the non-linear evolution of thin liquid films dewetting near soft elastomeric layers.Charitatos & Kumar (2020) followed Matar's work for droplet spreading on soft solid substrates, while Charitatos & Kumar (2021), extended their own work developing a one-sided model to examine droplet evaporation on viscoelastic substrates.
Building on the latter model, the present paper presents a detailed and comprehensive theoretical model for the investigation of the dynamics of a system of two evaporating droplets residing on a compliant solid substrate.The droplets may interact through the developed elastic stresses in the viscoelastic substrate which is modelled using the Kelvin-Voigt model.When the droplets are exposed in atmosphere, a further question is raised, concerning the effect of the local variations in vapour concentration on the evaporation rate of the droplets.To this end, we develop a two-sided evaporation model following a similar approach with earlier studies for rigid substrates (Sultan et al. 2005;Doumenc & Guerrier 2011).Thus, apart from the viscoelastic substrate, our model unravels the potential of communication of the droplets through the atmosphere, while also taking fully into consideration the effects of evaporative cooling and induced thermocapillary phenomena.To remove the stress singularity that arises at the moving contact line, we assume the presence of sufficiently thin precursor film.The precursor film is stabilised and evaporation therein is suppressed through the inclusion of a disjoining pressure.
The rest of the paper is organised as follows.The problem is formulated in section 2 and the equations governing the flow dynamics are discussed.The scaling and resulting evolution equations are presented in section 3 and section 4, respectively.Results are presented and discussed in section 5, followed by concluding remarks in section 6.

Description of the problem
We consider the behaviour of a single or a system of two two-dimensional sessile evaporating droplets, with initial cross-sectional area V , placed on an incompressible linear viscoelastic solid substrate, which is also referred to as soft substrate.At t = 0, the droplet is resting on the soft substrate and has an initial footprint half-width l0 and an initial height ĥ0 (Fig. 1a).The liquid-solid interface is originally located at ẑ = 0.The soft substrate is originally undistorted and attached to a rigid substrate at ẑ = − Ĥ; the rigid substrate is highly conductive with constant temperature Tb .The soft substrate is characterized by density ρs , viscosity ηs , thermal conductivity λw , shear modulus Ê and constant liquid-solid interfacial tension γ, which is independent of strain; the presence of an immiscible liquid solvent in the soft substrate is assumed.The droplet is assumed to have constant density ρl , viscosity ηl , thermal conductivity λ, specific heat capacity ĉp and saturation temperature Tsat .The liquid-gas interfacial tension, σ, is assumed to be a linear function of temperature where σ0 is the surface tension at the reference temperature Tref , and Ts is the local temperature at the liquid-gas interface.The reference temperature was considered to be equal to the bulk temperature of the gas Tb .At t > 0, the droplet, being in an unsaturated environment, evaporates causing the deformation of the soft solid substrate.The liquid-solid interface is located at ẑ = ξ(x, t), with an outward normal unit vector of ns = − while the liquidair interface is located at ẑ = ζ(x, t), with an outward normal unit vector of nl = for the liquid-solid and the liquid-air interface respectively.
We assume that the droplet is released into a thin precursor film; evaporation in the film is stabilised by the disjoining pressure which accounts for intermolecular van der Waals interactions.The inclusion of the precursor film removes the stress singularity that can arise at the moving contact line (see Fig. 1c).This approach also allows us to easily account for the evaporation of multiple droplets as well as their interactions; in Fig. 1b we depict a system of two evaporating droplets, of the same initial radius and height.Ahead of the contact line the dimensional precursor film thickness is denoted with β (see Fig. 1c) and the apparent contact angle is θ.Concerning the wetting ridge, its maximum height is denoted with ξmax .Moreover, the length of each droplet, that is the distance between the two contact lines, is noted as ∆x cl , while the length of the The soft solid deforms while the system of two droplets spreads and evaporates.(c) Magnified view of the contact line, where β is the precursor film thickness, θ is the apparent contact angle and ξmax denotes the maximum height of the wetting ridge.The local thickness of each droplet is given by ĥ(x, t) = ζ(x, t) − ξ(x, t).computational domain is noted as L x (i.e.0 < x < Lx ).In this domain, the global center of mass of the system is located at x = xcm,g and the center of mass of each droplet is located at x = xcm,l and x = xcm,r , respectively.Consequently, the distance between the two centers of mass is denoted as ∆x cm = xcm,r − xcm,l .The presented model geometry in Fig. 1 constitutes the reference layout for the rest of this paper.
In the present work, the droplets are assumed to be so thin that the droplet aspect ratio ϵ = ĥ0 / R0 is considered to be much smaller than unity; R0 is a characteristic length scale, defined as R0 = 2 l0 /3.Under this assumption, we will employ the lubrication approximation to derive a reduced set of governing evolution equations.Furthermore, gravitational forces are neglected, since the solid and liquid Bond numbers Bo s = ρs ĝ R2 0 /σ and Bo l = ρl ĝ R2 0 /σ are assumed to be less than unity; this condition typically holds for small droplets.A 2-D Cartesian coordinate system (x, ẑ) is used to model the velocity field, which is described by a function of v = (v x , vz ), whereas the solid displacement is given by û = (û x , ûz ).Our model can describe a typical system of water droplets drying on polydimethylsiloxane (PDMS) substrates and the physical properties of such a system are given in Table 1.

Liquid phase
The mass, momentum and energy conservation equations for the liquid are given by: Table 1: Properties of water and PDMS at 20 • C and 1 atm.
where pl is the liquid pressure, ∇ = (∂ x, ∂ ŷ ) is the gradient operator, T is the temperature and αl = λ ρl ĉp is the thermal diffusivity of the liquid.Along the free interface ẑ = ζ(x, t), the liquid velocity v = (v x , vz ) differs from the velocity of the interface vs = (v xs , vzs ).If the evaporative flux is denoted by Ĵ, then Ĵ = ρl (v − vs ) • nl .
(2.5) Furthermore, along the free interface, the local mass, energy and force balances are given by: Ĵ where ρg , λg , vg and Tg denote the density, the thermal conductivity, the velocity and the temperature of the gas phase respectively.Lv is the specific latent internal heat of vaporization, I is the identity tensor, κl is the mean curvature of the free interface, while ∇s is the surface gradient operator.In the above equations, κl = ∇s • nl and ∇s = (I − nl nl ) • ∇.Finally, Π stands for the disjoining pressure, which, taking into consideration the van der Waals interaction, equals to where Â1 = ÂHam / Â3 2 , is a constant that describes the intermolecular interactions between the liquid-gas and the liquid-solid interfaces, ÂHam the Hamaker constant, Â2 is a constant of the same order of magnitude as the precursor film thickness β. ĥ denotes the droplet thickness and n > c > 1.Moreover, the kinematic boundary condition along the moving interface ẑ = ζ(x, t), is described as: (2.10) 2.3.Gas phase The gas phase is assumed to comprise air and vapour, but it is not saturated by vapour.Typically, we may consider that Λ g = λ λg ≪ 1, and under this assumption the bulk temperature of the gas can be assumed to be constant and equal to Tb .Moreover, we assume that the viscosity of the gas, ηg , is much smaller than the viscosity of the liquid phase, i.e. ηg /η l ≪ 1 and therefore the gas can be considered as inviscid and passive with respect to the fluid.
The droplet evaporation is approached using the generalised diffusion-limited model developed by Sultan et al. (2005), in which evaporation is considered limited by the solvent vapour diffusion in the air; this model is able to capture the transition between the diffusion-limited and the one-sided model.The vapour concentration ρv in the gas phase is described by the Laplace's equation, due to the fact that the gas phase is considered to be at rest and the characteristic evaporation time is much larger than the respective diffusion time.As a result, we get ∇2 ρv = 0. (2.11) The vapour mass flux Ĵ is assumed to be limited by the rate of diffusion and thus where Dv the diffusion vapour coefficient.Considering also that the vapour mass flux Ĵ is proportional to the departure from equilibrium at the liquid-gas interface (Schrage 1953;Plesset & Prosperetti 1976;Moosman & Homsy 1980), the following linear constitutive equation, most commonly known as Hertz-Knudsen equation, for Ĵ can be used where α is the accommodation coefficient, usually considered equal to unity near equilibrium.R denotes the universal gas constant, Ts denotes the temperature of the liquid-gas interface, ρv is the local vapour concentration in the gas phase and ρve is the equilibrium vapour concentration.
In order to get a boundary condition for the vapour concentration ρv at ẑ = ζ, we can combine Eqs.(2.12) and (2.13), which leads to (2.14) Finally, following a similar procedure as described by (Moosman & Homsy 1980), the following equation for the equilibrium vapour concentration can be derived where ρv ref is the equilibrium vapour concentration at the reference temperature.
At the far-field boundary, the most natural choice would be to impose a constant vapour concentration at infinity.However, since it is known that for the diffusion-limited model there is no analytical solution in two-dimensional half-space (Schofield et al. 2020), we follow a similar approach to Schofield et al. (2020) considering a finite domain of the gas phase and the far-field condition is replaced by a similar Dirichlet condition at a distant, but finite, boundary.Thus, far from the droplet (ẑ = dg ), the vapour concentration is maintained at a constant initial vapour concentration ρvi : ρv | ẑ= dg = ρvi . (2.16)

Soft solid substrate
The mass, momentum and energy conservation equations for the soft solid substrate are given by: ∇ where ps is the pressure in the solid, αw is the thermal diffusivity of the solid, Tw is the temperature of the solid surface and Ts is the solid stress tensor.Following the work of Kumar & Matar (2004), Matar et al. (2005) and Charitatos & Kumar (2020, 2021), who modeled the soft elastomer layer as a linear viscoelastic material, we consider that the viscoelastic solid is described by the Kelvin-Voigt model and therefore the solid stress tensor is defined as (2.20) Finally, Eq. (2.18) gives ρs (2.21) At ẑ = − Ĥ, the application of the no-slip and no-displacement boundary condition yields: vx = vz = 0 and ûx = ûz = 0, while the temperature at the bottom of the solid substrate is considered to be equal to Tb , i.e.Tw | ẑ=− Ĥ = Tb .
Along the liquid-solid interface, at ẑ = ξ(x, t), we consider thermal equilibrium Tw | ξ = T | ξ and continuity of thermal flux: (2.22) In addition, the combination of the no-slip and the no-penetration boundary conditions at ẑ = ξ(x, t), that is the liquid-solid interface, form the continuity-of-velocity boundary condition: (2.23) The normal and tangential force balances on the liquid-solid interface lead to: where κs = ∇s • ns and Tl stands for the liquid stress tensor, defined as Tl = −p l I + ηl [ ∇v + ( ∇v) T ].

Scaling
In order to render the aforementioned equations and boundary conditions nondimensional, we use the scalings shown below: where ∆ T = ϵ 2 Tref .As Tref we consider the constant bulk temperature of the gas phase, Tb .As characteristic velocity we use Û = ϵ 3 σ0 /η l .Note that henceforth all the variables in the following equations are dimensionless.

Liquid phase
By substituting these scalings and taking into consideration that ϵ ≪ 1, the leading order equations for the liquid are: Along the free interface, i.e. z = ζ(x, t), we get for the mass, energy and force balances in the normal and tangential coordinate: which represents the ratio between the capillary time t c = R0 / Û and the evaporation time t e = ĥ2 0 ρl Lv λ∆ T , as it is derived from the scaling.In Eq.(3.8) the gas pressure has been set equal to zero (datum pressure) without loss of generality.
The kinematic equation, i.e.Eq. (2.10) in combination with Eq. (3.6) gives the evolution equation for the liquid-gas interface z = ζ(x, t): Finally, the scaled disjoining pressure is given by the following expression:

Soft solid substrate
Using the same scaling, the leading equations for the soft solid become: where m = ηs /η l .The ratio of elastic forces to interfacial tension forces is defined as G = ϵ −3 R0 / lec , where lec = σ0 / Ê denotes the elastocapillary length which sets the characteristic size of the deformation of the elastic substrate.
As far as the boundary conditions are concerned, at z = ξ(x, t) we get: The dimensionless continuity of thermal flux along z = ξ(x, t) is given by: where Λ w = λw / λ denotes the ratio of thermal conductivity between the solid and the liquid.
At the liquid-solid interface, i.e. z = ξ(x, t), both the continuity-of-velocity and the normal and tangential force balances render to the following form: where p cap,l = −C l −1 σ∂ 2 ζ/∂x 2 and p cap,s = −C s −1 ∂ 2 ξ/∂x 2 the capillary-like pressures in the liquid and the solid respectively, with C s = ηl Û ϵ 3 γ denoting the capillary number at the liquid-solid interface.

Gas phase
Since the gas phase in the atmosphere may extend to large distances above the liquid phase, the scaling in the z-direction shown in Eq. (3.1) is not appropriate and therefore we employ the same scaling with the x-direction (i.e.ẑ = R0 z ′ and ζ = R0 ζ ′ ).The dimensionless conservation equation for the vapour concentration is then given by: (3.21) The above equation is subjected to the following boundary conditions far from the droplet (z ′ = d g ) and along the liquid-gas interface (z ′ = ζ ′ (x, t)): RTs .In the above equation the equilibrium vapour concentration is given by: (3.25) In order to evaluate the precursor film thickness β, we can combine Eq. (3.24) with Eq. (3.25).Taking into account that far from the droplets the film is flat and at equilibrium with the environment, the following equation is derived An estimation of the order-of-magnitude of certain dimensionless parameters is depicted in Table 2.

Evolution equations
In order to derive the evolution equations, we make an approximation that the streamwise displacement follows a parabolic profile in z, since this is the simplest solution that satisfies Eq. (3.13) (Matar et al. 2005;Ghosh et al. 2016;Charitatos & Kumar 2020, 2021) for the soft solid: in which b 1 , b 2 and b 3 are functions of both space and time that will be determined later; a detailed derivation is given in the Appendix A.
Using the boundary condition at z = −H, i.e. u x = 0, we get for b 3 (x, t): Regarding the coefficient b 1 (x, t), introducing Eq. (4.1) into the x-component of the solid momentum, i.e.Eq. (3.13), yields the following expression: From Eq. (3.19), using Eqs.(4.1) and (3.13), as well as the expression of the xcomponent of the liquid velocity (i.e.Eq. (A 9) derived in the Appendix A), we conclude to an expression for the coefficient b 2 (x, t): In order to derive the evolution equation for ζ(x, t), we use the kinematic equation for the liquid, i.e.Eq. (3.10).Using the expressions of v x and v z (i.e.Eqs.(A 9) and (A 11) respectively derived in the Appendix A), and setting z = ζ, we get: Using the expressions of the material derivatives of ξ and u z (x, 0, t), (i.e.Eqs.(A 12) and (A 13) respectively derived in the Appendix A), leads to an evolution equation for Furthermore, we can easily result in an evolution equation for the droplet thickness h(x, t) = ζ(x, t) − ξ(x, t) by subtracting Eq. (4.6) from Eq. (4.5): where the liquid flowrate, q, is given by: By integrating the energy equation, i.e.Eq. (3.15), with respect to z, and using the continuity of thermal flux at z = ξ(x, t), i.e.Eq. (3.16) and the boundary condition, T w | −H = 0, the following evolution equation for the temperature in the soft solid substrate can be derived: Similarly, by integrating the energy equation, Eq. (3.5), with respect to z, and using the energy balance, Eq. (3.7), and the fact that T | ξ = T w | ξ , the following evolution equation for the temperature in the liquid phase can be derived: In summary, we solve numerically the evolution equations Eqs.(4.3), (4.4), (4.6), (4.7) and (3.24) on the domain 0 < x < L x and Eq.(3.21) on the 2D domain 0 < x < L x , 0 < z < L z .The latter equation is subjected to boundary conditions Eq. (3.22) and (3.23) in the z-direction and to the following condition in the x-direction: The numerical solution of the evolution equations Eqs.(4.3), (4.4), (4.6) and (4.7) is subjected to the following conditions at x = 0 and x = L x : where β = β/ ĥ0 is the dimensionless precursor film height.These conditions were concluded into, after the assumptions that both ζ and ξ are horizontal at x = 0 and x = L and that the dimensionless precursor film thickness equals to the distance between the two interfaces at these positions.Furthermore, the liquid flow rate and the solid displacement in the z-direction were considered equal to zero at all ends.
Concerning the initial conditions, we assumed a flat liquid-solid interface at t = 0: As far as the initial shape of the droplet thickness is concerned, we use a fourth order polynomial which satisfies ∂h ∂x = ∂ 3 h ∂x 3 = 0 at the droplet center (x = x cm,l or x cm,r ) and ∂h ∂x = 0 as well as h = β at distance l 0 from the droplet center, respectively.
The length and height of the computational domain was taken equal to L x = Lx / R0 = 16 and L z = Lz / R0 = 5, respectively.The dimensionless initial droplet footprint halfwidth, l 0 , was defined equal to l 0 = 1.5 and the dimensionless initial droplet crosssectional area was considered to be equal to V = 2/3.Moreover, the initial center of mass for each droplet was taken at x = x cm,l = 4 and at x = x cm,r = 12, hence the initial distance between the two centers of mass was given by ∆x cm = x cm,r − x cm,l = 8, while the center of mass of the system was initially located at x = x cm,g = 8.
The above set of equations is solved using the Finite Element Method and it has been implemented in COMSOL Multiphysics commercial software.We applied a fully implicit finite difference scheme (BDF) to solve the system of the evolution equations and we selected the PARDISO iterative solver for the intermediate time-stepping.Typically we use 10000 elements for the discretization of the system geometry and the moving mesh of the surrounding atmosphere was appropriately refined using free triangular cells; numerical checks showed that increasing the number of elements further led to negligible changes.The simulations stop when the system mass has decreased by 80%.

Results and discussion
Droplet evaporation on compliant substrates is a parametrically rich problem.We begin our study by examining the case of the evaporation of a single droplet on a soft substrate in section 5.1, while in section 5.2 we proceed with simulations for a system of two interacting droplets.Numerical solutions were obtained over a wide range of parameter values.The 'base' case, however, has broadly typical values of ϵ = 0.1, l 0 = 1.5, 5, unless noted otherwise in the text.In the figures that follow, we define a scaled time t ′ = t/t ev where t ev is defined as the time that the system mass has decreased by 80%.

Effect of thermocapillarity
To set the stage, we begin with the simplest configuration, i.e. the evaporation of a single droplet on a soft substrate.Fig. 2 depicts the typical time evolution of the liquidgas and the liquid-solid interfaces for a single sessile evaporating droplet, highlighting the contact line region in the inset of the same figure.Charitatos & Kumar (2021) considered a system similar to the present setup, albeit ignoring the effect of thermocapillarity and employing the one-sided model.In order to examine the effect of thermocapillary phenomena, we present in Fig. 2a the evolution for M a = 0 and in Fig. 2b for M a = 0.005.In the absence of thermocapillary stresses (Fig. 2a), in line with Charitatos & Kumar (2021), we notice a gradual decrease in the droplet footprint, which is accompanied by a small deformation of the soft substrate, due to the balance of the capillary forces along the liquid-solid interface and in the contact line region.Consequently, a wetting ridge is formed and as the droplet dries out, both the contact line and the wetting ridge retract as a result of the decrease in the droplet volume.
On the other hand, in the presence of thermocapillarity (see Fig. 2b), the Marangoni stresses drive liquid towards the colder region (i.e. at droplet apex, see Fig. 3b) causing a faster retraction of the droplet.The faster motion of the contact line results in significantly larger substrate deformation, since for example at t = 100 the maximum deformation of the wetting ridge (evaluated as the z-position of the contact line, see Fig. 1c), is ξ max = 0.038 and ξ max = 0.068, in Figs.2a and 2b, respectively.Furthermore, .This is due to the fact that the action of thermocapillary stresses leads to a considerably smaller droplet footprint with larger distance of the droplet apex from the rigid solid (at z = −H).The increased droplet height inhibits the supply of heat from the substrate (maintained at a constant temperature) to the interface, which is continuously being cooled due to the effect of latent heat.This consequently leads to lower temperature along the liquid-gas interface and in turn results in the overall decrease of the evaporation rate; the evolution of the local evaporation flux is presented in Fig. 3a.respectively.The rest of the system parameters are the same with the 'base' case.
To illustrate the vapour concentration field in the gas phase, we present the corresponding contour plot in Fig. 4, for the case of M a = 0.005.It is noted that the far-field boundary is taken to be very far from the droplet (i.e. at z = ϵ −1 z ′ = 50) and as a result the droplet is difficult to be seen in this figure.We notice, though, that in the neighbourhood of the droplet the vapour concentration is high and decreases moving away from the droplet as expected.

Effect of substrate elasticity and thickness
Here, we examine the effect of elasticity of the substrate by varying G = Ê R0 σ0ϵ 3 ; this parameter measures the ratio of elastic to liquid-gas interfacial tension forces.G is proportional to the shear modulus of the soft solid and therefore smaller values correspond to the case of softer substrates.By letting G → ∞, the case of the rigid substrate can be recovered.In Fig. 5a, we investigate the effect of substrate elasticity on the deformation of the soft solid, by plotting the evolution of the maximum deformation of the wetting ridge, ξ max , with time.Naturally, it can be seen that the softer substrates deform more easily.Figs.5b and 5c depict the time evolution of the contact radius and the apparent contact angle, respectively, of a single droplet evaporating on a rigid (G = 10 7 ) and on soft solid substrates with G = 1, 3, 10, 100.Following the work of Charitatos & Kumar (2021), the apparent contact angle is defined as the largest angle between the tangent of the liquid-air interface z = ζ(x, t) and z = 0. On the other hand, the contact radius is defined as the intersection point between the tangent of the largest angle and z = 0.
At early times and for all examined cases, the droplet contact radius quickly decreases (Fig. 5b), accompanied by an increase in the apparent contact angle (Fig. 5c), while in parallel the size of the wetting ridge grows (Fig. 5a).The initial droplet retraction, which is due to both droplet evaporation and the action of thermocapillary stresses, takes place faster for softer substrates.After the initial droplet retraction, the contact line remains apparently pinned for a significant amount of time (t ≈ 3−300) with a relatively constant droplet footprint, indicating the stick-phase of the droplet spreading (see inset of Fig. 5d).In fact, in the case of softer substrates (i.e.G = 1) the constant contact radius is maintained throughout evaporation accompanied with a continuous decrease of the contact angle (see Fig. 5c); the evaporation takes place in constant contact radius (CCR) mode.In contrast, for harder substrates (i.e.G ⩾ 10) the evaporation takes place in constant contact angle (CCA) mode, with a continuous slow decrease of the contact radius, in line with previous computational studies referring to droplet evaporation on rigid substrates (Pham & Kumar 2017).This CCR mode observed in softer substrates has been previously reported in experimental studies concerning the evaporation of water droplets on compliant PDMS substrates (Lopes & Bonaccurso 2012;Gerber et al. 2019).After de-pinning, the contact line retracts continuously until the droplet fully evaporates.
At the same time a non-monotonous behaviour of the contact angle is observed in line with previous experimental studies of Lopes & Bonaccurso (2012, 2013) and Yu et al. (2013).In Fig. 6a, we investigate the effect of substrate thickness on the deformation of the soft solid, by plotting the maximum deformation of the wetting ridge, ξ max , with time.It can be seen that the thicker substrates deform more easily than the thinner ones.Clearly, this is due to the fact that with decreasing thickness of the compliant substrate, less soft solid is available to deform, thereby increasing the resistance to the deformation of the substrate and leading to smaller wetting ridges.Consequently, making the substrate thinner can be seen as equivalent to making it more rigid, whereas thicker substrates behave similarly to softer ones.This effect is also reflected in the mode of evaporation.
As it can be seen in Fig. 6b, evaporation takes place in CCR mode for the thicker, hence softer, substrate (H = 10 −1 ), and in CCA mode for the thinner, hence harder, substrate (H = 10 −3 ), in line with the findings shown in Fig. 5.

Evaporation of a pair of droplets
Now that we have studied the basic characteristics of the flow for a single sessile evaporating droplet, we may proceed with the examination of a system of multiple volatile droplets.In particular, we will investigate the dynamics of a pair of droplets and focus on the effects of their interaction, either through the soft substrate or their atmosphere, on the dynamics of the drying process.
In Fig. 7, we depict the time evolution of a pair of droplets evaporating on a compliant substrate with G = 1.In these simulations, we fully take into account the effect of thermocapillarity and examine two cases with M a = 0.001 and M a = 0.005 in Figs.7a and 7b, respectively.An interesting observation is that in both cases the droplets appear to move away from each other as they dry out.Regarding the deformation of the liquid-solid interface near the two contact lines, we observe that for low values of M a the height of the left and the right wetting ridge of each droplet is nearly symmetric (see inset of Fig. 7a), whereas for higher values of M a the droplet deforms asymmetrically with the deformation of the soft solid in the inner region between the two droplets being somewhat smaller than the deformation in the outside region (see inset of Fig. 7b).These observations provide a clear indication of the interaction of two droplets which may communicate either through the gas phase or through the developed stresses in the underlying viscoelastic substrate.

Effect of the gas phase and thermocapillarity
In order to shed light on the physical mechanisms behind the observed dynamics, we will first focus on the gas phase and depict in Fig. 8 the vapour concentration in the  atmosphere of the two droplets.As shown in this figure, the vapour concentration is higher between the two droplets than in their periphery.Since the evaporation flux is limited by diffusion (see Eq. 3.23), the higher saturation of the gas phase with vapour in the region between the two droplets results in weaker evaporation in that region; the spatial dependence of the evaporation flux J is plotted in Fig. 9a.For all values of P e v that we have examined, J acquires an asymmetric profile along the liquid-gas interface of each droplet; this can be seen more clearly by plotting ∂J/∂x in the inset of the same figure for P e v = 0.1.
Due to the effect of the latent heat, the evaporation flux affects the local interfacial temperature, which is depicted in Fig. 9b; the liquid-gas interface is cooler than the rest of the drop and the interfacial temperature is lowest at the droplet apex.The presence of temperature gradients affects in turn the flow field inside the droplet due to the action of Marangoni stresses, the spatial dependence of which, is plotted in Fig. 9c; the Marangoni stress is proportional to h∂σ/∂x.Focusing first on each droplet, we notice that the Marangoni stresses, exhibiting opposite signs in the regions left and right from the droplet apex, act as a compressive force reducing the footprint of the droplet.The effect of thermocapillarity on the droplet footprint is shown very clearly in Fig. 10a where we plot the length of footprint of the left drop, ∆x cl (see also Fig. 1c), defined as the distance between the maxima of left and right wetting ridge.Additionally, the asymmetric profile of the evaporation flux along the liquid-gas interface also induces a symmetry breaking in the interfacial temperature profile; this is clearly shown in the inset of Fig. 9b where we plot the spatial dependence of ∂T s /∂x for the droplet on the left side of the domain for P e v = 0.1.As a result, the Marangoni stresses not only compress the droplet but also contribute to their repulsion.This is demonstrated in Fig. 10b where we plot the evolution of the distance between the centers of mass of the droplets on the left and the right side of the domain, ∆x cm = x cm,r − x cm,l (see also Fig. 1c), with time for different values of M a .A similar effect is also shown in Fig. 9d where enhanced repulsion is found for lower values of P e v ; increase of P e v corresponds to slower vapour diffusion enhancing the difference in the evaporation flux between the two sides of the droplets as shown in Fig. 9a.Regarding the droplet lifetime, thermocapillarity plays a dual role; on one hand enhancing the evaporation rate in the region between the two droplets, as they move away from each other, but at the same time reducing the overall evaporation due to the compressive action of Marangoni stresses on the droplets.As it can be seen in Fig. 10c, where we plot the total mass of the system, the droplet lifetime increases considerably with M a , thus indicating that the latter effect is more significant.

Effect of substrate elasticity in the absence of thermocapillarity
In order to investigate the effect of the substrate elasticity, we plot in Fig. 11 the time evolution of a pair of droplets evaporating on solid substrates with G = 10, 500, 10 5 ; here, we neglect the effect of Marangoni stresses, i.e.M a = 0. Interestingly, we find that in the case of soft substrates the droplets repulse as they dry out (see Fig. 11a), whereas in the case of stiffer substrates the droplets are attracted to each other (see Fig. 11b).We note that in the latter case the droplets approach each other but do not coalesce; this behaviour is found in a well-define range of G (i.e.300 ⩽ G ⩽ 2 × 10 4 for M a = 0 and 2 × 10 3 ⩽ G ⩽ 5 × 10 4 for M a = 10 −4 ).For very stiff substrates (see Fig. 11c for G = 10 5 ), the two droplets eventually coalesce, and the drying process continues as for a single droplet.The dynamics for these three cases are also presented in the form of space-time plots in Fig. 12 (see panels 12a, 12c and 12d for G = 10, 500, 10 5 , respectively).
Inspecting the space-time plots presented in Fig. 12 and comparing against the work of Henkel et al. (2021) for non-volatile droplets, we notice a significant difference.For all the cases of volatile droplets that we have examined, varying the softness of the substrate, we find that mass transfer from one drop to the other does not take place and hence there is no droplet coarsening through Ostwald ripening mode; the latter mode was found to be dominant in the study of Henkel et al. (2021).For substrates with G = 10 5 where droplets coalesce, the coarsening process takes place with a typical translation mode as shown in Fig. 12d.After the two contact lines of the neighboring droplets come in contact, there is fast droplet coalescence.This is in contrast to the case shown in Fig. 12c (G = 500) where after the two adjacent contact lines touch each other, the droplets do not coalesce.Furthermore, we note that in parallel to the results of Henkel et al. (2021), the translation mode which defines exclusively the coarsening behaviour in Fig. 12d leads to a symmetric movement of the droplets towards each other, without displacing the system center of mass (see also Fig. 11c and Fig. 14b).To check whether our model can capture the  Figure 14: Time evolution of (a) the distance between the two centers of mass ∆x cm and (b) the center of mass of the system x cm,g , for M a = 0.The rest of the system parameters are the same with the 'base' case.
emergence of Ostwald ripening in the case of non-volatile droplets, we further examine in Fig. 13a the limit of negligible evaporation (E = 10 −8 ) for a substrate with high rigidity (G = 10 7 ); the equilibrium contact angle is taken to be approximately equal to 44 o (A = 120) in order to consider a system with similar wetting characteristics as in Henkel et al. (2021).As it is clearly shown in this figure, in the limit of non-volatile droplets the Ostwald ripening emerges at late times, as expected, with mass transferring from one drop to the other until the smaller droplet vanishes, and all the mass is contained in the larger remaining droplet.On the other hand, in Fig. 13b we examine the same system with two volatile droplets (E = 10 −4 ).As shown, evaporation takes place on much faster time scales and coalescence eventually occurs with the translation mode.
As explained by Karpitschka et al. (2016), interaction of non-volatile droplets on the surface of elastic solids is determined by the balance of elasticity and capillary forces and the resulting local deformation of the soft solid in the contact line region.Depending on the stiffness (or thickness) of the substrate, the elastic meniscus in the contact line region between the two droplets rotates by an angle as compared to the meniscus of an isolated drop and the direction of the rotation determines whether the drop-drop interaction is attractive or repulsive.In the case of drying droplets, though, the shape of the wetting ridge is not determined merely by elastocapillary phenomena but can also be significantly affected by the local evaporation rate (see relevant discussion in Charitatos & Kumar (2021)).Given the fact that the evaporation mass flux between the two contact lines of each droplet differs, this will also contribute to the imbalance of forces between the inner and outer contact lines, thereby affecting the mode of droplet interaction.To examine in more detail the complex droplet dynamics of our system, we plot in Fig. 14a the evolution of the distance between the center of mass of the droplets, ∆x cm , with time.It can be seen that in the case of soft substrates (G ⩽ 100) the droplets repulse, since ∆x cm continuously increases throughout evaporation.For harder substrates (G > 100), though, the imbalance acts in the opposite direction pushing the droplets towards each other.In the case of nearly rigid substrates (G ⩾ 10 6 ), the deformation of the viscoelastic solid is so small that this imbalance does not play an important role and thus the distance between the two droplets does not change significantly during evaporation.
The different modes of the drying process affect also the lifetime of the droplets.As shown in Fig. 11d, the evaporation is faster in softer substrates, due to the increased distance between the two droplets and the fact that less amount of vapor is trapped amidst the repulsing droplets leading to enhanced evaporation fluxes.In contrast, the greater amount of vapour trapped amidst the droplets when they attract in stiffer substrates, retards the evaporation significantly.Moreover, we notice that although for soft substrates the symmetry of the system is preserved throughout the drying process, this is not the case for substrates with intermediate stiffness.In fact, as it can be seen in Fig. 11b (and Fig. 12c), the pair of droplets at late stages of evaporation starts moving to the left exhibiting a clear symmetry breaking; the mechanisms for this behaviour will be investigated in detail below.Similarly, as shown in Fig. 11c (and Fig. 12d) for G = 10 5 , the droplet that has emerged after the coalescence of the two droplets appears to move slightly to the right, also indicating a symmetry breaking of the system, albeit with a somewhat smaller droplet displacement from the system center of mass.
As noted above, the elasticity of the substrate affects not only the relative distance between the droplets but may also lead under conditions to a symmetry breaking with the center of mass of the system, x cm,g being displaced from its initial position, i.e. the droplets appear to be 'walking' along the viscoelastic substrate.In Fig. 14b, we depict the effect of G on the evolution of the position of the system center of mass, x cm,g .As shown in this figure, for very soft and very hard substrates (i.e.G = 10 and G ⩾ 10 6 ) the center of mass of the system remains at x cm,g = 8 and the symmetry is preserved throughout the drying process.This is not the case, though, for substrates with intermediate stiffness where symmetry breaking is found; we note that the system symmetry is considered broken when the center of mass of the system has moved ±10 % of its initial maximum height (i.e.0.1 dimensionless distance) from its initial position.It should be pointed that these asymmetric solutions are spontaneous and emerge due to disturbances of the numerical finite element scheme, while they appear to be stable with the increase in mesh resolution.To make sure that the symmetry breaking is not artificially introduced by the imposed boundary conditions, we varied the size of the domain or even applied periodic boundary conditions in the x-direction; these efforts are presented in detail in the Appendix B. As discussed therein, neither the type of imposed boundary conditions or the domain size qualitatively affect the observed droplet behaviour.It is important to note that a spontaneous symmetry breaking has been also a matter of interest in Next, we take into account the effect of thermocapillary stresses.In Fig. 15, we depict the droplet dynamics for M a = 0.005 and three different values of G = 1, 50, 500.Regardless of the elasticity strength of the substrate, it is shown that in all cases the droplets repulse from each other.This behaviour is markedly different from the one discussed previously in the absence of thermocapillary effects where the droplets are attracted to each other for substrates with intermediate or high values of G.As discussed in Fig. 10, the Marangoni stresses play a dual role both acting as to compress the droplet footprint and also contributing to droplet repulsion.As shown in Fig. 16a where we plot the distance ∆x cm between the two droplets for a wide range of G values for M a = 0.005, the latter contribution is dominant and thus always leading the droplets to repulse from each other at the early stages of the drying process.Nevertheless, we notice that at later stages and for substrates with intermediate values of G (i.e. for G = 50, 500, 2000, 10 4 ) the droplet distance eventually starts decreasing indicating that the droplets are attracted to each other.This can be attributed to the fact that, at these late stages of the drying process, the droplet distance has increased considerably allowing the vapour concentration to acquire a more uniform profile along the interface of each droplet, leading to a more uniform evaporation flux and in turn to smaller temperature gradients.As a result, the thermal Marangoni stresses are significantly reduced, while the capillary forces induced by the substrate elasticity become dominant and drive the droplets closer to each other.For harder substrates, the capillary forces, as explained above, are weaker due to the fact that the substrate is less susceptible to elastic deformations and therefore the droplets continue to repulse due to the action of Marangoni stresses throughout the drying process.
Interestingly, we also notice in Fig. 15c, i.e. for a substrate with intermediate stiffness (G = 500), that the droplets initially repulse, then they are attracted and eventually symmetry breaking takes place; in Fig. 15c, a dashed arrow is drawn to indicate the motion of each droplet.As shown in Fig. 16b where we plot the evolution of the global center of mass x cm,g with time, we find that the symmetry is preserved only for extremely soft and extremely stiff substrates, while for intermediate values of G the droplets appear to 'walk' along the substrate.It should be noted that the emergence of this symmetry breaking at late stages of evaporation takes place spontaneously (i.e. at no specific time instant) and there is no preferred direction; it is triggered by numerical disturbances and appears to be stable and persistent with the increase in mesh resolution.
In Fig. 16c, we make an effort to rationalize and elucidate the mechanisms responsible for the symmetry breaking, shown in Fig. 15c (i.e. for M a = 0.005 and G = 500).In order to examine the contribution of various forces, i.e. capillary, Marangoni and elastic, on the motion of the droplets, we evaluate their contributions to the mean velocity in the x-direction of each droplet, as follows (5.1) The terms v cap,l and v cap,s denote the contributions from the capillary forces along the liquid-gas and liquid-solid interfaces, respectively, while the term v el corresponds to the contribution of the elastic stresses and v M a corresponds to the contribution of the Marangoni stresses; the analytical expressions for the various contributions are given in the Appendix C. We note that for a number of different cases that we have examined the dominant contributions come from the capillary forces and the Marangoni stresses along the liquid-gas interface, evaluated by the terms v cap,l and v M a respectively; v cap,s and v el were typically found to be two orders of magnitude smaller than v cap,l and v M a , and thus neglected here.Nevertheless, it is important to note that despite the fact that v el is typically very small, substrate elasticity implicitly contributes to the effect of capillary forces of the liquid-gas interface through the induced deformation of the contact line region.The time evolution of v cap,l and v M a is depicted in Fig. 16c for both droplets; indexes 1 and 2 correspond to the droplet on the left and right, respectively.At early times (i.e.approximately for t < 800), the contribution of the capillary and Marangoni stresses have similar magnitudes in both droplets and the symmetry is preserved (see Fig. 16b).The capillary forces act antagonistically with the Marangoni stresses pushing the droplets in opposite directions.The droplets, however, repulse due to the slightly higher magnitude of the Marangoni contribution.At later times (i.e. for t > 800), a disturbance in the local deformation of the solid causes an imbalance between the two droplets (see Fig. 16c) leading to symmetry breaking and driving the system center of mass away from its initial position.From Fig. 16b, it becomes evident that whether some disturbance will lead to a symmetry breaking or not, is a matter of the substrate elasticity.On the one hand, when the substrate is soft, it is very flexible and its deformation is very large.The size of an arising disturbance is insignificant compared to the size of the total substrate deformation and, as a result, the system center of mass will remain constant and the symmetry will be preserved.On the other hand, in extremely stiff substrates, the deformation of the liquid-solid interface is very small, quickly damping any possible disturbance that could lead to an imbalance between the two droplets.However, at intermediate values of substrate elasticity, there can be a competition between this disturbance and the substrate deformation, which might eventually lead to an imbalance in the induced capillary and Marangoni stresses between the two droplets and thus to a symmetry breaking, if the size of the disturbance grows considerably as compared to the substrate deformation.

Conclusions
In this paper we have studied the two-dimensional dynamics of a system of one or two droplets evaporating on a viscoelastic solid substrate.Lubrication theory is used to simplify the equations of mass, momentum, energy and the force balances applied in the liquid and the solid phases, considering the Kelvin-Voigt model to account for substrate viscoelasticity.Our model takes into consideration the effect of thermal Marangoni stresses, as well as the droplet interaction through both the compliant substrate and the surrounding vapour.The model accounts for the presence of the vapour employing a two-sided approach and considering the diffusion-limited model.The contact line is modelled assuming a precursor film ahead of the droplet.
We have carried out a parametric study to investigate how the evaporation process, the flow dynamics and the interaction of droplets are affected by the physical properties of the compliant substrate (e.g.thickness, shear modulus) and vapour diffusion in the atmosphere affecting the local evaporation rate.In the case of a single droplet, it was found that for thinner substrates the elastic effects become decreasingly important and thus making the substrate thinner can be seen as equivalent to making it more rigid.Moreover, it is shown that on softer (or thicker) substrates the solid deforms affecting the wetting of the droplet and promoting evaporation in CCR mode, in line with experimental observations in the literature (Lopes & Bonaccurso 2012, 2013;Yu et al. 2013;Gerber et al. 2019); the CCA mode is observed for harder (or thinner) substrates.Lastly, the effect of evaporative cooling and the action of thermocapillary stresses lead to smaller droplet footprints, resulting in an overall decrease of the evaporation rate, capturing the trend observed in earlier studies (Talbot et al. 2012;Schofield et al. 2018;Dunn et al. 2009b) in the case of rigid substrates.On the other hand, in the case of a system of a pair of volatile droplets, it is shown that the droplets may communicate both through the viscoelastic substrate and the induced deformations of the liquid-solid interface and also through the vapour that diffuses in the atmosphere of the droplets.The delicate interplay between the elastic stresses in the substrate, the capillary pressure and the thermal Marangoni stresses determine the mode of droplet interaction.
To summarize the rich dynamics of this complex system, we produced the parametric map of the dynamic regimes depicted in Fig. 17, varying the values of Marangoni number, M a , and substrate elasticity G.In order to characterize the different regimes of the map, we consider that the system symmetry is broken when the center of mass of the system has moved ±10 % of the initial maximum droplet height (i.e.0.1 dimensionless distance) from its initial position.In the absence of thermocapillary stresses (i.e.M a = 0), the droplets repulse on soft substrates (region I) whereas they are attracted to each other (and even coalesce) on stiff substrates (region VI) with symmetry breaking arising in the case of substrates with intermediate stiffness (region V).Droplet coalescence that takes place for substrates with intermediate stiffness closely resemble the translation mode of droplet coarsening observed for intermediate elasticity in Henkel et al. (2021).The translation mode is mediated by elastic deformation recovering the inverted Cheerios effect (Karpitschka et al. 2016).On very stiff and very soft substrates, however, the dominance of the Ostwald ripening effect found by Henkel et al. (2021) in the case of non-volatile droplets is significantly suppressed, due to the effect of evaporation that suppresses mass transfer between the droplets.Themocapillarity, on the other hand, apart from being responsible for the asymmetry of the wetting ridges on the two sides Figure 17: Map of the dynamic regimes depending on the value of Marangoni number, M a , and substrate elasticity G.The rest of the system parameters are the same with the 'base' case.We note that the borders in this map have been added as a visual guide and are not precise.
of each droplet between the two contact lines of each droplet, has also a drastic effect on the dynamics of the two droplets causing the repulsion of the droplets at the early stages of the drying process, irrespective of the stiffness of the viscoelastic solid (regions I, II and III).For substrates with intermediate stiffness, though, the droplet repulsion may also be followed by a phase of droplet attraction at later stages of the evaporating process (regions II and III), and under conditions to a symmetry breaking (region III); the symmetry is preserved either for very soft or very stiff substrates (regions I, II).Spontaneous symmetry breaking was also found in the study of Leong & Le (2020) who examined the coalescence of inflating droplets on viscoelastic substrates, indicating that these systems are prone to instability as a result of the delicate interplay amongst elastic, capillary and Marangoni forces.
Our findings clearly indicate that the flow dynamics can be very interesting with important implications for the optimal design of soft substrates for controlled evaporation of droplets.Our comprehensive model can be easily extended to more realistic setups such as the simulation of multiple 3D droplets or more complex systems such as the evaporation of particle-laden droplets.We believe that the present work should be complemented in the future with detailed experimental studies.To the best of our knowledge, such studies are lacking, despite the vast experimental work that exists on evaporating droplets on rigid substrates.by setting z = ζ(x, t).Using the expression for ∂vx ∂z | ζ from the tangential stress balance in the liquid phase, i.e.Eq. (3.9), we get: ) is determined by setting z = ξ(x, t) on Eq. (A 5) and using Eqs.(3.17), (4.1) and (4.3): Consequently, the x-component of the liquid velocity, v x , taking into account Eqs.(A 7) and (A 8) in Eq. (A 5), equals to: f 3 (x, t) is determined by setting z = ξ(x, t) on Eq. (A 6) and using Eqs.(3.18), (A 4) and (A 9): (A 10) Substituting Eq. (A 10) in Eq. (A 6) and using Eq.(A 9) we finally get for v z : Taking the material derivative of both sides of ξ = u z (x, 0, t) allows us to derive an evolution equation for ξ(x, t).The material derivative of ξ is derived using Eqs.(3.16) and (4.1), while the material derivative of u z (x, 0, t) is derived using Eq.(A 3):

Appendix B. Effect of domain size and boundary conditions
Here, we examine the effect of the boundary conditions on the observed system dynamics.To this end, we vary the size of the domain both in x and z directions, as well as considering the application of periodic boundary conditions in the x-direction.In Fig. B1 we examine the effect of the position of the far-field boundary condition in z-direction.As explained in section 2.3, the most natural choice would be to impose the boundary condition of constant vapour concentration at z → ∞.Here, however, we follow a similar approach to Schofield et al. (2020) considering a finite domain of the gas phase and the far-field condition is replaced by a similar Dirichlet condition at a distant, but finite, boundary.To examine the effect of this simplification, we vary the distance, d g , where the Dirichlet condition is imposed.In Fig. B1 we depict the effect on the evolution of both the total droplet mass and the position of the system center of mass.As shown in Fig. B1a, decreasing d g results in faster evaporation due the higher concentration gradient near the liquid-gas interface that is implicitly imposed and the fact that evaporation rate depends on the rate of diffusion; we observe a rather slow convergence of the total evaporation time with increasing values of d g .
In order to examine the effect of d g on the dynamics of a pair of droplets, we focus on the case presented in Fig. 11c for M a = 0 and G = 10 5 and d g = 5.In Fig. B1c and B1d we depict the same case for d g = 10 and 15, respectively, at times that correspond to the same scaled time, t ′ ; we note the great similarity between the three cases, which is also reflected on the evolution of the position of the center of mass, presented in Fig. B1b.Clearly, the position of the far-field condition does not affect the qualitative characteristics of either droplet coalescence or the observed symmetry breaking.
To investigate the effect of the boundary condition in the x-direction we follow two different routes.The first is to simply examine the effect of the domain length, L x for the same case examined in Fig. 11c.In Fig. B2 we see that the length mildly affects the predicted evaporation time but does not have a significant impact on the the rest qualitative characteristics of the flow (e.g.see Fig. B2b).Secondly, we solve a case for M a = 0 and G = 300 imposing periodic boundary conditions at the edges of the domain to check whether in cases where a symmetry breaking appears this might be due to the symmetry being already broken by the boundary conditions; the remaining parameters are kept the same with the 'base' case.As shown in Fig. B3, imposing either open or closed boundary conditions does not affect significantly the qualitative characteristics of the system dynamics and the symmetry breaking persists irrespective of the applied boundary conditions.

Appendix C. Mean velocity
In order to compute the mean velocity of each droplet we first computed the average x-velocity, v x,ave (using Eq.(A 9)), as follows (C 1) The total mean velocity in x-direction of each droplet can be evaluated by the following expression:  (C 2) The above integrals are solved between the two contact lines of each droplet.The total mean velocity of each one of the two droplets is the sum of the contribution of the capillary forces, the Marangoni stresses and the forces due to the elastic substrate.The contribution of the solid substrate, both in terms of capillary forces and in terms of elasticity is O(10 −4 ).More specifically, using Eq.(4.3) and Eq.(4.4), we get:

Figure 1 :
Figure 1: Schematic diagram of model geometry.(a) Initial configuration of a droplet with initial half-width l0 and initial height ĥ0 resting on an undeformed compliant substrate at ẑ = 0, which is attached to a rigid substrate at ẑ = − Ĥ.(b) The soft solid deforms while the system of two droplets spreads and evaporates.(c) Magnified view of the contact line, where β is the precursor film thickness, θ is the apparent contact angle and ξmax denotes the maximum height of the wetting ridge.The local thickness of each droplet is given by ĥ(x, t) = ζ(x, t) − ξ(x, t).
.23) where H = ρvi /ρ v ref denotes the relative humidity and P e v = λ∆ T R0 ĥ0 Dv ρv ref Lv .Solving Eqs.(3.21)-(3.23)we evaluate the vapour concentration in the gas phase, therefore making it possible to compute the vapour mass flux using the following dimensionless constitutive equation KJ = ρ ve − ρ v | ζ ′ .

Figure 2 :Figure 3 :
Figure 2: Time evolution of the liquid-air (ζ) and the liquid-solid (ξ) interfaces for a single droplet for (a) M a = 0 (t ev = 4242) and (b) M a = 5 • 10 −3 (t ev = 5366) respectively, for G = 3.The inset is an enlargement of the contact line region.The rest of the system parameters are the same with the 'base' case.

Figure 4 :
Figure 4: Gas phase concentration profiles at different time instants for M a = 5•10 −3 and G = 3.From top to bottom: t'=0, t'=0.02(left column), t'=0.4,t'=0.8 (right column) respectively.The rest of the system parameters are the same with the 'base' case.

Figure 5 :
Figure 5: Time evolution of (a) the point of maximum deformation of the wetting ridge ξ max , (b) the contact radius r and (c) the apparent contact angle θ for a single droplet, varying substrate elasticity G and for M a = 0.005.(d) Space-time plot of the droplet profiles at a soft substrate with G = 1 and for M a = 0.005.The inset is a magnified view of the wetting ridge profiles during droplet spreading.The rest of the system parameters are the same with the 'base' case.

Figure 6 :
Figure 6: Time evolution of (a) the point of maximum deformation of the wetting ridge ξ max and (b) the apparent contact angle θ for a single droplet, varying substrate thickness H and for G = 1, M a = 0.005.The rest of the system parameters are the same with the 'base' case.

Figure 7 :
Figure 7: Time evolution of the liquid-air (ζ) and the liquid-solid (ξ) interfaces for 2 droplets drying on a soft substrate with G = 1 and for (a) M a = 0.001 (t ev = 4667) and (b) M a = 0.005 (t ev = 5889), respectively.The inset is an enlargement of the height-range of the contact line region of the left drop at t ′ = 0.4.The rest of the system parameters are the same with the 'base' case.

Figure 9 :
Figure 9: Effect of P e v on the spatial profile of (a) the evaporation rate J, (b) the interfacial temperature T s , (c) the Marangoni stresses, h ∂σ ∂x , at t ′ = 0.5.(d) The time evolution of the distance between the droplets' centers of mass, ∆x cm = x cm,r − x cm,l for G = 1, M a = 0.005 and the rest of the system parameters are the same with the 'base' case.The insets in panels (a) and (b) depict the spatial profiles for P e v = 0.1 of ∂J/∂x and ∂T s /∂x, respectively.

Figure 10 :
Figure 10: Time evolution of (a) the length of footprint of the left drop, ∆x cl , (b) the distance between the two centers of mass, ∆x cm , and (c) the system mass, for different values of M a and for G = 1.The rest of the system parameters are the same with the 'base' case.

Figure 11 :
Figure 11: Time evolution of the liquid-air (ζ) and the liquid-solid (ξ) interfaces for 2 droplets with (a) G = 10 (t ev = 4537), (b) G = 500 (t ev = 4934) and (c) G = 10 5 (t ev = 5048) respectively, for M a = 0.The inset is an enlargement of the height-range of the contact line region of the left drop at t ′ = 0.4.(d) Time evolution of the system mass, varying substrate elasticity G.The rest of the system parameters are the same with the 'base' case.

Figure 13 :
Figure 13: Space-time plots of the droplet profiles for (a) E = 10 −8 and (b) E = 10 −4 , for M a = 0 and G = 10 7 .The domain length is L = 48, A = 120 and the rest of the system parameters are the same with the 'base' case.

Figure 15 :
Figure 15: Time evolution of the liquid-air (ζ) and the liquid-solid (ξ) interfaces for 2 droplets with (a) G = 1 (t ev = 5889), (b) G = 50 (t ev = 5329) and (c) G = 500 (t ev = 5222) respectively, for M a = 0.005.The inset is an enlargement of the heightrange of the contact line region of the left drop at t ′ = 0.4.(d) Time evolution of the system mass varying substrate elasticity G.The rest of the system parameters are the same with the 'base' case.

Figure 16 :
Figure16: Time evolution of (a) the distance between the two centers of mass ∆x cm and (b) the center of mass of the system x cm,g , for M a = 0.005.(c) The contribution of the Marangoni stresses and capillary forces in the average x-velocity of droplet 1 and 2 for G = 500 and M a = 0.005.The rest of the system parameters are the same with the 'base' case.

Figure B1 :
Figure B1: Time evolution of (a) the system mass and (b) the center of mass of the system x cm,g varying the height d g of the gas phase.Time evolution of the liquid-air (ζ) and the liquid-solid (ξ) interfaces for 2 droplets with (c) dg = 10 (t ev = 6810) and (d) d g = 15 (t ev = 7993) respectively.In all panels M a = 0 and G = 10 5 .The rest of the system parameters are the same with the 'base' case.

Figure B2 :
FigureB2: Time evolution of (a) the system mass and (b) the center of mass of the system x cm,g varying the length L of the solid substrate, for M a = 0 and G = 10 5 .The initial value of x cm,g has been moved to L = 8 for all cases for presentational purposes.The rest of the system parameters are the same with the 'base'

Figure B3 :
Figure B3: Time evolution of (a) the system mass and (b) the center of mass of the system x cm,g applying open and closed boundary conditions.Time evolution of the liquid-air (ζ) and the liquid-solid (ξ) interfaces for 2 droplets with (c) open boundary conditions (t ev = 4882) and (d) closed boundary conditions (t ev = 4716) respectively.In all panels M a = 0 and G = 300.The rest of the system parameters are the same with the 'base' case.
cap,l + v cap,s + v M a + v el .