1 Introduction

The evaporation of sessile droplets has been studied extensively in recent years [1]. Topics of particular interest include particle transport during drying [2,3,4] and the resulting instabilities [5], nano-scale phenomena [6], and droplet lifetimes [7,8,9].

A physically realistic model of droplet evaporation must describe vapour transport in the surrounding atmosphere. In the simplest setting, this requires us to solve Laplace’s equation for the vapour concentration \(\hat{c}\) in the atmosphere, subject to a saturation condition \(\hat{c}=\hat{c}_{\text {sat}}\) on the surface of the droplet and a no-flux condition on the substrate [10, 11]. The diffusive mass flux \(\hat{J}\) 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 \(\hat{c}\). Thus, the mathematical problem typically becomes that of solving for \(\hat{c}\) 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).

2 One-droplet problem

2.1 Model

Consider a thin two-dimensional sessile droplet with constant surface tension \(\hat{\sigma }\) and density \(\hat{\rho }\), 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-width \(\hat{R}(\hat{t})\), contact angle \(\hat{\theta }(\hat{t})\) and cross-sectional area \(\hat{A}(\hat{t})\). Using Cartesian co-ordinates \((\hat{x},\hat{y})\) with origin at the centre of the base of the droplet, the droplet evaporates into a surrounding atmosphere with constant coefficient of vapour diffusion \(\hat{D}\), vapour saturation concentration \(\hat{c}=\hat{c}_\text {sat}\), and ambient vapour concentration \(\hat{c}=\hat{c}_\infty \, (<\hat{c}_\text {sat})\). The vapour concentration in the atmosphere is denoted by \(\hat{c}(\hat{x},\hat{y},\hat{t})\), and the diffusive mass flux from the surface of the droplet by \(\hat{J}(\hat{x},\hat{t})\).

Following the approach of [9, 12, 13], we nondimensionalise and scale according to

$$\begin{aligned} \begin{aligned}&(\hat{x},\hat{y})=\hat{R}_0(x,y), \quad \hat{R}=\hat{R}_0R, \quad \hat{\theta }=\hat{\theta }_0\theta , \quad \hat{A}=\hat{R}_0^2\hat{\theta }_0A, \\&\hat{c}=\hat{c}_\infty +(\hat{c}_\text {sat}-\hat{c}_\infty )c, \quad \hat{J}=\frac{\hat{D}\left( \hat{c}_\text {sat}-\hat{c}_\infty \right) }{\hat{R}_0}J, \quad \hat{t}=\frac{\hat{\rho }\hat{\theta }_0\hat{R}^2_0}{\hat{D}\left( \hat{c}_\text {sat}-\hat{c}_\infty \right) }t, \end{aligned} \end{aligned}$$
(1)

where \(\hat{R}_0=\hat{R}(0)\) and \(\hat{\theta }_0=\hat{\theta }(0)\).

The vapour concentration is assumed to be quasi-steady, and so c satisfies Laplace’s equation

$$\begin{aligned} \nabla ^2c = 0, \end{aligned}$$
(2)

throughout the atmosphere.

Assuming that the droplet is sufficiently small, the Eötvös–Bond number \(Eo = \hat{\rho }\hat{g}\hat{R}^2_0/\hat{\sigma }\) will be small; under these conditions the free surface of the droplet is approximately parabolic and its cross-sectional area is given by

$$\begin{aligned} A=\frac{2}{3}R^2\theta . \end{aligned}$$
(3)

The flux from the droplet is given by

$$\begin{aligned} J=- \frac{\partial c}{\partial y} \quad \text {for} \quad |x| \le R, \end{aligned}$$
(4)

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

$$\begin{aligned} c(x,0)=1 \quad \text {for} \quad |x| < R, \qquad \frac{\partial c}{\partial y} (x,0)=0 \quad \text {for} \quad |x| > R, \end{aligned}$$
(5)

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.

2.2 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

$$\begin{aligned} c\rightarrow 0 \qquad \text {as} \qquad x^2+y^2\rightarrow \infty \quad \text {in} \quad y > 0, \end{aligned}$$
(6)

as well as to a mixed boundary condition on \(y=0\) of the form

$$\begin{aligned} c(x,0)=f(x) \, (>0) \quad \text {for} \quad |x| < R, \qquad \frac{\partial c}{\partial y} (x,0)=0 \quad \text {for} \quad |x| > R. \end{aligned}$$
(7)

Applying a cosine transform to (2) and imposing the far-field condition (6) leads to a solution of the form

$$\begin{aligned} c=\int _0^{\infty } u^{-1}A(u)\text {e}^{-uy}\cos (ux)\,\text {d}u, \end{aligned}$$
(8)

where the function A(u) is to be determined. Imposing the boundary condition (7) requires that

$$\begin{aligned}&\int _0^\infty u^{-1}A(u)\cos (xu)\,\text {d}u =f(x) \quad \text {for} \quad |x|<R, \end{aligned}$$
(9)
$$\begin{aligned}&\int _0^\infty A(u)\cos (xu)\,\text {d}u =0 \quad \text {for} \quad |x|>R. \end{aligned}$$
(10)

The work of Sneddon [14, Sect. 4.5] shows that requiring regularity of c at the contact line \(x=R\) imposes the condition

$$\begin{aligned} \int _0^R\frac{f(x)}{\sqrt{R^2-x^2}}\,\text {d}x=0, \end{aligned}$$
(11)

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].

2.3 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

$$\begin{aligned} \nabla ^2 c=0 \quad \text {in} \quad y>0, \quad x^2+y^2<\gamma ^2, \end{aligned}$$
(12)

subject to the standard boundary conditions on \(y=0\),

$$\begin{aligned} c(x,0)=1 \quad \text {for} \quad |x|< R, \qquad \frac{\partial c}{\partial y} (x,0)=0 \quad \text {for} \quad R<|x|<\gamma , \end{aligned}$$
(13)

and the relaxed boundary condition

