The shielding effect extends the lifetimes of two-dimensional sessile droplets

We consider the diffusion-limited evaporation of thin two-dimensional sessile droplets either singly or in a pair. A conformal-mapping technique is used to calculate the vapour concentrations in the surrounding atmosphere, and thus to obtain closed-form solutions for the evolution and the lifetimes of the droplets in various modes of evaporation. These solutions demonstrate that, in contrast to in three dimensions, in large domains the lifetimes of the droplets depend logarithmically on the size of the domain, and more weakly on the mode of evaporation and the separation between the droplets. In particular, they allow us to quantify the shielding effect that the droplets have on each other, and how it extends the lifetimes of the droplets.

. The diffusive mass fluxĴ from the free surface of the droplet then controls the evolution, and hence the lifetime, of the droplet.
In the limit in which the droplet is thin [9,12,13], the problem simplifies further because the profile of the droplet may be neglected when imposing the boundary conditions onĉ. Thus, the mathematical problem typically becomes that of solving forĉ in a half-space or other large domain, subject to appropriate mixed boundary conditions. Similar mixed boundary-value problems occur in physical contexts including elastostatics [14], electrostatics [15], thermostatics [16], and hydrodynamics [17].
A range of mathematical techniques can be deployed to solve such problems [14,18]; contributions go back at least as far as the work of Weber [19], who presented what is effectively the vapour concentration field induced by a thin circular droplet. Subsequent work has employed methods including separation of variables [20], orthogonal polynomial expansions [21], Fourier or Hankel transforms [22,23], and Green's functions [24].
In two dimensions, additional techniques become available, notably conformal mapping [16,25]. This makes two-dimensional analogues of droplet evaporation problems appealing from the modeller's point of view: although two-dimensional problems may be somewhat artificial, their greater tractability allows more thorough analysis to be carried out. However, in two dimensions there is a fundamental difficulty concerning the specification of appropriate boundary conditions [14], which we will overcome, in the spirit of the work of Yarin et al. [26], by considering a suitably relaxed boundary condition.
In practice, droplets rarely occur in isolation, and so it is important to understand how droplets evaporate in the presence of other evaporating droplets. Previous studies of the evaporation of multiple sessile droplets have employed a variety of experimental, numerical and analytical approaches [27][28][29][30][31][32][33][34][35][36][37][38]. The critical difference between the evaporation of single and of multiple droplets is the occurrence of the shielding effect, namely that the presence of other evaporating droplets increases the local vapour concentration, and so each droplet evaporates more slowly than it would in isolation.
Again, analogous problems have been studied in other physical contexts. The first relevant study by Greenwood [39] examined the interaction of large numbers of microcontacts in electric contact theory, treating them as independent at leading order, and introducing an interaction term at higher order. Similar approaches have since been applied to elastic punches [40] and flow through pores [41], and have been put on a more rigorous asymptotic basis [42][43][44]. All these studies essentially considered the equivalent of thin circular droplets in three dimensions; recent work has used a variety of approaches to investigate the closely related problem of the dissolution of immersed nanobubbles and nanodroplets [45][46][47][48].
In this study we consider the evaporation of thin two-dimensional sessile droplets. In Sect. 2 we consider the one-droplet problem. We present the governing equations (Sect. 2.1), show that the most apparently natural problem does not have a solution (Sect. 2.2), and then show that by considering a suitably relaxed boundary condition we can obtain a physically acceptable solution via a conformal-mapping technique (Sect. 2.3). We validate this solution against numerical simulations (Sect. 2.4), and use it to obtain closed-form solutions for the evolution and lifetimes of the droplet in various modes of evaporation (Sect. 2.5). We then develop asymptotic expressions for these lifetimes in a large domain (Sect. 2.5.4). In Sect. 3 we consider the two-droplet problem. We obtain a solution to this problem (Sect. 3.1), which we again validate against numerical simulations (Sect. 3.2), before using it to obtain closed-form solutions for the evolution and lifetimes of the droplets (Sect. 3.3). We develop asymptotic expressions for these lifetimes (Sect. 3.3.3), and use these expressions to compare the lifetimes of a single droplet and a pair of droplets in dimensional terms (Sect. 4).

