Contact-line deposits from multiple evaporating droplets

Building on the recent theoretical work of Wray, Duffy and Wilson [J. Fluid Mech. 884, A45 (2020)] concerning the competitive diffusion-limited evaporation of multiple thin sessile droplets in proximity to each other, we obtain theoretical predictions for the spatially non-uniform densities of the contact-line deposits (often referred to as"coffee stains"or"ring stains") left on the substrate after such droplets containing suspended solid particles have completely evaporated. Neighbouring droplets interact via their vapour fields, which results in a spatially non-uniform"shielding"effect. We give predictions for the deposits from a pair of identical droplets, which show that the deposit is reduced the most where the droplets are closest together, and demonstrate excellent quantitative agreement with experimental results of Pradhan and Panigrahi [Coll. Surf. A 482, 562-567 (2015)]. We also give corresponding predictions for a triplet of identical droplets arranged in an equilateral triangle, which show that the effect of shielding on the deposit is more subtle in this case.


I. INTRODUCTION
The evaporation of sessile droplets has been the subject of extensive experimental, numerical and analytical investigation in recent years (see, for example, [1][2][3][4][5], and the references therein), partly motivated by the wide range of everyday and industrial situations, such as protein crystallography [6], surface patterning [7], ink-jet printing, including that of OLED displays [8], and agrochemical spraying of plants [9], in which it occurs. Particular attention has been paid to the so-called "coffee-stain" or "ring-stain" effect: when a droplet of coffee (or indeed a droplet of any fluid containing suspended solid particles) with a pinned (i.e. a fixed) contact line evaporates it tends to deposit the majority of the particles close to the location of its contact line, even if the particles were initially distributed uniformly throughout the bulk of the droplet. The explanation of this phenomenon, as first given by Deegan et al. [10], is that as the droplet evaporates its free surface adjusts quasistatically under the effect of capillarity, inducing a flow within the droplet that advects the particles suspended within it towards its contact line, resulting in a characteristic ring-like contact-line deposit on the substrate after the droplet has completely evaporated. Since the seminal work by Deegan et al. [10], many aspects of this phenomenon have been investigated in considerable detail (see, for example, [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] and the reviews [2,26,27]). Note that although Deegan et al. [10] considered the most commonly studied situation of diffusionlimited evaporation into a quiescent atmosphere with a uniform far-field concentration of vapour, which has a large (theoretically singular) evaporative flux close to the contact line, the effect is quite robust, and even, for example, a spatially uniform evaporative flux will lead to advection of particles towards the contact line (see, for example, [12,14,23]).
The vast majority of the previous work on deposition from evaporating sessile droplets has, for obvious reasons, focused on axisymmetric deposits from axisymmetric droplets.
There has, however, been some work on non-axisymmetric deposits from non-axisymmetric droplets (see, for example, [9,10,[28][29][30][31][32]). In particular, Du and Deegan [29] examined a twodimensional droplet on an inclined substrate numerically, and found that, depending on the initial volume of the droplet and the angle of inclination of the substrate, the larger deposit can occur at either the upper or the lower contact line, while Sáenz et al. [31] investigated a variety of non-axisymmetric droplets both experimentally and numerically, and found that larger deposits occur where the contact line has the largest curvature (e.g. near the tips of a droplet with a triangular contact line). However, their theoretical modelling of the density of the deposit was essentially phenomenological.
Non-axisymmetric deposits also occur as a result of the non-axisymmetric evaporation of multiple droplets in proximity to each other, a situation that occurs much more commonly in practice than single droplets in isolation [8]. Specifically, neighbouring droplets undergoing diffusion-limited evaporation interact via their vapour fields, which results in a spatially non-uniform "shielding" effect that reduces the evaporation rate. While there have been some analytical studies of mathematically analogous situations concerning clusters of microcontacts and nanobubbles (see, for example, [33,34]), analytical work on the evaporation of multiple droplets is rather limited. To a large extent this is explained by the inherent difficulty of analysing such situations, and while the evaporation of multiple droplets in various configurations has been the subject of growing recent interest (see, for example, [12,[35][36][37][38][39][40][41][42][43][44][45][46][47][48]), the previous studies have been predominantly numerical or experimental. Two notable exceptions are the recent work of Wray et al. [48], who, building on the earlier work of Fabrikant [49] concerning a model for diffusion through a porous membrane, analysed the spatially non-uniform shielding that occurs in arbitrary configurations of thin droplets with circular contact lines, and that of Schofield et al. [47], who used conformal-mapping techniques to analyse the analogous spatially non-uniform shielding that occurs in the closely related two-dimensional situation of a pair of evaporating ridges. In particular, Wray et al. [48] gave explicit formulae for the evaporative flux of arbitrary configurations of droplets that were found to be remarkably accurate up to and including the limit of touching droplets, and led to theoretical predictions for the evolution of an arrangement of seven droplets that were found to be in excellent agreement with experimental results of Khilifi et al. [45].
In the present contribution we build on the work of Wray et al. [48] in order to analyse the spatially non-uniform densities of the deposits left on the substrate by the diffusion-limited evaporation of multiple thin droplets with pinned circular contact lines in proximity to each other. Specifically, in Secs. II and III we formulate and solve the evaporation, hydrodynamic, and particle-transport problems. In Sec. IV we give theoretical predictions for the densities of the deposits from a pair of identical droplets, and demonstrate excellent quantitative agreement with experimental results of Pradhan and Panigrahi [39]. In Sec. V we also give corresponding predictions for a triplet of identical droplets arranged in an equilateral triangle. Finally, we summarise our conclusions in Sec. VI. The present analysis is for the most commonly studied case of small droplets, in which capillary effects dominate over gravitational effects, corresponding to the limit of small Bond number. In Appendix A we describe the corresponding analysis for the less commonly studied case of large droplets, corresponding to the limit of large Bond number, in which even greater analytical progress is possible.

II. PROBLEM FORMULATION
Consider N (N = 1, 2, 3, . . .) thin axisymmetric sessile droplets with pinned circular contact lines with constant radii a k and fixed centres at (x k , y k ) for k = 1, 2, . . . , N on a planar solid substrate z = 0, as shown in Fig. 1(a). The droplets undergo quasi-static diffusion-limited evaporation, which, as described in Sec. I, induces flows within the droplets that advect the solid particles suspended within them towards their contact lines. The goal of the present work is to determine the spatially non-uniform densities of the deposits left on the substrate after the droplets have completely evaporated. Three coupled problems must therefore be solved: the evaporation problem for the concentration of vapour in the atmosphere (which determines the rates of evaporation of the droplets), the hydrodynamic problem for the fluid flow that is induced in each droplet, and the advection problem for the motion of the particles suspended within each droplet, as shown in Fig. 1(b). We now discuss each of these problems in turn.

A. The evaporation problem
According to the diffusion-limited model, the quasi-static concentration of vapour in the atmosphere, denoted by c (v) = c (v) (r, φ, z), satisfies Laplace's equation ∇ 2 c (v) = 0 subject to conditions of complete saturation at the free surfaces of the droplets and of no flux of vapour through the unwetted part of the substrate. Following Wray et al. [48], we scale and nondimensionalise variables appropriately for the atmosphere according to  ∞ are the constant saturation concentration and far-field concentration of vapour, and J k = J k (r, φ) and F k are the local evaporative flux and the integral evaporative flux from the k th droplet, respectively, which are related by where S k denotes the free surface of the k th droplet. Since the contact lines of the droplets are pinned, J k and F k are independent of time except for discontinuous jumps when any droplet completely evaporates; in particular, J k and F k jump instantaneously to zero when the k th droplet completely evaporates.
For clarity, we immediately drop the star superscripts on non-dimensional quantities, and so the boundary conditions on c (v) become c (v) = 1 on z = h k for k = 1, 2, . . . , N, ∂c (v) /∂z = 0 on the unwetted part of the substrate z = 0, and the far-field condition As Wray et al. [48] described, the earlier work of Fabrikant [49] on diffusion through a porous membrane, when interpreted in terms of the evaporation of multiple thin sessile droplets, shows that the integral evaporative flux F k is given, to a high degree of accuracy, by the solution of the linear system where r k,n (≥ a k + a n ) is the distance between the centres of the k th and the n th droplets shown in Fig. 1(a) and given by Wray et al. [48] also showed that, to the same high degree of accuracy, the local evaporative flux J k is given by where is the local evaporative flux from the k th droplet in isolation, and ψ k,n is the angle between the x axis and the line joining the centres of the k th and the n th droplets, also shown in Fig. 1(a) and given by B. The hydrodynamic problem The velocity and pressure within the k th droplet, denoted by u k = u k (r, φ, z, t) = u k e r + v k e φ + w k e z and p k = p k (r, φ, z, t), where t denotes time, satisfy the usual mass-conservation and Stokes equations subject to the usual boundary conditions, and the free surface, contact angle and volume of the k th droplet are denoted by z = h k = h k (r, t), θ k = θ k (t) (≪ 1) and We scale and nondimensionalise variables appropriately for the droplet according to respectively, where ρ is the constant density of the fluid.
At leading order in θ ref ≪ 1 the governing equations for the k th droplet are, with the hats dropped for clarity, where Ca is an appropriate capillary number, defined by where µ is the constant viscosity of the fluid. Equation (10) is to be solved subject to zero velocity at the substrate, balances of normal and tangential stress at the free surface of the droplet, and the kinematic condition, where Q are the local radial and azimuthal volume fluxes of fluid within the droplet.
We consider the situation in which capillary effects are strong, corresponding to small values of the capillary number Ca, and so we seek asymptotic solutions of the form in the limit Ca → 0. As we shall see, for the analysis of particle transport presented in Sec. II C, we require only the leading order velocity components u k0 and v k0 , which in turn require p k0 and p k1 .
At leading order in Ca ≪ 1, equations (10) and (13) show that the leading-order pressure is independent of r, φ and z, and is given by p k0 = p k0 (t) = 2θ k /a k , and the leading-order free surface z = h k (r, t) takes the familiar paraboloidal form The leading-order volume of the droplet is therefore given by At first order in Ca ≪ 1, equations (10), (13) and (12) lead to and so the leading-order local fluid fluxes (15) are Dropping the subscript "1" on p k1 henceforth for clarity, the kinematic condition (14) there- We may use this condition to obtain the differential equation satisfied by p k by noting that h k = h k (r, t) is independent of the azimuthal coordinate φ, and that, by global mass conservation, the volume V k satisfies so that, with (17) and (18), ∂h k /∂t may be written as Finally, therefore, the kinematic condition (21) may be expressed in the form which is a partial differential equation for p k , with all of the other quantities in (24) being known. Once p k is determined from (24), the local fluid fluxes Q (r) are given by (20).
The depth-averaged radial and azimuthal velocities, denoted byū k =ū k (r, φ, t) and v k =v k (r, φ, t), are defined bȳ respectively. For future reference, note that the streamlines of the depth-averaged flow are determined by solving Sinceū k andv k have the same functional dependence on h, the θ k has cancelled out of (26), and so the streamlines depend on time only via changes in the flux J k . This means that in certain situations the computation of the streamlines is simplified somewhat by the fact (mentioned earlier) that J k and hence F k are independent of t except for discontinuous jumps when any droplet completely evaporates. Specifically, if the droplets are arranged in such a way that all of them completely evaporate at exactly the same time then their streamlines remain unchanged throughout the evaporation. In this case the determination of the density of the deposit reduces to performing a single integral, as will be described in Sec. IV B below.

C. The particle-transport problem
The motion of the particles suspended within each droplet is due to a combination of advection by the flow and diffusion, and so the concentration of particles in the k th droplet, denoted by c k = c k (r, φ, z, t), satisfies the (scaled) advection-diffusion equation in which c k has been nondimensionalised according to is an appropriate Péclet number, and D (p) is the constant diffusion coefficient for the particles in the fluid, and where the star subscript has again been dropped for clarity. Equation (27) is subject to conditions of no flux of particles through either the free surface of the droplet or the substrate, where n k denotes the outward unit normal to the free surface of the k th droplet.
As is well known (see, for example, [22]), when the Péclet number Pe is such that where the depth-averaged radial and azimuthal velocitiesū k andv k are given by (25). Equation (29) may be solved by the method of characteristics: subject to a prescribed initial condition c k = c k (r, φ, 0) at t = 0. For simplicity, in all of the results presented below we assume that the initial concentration of particles takes the same uniform value in all of the droplets, which we may, without loss of generality, take to be unity.

III. SOLUTION FOR THE PRESSURE p k
In general, the expression for the evaporative flux J k given in (5) is rather complicated, and precludes solving (24) for the pressure p k in closed form; however, we may determine p k to arbitrary accuracy as follows.
For r ≤ a k (< r k,n ) we expand J k given by (5) as the convergent series which we may rearrange in the form of a truncated Fourier series for a chosen number of modes M (M = 1, 2, 3, . . .), with known functions j 0 = j 0 (r) and j n,m = j n,m (r). To determine p k we decompose it into a corresponding form, namely substitution of which into (24) leads to a sequence of differential equations for p Specifically, the homogeneous version of (35), namely which satisfies and which satisfies where M ± = m ± √ m 2 + 9, and so the solution of the inhomogeneous equation (35) is where is the Wronskian, which can be evaluated to give Note that we have chosen the forms of the homogeneous solutions and imposed the boundary conditions by selecting the lower limits of the integrals in equation (41) so as to ensure regularity at the origin (first integral) and at the contact line (second integral). In general, the integrals in equation (41) must be evaluated numerically to obtain p (n,m) k (r). Note, however, that, as described in Appendix A, even greater analytical progress is possible for the corresponding problem in the limit of large Bond number.

IV. A PAIR OF IDENTICAL DROPLETS
In this Section we apply the general methodology developed in Secs. II and III to determine the densities of the deposits from a pair of identical droplets, a situation for which the predictions of the present asymptotic theory for the local evaporative flux J k and the inte-  Consider the evaporation of a pair of identical droplets, which we may, without loss of generality, take to be of unit radius a 1 = a 2 = 1, with their centres located at (±b/2, 0), i.e.
with their centres a distance r 1,2 = r 2,1 = b (> 2) apart. As Wray et al. [48] showed, the two droplets have the same integral evaporative flux given by By symmetry, it is sufficient to consider only the left-hand droplet with its centre located at (−b/2, 0), corresponding to k = 1. The local evaporative flux from the surface of the droplet is given by which may be expanded as where and is the flux from the same droplet in isolation. The expansion (46) with (47)  In particular, Fig. 3(b) shows that the shielding effect described in Sec. IV A and shown in Fig. 2 leads to a spatially non-uniform deposit with the smallest density where the shielding effect is strongest (i.e. at φ = 0) and the largest density where it is weakest (i.e. at φ = π). Note that, by conservation of mass, the total mass of the deposit is the same for all of the values of b used in Fig. 3(b).

C. Comparison with experimental results
In general, comparing theoretical predictions for the density of a deposit with experimentally obtained images of deposition patterns is a challenging task. In particular, because the depth of the deposit cannot usually be readily determined from images taken from above the droplet, obtaining a quantitative measure of the amount of deposit from experimental images may often not be possible due to saturation and non-linearity of the data.  In order to compare the results of Pradhan and Panigrahi [39] with the present theoretical predictions, it is first necessary to convert the information contained within the image into a form proportional to the density distribution D. The raw data extracted from the image corresponds to the reflectance of the deposit; this must be converted into its absorbance, which can then be related to its concentration via the Beer-Lambert law [50], where A is the absorbance, P is the radiant power (reflectance) or intensity of the light as measured at each pixel in the image, P 0 is the initial radiant power of the light before absorbance, ǫ is the molar absorptivity of the deposit (which may reasonably be assumed to be constant), l is the path length of the light in the deposit, and C is the concentration of the absorbing species. The values for P are known, as this is the light gathered by each pixel in the sensor of the camera to generate the image; however, the value of P 0 is unknown.
(Ideally P 0 would have been determined through the collection of a reference image in which the light reflected from a calibrated sample, such as a 99% reflectance standard, was captured.) Therefore, it is necessary to make an estimate of P 0 , and to do this we used the maximum possible brightness value of 255 as the reference value for all of the pixels (consistent with the brightness of the image outside the footprint of the droplets). In order to gather only density data from the deposit near the contact line of each droplet, and not the density of any residual deposit left within it, the data was taken from an annular region around the edge of the footprint of each droplet. Slightly unfortunately, as Fig. 4 shows, in the published image the deposit from the left-hand droplet is overlaid by the head of an arrow that the authors added to indicate the region of weakest deposition. In order to reduce artifacts associated with this arrowhead, it was removed from the data by interpolating from the neighbouring pixels. The procedure for the extraction of the data for each droplet was implemented in Python [51], and is detailed below: 1. Convert the image to binary, and determine its centre of mass.
2. Trace out the outer perimeter of the deposit.
3. Define a second inner perimeter at 80% of the radius of the outer perimeter measured relative to the centre of mass.

Divide the annular region between the inner and outer perimeters into N sectors
subtending equal angles at the centre of mass. In practice, N = 300 sectors were used.
5. Use the Beer-Lambert law (51) to calculate the absorbance of each pixel, which is proportional to the mass of residue per unit area.
6. Integrate the mass per unit area numerically over each sector, and divide by the angle subtended, to determine the corresponding proportional mass per unit length of the contact line (i.e. the density to be compared with D). The constant of proportionality for comparison with D is determined by imposing the condition that the integral of the density is π/4, in line with the nondimensionalisation in Sec. II C. Figure 5 shows a comparison between the present theoretical predictions and the experimental results extracted from Fig. 2(b) of Pradhan and Panigrahi [39] using the procedure described above. Specifically, Fig. 5(a) compares the density of the deposit predicted by (50) (shown with the black solid curve) with the experimental results for the left-hand droplet (shown with the blue dashed curve) and the right-hand droplet (shown with the red dotted curve). Given the approximations made in both the mathematical model and in the processing of the experimental image, the agreement between theory and experiment shown in Fig. 5(a) is remarkably good, especially when it is noted that no fitting parameters have been used. The only scaling used is to ensure that the total mass (not given by Pradhan and Panigrahi [39]) is the same for theory and experiment.
While the agreement between theory and experiment shown in Fig. 5(a) is already good, with an integral absolute relative error of around 7% for the left-land droplet and 9% for the right-hand droplet, the experimental results are inevitably rather noisy. Figure 5(b) compares the density of the deposit predicted by (50) (shown again with the black solid curve) with the experimental results averaged across both droplets, as well as about φ = 0 (shown with the black dashed curve). In particular, Fig. 5(b) shows that, as expected, averaging reduces the noise in the experimental results and, rather pleasingly, leads to even better agreement between theory and experiment than that shown in Fig. 5(a), with an integral absolute relative error of around 3%. Figure 5(b) also includes the corresponding phenomenological estimate of the radially-integrated fluid flux (now interpreted as the density of the deposit) derived by Wray et al. [48] by using the approach of Sáenz et al. [31], denoted here by R 1 (φ, b), given by which can be evaluated to yield where k = 1/b (< 1) and F is given by (44). In particular, Fig. 5(b) shows that the estimate (53) is reasonably accurate, with an integral absolute relative error of around 9%, and captures the experimental results qualitatively but not quite quantitatively.

V. A TRIPLET OF IDENTICAL DROPLETS
In this Section we use the same approach as that described in Sec. IV to determine the densities of the deposits from a triplet of identical droplets of unit radius with their centres located at (−b/2, 0) and (b/2)(cos(π/3), ± sin(π/3)), i.e. with their centres a distance r 1,2 = r 2,3 = r 3,1 = b sin(π/3) (> 2) apart. Note that, in the terminology used by Wray et al.  Fig. 7(a)), and the density of the deposit develops local maxima and mimima (as shown by the curve for b sin(π/3) = 2.1 in Fig. 7(b)). Note that, as in Fig. 3(b), by conservation of mass, the total mass of the deposit is the same for all of the values of b used in Fig. 7(b).

VI. CONCLUSIONS
In the present work we obtained theoretical predictions for the spatially non-uniform densities of the contact-line deposits left on the substrate after the competitive diffusionlimited evaporation of multiple thin axisymmetric sessile droplets in proximity to each other.
In particular, we gave predictions for the deposits from a pair of identical droplets, which showed that the deposit is reduced the most where the droplets are closest together, and demonstrated excellent quantitative agreement with experimental results of Pradhan and Panigrahi [39]. We also gave corresponding predictions for a triplet of identical droplets arranged in an equilateral triangle, which showed that the effect of shielding on the deposit is more subtle in this case.
We note that, while the present analysis is formally restricted to thin droplets, the fact that much of the deposition from a non-thin droplet occurs towards the end of its lifetime when the contact angle is small and the velocity within the droplet is large (sometime referred to as "the rush hour", see, for example, [16,18]) means that the results of the present analysis are also expected to provide useful predictions for the contact-line deposits Finally, we note that the approach described in the present work is rather general and can, in principle, be applied to any arrangement of any number of thin droplets with pinned or unpinned contact lines.
progress is possible.
In the limit of large Bond number, the free surface of the droplet is flat, i.e. h k = h k (t), except in a narrow region near the contact line which we may neglect [52], and so the volume of the (nearly cylindrical) droplet is now given by V k = πa 2 k h k . The local fluid fluxes are again given by (15), while the kinematic condition (21) simplifies to and hence (23) becomes and so the partial differential equation for p k given by (24) becomes In principle, the same approach as that used in the main body of the present work can be used to solve the corresponding problem for large droplets. However, for brevity, in this Appendix we simply show how to obtain explicit asymptotic expressions for the pressure, and hence for the fluid fluxes and the density of the deposit, for a well-separated pair of identical droplets.
Adopting the same notation as in Sec. IV, and, without loss of generality, taking a 1 = a 2 = 1 and h 1 (0) = h 2 (0) = 1, the evaporative flux from the left-hand droplet is again given by equation (45), which can be expanded as in the limit of well-separated droplets, b → ∞. Hence, the equation for the pressure in the left-hand droplet, obtained by setting k = 1 in (A.3), can be expanded as (A.8) We now seek to determine the density of the deposit at a point on the contact line with polar angle φ = φ C , i.e. D(φ C ), where D is given by (50) with M given by (49). Note that M(2π) = π for this cylindrical droplet. In order to do this we must first locate the source from which all of the streamlines emanate. By symmetry this must lie on the line of symmetry, φ = 0, and, by expanding Q The deposit at φ = φ C is then given by Higher-order corrections to D(φ C ) may also be obtained, but rapidly become more complicated. Table I shows a comparison between the coefficients of the terms b −2 cos φ C and b −3 cos(2φ C ) obtained from the analytical prediction for D given by (A.14) and the values of D obtained from (49) and (50) Table I are due to the omitted higher-order corrections. In particular, analysis of these corrections indicates that the expression for b −2 cos φ C has relative error O (b −2 ), while the expression for b −3 cos(2φ C ) has relative error O (b −1 ), explaining the superior performance of the former.