$$\begin{aligned} c=0 \quad \text {for} \quad y>0, \quad x^2+y^2=\gamma ^2. \end{aligned}$$
(14)

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

$$\begin{aligned} z=x+\mathrm{i}y, \quad w=u+\mathrm{i}v. \end{aligned}$$
(15)

Then the mapping

$$\begin{aligned} z=g(w)=-R\,\cos \left( \frac{\pi }{2}(w+1)\right) \end{aligned}$$
(16)

maps the semi-infinite strip \((u,v)\in (-1,1)\times (0,\infty )\) in the w-plane to the upper half of the z-plane. In particular, the rectangle \((-1,1)\times (0,S)\) shown in Fig. 1a is mapped to the semi-ellipse with semi-major axis length \(\sqrt{\varPsi ^2+R^2}\) and semi-minor axis length \(\varPsi \) shown in Fig. 1b and given by

$$\begin{aligned} z=\sqrt{\varPsi ^2+R^2}\cos (s)+\mathrm{i}\varPsi \sin (s) \quad \text {for} \quad 0\le s\le \pi , \end{aligned}$$
(17)

where

$$\begin{aligned} \varPsi =R\sinh \left( \frac{\pi S}{2}\right) . \end{aligned}$$
(18)

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 \(\varPsi \). 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, \(\varPsi \gg R\), in which the domain is large, Eq. (17) gives

$$\begin{aligned} z = \varPsi \text {e}^{\mathrm{i}s}\left[ 1 + O(\varPsi ^{-2})\right] , \end{aligned}$$
(19)

and so the domain is semi-circular with radius \(\varPsi \) and independent of R up to \(O(\varPsi ^{-2}) \ll 1\).

Fig. 1
figure 1

a The rectangular domain in the w-plane and b the semi-elliptical domain in the z-plane for the one-droplet problem

In the rectangular domain in the w-plane we seek a harmonic function \(\varPhi (u,v)\) satisfying

$$\begin{aligned} \varPhi (u,0) = 1, \quad \dfrac{\partial \varPhi }{\partial u}(-1,v) = 0 = \dfrac{\partial \varPhi }{\partial u}(1,v), \quad \varPhi (u,S) = 0. \end{aligned}$$
(20)

Solving the problem for \(\varPhi \) 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

$$\begin{aligned} \varPhi (u,v)=1-\frac{v}{S}=1-\frac{\mathfrak {I}(w)}{S}, \end{aligned}$$
(21)

so

$$\begin{aligned} c(x,y)=1-\frac{\mathfrak {I}\left( g^{-1}(z)\right) }{S}=1-\frac{1}{{{\,\mathrm{arcsinh}\,}}\left( \varPsi /R\right) }\mathfrak {I}\left[ \arccos \left( -\dfrac{z}{R}\right) \right] , \end{aligned}$$
(22)

and the flux is given by

$$\begin{aligned} J(x)=- \frac{\partial c}{\partial y} (x,0)=\frac{1}{{{\,\mathrm{arcsinh}\,}}\left( \varPsi /R\right) }\frac{1}{\sqrt{R^2-x^2}} \quad \text {for} \quad |x|<R. \end{aligned}$$
(23)

In particular, the flux satisfies

$$\begin{aligned} J \sim \frac{1}{\sqrt{2R}{{\,\mathrm{arcsinh}\,}}\left( \varPsi /R\right) }\frac{1}{\sqrt{R-x}} \quad \text {as} \quad x\rightarrow R^-, \end{aligned}$$
(24)

and so it has the same (integrable) square-root singularity at the contact line \(x=R\) as in the corresponding three-dimensional problem [10].

Fig. 2
figure 2

The vapour concentration and flux for the one-droplet problem for \(R=1\) and \(\varPsi =2\) [(a), (d), (g)], \(\varPsi =10\) [(b), (e), (h)], and \(\varPsi =100\) [(c), (f), (i)]. ac Solutions for the vapour concentration along the x-axis, c(x, 0); df Solutions for the vapour concentration along the y-axis, c(0, y); gi Solutions for the flux, J(x). Solid lines denote the analytical solutions in the semi-elliptical domain given by (22) and (23), and dashed lines denote the corresponding numerical solutions in the semi-circular domain

2.4 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 2a–c shows solutions for the vapour concentration along the x-axis, c(x, 0); d–f shows solutions for the vapour concentration along the y-axis, c(0, y); g–i shows solutions for the flux, J(x). The first column [(a), (d), (g)] shows results for \(\varPsi =2\), the second column [(b), (e), (h)] shows results for \(\varPsi =10\), and the third column [(c), (f), (i)] shows results for \(\varPsi =100\). In all cases the solid lines denote the analytical solutions in the semi-elliptical domain, and the dashed lines denote the corresponding numerical solutions in the semi-circular domain. Figure 2a–f shows that the vapour concentration c takes its saturation value \(c=1\) on the surface of the droplet and decreases monotonically to its ambient value \(c=0\) at the edge of the domain, and Fig. 2g–i shows that the flux J increases monotonically with distance from a minimum value at the centre of the droplet and is singular at the contact line \(x=R\). Figure 2 also shows that the analytical solutions accurately capture the behaviour of the numerical solutions provided that \(\varPsi \) 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.

2.5 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,