Model
Consider a thin two-dimensional sessile droplet with constant surface tensionσ and densityρ, evaporating in the diffusion-limited regime. (For simplicity, we shall refer to the fluid throughout as a droplet; viewed in three dimensions it is more accurately described as a ridge or line.) Let it have semi-widthR(t), contact angleθ(t) and cross-sectional areaÂ(t). Using Cartesian co-ordinates (x,ŷ) with origin at the centre of the base of the droplet, the droplet evaporates into a surrounding atmosphere with constant coefficient of vapour diffusionD, vapour saturation concentrationĉ =ĉ sat , and ambient vapour concentrationĉ =ĉ ∞ (<ĉ sat ). The vapour concentration in the atmosphere is denoted byĉ(x,ŷ,t), and the diffusive mass flux from the surface of the droplet byĴ (x,t).
Assuming that the droplet is sufficiently small, the Eötvös-Bond number Eo =ρĝR 2 0 /σ will be small; under these conditions the free surface of the droplet is approximately parabolic and its cross-sectional area is given by The flux from the droplet is given by which may be evaluated at y = 0 due to the thinness of the droplet. Similarly, the saturation condition, c = 1, on the surface of the droplet may also be imposed on y = 0. The saturation condition on the droplet and the no-flux condition on the substrate thus become respectively. To complete the problem we require a suitable boundary condition to be imposed in the "far field"; this turns out to be non-trivial to specify.

Absence of a solution in an infinite half-space
The simplest problem to specify is evaporation into an infinite half-space, so we aim to solve (2) subject to the far-field condition c → 0 as as well as to a mixed boundary condition on y = 0 of the form Applying a cosine transform to (2) and imposing the far-field condition (6) leads to a solution of the form where the function A(u) is to be determined. Imposing the boundary condition (7) requires that The work of Sneddon [14,Sect. 4.5] shows that requiring regularity of c at the contact line x = R imposes the condition so specifying that the function f (x) is any positive constant is not an admissible boundary condition, and so, as could have been anticipated from the behaviour of the fundamental solution of Laplace's equation in two dimensions, the problem specified by (2), (5) and (6) has no solution. We note that (11) precludes not only solutions to the simplest problem in which the saturation concentration is constant on the droplet, but also solutions to more general problems in which it varies along the droplet surface due, for example, to changes in temperature [11,49].

Solution in a finite domain via conformal mapping
Since the most apparently natural problem does not have a solution, we instead look for a closely related analogue that does. We therefore consider a slightly modified problem in which the far-field condition (6) is replaced by a similar Dirichlet condition at a distant, but finite, boundary. We therefore aim to solve subject to the standard boundary conditions on y = 0, and the relaxed boundary condition While it is difficult to find an analytical solution in a domain that is exactly semi-circular, we can obtain a solution in a semi-elliptical domain that approaches a semi-circular shape as it becomes large.
We proceed using conformal mapping. Let maps the semi-infinite strip (u, v) ∈ (−1, 1) × (0, ∞) in the w-plane to the upper half of the z-plane. In particular, the rectangle (−1, 1)×(0, S) shown in Fig. 1a is mapped to the semi-ellipse with semi-major axis length √ Ψ 2 + R 2 and semi-minor axis length Ψ shown in Fig. 1b and given by where An important point to note is that the shape of the semi-elliptical domain in the z-plane given by (17) depends on R as well as on Ψ . Thus, in general, for a droplet whose semi-width changes as it evolves, the shape of the domain also changes. However, in the regime of most interest, Ψ R, in which the domain is large, Eq. (17) gives and so the domain is semi-circular with radius Ψ and independent of R up to O(Ψ −2 ) 1. In the rectangular domain in the w-plane we seek a harmonic function Φ(u, v) satisfying Solving the problem for Φ in the rectangular domain immediately gives the corresponding solution for c in the semi-elliptical domain. By inspection, the solution in the rectangular domain is given by so and the flux is given by In particular, the flux satisfies and so it has the same (integrable) square-root singularity at the contact line x = R as in the corresponding three-dimensional problem [10].

Numerical validation
In order to validate the solution obtained in Sect. 2.3 (i.e. in order to assess the accuracy of the solution and to quantify the effect of the non-circularity of the domain), we solved the problem in the semi-circular domain using COMSOL Multiphysics [50]. In Fig. 2 we compare these numerical solutions to the analytical solutions in the semi-elliptical domain given by (22) and (23).   Figure 2 also shows that the analytical solutions accurately capture the behaviour of the numerical solutions provided that Ψ is sufficiently large, which is exactly as expected since it is for smaller domains that the semi-circular and semi-elliptical domains are most different.

Evolution and lifetime of the droplet
The rate of change of the cross-sectional area (3) is given by the flux (23) integrated over the surface of the droplet, In order to integrate (25) to determine the evolution and lifetime of the droplet, we require additional information about the behaviour of R(t) and θ(t), i.e. we need to specify the mode in which the droplet is evaporating. The works by Stauber et al. [7,8] and Schofield et al. [9] describe various modes of evaporation for axisymmetric droplets. In the present work, we consider the two-dimensional analogues of three of these modes: the constant-radius (CR) mode, the constant-angle (CA) mode, and the stick-slide (SS) mode. (Throughout, for consistency with the threedimensional problem, we will refer to modes in which R is fixed as "constant-radius" modes, although strictly R is not the radius but the semi-width of the two-dimensional droplet.)  (22) and (23), and dashed lines denote the corresponding numerical solutions in the semi-circular domain

Constant-radius (CR) mode
In the constant-radius (CR) mode, R(t) ≡ R 0 = 1. Noting that θ(0) = θ 0 = 1, we may immediately integrate (25) to obtain Fig. 3 Evolution and lifetime of a single droplet evaporating in the CR and CA modes: a contact angle θ(t) in the CR mode given by (26), b semi-width R(t) in the CA mode given by (28), and c areas A(t) given by (26) and (28), plotted as functions of t for Ψ = 10, 100 and 1000 with the arrows indicating the direction of increasing Ψ ; d lifetimes t CR and t CA given by (27) and (29)  Hence the lifetime of a single droplet evaporating in the CR mode is Figure 3a, c, d shows the evolution and lifetime of a single droplet evaporating in the CR mode. The contact angle θ and the area A both decrease linearly with time t (Fig. 3a, c). As the size of the domain Ψ increases, the contact angle θ and the area A decrease more slowly, and so the lifetime t CR increases monotonically with Ψ (Fig. 3d). This is because in two dimensions the distance to the outer boundary sets the distance over which the concentration decays to zero, and thus controls the concentration gradient close to the droplet, as seen in (23). This strong role of the outer boundary is a fundamental difference from the corresponding three-dimensional problem, in which the distance to the outer boundary becomes irrelevant for a sufficiently large domain, and so a far-field boundary condition can be safely imposed "at infinity".

Constant-angle (CA) mode
In the constant-angle (CA) mode, θ(t) ≡ θ 0 = 1. Noting that R(0) = R 0 = 1, we may integrate (25) implicitly to obtain Hence the lifetime of a single droplet evaporating in the CA mode is which can be re-written as Figure 3b-d shows the evolution and lifetime of a single droplet evaporating in the CA mode. In contrast to the CR mode, the semi-width R and the area A now both decrease nonlinearly with time t (Fig. 3b, c). However, as in the CR mode, the lifetime t CA increases monotonically with Ψ (Fig. 3d). Figure 3d also illustrates, as (30) also shows, that due to its pinned contact lines, a droplet evaporating in the CR mode always has a larger surface area, and hence a larger total flux and thus a shorter lifetime, than the same droplet evaporating in the CA mode, i.e. t CR ≤ t CA for all Ψ .