$$\begin{aligned} \frac{\text {d} A}{\text {d} t} =\frac{2}{3} \frac{\text {d} }{\text {d} t} \left[ R^2(t)\theta (t)\right] =-\int _{-R}^RJ(x)\,\text {d}x=-\frac{\pi }{{{\,\mathrm{arcsinh}\,}}(\varPsi /R)}. \end{aligned}$$
(25)

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 \(\theta (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 three-dimensional 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.)

Fig. 3
figure 3

Evolution and lifetime of a single droplet evaporating in the CR and CA modes: a contact angle \(\theta (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 \(\varPsi =10\), 100 and 1000 with the arrows indicating the direction of increasing \(\varPsi \); d lifetimes \(t_\text {CR}\) and \(t_\text {CA}\) given by (27) and (29) plotted as functions of \(\varPsi \). In c and d solid lines denote the CA mode and dashed lines denote the CR mode

2.5.1 Constant-radius (CR) mode

In the constant-radius (CR) mode, \(R(t)\equiv R_0=1\). Noting that \(\theta (0)=\theta _0=1\), we may immediately integrate (25) to obtain

$$\begin{aligned} \theta (t)=1-\frac{3\pi }{2{{\,\mathrm{arcsinh}\,}}\varPsi }t, \quad A(t)=\frac{2}{3}\left[ 1-\frac{3\pi }{2{{\,\mathrm{arcsinh}\,}}\varPsi }t\right] . \end{aligned}$$
(26)

Hence the lifetime of a single droplet evaporating in the CR mode is

$$\begin{aligned} t_\text {CR}=\frac{2}{3\pi }{{\,\mathrm{arcsinh}\,}}\varPsi . \end{aligned}$$
(27)

Figure 3a, c, d shows the evolution and lifetime of a single droplet evaporating in the CR mode. The contact angle \(\theta \) and the area \(A\) both decrease linearly with time t (Fig. 3a, c). As the size of the domain \(\varPsi \) increases, the contact angle \(\theta \) and the area \(A\) decrease more slowly, and so the lifetime \(t_\text {CR}\) increases monotonically with \(\varPsi \) (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”.

2.5.2 Constant-angle (CA) mode

In the constant-angle (CA) mode, \(\theta (t)\equiv \theta _0=1\). Noting that \(R(0)=R_0=1\), we may integrate (25) implicitly to obtain

$$\begin{aligned}&t=\frac{2}{3\pi }\left[ {{\,\mathrm{arcsinh}\,}}\varPsi -R^2(t){{\,\mathrm{arcsinh}\,}}\left( \frac{\varPsi }{R(t)}\right) +\,\varPsi \left( \sqrt{\varPsi ^2+1}-\sqrt{\varPsi ^2+R^2(t)}\right) \right] ,\nonumber \\&A(t)=\frac{2R^2(t)}{3}. \end{aligned}$$
(28)

Hence the lifetime of a single droplet evaporating in the CA mode is

$$\begin{aligned} t_\text {CA}=\frac{2}{3\pi }\left[ {{\,\mathrm{arcsinh}\,}}\varPsi +\varPsi \left( \sqrt{\varPsi ^2+1}-\varPsi \right) \right] , \end{aligned}$$
(29)

which can be re-written as

$$\begin{aligned} t_\text {CA}=t_\text {CR}+\frac{2\varPsi }{3\pi }\left( \sqrt{\varPsi ^2+1}-\varPsi \right) . \end{aligned}$$
(30)

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_\text {CA}\) increases monotonically with \(\varPsi \) (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_\text {CR}\le t_\text {CA}\) for all \(\varPsi \).

2.5.3 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 \(\theta =\theta ^*\) (\(0\le \theta ^*\le 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

$$\begin{aligned} R(t)\equiv 1 \quad \text {for}\quad 0<t<t^*,\quad \theta (t)\equiv \theta ^*\quad \text {for}\quad t^*<t<t_\text {SS}. \end{aligned}$$
(31)

In the pinned (i.e. the CR) phase, \(0<t<t^*\), the droplet evolves according to (26), so that

$$\begin{aligned} t^*=\frac{2(1-\theta ^*){{\,\mathrm{arcsinh}\,}}\varPsi }{3\pi }, \end{aligned}$$
(32)

while in the de-pinned (i.e. the CA) phase, \(t^*<t<t_\text {SS}\), the droplet evolves according to

$$\begin{aligned} t=t^*+\frac{2\theta ^*}{3\pi }\left[ {{\,\mathrm{arcsinh}\,}}\varPsi -R^2(t){{\,\mathrm{arcsinh}\,}}\left( \frac{\varPsi }{R(t)}\right) +\varPsi \left( \sqrt{\varPsi ^2+1}-\sqrt{\varPsi ^2+R^2(t)}\right) \right] . \end{aligned}$$
(33)

Combining (32) and (33), the lifetime of a single droplet evaporating in the SS mode is

$$\begin{aligned} t_\text {SS}=\frac{2}{3\pi }\left[ {{\,\mathrm{arcsinh}\,}}\varPsi +\theta ^*\varPsi \left( \sqrt{\varPsi ^2+1}-\varPsi \right) \right] , \end{aligned}$$
(34)

which can be re-written as

$$\begin{aligned} t_\text {SS}=t_\text {CR}+\frac{2\theta ^*\varPsi }{3\pi }\left( \sqrt{\varPsi ^2+1}-\varPsi \right) . \end{aligned}$$
(35)
Fig. 4
figure 4

Evolution and lifetime of a single droplet evaporating in the SS mode: a semi-width R(t), b contact angle \(\theta (t)\), and c area \(A(t)\) given by (31)–(33), plotted as functions of t for \(\theta ^{*}=0\), \(\theta ^*=1/4\), \(\theta ^*=1/2\), \(\theta ^*=3/4\) and \(\theta ^*=1\) with the arrows indicating the direction of increasing \(\theta ^*\); d de-pinning time \(t^*\) given by (32) (dashed line) and lifetime \(t_\text {SS}\) given by (34) (solid line) plotted as functions of \(\theta ^{*}\). In all cases \(\varPsi =100\)

Figure 4 shows the evolution and lifetime of a single droplet evaporating in the SS mode. The de-pinning time \(t^*\) decreases linearly with \(\theta ^*\), while the lifetime \(t_\text {SS}\) increases linearly with \(\theta ^*\) (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_\text {CR}\le t_\text {SS}\le t_\text {CA}\) for all \(\varPsi \) and \(\theta ^*\). In the limit \(\theta ^*\rightarrow 1^-\) the SS mode approaches the CA mode and thus \(t_\text {SS}\rightarrow {t_\text {CA}}^-\), while in the limit \(\theta ^*\rightarrow 0^+\) the SS mode approaches the CR mode and thus \(t_\text {SS}\rightarrow t_\text {CR}{^+}\).

We note that in Fig. 4a all of the curves for which \(\theta ^*\ne 0\) intersect at \(t=t_\text {CR}\), and from (27), (32) and (33), the semi-width of the droplet at this time, \(R(t_\text {CR})\), satisfies

$$\begin{aligned} R^2(t_\text {CR}){{\,\mathrm{arcsinh}\,}}\left( \frac{\varPsi }{R(t_\text {CR})}\right) =\varPsi \left( \sqrt{\varPsi ^2+1}-\sqrt{\varPsi ^2+R^2(t_\text {CR})}\right) . \end{aligned}$$
(36)

Note that \(R(t_\text {CR})\) is a monotonically decreasing function of \(\varPsi \) which takes its maximum value \(R(t_\text {CR}) = 1/2\) in the limit \(\varPsi \rightarrow 0^+\) and satisfies \(R(t_\text {CR}) \rightarrow 0^+\) as \(\varPsi \rightarrow \infty \).

2.5.4 Asymptotic behaviour of the lifetimes in a large domain, \(\varPsi \gg R\)

Consider the regime of most interest, \(\varPsi \gg 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

$$\begin{aligned} t_\text {CR}&=\frac{2}{3\pi }\ln (2\varPsi )+O\left( \frac{1}{\varPsi ^2}\right) , \end{aligned}$$
(37)
$$\begin{aligned} t_\text {CA}&=\frac{2}{3\pi }\ln (2\varPsi )+\frac{1}{3\pi }+O\left( \frac{1}{\varPsi ^2}\right) , \end{aligned}$$
(38)
$$\begin{aligned} t_\text {SS}&=\frac{2}{3\pi }\ln (2\varPsi )+\frac{\theta ^*}{3\pi }+O\left( \frac{1}{\varPsi ^2}\right) , \end{aligned}$$
(39)

respectively. Equations (37)–(39) show that the lifetimes of the droplets all depend logarithmically on the size of the domain, and differ by an amount of \(O(1)\) which depends on the mode of evaporation. The corrections at \(O(\varPsi ^{-2})\) are of the same order as the deviation of the domain from circularity.

3 Two-droplet problem

We now consider the analogous two-droplet problem in the \(\zeta \)-plane, where \(\zeta =\eta +\mathrm{i}\xi \). We assume that the droplets are identical, and use the initial semi-width of the droplets as the characteristic length scale in the non-dimensionalisation. The droplets are located so that they have inner contact lines at \(\eta =\pm \,I\) and outer contact lines at \(\eta =\pm \varOmega \), where \(\varOmega >I\), as shown in Fig. 5. The cross-sectional area of each droplet is then given by

$$\begin{aligned} A=\frac{\left( \varOmega -I\right) ^2\theta }{6}. \end{aligned}$$
(40)
Fig. 5
figure 5

The quasi-semi-elliptical domain in the \(\zeta \)-plane for the two-droplet problem

3.1 Solution in a finite domain via conformal mapping

Consider the conformal map

$$\begin{aligned} \zeta =\varGamma (z)=\sqrt{I^2+z^2} \end{aligned}$$
(41)

from the z-plane to the \(\zeta \)-plane. This maps the real interval (0, R) in the z-plane to the real interval \((I,\varOmega )\) where \(\varOmega =\sqrt{I^2+R^2}\) in the \(\zeta \)-plane, preserving the saturation condition on the droplet. It maps the real interval \((R,\sqrt{\varPsi ^2+R^2})\) in the z-plane to the real interval \((\varOmega ,\sqrt{\varPsi ^2+\varOmega ^2})\) in the \(\zeta \)-plane, preserving the no-flux condition on the substrate. It maps the imaginary interval \((0,\mathrm{i}I)\) in the z-plane to the real interval \((0,I)\) in the \(\zeta \)-plane: the symmetry condition on the imaginary axis in the z-plane now becomes a no-flux condition on the real interval \((0,I)\) in the \(\zeta \)-plane. With a suitable choice of branch cut, the equivalent regions in the left half of the z-plane are mapped into corresponding regions in the left half of the \(\zeta \)-plane. The outer boundary of the domain, given by the rectangle in the z-plane and the semi-ellipse (17) in the z-plane, is mapped to the quasi-semi-elliptical curve shown in Fig. 5 and given by

$$\begin{aligned} \zeta = \left[ I^2+\left( \sqrt{\varPsi ^2+\varOmega ^2-I^2}\,\cos (s)+\mathrm{i}\varPsi \sin (s)\right) ^2\right] ^{1/2} \quad \text {for} \quad 0 \le s \le \pi . \end{aligned}$$
(42)

However, as in the one-droplet problem, in the regime of most interest, \(\varPsi \gg I\), in which the domain is large, Eq. (42) gives

$$\begin{aligned} \zeta = \varPsi \text {e}^{\mathrm{i}s}\left[ 1+O(\varPsi ^{-2})\right] , \end{aligned}$$
(43)

and so the domain is again semi-circular with radius \(\varPsi \) and independent of \(I\) and \(\varOmega \) up to \(O(\varPsi ^{-2}) \ll 1\).

The solution in the quasi-semi-elliptical domain is given by

$$\begin{aligned} c(\eta ,\xi ) = c(\zeta )=c(w(z(\zeta ))) = 1-\frac{1}{{{\,\mathrm{arcsinh}\,}}\left( \varPsi /\sqrt{\varOmega ^2-I^2}\right) }\mathfrak {I}\left[ \arccos \left( -\sqrt{\frac{\zeta ^2-I^2}{\varOmega ^2-I^2}}\right) \right] . \end{aligned}$$
(44)
Fig. 6
figure 6

Contours of the vapour concentration ac(xy) for the one-droplet problem given by (22) when \(R=1\), and b\(c(\eta ,\xi )\) for the two-droplet problem given by (44) when \(I=1\) and \(\varOmega =3\). In both cases \(\varPsi =100\) and the contours are shown with increments of 0.05

Figure 6 shows the contours of the vapour concentration \(c(\eta ,\xi )\) for the two-droplet problem given by (44) and the corresponding contours of c(xy) 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

$$\begin{aligned} J(\eta )=- \frac{\partial c}{\partial \xi } (\eta ,0)=\frac{1}{{{\,\mathrm{arcsinh}\,}}\left( \varPsi /\sqrt{\varOmega ^2-I^2}\right) }\frac{\eta }{\sqrt{\varOmega ^2-\eta ^2}\sqrt{\eta ^2-I^2}}. \end{aligned}$$
(45)

In particular, the flux satisfies

$$\begin{aligned} J \sim \frac{1}{\sqrt{2(\varOmega ^2-I^2)}{{\,\mathrm{arcsinh}\,}}\left( \varPsi /\sqrt{\varOmega ^2-I^2}\right) } \times {\left\{ \begin{array}{ll} \displaystyle \sqrt{\frac{I}{\eta -I}} &{} \text {as} \quad \eta \rightarrow I^{+}, \\ \displaystyle \sqrt{\frac{\varOmega }{\varOmega -\eta }} &{} \text {as} \quad \eta \rightarrow \varOmega ^{-}, \end{array}\right. } \end{aligned}$$
(46)

confirming that it again has square-root singularities at both contact lines.

3.2 Numerical validation

Fig. 7
figure 7

The vapour concentration and flux for the two-droplet problem for \(I=1\), \(\varOmega =3\) and \(\varPsi =4\) [(a), (d), (g)], \(\varPsi =10\) [(b), (e), (h)], and \(\varPsi =100\) [(c), (f), (i)]. ac Solutions for the vapour concentration along the \(\eta \)-axis, \(c(\eta ,0)\); df Solutions for the vapour concentration along the \(\xi \)-axis, \(c(0,\xi )\); gi Solutions for the flux, \(J(\eta )\). Solid lines denote the analytical solutions in the quasi-semi-elliptical domain given by (44) and (45), and dashed lines denote the corresponding numerical solutions in the semi-circular domain

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).

Figure 7a–c shows solutions for the vapour concentration along the \(\eta \)-axis, \(c(\eta ,0)\); d–f shows solutions for the vapour concentration along the \(\xi \)-axis, \(c(0,\xi )\); g–i shows solutions for the flux, \(J(\eta )\). The first column [(a), (d), (g)] shows results for \(\varPsi =4\), the second column [(b), (e), (h)] shows the corresponding results for \(\varPsi =10\), and the final column [(c), (f), (i)] shows the corresponding results for \(\varPsi =100\). In all cases the solid lines denote the analytical solutions in the quasi-semi-elliptical domain, and the dashed lines denote the corresponding numerical solutions in the semi-circular domain. As in the one-droplet problem, 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=\varOmega \), and that the analytical solutions accurately capture the behaviour of the numerical solutions provided that \(\varPsi \) 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 \(\varPsi \): 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).

3.3 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

$$\begin{aligned} \frac{\text {d} A}{\text {d} t} =\frac{1}{6} \frac{\text {d} }{\text {d} t} \left[ \left( \varOmega (t)-I(t)\right) ^2\theta (t)\right] =-\int _{I}^{\varOmega }J(\eta )\,\text {d}\eta =-\frac{\pi }{2{{\,\mathrm{arcsinh}\,}}\left( \varPsi /\sqrt{\varOmega ^2-I^2} \right) }. \end{aligned}$$
(47)

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 \(\xi \)-axis. Specifically, we consider the following modes of evaporation:

  1. 1.

    The constant-inner-and-outer-contact-line (CIO) mode: the inner and outer contact lines of both droplets are pinned at \(\eta =\pm I(0)=\pm I_0\) and \(\eta =\pm \varOmega (0)=\pm \varOmega _0\) as the droplets evaporate.

  2. 2.

    The constant-angle-centred (CAC) mode: both droplets evaporate with constant contact angle and remain centred at \(\eta =\pm (I+\varOmega )/2=\pm (I_0+\varOmega _0)/2\).

  3. 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 \(\eta =\pm I_0\).

  4. 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 \(\eta =\pm \varOmega _0\).

3.3.1 Constant-inner-and-outer-contact-line (CIO) mode

In the CIO mode, the inner and outer contact lines of both droplets are pinned, \(I\equiv I_0\) and \(\varOmega \equiv \varOmega _0=I_0+2\). We may then immediately integrate (47) to obtain

$$\begin{aligned} \theta (t)=1-\frac{3\pi }{4{{\,\mathrm{arcsinh}\,}}\left( \varPsi /\sqrt{\varOmega _0^2-I_0^2}\right) }t,\quad A=\frac{(\varOmega _0-I_0)^2}{6}\left[ 1-\frac{3\pi }{4{{\,\mathrm{arcsinh}\,}}\left( \varPsi /\sqrt{\varOmega _0^2-I_0^2}\right) }t\right] . \end{aligned}$$
(48)

Hence the lifetime of a pair of droplets evaporating in the CIO mode is

$$\begin{aligned} t_\text {CIO}=\frac{4{{\,\mathrm{arcsinh}\,}}\left( \varPsi /\sqrt{\varOmega _0^2-I_0^2}\right) }{3\pi }. \end{aligned}$$
(49)
Fig. 8
figure 8

Evolution and the lifetime of a pair of droplets evaporating in the CIO mode: a contact angle \(\theta (t)\) and b area \(A(t)\) given by (48) plotted as functions of t for \(\varPsi =10\), \(\varPsi =100\) and \(\varPsi =1000\) with the arrows indicating the direction of increasing \(\varPsi \); lifetime \(t_\text {CIO}\) given by (49) plotted c as a function of \(\varPsi \) and d as a function of \(I_0\) when \(\varPsi =100\). In ac all curves are for \(I_0=1\) and \(\varOmega _0=3\)

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 \(\theta \) and the area \(A\) both decrease linearly with time t (Fig. 8a, b) and the lifetime \(t_\text {CIO}\) increases monotonically with \(\varPsi \) (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_\text {CIO}\) decreases monotonically with the separation between the droplets, \(2I_0\) (Fig. 8d).

3.3.2 Constant-angle (CAC, CAI, CAO) modes

In the CAC, CAI and CAO modes, the contact angle remains constant, \(\theta (t)\equiv \theta _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 \(\eta =\pm (I+\varOmega )/2=\pm (I_0+\varOmega _0)/2\). It is therefore convenient to write

$$\begin{aligned} I(t)=\varGamma -\varDelta (t), \quad \varOmega (t)=\varGamma +\varDelta (t), \end{aligned}$$
(50)

where \(\varGamma =(I_0+\varOmega _0)/2\) is the position of the centre of the right-hand droplet and \(\varDelta (t)=(\varOmega -I)/2\) is its semi-width. We may then integrate (47) implicitly to obtain

$$\begin{aligned} 18\pi \varGamma ^2t&= 24\varGamma ^2\left( {{\,\mathrm{arcsinh}\,}}\left( \frac{\varPsi }{2\sqrt{\varGamma }}\right) -\varDelta ^2{{\,\mathrm{arcsinh}\,}}\left( \frac{\varPsi }{2\sqrt{\varGamma \varDelta }}\right) \right) \nonumber \\&\quad +\sqrt{4\varGamma \varDelta +\varPsi ^2}\left( \varPsi ^2-2\varGamma \varDelta \right) \varPsi -\sqrt{4\varGamma +\varPsi ^2}\left( \varPsi ^2-2\varGamma \right) \varPsi . \end{aligned}$$
(51)

Hence the lifetime of a pair of droplets evaporating in the CAC mode is

$$\begin{aligned} t_\text {CAC}= \frac{1}{18\pi \varGamma ^2}\left[ 24\varGamma ^2{{\,\mathrm{arcsinh}\,}}\left( \frac{\varPsi }{2\sqrt{\varGamma }}\right) +\varPsi ^4-\sqrt{4\varGamma +\varPsi ^2}\left( \varPsi ^2-2\varGamma \right) \varPsi \right] . \end{aligned}$$
(52)

In the constant-angle-and-inner-contact-line (CAI) mode, the inner contact line is pinned, \(I\equiv I_0\). We may then integrate (47) implicitly to obtain

$$\begin{aligned} F_\text {CAI}(\varOmega )=F_\text {CAI}(\varOmega _0)+\frac{3\pi }{2}t, \end{aligned}$$
(53)

where

$$\begin{aligned} F_\text {CAI}(\varOmega )&= \frac{I_0^2}{4}\left[ 3{{\,\mathrm{arctanh}\,}}\left( \frac{\varPsi ^2-\varOmega I_0-I_0^2}{\varPsi \sqrt{\varPsi ^2+\varOmega ^2-I_0^2}}\right) -{{\,\mathrm{arctanh}\,}}\left( \frac{\varPsi ^2+\varOmega I_0-I_0^2}{\varPsi \sqrt{\varPsi ^2+\varOmega ^2-I_0^2}}\right) \right] \nonumber \\&\quad +\frac{\varOmega }{2}\left( 2I_0-\varOmega \right) {{\,\mathrm{arcsinh}\,}}\left( \frac{\varPsi }{\sqrt{\varOmega ^2-I_0^2}}\right) -\frac{\varPsi }{2}\sqrt{\varPsi ^2+\varOmega ^2-I_0^2} \nonumber \\&\quad +\varPsi I_0\ln \left( \varOmega +\sqrt{\varPsi ^2+\varOmega ^2-I_0^2}\right) . \end{aligned}$$
(54)

Hence the lifetime of a pair of droplets evaporating in the CAI mode is

$$\begin{aligned} t_\text {CAI}=\frac{2}{3\pi }\left[ F_\text {CAI}(I_0)-F_\text {CAI}(\varOmega _0)\right] . \end{aligned}$$
(55)

In the constant-angle-and-outer-contact-line (CAO) mode, the outer contact line is pinned, \(\varOmega \equiv \varOmega _0\). We may then integrate (47) implicitly to obtain

$$\begin{aligned} F_\text {CAO}(I)=F_\text {CAO}(I_0)+\frac{3\pi }{2}t, \end{aligned}$$
(56)

whereFootnote 1

$$\begin{aligned} F_\text {CAO}(I)&= \frac{\varOmega _0^2}{4}\left[ 3{{\,\mathrm{arctanh}\,}}\left( \frac{\varPsi ^2+\varOmega _0I+\varOmega _0^2}{\varPsi \sqrt{\varPsi ^2+\varOmega _0^2-I^2}}\right) -{{\,\mathrm{arctanh}\,}}\left( \frac{\varPsi ^2-\varOmega _0I+\varOmega _0^2}{\varPsi \sqrt{\varPsi ^2+\varOmega _0^2-I^2}}\right) \right] \nonumber \\&\quad +\frac{I}{2}\left( 2\varOmega _0-I\right) {{\,\mathrm{arcsinh}\,}}\left( \frac{\varPsi }{\sqrt{\varOmega _0^2-I^2}}\right) +\frac{\varPsi }{2}\sqrt{\varPsi ^2+\varOmega _0^2-I^2} \nonumber \\&\quad +\varPsi \varOmega _0\arctan \left( \frac{I}{\sqrt{\varPsi ^2+\varOmega _0^2-I^2}}\right) . \end{aligned}$$
(57)

Hence the lifetime of a pair of droplets evaporating in the CAO mode is

$$\begin{aligned} t_\text {CAO}=\frac{2}{3\pi }\left[ F_\text {CAO}(\varOmega _0)-F_\text {CAO}(I_0)\right] . \end{aligned}$$
(58)
Fig. 9
figure 9

Evolution and lifetime of a pair of droplets evaporating in the CAC (solid lines), CAI (dotted lines), and CAO (dashed lines) modes: a contact line positions \(I(t)\), \(\varOmega (t)\) and b area \(A(t)\) given by (40), (51), (53) and (56) plotted as functions of t for \(\varPsi =10\), \(\varPsi =100\) and \(\varPsi =1000\) with the arrows indicating the direction of increasing \(\varPsi \); lifetimes \(t_\text {CAC}\), \(t_\text {CAI}\) and \(t_\text {CAO}\) given by (52), (55) and (58) plotted c as functions of \(\varPsi \) and d as functions of \(I_0\) for \(\varPsi =100\). In ac all curves are for \(I_0=1\) and \(\varOmega _0=3\), while d also includes the leading-order behaviour in the asymptotic regime \(1 \ll I_0 \ll \varPsi \) given by \(4/(3\pi )\ln (\varPsi /\sqrt{I_0})\) (dot-dashed line)

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_\text {CAC}\), \(t_\text {CAI}\) and \(t_\text {CAO}\) increase monotonically with \(\varPsi \) (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_\text {CAC}\), \(t_\text {CAI}\) and \(t_\text {CAO}\) becomes negligible when \(I_0 > rsim 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_\text {CAI}\), \(t_\text {CAC}\), \(t_\text {CAO}\) and \(t_\text {CIO}\) plotted as functions of \(\varPsi \). 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.

Fig. 10
figure 10

Lifetimes of a pair of droplets evaporating in the CAI (top curve in each set), CAC, CAO and CIO (bottom curve in each set) modes given by (49), (52), (55) and (58) plotted as functions of \(\varPsi \) when \(I_0=0\) (top set of curves), \(I_0=10\) and \(I_0=100\) (bottom set of curves)

3.3.3 Asymptotic behaviour of the lifetimes in a large domain, \(\varPsi \gg I\)

The results shown in Fig. 10 motivate us to derive asymptotic expressions for the lifetimes of the droplets when \(\varPsi \gg I\). Noting the difference between closely-spaced and widely-separated droplets, we consider the regimes \(I_0 \ll 1\) and \(I_0 \gg 1\) separately.

In the regime \(I_0 \ll 1 \ll \varPsi \), the initial droplet separation is much less than the initial droplet semi-width. Equations (49), (52), (55) and (58) then yield

$$\begin{aligned} t_\text {CIO}&=\frac{4}{3\pi }\ln (\varPsi ) -\frac{2I_0}{3\pi }+O\left( {I_0}^2, \, \frac{1}{\varPsi ^2}\right) , \end{aligned}$$
(59)
$$\begin{aligned} t_\text {CAC}&=\frac{4}{3\pi }\ln (\varPsi )+\frac{1}{3\pi } -\frac{2I_0}{3\pi }+O\left( {I_0}^2, \, \frac{1}{\varPsi ^2}\right) , \end{aligned}$$
(60)
$$\begin{aligned} t_\text {CAI}&=\frac{4}{3\pi }\ln (\varPsi )+\frac{2}{3\pi } -\frac{4I_0}{3\pi }+O\left( {I_0}^2\log I_0, \, \frac{1}{\varPsi ^2}\right) , \end{aligned}$$
(61)
$$\begin{aligned} t_\text {CAO}&=\frac{4}{3\pi }\ln (\varPsi )+\frac{2}{\pi }\left( 1-\frac{4}{3}\ln {2}\right) +\frac{4I_0}{3\pi }\left( 1-2\ln {2}\right) +O\left( {I_0}^2, \, \frac{1}{\varPsi ^2}\right) , \end{aligned}$$
(62)

respectively.

On the other hand, in the regime \(1 \ll I_0 \ll \varPsi \), the initial droplet separation is much greater than the initial droplet semi-width. Equations (49), (52), (55) and (58) then yield

$$\begin{aligned} t_\text {CIO}&=\frac{4}{3\pi }\ln \left( \frac{\varPsi }{\sqrt{I_0}}\right) -\frac{2}{3\pi I_0}+O\left( \frac{1}{I_0^2}, \, \frac{I_0}{\varPsi ^2}\right) , \end{aligned}$$
(63)
$$\begin{aligned} t_\text {CAC}&=\frac{4}{3\pi }\ln \left( \frac{\varPsi }{\sqrt{I_0}}\right) +\frac{1}{3\pi }-\frac{2}{3\pi I_0}+O\left( \frac{1}{I_0^2}, \, \frac{I_0}{\varPsi ^2}\right) , \end{aligned}$$
(64)
$$\begin{aligned} t_\text {CAI}&=\frac{4}{3\pi }\ln \left( \frac{\varPsi }{\sqrt{I_0}}\right) +\frac{1}{3\pi }-\frac{4}{9\pi I_0}+O\left( \frac{1}{I_0^2}, \, \frac{I_0}{\varPsi ^2}\right) ,\end{aligned}$$
(65)
$$\begin{aligned} t_\text {CAO}&=\frac{4}{3\pi }\ln \left( \frac{\varPsi }{\sqrt{I_0}}\right) +\frac{1}{3\pi }-\frac{8}{9\pi I_0}+O\left( \frac{1}{I_0^2}, \, \frac{I_0}{\varPsi ^2}\right) , \end{aligned}$$
(66)

respectively.

Equations (59)–(66) show that, as in the one-droplet problem, in the regime \(\varPsi \gg I\), the lifetimes of all four modes depend logarithmically on the size of the domain.

When \(I_0 \ll 1\), Eqs. (59)–(62) show that the lifetimes depend on the mode of evaporation at \(O(1)\). The lifetime of the CIO mode is the shortest because, due to their pinned contact lines, the droplets in this mode have the largest surface area, and hence evaporate the fastest. Of the three constant-angle modes, the CAO mode has the shortest lifetime and the CAI mode the longest lifetime. This is because when the inner contact lines are pinned the droplets remain closer together and hence evaporate more slowly than in the CAC mode due to a stronger shielding effect, whereas when the outer contact lines are pinned the droplets move further apart and hence evaporate more quickly than in the CAC mode due to a weaker shielding effect.

On the other hand, when \(I_0 \gg 1\), Eqs. (63)–(66) show that the influence of the different behaviours of the contact lines on the three constant-angle modes is very weak and at \(O(1)\) the lifetimes of the CAC, CAI and CAO modes all coincide with each other, but differ from the lifetime of the CIO mode by \(1/(3\pi )\).

4 Comparison between the lifetimes of a single droplet and a pair of droplets

To provide further insight into the shielding effect the droplets have on each other, we compare the lifetimes of a single droplet and a pair of droplets in dimensional terms. For simplicity, we consider only the leading-order behaviour of the lifetimes in the regime of most interest in which the domain is large and approximately semi-circular with radius \(\hat{\varPsi }\).

Our reference point is the lifetime of a single droplet with initial semi-width \(\hat{R}_0\), which from (1), (37) and (38) is given by

$$\begin{aligned} \hat{t}_{\text {single}} \sim \dfrac{2}{3\pi }\ln \left( \dfrac{2\hat{\varPsi }}{\hat{R}_0}\right) \hat{T}, \qquad \text {where} \qquad \hat{T} = \dfrac{\hat{\rho }\hat{\theta }_0\hat{R}_0^2}{\hat{D}(\hat{c}_{\text {sat}}-\hat{c}_{\infty })}. \end{aligned}$$
(67)

A first natural comparison is with a pair of droplets, each with initial semi-width \(\hat{R}_0/2\), which together have the same total surface area as the single droplet (i.e. their initial separation is \(2\hat{I}_0=0\)). In this case the vapour concentration and flux are identical in the two problems. However, the lifetimes are not the same, because the cross-sectional area of the single droplet is twice the total cross-sectional area of the two droplets. Specifically, from (59)–(62) we obtain

$$\begin{aligned} \hat{t}_{\text {area}} \sim \dfrac{4}{3\pi }\dfrac{\hat{\rho }\hat{\theta }_0(\hat{R}_0/2)^2}{\hat{D}(\hat{c}_{\text {sat}}-\hat{c}_{\infty })} \ln \left( \dfrac{\hat{\varPsi }}{\hat{R}_0/2}\right) = \dfrac{1}{3\pi }\ln \left( \dfrac{2\hat{\varPsi }}{\hat{R}_0}\right) \hat{T}, \end{aligned}$$
(68)

so that, as expected, the lifetime of the pair of droplets is half that of the single droplet, i.e. \(\hat{t}_{\text {area}}\sim \hat{t}_{\text {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-width \(\hat{R}_0/\sqrt{2}\). If the droplets are closely spaced then from (59)–(62) we obtain

$$\begin{aligned} \hat{t}_{\text {close}} \sim \dfrac{4}{3\pi }\dfrac{\hat{\rho }\hat{\theta }_0(\hat{R}_0/\sqrt{2})^2}{\hat{D}(\hat{c}_{\text {sat}}-\hat{c}_{\infty })} \ln \left( \dfrac{\hat{\varPsi }}{\hat{R}_0/\sqrt{2}}\right) = \dfrac{2}{3\pi }\left[ \ln \left( \dfrac{2\hat{\varPsi }}{\hat{R}_0}\right) -\dfrac{1}{2}\ln {2}\right] \hat{T}. \end{aligned}$$
(69)

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

$$\begin{aligned} \hat{t}_{\text {wide}} \sim \dfrac{4}{3\pi }\dfrac{\hat{\rho }\hat{\theta }_0(\hat{R}_0/\sqrt{2})^2}{\hat{D}(\hat{c}_{\text {sat}}-\hat{c}_{\infty })} \ln \left( \dfrac{\hat{\varPsi }}{\hat{R}_0/\sqrt{2}}\sqrt{\dfrac{\hat{R}_0/\sqrt{2}}{\hat{I}_0}}\right) = \dfrac{2}{3\pi }\left[ \ln \left( \dfrac{2\hat{\varPsi }}{\hat{R}_0}\right) - \dfrac{1}{2}\ln \left( 2^{3/2}\dfrac{\hat{I}_0}{\hat{R}_0}\right) \right] \hat{T}. \end{aligned}$$
(70)

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 (\hat{I}_0/\hat{R}_0))\) correction due to a weaker shielding effect when the droplets are widely separated.

To illustrate these results we take the representative parameter values \(\hat{\rho }=998\) kg m\(^{-3}\), \(\hat{D} = 2.44\times 10^{-5}\) m\(^2\) s\(^{-1}\), \(\hat{c}_{\text {sat}} = 1.94\times 10^{-2}\) kg m\(^{-3}\) and \(\hat{c}_{\infty } = 7.76\times 10^{-3}\) kg m\(^{-3}\), corresponding to water at 295 K, evaporating into an environment with ambient vapour concentration \(\hat{c}_{\infty }=0.4\,\hat{c}_{\text {sat}}\) [11]. We further take \(\hat{\varPsi }=1\) m together with \(\hat{\theta }_0=0.1\) and \(\hat{R}_0 = 1\) mm, so that the droplet has a cross-sectional area of approximately \(6.7\times 10^{-8}\) m\(^2\).

With these parameter values, the timescale \(\hat{T} \approx 351\) s and the lifetime of a single droplet is \(\hat{t}_{\text {single}} \approx 567\) s. The lifetime of a pair of droplets with the same total surface area as the single droplet is \(\hat{t}_{\text {area}} \approx 283\) s. The lifetime of two closely-spaced droplets with the same total cross-sectional area as the single droplet is \(\hat{t}_{\text {close}} \approx 541\) s, whereas the lifetime of two widely-separated droplets with the same total cross-sectional area as the single droplet is \(\hat{t}_{\text {wide}} \approx 442\) s if the droplets are separated by \(2\hat{I}_0 = 2\) cm, and \(\hat{t}_{\text {wide}} \approx 356\) s if the droplets are separated by \(2\hat{I}_0 = 20 \hbox { cm}\).

5 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.