Stick-slide (SS) mode
In the stick-slide (SS) mode, the contact line is initially pinned, while the contact angle decreases until it reaches its critical de-pinning (receding) value θ = θ * (0 ≤ θ * ≤ 1) at the de-pinning time t = t * . After de-pinning, the contact angle remains at its critical value, while the semi-width decreases until it reaches zero. Thus we have In the pinned (i.e. the CR) phase, 0 < t < t * , the droplet evolves according to (26), so that while in the de-pinned (i.e. the CA) phase, t * < t < t SS , the droplet evolves according to Combining (32) and (33), the lifetime of a single droplet evaporating in the SS mode is which can be re-written as Figure 4 shows the evolution and lifetime of a single droplet evaporating in the SS mode. The de-pinning time t * decreases linearly with θ * , while the lifetime t SS increases linearly with θ * (Fig. 4d). Comparing (27), (30) and  (35) shows that, as might have been anticipated, the lifetime of a droplet evaporating in the SS mode always lies between those of the same droplet in the CR and CA modes, i.e. t CR ≤ t SS ≤ t CA for all Ψ and θ * . In the limit θ * → 1 − the SS mode approaches the CA mode and thus t SS → t CA − , while in the limit θ * → 0 + the SS mode approaches the CR mode and thus t SS → t CR + .
We note that in Fig. 4a all of the curves for which θ * = 0 intersect at t = t CR , and from (27), (32) and (33), the semi-width of the droplet at this time, R(t CR ), satisfies Note that R(t CR ) is a monotonically decreasing function of Ψ which takes its maximum value R(t CR ) = 1/2 in the limit Ψ → 0 + and satisfies R(t CR ) → 0 + as Ψ → ∞.

Asymptotic behaviour of the lifetimes in a large domain, Ψ R
Consider the regime of most interest, Ψ R, in which the domain is large and approximately semi-circular, and the condition at the outer boundary corresponds most closely to a far-field condition. From (27), (29) and (34) we obtain Fig. 5 The quasi-semi-elliptical domain in the ζ -plane for the two-droplet problem respectively. Equations (37)

Two-droplet problem
We now consider the analogous two-droplet problem in the ζ -plane, where ζ = η + iξ . We assume that the droplets are identical, and use the initial semi-width of the droplets as the characteristic length scale in the nondimensionalisation. The droplets are located so that they have inner contact lines at η = ± I and outer contact lines at η = ±Ω, where Ω > I , as shown in Fig. 5. The cross-sectional area of each droplet is then given by

Solution in a finite domain via conformal mapping
Consider the conformal map from the z-plane to the ζ -plane. This maps the real interval (0, R) in the z-plane to the real interval (I, Ω) where Ω = √ I 2 + R 2 in the ζ -plane, preserving the saturation condition on the droplet. It maps the real interval (R,  c(x, y) for the one-droplet problem given by (22) when R = 1, and b c(η, ξ ) for the two-droplet problem given by (44) when I = 1 and Ω = 3. In both cases Ψ = 100 and the contours are shown with increments of 0.05 However, as in the one-droplet problem, in the regime of most interest, Ψ I , in which the domain is large, Eq.
and so the domain is again semi-circular with radius Ψ and independent of I and Ω up to O(Ψ −2 ) 1. The solution in the quasi-semi-elliptical domain is given by Figure 6 shows the contours of the vapour concentration c(η, ξ ) for the two-droplet problem given by (44) and the corresponding contours of c(x, y) for the one-droplet problem given by (22). In both cases, far from the droplet(s) the contours approach the (semi-elliptical or quasi-semi-elliptical) shape of the outer boundary, and near to the droplet(s) the contours approach the flat shape(s) of the droplet(s). For the two-droplet problem the concentration in the region between the droplets is increased relative to that in the one-droplet problem, and near to the droplets the concentration falls away more gradually than it does in the one-droplet problem, resulting in the shielding effect described in Sect. 1.
The flux is given by In particular, the flux satisfies confirming that it again has square-root singularities at both contact lines.

Numerical validation
As we did in the one-droplet problem, in order to validate the solution obtained in Sect. 3.1, we solved the two-droplet problem in the semi-circular domain using COMSOL Multiphysics [50]. In Fig. 7 we compare these numerical solutions to the analytical solutions in the quasi-semi-elliptical domain given by (44) and (45).  Fig. 7 shows that c takes its saturation value on the surface of the droplets and decreases monotonically to its ambient value at the edge of the domain, that J is singular at the contact lines x = I and x = Ω, and that the analytical solutions accurately capture the behaviour of the numerical solutions provided that Ψ is sufficiently large.
However, Fig. 7a-f also shows that c decreases monotonically to an (unsaturated) minimum value between the droplets, and that this value is an increasing function of Ψ : this latter behaviour reflects the smaller concentration gradients, and thus the higher concentrations, which occur near to the droplets in larger domains. In addition, Fig.  7g-i clearly illustrate the shielding effect that the droplets have on each other. In particular, as (46) shows, the flux near to the outer contact line is suppressed less by the presence of the other droplet, and so remains larger than the flux near to the inner contact line, resulting in the skewed flux profiles shown in Fig. 7g-i. In particular, the minimum value of the flux no longer occurs at the centre of each droplet (as it does in the one-droplet problem).

Evolution and lifetime of the droplets
Using the solution for the flux given by (45), the evolution and lifetime of the droplets are determined by In the one-droplet problem we investigated three modes of evaporation (namely the CR, CA and SS modes), but in the two-droplet problem there is a much richer variety of possible behaviours because any of the four contact lines may, in principle, move independently of the other three. In the present work, we consider four canonical behaviours, in each of which the droplets remain symmetric about the ξ -axis. Specifically, we consider the following modes of evaporation: 1. The constant-inner-and-outer-contact-line (CIO) mode: the inner and outer contact lines of both droplets are pinned at η = ± I (0) = ± I 0 and η = ± Ω(0) = ± Ω 0 as the droplets evaporate. 2. The constant-angle-centred (CAC) mode: both droplets evaporate with constant contact angle and remain centred at η = ± (I + Ω)/2 = ± (I 0 + Ω 0 )/2. 3. The constant-angle and constant-inner-contact-lines (CAI) mode: both droplets again evaporate with constant contact angle, but now their inner contact lines are pinned at η = ± I 0 . 4. The constant-angle and constant-outer-contact-line (CAO) mode: both droplets again evaporate with constant contact angle, but now their outer contact lines are pinned at η = ± Ω 0 .

Constant-inner-and-outer-contact-line (CIO) mode
In the CIO mode, the inner and outer contact lines of both droplets are pinned, I ≡ I 0 and Ω ≡ Ω 0 = I 0 + 2. We may then immediately integrate (47) to obtain Hence the lifetime of a pair of droplets evaporating in the CIO mode is  Figure 8 shows the evolution and the lifetime of a pair of droplets evaporating in the CIO mode. As for a single droplet in the CR mode, the contact angle θ and the area A both decrease linearly with time t (Fig. 8a, b) and the lifetime t CIO increases monotonically with Ψ (Fig. 8c). In addition, since the shielding effect is weaker, and hence evaporation is faster, when the droplets are more widely separated, the lifetime t CIO decreases monotonically with the separation between the droplets, 2I 0 (Fig. 8d).

Constant-angle (CAC, CAI, CAO) modes
In the CAC, CAI and CAO modes, the contact angle remains constant, θ(t) ≡ θ 0 = 1. The three modes are distinguished by the different behaviours of the contact lines.
In the constant-angle-centred (CAC) mode, the droplets remain centred at η = ± (I + Ω)/2 = ± (I 0 + Ω 0 )/2. It is therefore convenient to write where Γ = (I 0 + Ω 0 )/2 is the position of the centre of the right-hand droplet and Δ(t) = (Ω − I )/2 is its semi-width. We may then integrate (47) implicitly to obtain Hence the lifetime of a pair of droplets evaporating in the CAC mode is In the constant-angle-and-inner-contact-line (CAI) mode, the inner contact line is pinned, I ≡ I 0 . We may then integrate (47) implicitly to obtain where Hence the lifetime of a pair of droplets evaporating in the CAI mode is In the constant-angle-and-outer-contact-line (CAO) mode, the outer contact line is pinned, Ω ≡ Ω 0 . We may then integrate (47) implicitly to obtain Hence the lifetime of a pair of droplets evaporating in the CAO mode is (58) Figure 9 shows the evolution and the lifetime of a pair of droplets evaporating in the three constant-angle modes. The difference between the modes is most clearly seen in Fig. 9a, which shows the inner and outer contact lines moving towards the centre of the droplet in the CAC mode, the outer contact line moving inward in the CAI mode, and the inner contact line moving outward in the CAO mode. Despite these differences, the evolution of the area A, which decreases nonlinearly with t, is similar for all three modes (Fig. 9b). As in the CAI mode, the lifetimes t CAC , t CAI and t CAO increase monotonically with Ψ (Fig. 9c) and decrease monotonically with the separation between the droplets, 2I 0 (Fig. 9d).
As Fig. 9c, d shows, the lifetimes of the three constant-angle modes are very similar, and it is only when the separation between the droplets is small that the difference between them becomes important. Specifically, Fig.  9d shows that the difference between t CAC , t CAI and t CAO becomes negligible when I 0 5 (i.e. when the droplet separation is several times the width of the droplets). As the droplets evaporate, the droplet separation is smallest in the CAI mode and largest in the CAO mode, resulting in the slowest evaporation, and hence the longest lifetime, in the CAI mode and the fastest evaporation, and hence the shortest lifetime, in the CAO mode. This point is further illustrated in Fig. 10, which shows the lifetimes t CAI , t CAC , t CAO and t CIO plotted as functions of Ψ . In particular, Fig. 10 shows that as the droplet separation increases the lifetimes of the three constant-angle modes (but not that of CIO mode) become virtually indistinguishable. We will discuss the latter behaviour in more detail in Sect. 3.3.3 below.

Asymptotic behaviour of the lifetimes in a large domain, Ψ I
The results shown in Fig. 10 motivate us to derive asymptotic expressions for the lifetimes of the droplets when Ψ I . Noting the difference between closely-spaced and widely-separated droplets, we consider the regimes I 0 1 and I 0 1 separately. In the regime I 0 1 Ψ , the initial droplet separation is much less than the initial droplet semi-width. Equations (49), (52), (55) and (58) then yield respectively.
On the other hand, in the regime 1 (59)-(62) we obtain so that, as expected, the lifetime of the pair of droplets is half that of the single droplet, i.e.t area ∼t single /2. Alternatively, we can consider the same total cross-sectional area of fluid, arranged either as two closely-spaced or two widely-separated droplets. In both cases the droplets have initial semi-widthR 0 / √ 2. If the droplets are closely spaced then from (59)-(62) we obtain At leading order the lifetime of the pair of droplets is the same as that of the single droplet, but there is a negative O(1) correction because the two droplets have a larger total surface area than the single droplet. On the other hand, if the droplets are widely separated then from (63)-(66) we obtain At leading order the lifetime of the pair of droplets is again the same as that of the single droplet, but now there is a larger negative O(ln(Î 0 /R 0 )) correction due to a weaker shielding effect when the droplets are widely separated.
With these parameter values, the timescaleT ≈ 351 s and the lifetime of a single droplet ist single ≈ 567 s. The lifetime of a pair of droplets with the same total surface area as the single droplet ist area ≈ 283 s. The lifetime of two closely-spaced droplets with the same total cross-sectional area as the single droplet ist close ≈ 541 s, whereas the lifetime of two widely-separated droplets with the same total cross-sectional area as the single droplet ist wide ≈ 442 s if the droplets are separated by 2Î 0 = 2 cm, andt wide ≈ 356 s if the droplets are separated by 2Î 0 = 20 cm.

Discussion and conclusion
In this contribution, we considered the diffusion-limited evaporation of thin two-dimensional sessile droplets either singly or in a pair. This two-dimensional problem is qualitatively different from the corresponding three-dimensional problem because, in contrast to in three dimensions, in two dimensions the size of the domain remains important even when it is much larger than the width of the droplets; it is therefore not possible to obtain a solution to the two-dimensional problem with a far-field boundary condition imposed "at infinity". We therefore formulated a slightly modified problem in which the far-field condition was replaced by a relaxed condition at a distant, but finite, boundary. We then showed how a conformal-mapping technique may be used to calculate the vapour concentrations, and hence obtain closed-form solutions for the evolution and the lifetimes of the droplets in various modes of evaporation. These solutions demonstrate that in large domains the lifetimes of the droplets depend logarithmically on the size of the domain, and more weakly on the mode of evaporation and the separation between the droplets. In particular, they allowed us to quantify the shielding effect that the droplets have on each other, and how it extends the lifetimes of the droplets.
Although the present two-dimensional problem may be somewhat artificial, it has direct practical applications, for example to the inkjet printing of circuits [26], and may be realisable experimentally using a Hele-Shaw cell geometry. More fundamentally, it provides a rare opportunity to obtain a closed-form description of the behaviour of interacting droplets and to quantify the shielding effect. It therefore permits asymptotic and analytical insight into a class of problems of increasing scientific and industrial interest.