Kinetic effects regularize the mass-flux singularity at the contact line of a thin evaporating drop

We consider the transport of vapour caused by the evaporation of a thin, axisymmetric, partially wetting drop into an inert gas. We take kinetic effects into account through a linear constitutive law that states that the mass flux through the drop surface is proportional to the difference between the vapour concentration in equilibrium and that at the interface. Provided that the vapour concentration is finite, our model leads to a finite mass flux in contrast to the contact-line singularity in the mass flux that is observed in more standard models that neglect kinetic effects. We perform a local analysis near the contact line to investigate the way in which kinetic effects regularize the mass-flux singularity at the contact line. An explicit expression is derived for the mass flux through the free surface of the drop. A matched-asymptotic analysis is used to further investigate the regularization of the mass-flux singularity in the physically relevant regime in which the kinetic timescale is much smaller than the diffusive one. We find that the effect of kinetics is limited to an inner region near the contact line, in which kinetic effects enter at leading order and regularize the mass-flux singularity. The inner problem is solved explicitly using the Wiener–Hopf method and a uniformly valid composite expansion is derived for the mass flux in this asymptotic limit.


Introduction
The evaporation of a liquid drop on a solid substrate has many important biomedical, geophysical, and industrial applications. Such applications include DNA mapping and gene-expression analysis, the water cycle, and the manufacture of semiconductor and micro-fluidic devices (see, for example, [1][2][3][4][5][6][7] and references therein). Modelling mass transfer from a partially wetting liquid drop is complicated because one must consider the transport of mass, momentum, and energy within and between three phases: the solid substrate, the liquid, and the surrounding atmosphere (assumed here to be a mixture of the liquid vapour and an inert gas). A key ingredient of any such model is an expression for the mass flux across the liquid-gas interface.
A commonly used model of a drop evaporating into an inert gas is the 'lens' model [2,5,[8][9][10][11][12]. The lens model is based on the assumptions that the drop is axisymmetric, the vapour concentration field is stationary, and the vapour immediately above the liquid-gas interface is at thermodynamic equilibrium, with the equilibrium vapour concentration being constant. These assumptions imply that evaporation is limited by the diffusion of vapour away from the interface. Notably, however, the lens model is thought not to apply to water [2,11].
The 'lens' model is so-called because the mixed-boundary-value problem for the vapour concentration is mathematically equivalent to that of finding the electric potential around a lens-shaped conductor [10,13]. Furthermore, if the drop is thin, this problem reduces to one equivalent to that of finding the electric potential around a disc charged to a uniform potential. The analytical solution of this electrostatic problem [14], translated to the evaporation problem, shows that the mass flux E * per unit area per unit time has the form where R is the radius of the circular contact set and r * is the distance from the axis of symmetry of the thin drop. The expression (1) for the mass flux has an inverse-square-root singularity at the contact line. Since this singularity is integrable, the total mass flux out of the drop is not singular, and physically reasonable predictions for the evolution of the drop volume are obtained even without regularization of the mass-flux singularity [10,12]. However, the need to supply a diverging mass flux means that there is a singularity in the depth-averaged radial velocity of the liquid flow within the drop [10,12]. Such a divergent velocity is clearly unphysical. In reality the mass flux at the contact line must be finite. Relaxing the assumption that the vapour concentration is stationary affects only the coefficient of the singularity. Instead, the assumption that the vapour immediately above the liquid-gas interface is at equilibrium must be invalid in the vicinity of the contact line. If the gas phase surrounding the drop instead consists of its vapour only (and no inert gas), an alternative boundary condition to apply on the liquid-gas interface is the Hertz-Knudsen relation, derived from the kinetic theory of gases [15]. The Hertz-Knudsen relation states that the mass flux across the drop surface per unit area per unit time is proportional to the difference between the equilibrium vapour density and the density of the vapour immediately above the drop. Formulated in terms of the vapour concentration (rather than the vapour density), on the free surface of the drop, we have where M is the molar mass of the liquid vapour, v k is a typical kinetic velocity (which we define later in the paper), c * e is the equilibrium vapour concentration, and c * is the vapour concentration at the interface. It is immediately apparent from the expression (2) that, provided the vapour concentration c * is finite, the mass flux is non-singular. The Hertz-Knudsen relation or the modified versions formulated in terms of vapour pressure, density, or temperature, have previously been used to model the evaporation of thin films [16], vapour bubbles in microchannels [17], and droplet evaporation on a precursor film [18]. While the assumptions required to derive the Hertz-Knudsen relation are not strictly satisfied when an inert gas is present, there is some experimental evidence that the Hertz-Knudsen relation is valid in such situations [19]. A possible explanation for this is that immediately above the drop, the gas phase is almost entirely vapour. It may therefore be reasonable to use the Hertz-Knudsen relation to model evaporation into an inert gas [20,21]. To close a model based upon the Hertz-Knudsen relation (2), it is necessary to prescribe a constitutive law for the equilibrium vapour concentration c * e (of course, such a constitutive law is also necessary if one makes the equilibrium assumption that c * = c * e on the liquid-gas interface). The simplest choice of constitutive law is to assume that the equilibrium vapour concentration is constant (as in the lens model). For a constant equilibrium vapour concentration, a kinetics-based model has the major advantage that, to leading order in the thin-film limit, the vapour transport problem depends on the liquid flow solely through the geometry of the contact set (and not through the drop thickness). This means that the vapour transport problem may be solved independently of the liquid problem. In this study, we shall exploit the simplicity of a kinetics-based model with a constant equilibrium vapour concentration to perform a mathematical analysis of the model and investigate the way in which kinetic effects regularize the mass-flux singularity.
Another possible constitutive law for the equilibrium vapour concentration is Kelvin's equation; this takes into account the variation in vapour pressure due to the curvature of the liquid-gas interface [22]. This approach has been used to model the evaporation of liquid drops in the presence of an ultra-thin precursor film that wets the substrate ahead of the drop [8,23]. In the bulk of the drop (away from the contact line), the dominant term in a linearized version of Kelvin's equation is independent of the drop thickness. As a result, in an outer region away from the contact line, a constant vapour concentration is prescribed on the liquid-gas interface and the mass flux appears to have a singularity at the contact line [23]. This singularity is in fact regularized in an inner region in the vicinity of the contact line, in which the other terms in Kelvin's equation become important [24]. In problems with a moving contact line, this evaporation model has the significant advantage that it also regularizes the stress singularity at the contact line [25,26]. Another advantage is the compatibility of the model with a precursor film; there is experimental evidence that such films exist in at least some parameter regimes [27,28]. We shall neglect the Kelvin effect in this paper, and establish a posteriori the regimes in which it is appropriate to do so (see Appendix 6).
In this paper, we adopt a linear, kinetics-based constitutive law for the mass flux across the liquid-gas interface, inspired by the Hertz-Knudsen relation (2); we assume that the equilibrium vapour concentration is constant. We will have two main goals. The first is to investigate the way in which kinetic effects regularize the mass-flux singularity at the contact line. The second is to derive an explicit expression for the evaporation rate. In Sect. 2, we formulate and non-dimensionalize the mixed-boundary-value problem for the vapour concentration. In Sect. 3, we perform a local analysis of both the lens evaporation model and the kinetics-based model to investigate the regularization of the mass-flux singularity at the contact line. In Sect. 4, we solve the mixed-boundary-value problem formulated in Sect. 2 to obtain an explicit expression for the evaporation rate. In Sect. 5, we perform an asymptotic analysis in the physically relevant limit in which the timescale of vapour diffusion is much longer than the timescale of kinetic effects to gain further insight into how kinetic effects regularize the mass-flux singularity. We find that there is an outer region away from the contact line where the equilibrium assumption (which leads to the mass-flux singularity) is recovered from our constitutive law and an inner region near the contact line where kinetic effects regularize the mass-flux singularity. The inner problem is solved explicitly using the Wiener-Hopf method, allowing us to derive a uniformly valid composite expansion for the mass flux in this asymptotic limit. In Sect. 6, we summarize our results and outline some possible directions for future work.

Formulation
We consider a three-dimensional, axisymmetric drop on a rigid, flat, impermeable substrate. We introduce cylindrical polar coordinates (r * , z * ) measuring the radial distance from the axis of symmetry of the drop and the normal distance from the substrate, respectively (here and hereafter, starred variables denote dimensional quantities). The contact set of the drop is 0 ≤ r * < R, so that (r * , z * ) = (R, 0) is the location of the contact line (at which the drop thickness vanishes). A mixture of liquid vapour and an inert gas occupies the region above the drop and substrate. A definition sketch is shown in Fig. 1. We assume that the drop is thin: the slope everywhere is comparable to the microscopic contact angle, Φ 1. Thus, the vertical extent of the drop is much smaller than the radius of the Cylindrical polar coordinates (r * , z * ) measure the radial distance from the axis of symmetry of the drop and the normal distance from the substrate, respectively. The location of the contact line is (r * , z * ) = (R, 0). circular contact set of the drop; since the latter is the relevant lengthscale for the transport of liquid vapour, the gas phase occupies the region z * > 0 to leading order in the limit of a thin drop. We assume that the dynamics of the vapour may be reduced to a diffusion equation for the vapour concentration c * , with constant diffusion coefficient D. We further assume that the timescale of vapour diffusion is much shorter than the timescale of the liquid flow (a common assumption in the literature [9,10,12]). Thus, transport of the vapour is governed to leading order in the thin-film limit by Laplace's equation, with We assume that the vapour concentration in the far field takes a constant value c ∞ , so that The inert gas is assumed to be insoluble in the liquid, so that the mass flux E * across the interface per unit area per unit time is entirely accounted for by the mass flux of liquid vapour. Since the substrate is impermeable, we have a condition of no flux of vapour through the substrate. After linearizing the boundary condition on the surface of the drop onto z * = 0, we obtain, to leading order in the thin-film limit, the boundary conditions where M is the molar mass of the liquid vapour. We assume that the mass flux out of the drop is governed by a linear constitutive law, given by where the equilibrium vapour concentration c e is a constant. The constitutive law (7) is inspired by the Hertz-Knudsen relation [15]. As discussed in Sect. 1, the Hertz-Knudsen relation is strictly only valid when the gas phase consists of pure vapour. However, there is experimental evidence that it may be valid for a vapour-inert gas mixture [19], and it has previously been used to model such situations [20,21]. The constant v k is a typical kinetic velocity, given by where R u is the universal gas constant and T in is the interfacial temperature. The (dimensionless) evaporation coefficient σ e is the fraction of the maximum possible evaporating flow rate that actually occurs [15]. One disadvantage of the constitutive law (7) is that the evaporation coefficient σ e is difficult to estimate; although a value of unity has been reported for many standard liquids, smaller values (anywhere between about 10 −4 and 1) have been reported in other cases. A quantity of interest is the surface-integrated flux out of the drop Q * , given by The quantity Q * is needed to determine the evolution of the volume of the drop and thus the extinction time (at which the drop volume vanishes), even in models that do not consider the detailed hydrodynamics of motion [29,30]. We see that if the contact line is pinned (so that the contact-set radius R is constant) the model (3)-(7) is independent of time-i.e. the problem is steady. If instead the contact line is allowed to move (so that R depends on time), then the problem is quasi-steady; the time dependence would become important if the expression that we ultimately derive for the mass flux were to be used as an input for a model for the evolution of the liquid drop. We shall use the contact-set radius R as a typical lengthscale on which to non-dimensionalize, suppressing the dependence of R on time in the case that the contact line is allowed to move. Thus, the expression that we shall ultimately derive for the evaporation rate will be valid for drops with either pinned or moving contact lines.
We non-dimensionalize (3)-(7) by scaling r * = Rr, z * = Rz, c * = c ∞ + (c e − c ∞ )c, and E * = DM(c e − c ∞ )E/R. We obtain thereby the following mixed-boundary-value problem for the dimensionless vapour concentration c(r, z): where non-dimensionalization has introduced a dimensionless parameter, namely the kinetic Péclet number, The kinetic Péclet number is the ratio of the timescales of diffusive and kinetic effects (over the radius of the circular contact set of the drop: R 2 /D and R/v k , respectively) and is the only parameter remaining in the problem following non-dimensionalization. We note the physical significance of two extreme cases: Pe k = 0 corresponds to the case of no mass transfer, while Pe k = ∞ corresponds to the case in which the vapour immediately above the free surface is at thermodynamic equilibrium, so that c = 1 on z = 0, 0 ≤ r < 1. Since this is the limit used in the lens model, we expect to obtain a diverging mass flux at the contact line as Pe k → ∞ (as will be discussed in Sect. 3.1). In Table 1, we give typical values of the relevant physical parameters for various liquids and various drop radii. We see that the kinetic Péclet number may take a wide range of values, but that it is at least moderately large for all but very small drops.
The key quantity of interest, the dimensionless evaporation rate E(r ), is given by Pe k ,R = 10 µm (-) 220 84 46 The equilibrium vapour concentration c e is evaluated using the saturation vapour pressure. In calculating the typical kinetic velocity v k from (8), we assume that the evaporation coefficient σ e = 1 and that the interfacial temperature T in is constant at 25 • C. We assume that c ∞ = 0 for each of the liquids in the table. The kinetic Péclet number Pe k = Rv k /D is given for (thin) drops with contact-set radii R = 1 mm and R = 10 µm A related quantity of interest, and a useful proxy, is the evaporation rate at the contact line, E(1 − ); the liquid motion has a strong dependence upon the size of this quantity [33]. We note that with Pe k = ∞, where the total (dimensionless) flux out of the drop Q is given by We emphasize that the three quantities E(r ), E(1 − ), and Q are all functions of the kinetic Péclet number Pe k . They therefore depend on the contact-set radius R (but not, in the thin-film limit, on the drop thickness).

Local analysis near the contact line
In this section, we perform a local analysis near the contact line of both the lens model and the kinetics-based model (considering the former puts the latter into context). This will demonstrate explicitly that the lens model has a mass-flux singularity at the contact line, while the kinetics-based model does not. Comparing the local expansions for the two models should also give us some insight into the way in which the kinetics-based model regularizes the mass-flux singularity.

Lens model
For the lens model, the boundary condition (12) is replaced by As noted earlier, this may be viewed as a special case of (12) with Pe k = ∞. Recall that the lens model (10), (11), (13), and (17) is mathematically equivalent to the problem of finding the electric potential around a disc charged to a uniform potential [14]. Assuming continuity of c at r = 1, this electrostatic problem has an exact solution [34,35], given by We deduce from (18) that the evaporation rate is given by We note from (19) that the total flux, Q = 4, is finite. From the exact solution (18), we deduce that the local expansion of the solution near the contact line is given by as ρ → 0 + , 0 ≤ θ < π, where (ρ, θ ) are local polar coordinates defined by r = 1 + ρ cos θ , z = ρ sin θ . The corresponding evaporation rate near the contact line has the local expansion Thus, we see clearly that there is an inverse-square-root singularity in the evaporation rate at the contact line, r = 1.
In Appendix 1, we show how this singularity leads to a singularity in the depth-averaged radial velocity of the liquid drop, which is unphysical.

Kinetics-based model
We now return to the mixed-boundary-value problem (10)-(13) for finite Pe k . We assume that c is continuous at the contact line and takes the value c L (Pe k ) there, with c L (Pe k ) not equal to 0 or 1. Under these assumptions, a local analysis near the contact line implies that is a degree of freedom. We then use (15) to find that the local expansion for the evaporation rate E(r ) near the contact line is given by In particular, this implies that the evaporation rate at the contact line E(1 − ) is given by Thus, the evaporation rate at the contact line (and everywhere else) is finite. In Appendix 1, we show that the depth-averaged radial velocity of the liquid drop is also finite. We recall that the lens model is a special case of the kinetics-based model with Pe k = ∞. Thus, for the local expansions (20) and (22) to be in agreement, it must be the case that but with c L < 1 for finite Pe k . Hence, we will be interested in determining the degree of freedom c L (Pe k ) by solving the mixed-boundary-value problem (10)-(13).

Explicit expression for the evaporation rate
We shall now solve the mixed-boundary-value problem (10)- (13). An important aim of this calculation is to determine the degree of freedom c L (Pe k ), appearing in (22), which will put the results of Sect. 3 in context. We will also obtain an explicit expression for the evaporation rate; this expression would be a key ingredient in investigations of the evolution of the drop.

Solution of the mixed-boundary-value problem
We note that the mixed-boundary-value problem (10)- (13) is mathematically equivalent to that of finding the temperature around a partially thermally insulated disc whose exterior is completely insulated; this problem was solved by Gladwell et al. [36] using Hankel, Fourier cosine, and Abel transforms, as well as properties of Legendre polynomials. The solution is given by where J 0 (kr) is the Bessel function of first kind of order zero, and the function f (x) satisfies the Abel integral equation given by By writing f (x) = ∞ n=0 a n sin[(2n + 1) cos −1 (x)] and expanding (27) in Legendre polynomials [36], we obtain where the coefficients a n (Pe k ) satisfy a system of infinitely many linear algebraic equations, given by where b mn = 1 2 and δ 0n is the Kronecker delta.
Using (15) we deduce that, for 0 ≤ r < 1, the evaporation rate is given by We integrate by parts once with respect to x and then change the order of integration. The resulting integral with respect to k may be evaluated explicitly, yielding for 0 ≤ r < 1. From this expression it is not clear, without further analysis, how E behaves as the contact line is approached, i.e. as r → 1 − . In Appendix 2, we analyse (32) as r → 1 − to find that the evaporation rate at the contact line E(1 − ) is given by By comparing the expression (33) for the evaporation rate at the contact line to the earlier expression (24) for the same quantity in terms of the concentration c L (Pe k ) at the contact line, we deduce that In practical applications, we may be interested in the total flux out of the drop, Q, given by We switch the order of integration since the integral with respect to r can be evaluated analytically. We then evaluate the remaining integral via the substitution x = cos(θ ) to obtain which is finite (as is also the case for infinite kinetic Péclet number).

Computing the evaporation rate
We have now deduced expressions for the evaporation rate E(r ) for 0 ≤ r < 1, the concentration c L (Pe k ) at the contact line, the evaporation rate at the contact line E(1 − ), and the total flux out of the drop Q in terms of a set of coefficients a n (Pe k ) that satisfy a system of infinitely many linear algebraic equations (29). We shall now describe how to solve numerically this algebraic system and thus how to compute the evaporation rate in practice. Previous work has shown that the system is regular [37] (in the sense that a n+1 a n as n → ∞) and may therefore be solved by truncation. In Fig. 2a, we plot a n (Pe k ) as a function of n for several values of Pe k . We observe that a n = O(n −4 ) as n → ∞; this rapid decay confirms that truncating the system (at a suitably large value of n) is appropriate.
It remains to determine a suitable value of n at which to truncate the system (29). We define the truncation error T M (Pe k ) in the evaporation rate at the contact line by where the coefficients a n satisfy the system (29) truncated at n = M. We define M * (Pe k ) to be the smallest value of M for which T M (Pe k ) ≤ 10 −4 . We calculate M * for a range of values of Pe k to create a lookup table, and then the value of M * for general Pe k is determined by spline interpolation (rounding up to the nearest integer). We plot M * as a function of Pe k in Fig. 2b. Thus, to compute the coefficients a n (Pe k ) in practice, we first use a lookup table and spline interpolation to determine a suitable value n = M * (Pe k ) at which to truncate the system (29). The resulting finite linear algebraic system is then solved using Matlab's backslash command (since the system is symmetric positive definite, this uses Cholesky factorization). Once the coefficients a n (Pe k ) have been determined numerically, the evaporation rate E(r ) is approximated by (32) with the sum truncated at n = M * (Pe k ). The integral in (32) is evaluated numerically using the integral command in Matlab. We check convergence in the usual way by reducing the error tolerances. We plot a scaled evaporation rate Pe k −1/2 E(r ) as a function of r for several values of Pe k in Fig. 2c. We see that the evaporation rate is everywhere finite for the values of Pe k plotted (which we note from Table 1 covers physically realistic values). We note from Fig. 2c that for large values of Pe k there appears to be a boundary layer near to the contact line in which the evaporation rate is much larger. We also observe from Fig. 2c that there appears to be a large-Pe k asymptote for the evaporation rate at the contact line of the form for some constant α ≈ 0.798 (with this asymptote presented as the dashed line in Fig. 2c). We deduce from (34) that c L (Pe k ) < 1 for finite Pe k and that c L → 1 − as Pe k → ∞, in agreement with our local analysis. Together with the fact that Pe k is typically large in practice (see Table 1), this motivates us to undertake an asymptotic analysis of the limit Pe k → ∞. It is not obvious how to find the coefficients a n (Pe k ) as Pe k → ∞ in the algebraic system (29), nor is it obvious how to analyse the integral equation (27) as Pe k → ∞, so we instead proceed by analysing the mixed-boundary-value problem (10)-(13) rather than the exact solution (32).

Asymptotic analysis in the limit of large kinetic Péclet number
In this section, we perform a matched-asymptotic analysis of the limit Pe k → ∞ to gain further insight into the way in which kinetic effects regularize the mass-flux singularity at the contact line. This is a singular perturbation problem; the asymptotic structure consists of an outer region in which |1 − r |, z = O(1) as Pe k → ∞, and an inner region near the contact line in which there is a full balance of terms in the boundary condition (12) on the free surface of the drop. We see that this happens when z = O(Pe k −1 ) and that to keep a full balance of terms in

Outer region
We expand c ∼ c 0 as Pe k → ∞. We find that the leading-order vapour concentration c 0 (r, z) satisfies (10), (11), and (13), but the boundary condition (12) is replaced by The leading-order vapour concentration therefore satisfies the mixed-boundary-value problem considered in Sect.
3.1 and we deduce that as Pe k → ∞ with (1 − r ) = O(1), We see that this outer evaporation rate has an inverse-square-root singularity as r → 1 − ; we expect this singularity to be regularized in an inner region close to r = 1.

The leading-order-inner problem
In an inner region near the contact line, we set r = 1 + Pe k −1 X , z = Pe k −1 Y , and expand c(r, z) ∼ 1 − Pe k −1/2 C(X, Y ) as Pe k → ∞. To leading order, the vapour transport equation (10) and the mixed-boundary conditions (12) and (13) become Finally, matching with the leading-order-outer solution (18) gives the conditions as ρ → ∞, where (ρ, θ ) are now plane polar coordinates related to (X, Y ) by X = ρ cos θ , Y = ρ sin θ . A local analysis of (41) subject to (42) and (43), assuming C to be continuous and non-zero at the contact line, implies that The value of the leading-order-inner solution at the contact line, C O := C(0, 0), is a degree of freedom in this expansion and we note that it is related to the degree of freedom c L (Pe k ) in the local expansion (23) of the full problem by The expression (24) for the evaporation rate at the contact line in terms of c L (Pe k ) then tells us that We shall solve the mixed-boundary-value problem (41)-(44) using the Wiener-Hopf method. The methodology employed is analogous to that used by Thompson [38] to solve a similar problem (known as the 'dock problem') consisting of (41) and (43), but with a sign change to the right-hand side of (42) and with different far-field behaviours.

Regularized inner problem
We begin by defining the functions with corresponding one-sided Fourier transforms C ± (k) given by We shall assume (and verify a posteriori) that C(X, 0) is infinitely differentiable on (−∞, 0) and (0, ∞). Then, using the far-field behaviour (44) and the local expansion (45), the Abelian Theorem in Appendix 3 tells us that C + (k) is holomorphic in Im(k) > 0, with and C − (k) is holomorphic in Im(k) < 0, with Moreover, a standard asymptotic analysis implies that the behaviour of C ± (k) as k → 0 is dominated by the behaviour of C(X, 0) as X → ±∞, with where k 3/2 + and k 1/2 − are defined as follows: Here, k The Abelian Theorem tells us that there is no value of k for which both C + (k) and C − (k) exist, so we are unable to apply the Wiener-Hopf method to the problem as it stands. Instead, we consider in the usual way [39,40] the regularized problem for the function C ε (X, Y ), given by We shall subsequently take the limit ε → 0 + to recover the leading-order-inner solution C(X, Y ) = lim ε→0 + C ε (X, Y ). A local analysis of (56) subject to (57) and (58), assuming C ε to be continuous and non-zero at the contact line, implies that C ε (X, Y ) has the same local expansion (45) at the origin as C(X, Y ), but with C O replaced by C ε O := C ε (0, 0). A far-field analysis, admitting only exponentially decaying separable solutions, implies that we require as ρ → ∞, where in order to recover (44) in the limit ε → 0 + , it is necessary for the constant A ε to satisfy the condition lim ε→0 + We now define F ε (X ) = ∂C ε /∂ X (X, 0). Using (59), we deduce from the Abelian Theorem in Appendix 3 that F ε where the functions C ε ± (X ) and their Fourier transforms C ε ± (k) are defined analogously to (48) and (49). By applying analytic continuation, we deduce that C ε + (k) is holomorphic in Im(k) > −ε except for a simple pole at k = 0 and C ε − (k) is holomorphic in Im(k) < ε except for a simple pole at k = 0. The presence of a simple pole at the origin in both C ε + (k) and C ε − (k) is consistent with the constants a ± being non-zero in the far-field expansion which follows from (59). We shall therefore apply the Wiener-Hopf method to the functions F ε ± (k). This is equivalent to applying it to the functions kC ± (k) due to (61) and (62) is holomorphic in Im(k) < ε, so that these functions are both holomorphic in the overlap strip −ε < Im(k) < ε. Before proceeding with the Wiener-Hopf method in the next section, we note that the Abelian Theorem in Appendix 3, together with the identities (61) and (62) (extended to Im(k) > −ε and Im(k) < ε, respectively), gives the far-field behaviour

Wiener-Hopf method
We begin by defining branches of the square roots (k ± iε) 1/2 : Thus, the square root (k + iε) 1/2 has branch cut S − = {k ∈ C : Re(k) = 0, Im(k) ≤ −ε}, while (k − iε) 1/2 has branch cut S + = {k ∈ C : Re(k) = 0, Im(k) ≥ ε}. We then define which has positive real part everywhere on the cut plane C \ (S + ∪ S − ). Now we define the Fourier transform in X of C ε (X, Y ) by Taking a Fourier transform in X of (56), we find that We therefore expect C ε (k, Y ) to be holomorphic in the strip −ε < Im(k) < ε except for a simple pole at the origin.
The boundary conditions (57) and (58) imply that so that eliminating B(k) between (70) and (71), and using (62) and (62) gives the following Wiener-Hopf equation for the functions F ε ± (k): In order to apply the Wiener-Hopf method to (72), we must find a product factorization of the function 1 + (k 2 + ε 2 ) −1/2 , namely where P ε + (k) is holomorphic in some upper half-plane Im(k) > γ + , and P ε − (k) is holomorphic in some lower half-plane Im(k) < γ − , with −ε ≤ γ + < γ − ≤ ε. The details of this standard factorization are given in Appendix 4 and reveal that suitable P ε ± (k) may be found with P ε + (k) holomorphic in C \ S − and P ε − (k) holomorphic in C \ S + . Given the product factorization (73), we may rewrite the Wiener-Hopf equation (72) as Since both sides of (74) are equal in the overlap strip −ε < Im(k) < ε, we deduce from the identity theorem that the right-hand side is the analytic continuation of the left-hand side into the upper half-plane. In the usual way, this allows us to define an entire G(k), given by Using the large-k behaviour (64) and (65) of F ε ± (k) and the fact that, by construction, P ε ± (k) → 1 as k → ∞ (see Appendix 4), we deduce that the large-k behaviour of G(k) is given by G(k) ∼ C ε O as k → ∞. Then applying Liouville's theorem -that a bounded, entire function is constant -to G(k) tells us that G(k) ≡ C ε O and we deduce that Solving for C ε ± (k) using (61) and (62), and taking the limit ε → 0 + , we obtain where P ± (k) := lim ε→0 + P ε ± (k) and C O = lim ε→0 + C ε O . We use the behaviour (120) and (121) of P ± (k) near the origin that we derive in Appendix 4 to deduce the behaviour of C ± (k) near the origin which is given by with k 3/2 + and k 1/2 − as defined in (54) and (55). We compare (78) and (79) to the asymptotic results (52) and (53) to deduce that the degree of freedom C O is given by We note that (80) may also be derived by, for example, inverting F ε + (k) to find F ε + (X ) for X > 0 (cf. Sect. 5.2.4) and then using Laplace's method to deduce that We may then apply (59) and (60), together with the fact that P ε − (−iε) ∼ (2ε) 1/2 as ε → 0 + , to deduce that

Inversion to find the inner mass flux
To find the mass flux in the inner region, we see from (57) that it is sufficient to find C(X, 0) for X < 0. (The full solution C(X, Y ) of the leading-order-inner problem is given for completeness in Appendix 5). Since C(X, 0) = C − (X ) for X < 0, we will invert C ε − (k) to find C ε − (X ) and take the limit ε → 0 + . We have The inversion contour lies below the singularities of C ε − (k) (namely, the branch cut S + and the pole at k = 0), so that for X > 0 we may close in the lower half-plane, where Re(−ik X) < 0, and use Cauchy's Theorem to obtain For X < 0, we deform into the upper half-plane, where Re(−ik X) < 0, with a 'keyhole' incision around S + , We note that this encloses the pole at k = 0. We obtain thereby, for X < 0, We take the limit ε → 0 + and use the expression (119) for P + (it) derived in Appendix 4, as well as the fact that P ε + (0) ∼ ε 1/2 as ε → 0 + , to deduce that where the function I (t) is given by Thus, as Pe k → ∞ with X = Pe k (r − 1) = O(1), X < 0, the inner mass flux is given by

Conclusions from the matched-asymptotic analysis
The evaporation rate (40)  We recall that the degree of freedom c L (Pe k ) belonging to the finite-Pe k mixed-boundary-value problem (10)-(13) is related to the degree of freedom C O of the leading-order-inner problem (41)-(44) by the expression (46). Using the expression (80) for C O , obtained from our matched-asymptotic analysis, we find that This result is in agreement with the conclusion (25) (which we made after performing a local analysis of the lens and kinetics-based models) about the way in which kinetic effects regularize the mass-flux singularity. In particular, this tells us, via (47), that the evaporation rate at the contact line E(1 − ) is given by Thus our matched-asymptotic expansion is in agreement with the numerics for the exact solution; in our prediction (38) for the large-Pe k behaviour of E(1 − ), we have α = (2/π ) 1/2 ≈ 0.798, which is presented as the horizontal dashed line in Fig. 2c. From the expressions (40) and (88) for the evaporation rate in the outer and inner regions, respectively, we deduce that a leading-order additive composite expansion for the evaporation rate E(r ), uniformly valid for 0 ≤ r < 1 as Pe k → ∞, is given by

Validation of asymptotic results
We shall now validate our leading-order asymptotic predictions against the finite-Pe k solutions that we obtained in Sect. 4. We shall consider the predictions for the total flux out of the drop Q, the evaporation rate at the contact line E(1 − ), and the evaporation rate E(r ) as a function of r . In Fig. 3a, we take the finite-Pe k solution for the total flux Q and plot Pe k (4 − Q) as a function of log 10 (Pe k ). For large values of Pe k (between 10 2 and 10 4 ) we fit a linear relationship, which we plot on the same axes. We see that for the physically realistic values of Pe k (40 and higher; see Table 1), there is very good agreement between the fit and the data. This gives us confidence that the leading-order asymptotic prediction (89) is correct, but with where we find numerically that A ≈ 1.28 and B ≈ 2.85. We do not investigate further in this paper such higher-order terms.
We plot the finite-Pe k solution for the evaporation rate at the contact line E(1 − ), given by (33), as a function of Pe k in Fig. 3b. On the same axes we plot the leading-order asymptotic prediction (91). We see that there is good agreement between the solutions even for moderately large values of Pe k . We note that both the form of the asymptote (91) and its validity for moderately large kinetic Péclet numbers are consistent with the observations that we made about the finite-Pe k solution following Fig. 2c.
To evaluate numerically the leading-order composite expansion for the evaporation rate (92), we first write the function I (t) as so that the integrand in (94) is bounded at the endpoints of the integration range. We then make the substitution t = τ 2 in order to remove the integrable singularity in the integrand of the second term in (92); we obtain, as Pe k → ∞, The integrals in (95) are computed in Matlab with the same methods used in the evaluation of (32). We plot the composite evaporation rate (95) as a function of r for Pe k = 10 1 , 10 2 , 10 3 , 10 4 in Fig. 3c. On the same axes we plot the finite-Pe k solutions (32); we see good agreement between the two solutions even for only moderately large values of Pe k . In Fig. 3d, we plot the relative error in E(1/2) between the finite-Pe k solution (32) and the leadingorder asymptotic prediction (95). The sharp dip in Fig. 3d is because for Pe k = O(1), the asymptotic prediction is an overestimate, while for large Pe k , it is an underestimate (i.e. the correction changes sign). For physically realistic values of the kinetic Péclet number (see Table 1), the relative error in E(1/2) is below 2% and is a decreasing function of Pe k , illustrating very good agreement between the two solutions (32) and (95).

Discussion
Our first aim in this paper was to investigate how the mass-flux singularity at the contact line of a thin, evaporating drop is regularized by applying a linear constitutive law on the liquid-gas interface that takes kinetic effects into account. Our second aim was to derive an explicit expression for the evaporation rate. In Sect. 2, we formulated a model for the transport of liquid vapour within the gas phase, assuming that the vapour concentration is steady, there is no flux of vapour through the solid substrate, the mass flux through the liquid-gas interface is governed by a linear, kinetics-based constitutive law, and the diffusion coefficient and the equilibrium and far-field vapour concentrations are all constant. The model was non-dimensionalized, leaving us with a single dimensionless parameter, the kinetic Péclet number Pe k (the ratio of the timescales of diffusive and kinetic effects). We tabulated the values of the physical parameters for hexane, isopropanol, and HFE-7100 and saw that Pe k was typically large for all but the smallest drops.
In Sect. 3, we performed a local analysis in the vicinity of the contact line on the kinetics-based model and also on the more standard lens evaporation model (which leads to a mass-flux singularity at the contact line). This demonstrated that the vapour concentration at the contact line c L (Pe k ) in the kinetics-based model was key to how the mass-flux singularity is regularized, with c L → 1 − as Pe k → ∞, but with c L < 1 for finite Pe k . This motivated the need to solve the mixed-boundary-value problem formulated in Sect. 2 and determine the degree of freedom c L .
In Sect. 4, we solved the mixed-boundary-value problem and deduced an expression for the mass flux in terms of a set of coefficients that satisfy a system of infinitely many linear algebraic equations. Analysis of the expression for the mass flux confirmed the hypotheses made in Sect. 3 about the degree of freedom c L and how the mass-flux singularity is regularized by kinetic effects. Our numerical simulations suggested that there was a boundary layer close to the contact line in which the evaporation rate was of size O(Pe k 1/2 ) as Pe k → ∞. This motivated us to further analyse the physically relevant limit of large kinetic Péclet number.
In Sect. 5, we performed a matched-asymptotic analysis of our model in the physically relevant regime of large kinetic Péclet number. We found that the asymptotic structure of the problem consists of an outer region away from the contact line, in which the vapour immediately above the liquid-gas interface is at equilibrium to leading order (as is assumed in the lens model). However, there is also an inner region near the contact line, in which kinetic effects enter at leading order. The leading-order-outer problem is equivalent to the lens model, while the leading-order-inner problem was solved readily using the Wiener-Hopf method. We found that the assumption that the vapour immediately above the drop surface is at thermodynamic equilibrium is valid in the outer region, with the mass-flux singularity being regularized in the inner region. We deduced from our leading-order asymptotic solution that c L ∼ 1 − (2/π ) 1/2 Pe k −1/2 as Pe k → ∞, quantifying the way in which kinetic effects regularize the mass-flux singularity. We also constructed a leading-order additive composite expansion and validated this asymptotic prediction by comparison with the solution found in Sect. 4; we found good agreement for physically realistic values of the kinetic Péclet number. Thus, for such values of the kinetic Péclet number, either solution for the mass flux may be used as an input to a model for the evolution of a liquid drop. The most important direction for future work is to incorporate our expression for the mass flux into a model for the evolution of the liquid drop. This would allow us to obtain predictions for the evaporation time, the evolution of the drop volume (or, equivalently, the dynamic contact angle or drop thickness), and, in the case of a moving contact line, the evolution of the contact-set radius within this model. Previous theoretical work has obtained such predictions for the lens evaporation model (with a mass-flux singularity at the contact line) [11,12,41,42] and for other evaporation models [7,11,18,[43][44][45][46][47][48][49]. In particular, it would be informative to compare the predictions of this previous work to the corresponding predictions for the model considered here. This comparison would give us some indication of what net result the inclusion of kinetic effects has on the liquid motion beyond regularizing the mass-flux singularity.
For a pinned drop, the evolution of the drop volume is fully described by the global conservation of mass equation (105). We have seen that, in the physically relevant limit when kinetic effects are weak compared to diffusive effects, the leading-order total flux out of the drop per unit time is the same for the kinetics-based model and the lens model.
For a drop with a moving contact line, we expect an important factor in determining the effect of kinetics to be the relative widths of the inner region in which kinetic effects come into play and the region in which the force singularity at a moving contact line is regularized. If the kinetic region is smaller, presumably the only noticeable effect of kinetics is to regularize the mass-flux singularity, while the remainder of the drop dynamics is the same as for the lens model (which we have shown is the leading-order approximation to the kinetics-based model away from the contact line when kinetic effects are weak compared to diffusive effects). On the other hand, if the kinetic region is at least as large as the region in which the force singularity is regularized, we expect that kinetics will have a more significant effect on the drop dynamics. Analysis of the drop dynamics for the lens model [12] suggests that this effect may be through an effective microscopic contact angle (different to both the true microscopic contact angle and the effective one for the lens model) that appears in the contact-line law.
Our analysis assumed that the timescale of vapour diffusion was much shorter than the timescale of interest (set by the liquid evolution). However, there are some situations in which the timescale of diffusion is comparable to the shortest timescale on which mass loss is important [12]. In such cases, Laplace's equation must be replaced by the unsteady diffusion equation. The resulting problem for the vapour concentration may be solved analytically [50]; we expect the solution on the timescale of vapour diffusion to converge in the long-time limit to the solution of the steady problem. A more thorough study of vapour transport would therefore be an interesting direction for future work. This point is particularly relevant for water, for which it is thought that the effect of the atmosphere may be important [51][52][53].
We made the assumption that the equilibrium vapour concentration is constant. However, there are many experimentally relevant scenarios in which it is more reasonable to assume that the equilibrium vapour concentration varies with temperature [21,45,54] or with the curvature of the interface [20,24,26]. In these cases, the appropriate modification of the mass flux is not independent of the drop thickness. A more thorough investigation of these scenarios would be of interest. In Appendix 6, we use the analysis of this paper to determine the range of lengthscales over which it is appropriate to neglect the effect of variations in the equilibrium concentration due to curvature (i.e. the Kelvin effect) compared to kinetic effects.
We also assumed that the problem is axisymmetric. In the non-axisymmetric case, in the large-Pe k limit, we expect that the details of the inner region would be the same in each plane perpendicular to the contact line, provided that the contact line is smooth. It would be interesting to investigate this point further and compare the results to previous work on non-axisymmetric drops [48].
The analysis presented in this paper pertains to thin drops, with a small microscopic contact angle Φ, and is only valid to leading order in the thin-film limit. Since we linearized the boundary condition on the free surface of the drop onto the substrate, the corrections to our analysis are of size O(Φ). While the leading-order prediction is independent of the drop profile, the O(Φ)-corrections would depend on the shape of the drop. Dependence on the drop profile is an ingredient in different mass-transfer models, such as those utilizing the Kelvin effect [23,26].
The expression (2) for the mass flux suggests that the inclusion of kinetic effects also ensures a finite mass flux for thick drops (where the aspect ratio is of order unity). A local analysis of the lens model near the contact line for 0 < Φ < π [2,55], assuming c to be continuous at the contact line, implies that, as r → 1 − , Thus, there is a mass-flux singularity at the contact line for 0 < Φ < π/2. The expression (96) is consistent with the corresponding expression for a thin drop (21) in the limit Φ → 0. On the other hand, a local analysis of the kinetics-based model near the contact line for 0 < Φ < π reveals that, as r → 1 − , where c L (Φ, Pe k ) and β(Φ, Pe k ) are degrees of freedom (the '. . .' indicating that other terms may impinge between those given). Thus the mass flux at the contact line is finite for 0 < Φ < π. The expression (97) is consistent with the corresponding expression for a thin drop (23) in the limit Φ → 0, provided that It would be interesting to investigate more thoroughly how the mass-flux singularity for 0 < Φ < π/2, Φ = O(1) is regularized by kinetic effects.

Appendix 1: The liquid phase
In Sect. 1, we noted that the lens model leads to a singularity in the liquid flow at the contact line. In this appendix, we give some brief details about the typical mathematical model for the liquid phase, assuming the flow to be axisymmetric. We use this model to show explicitly that the lens model leads to a singularity in the liquid velocity, while no such singularity is present for the kinetics-based model.

Formulation
Conservation of mass implies that, in the thin-film limit, the dimensional drop thickness h * (r * , t * ) is governed by [10,12,56] where t * is time, u * (r * , t * ) is the depth-averaged radial velocity of the liquid flow, and ρ is the density of the liquid (assumed to be constant). We assume that there are no body forces, the surface-tension γ of the liquid-gas interface is constant, and the liquid slips on the substrate according to a Navier slip law [57,58]. Under these assumptions, an expression for u * is given by [12] u * = γ 3μ where μ is the viscosity of the liquid and is the slip length (both assumed to be constant). A typical radial lengthscale is given by the initial contact-set radius R 0 (R 0 = R if the contact line is pinned). A typical timescale τ of capillary action may be identified from a balance of the two terms on the left-hand side of the thin-film equation (99) [12]. The thin-film approximation required to derive (99) is valid when the microscopic contact angle Φ 1 and the reduced Reynolds number Φ 2 ρ R 2 0 /μτ 1. We non-dimensionalize by setting r * = R 0 r , t * = τ t, R = R 0 s, h * = Φ R 0 h, u * = R 0 u/τ , and E * = DM(c e − c ∞ ) E/R 0 (so that r = sr and E = E/s). We obtain thereby the dimensionless thin-film equation given by ∂h ∂t with Non-dimensionalization has introduced two dimensionless parameters: the ratio α of the timescales of capillary action and mass loss, and the slip coefficient λ that measures the ratio of the drop thickness to the slip length. These dimensionless parameters are given by Appropriate boundary conditions subject to which to solve the thin-film equation (101) are given by The two boundary conditions at r = 0 (104a, b) are symmetry conditions. The third boundary condition (104c) states that the drop thickness vanishes at the contact line. The fourth boundary condition (104d) states that the dimensionless (small) microscopic contact angle is 1. We note that a local analysis of the thin-film equation (101) and (102) subject to the contact-line boundary conditions (104c, d) implies that there is no flux of liquid through the contact line [12]. We deduce from the thin-film equation (101) and the no-flux boundary conditions that the expression representing global conservation of mass of the drop is given by where V is the (dimensionless) volume of the drop and Q = s Q is the total (dimensionless) mass flux out of the drop per unit time.

Local analysis
To put our analysis of the lens and kinetics-based models into context, let us first consider the case of no evaporation ( E = 0). A local analysis of the thin-film equation (101) subject to the boundary conditions (104c, d) reveals that, for a moving contact line, whereṡ = ds/dt. Let us next consider the lens evaporation model. We write the local expansion (21) for the evaporation rate near the contact line in terms of liquid variables: A local analysis of the thin-film equation (101) subject to the boundary conditions (104c, d) therefore reveals that u ∼ 2 1/2 α s 1/2 (s − r ) 1/2 as r → s − ; there is an inverse-square-root singularity in the depth-averaged radial velocity at the contact line. Let us now consider the kinetics-based evaporation model. We write the local expansion (23) for the evaporation rate near the contact line in terms of liquid variables: A local analysis of the thin-film Eq. (101) subject to the boundary conditions (104c, d) therefore reveals that Thus, for the kinetics-based model, there is no singularity in the depth-averaged radial velocity at the contact line. Using the expression for u (102) and the three local expansions (106), (108), and (110), we may make the following deductions. First, the stress singularity at the contact line in the kinetics-based model has the same strength as the one for moving contact lines in the absence of mass transfer (see [59] for further details about the form of the singularity in the latter case). Moreover, this singularity is present for both moving and pinned contact lines (with a different coefficient in each of the two cases). On the other hand, the lens model has a stress singularity at the contact line that is stronger than the classical one for a moving contact line in the absence of evaporation, and this singularity is present even when the contact line is pinned (see [12] for details of the resulting local expansions for h at the contact line).

Appendix 4: Wiener-Hopf product factorization
We now outline the details of the Wiener-Hopf product factorization (73). With our choice of branch for (k 2 +ε 2 ) 1/2 , defined by (66)-(68), the left-hand side of (73) is holomorphic and non-zero in the strip −ε < Im(k) < ε, and tends to unity as k → ∞ in this strip. Moreover, the image of the strip under the principal branch of the logarithm, which we denote by log, does not encircle the branch point at the origin. We therefore apply the general method of Carrier et al. [39] and take logarithms to change the problem from a product decomposition to an additive decomposition. We find that where ± = {ζ ∈ C : ζ = x + iγ ± , x ∈ R}, with the real numbers γ ± chosen such that −ε < γ + < γ − < ε. We note that, by construction, P ε ± (k) → 1 as k → ∞. Following [34], we deform + into the lower half-plane with a 'keyhole' incision around the branch cut S − and − into the upper half-plane with an incision around the branch cut S + . We deduce thereby that P ε ± (k) is holomorphic for k ∈ C \ S ∓ , with It follows from (118) that P ± (±it) = t ∓1/2 (1 + t 2 ) ±1/4 exp ∓ 1 π t 0 log(s) 1 + s 2 ds for t > 0.

Appendix 5: Inversion to find C(X, Y )
In this appendix, we shall invert the Fourier transform in X of ∂C ε /∂ X (X, Y ), which by (76) is given by Provided that the contact-set radius R is much larger than R 1 , the Kelvin effect may be neglected in the outer region. Provided that R is much smaller than R 2 , the Kelvin effect may be neglected in the inner region where γ is the surface tension of the liquid-gas interface and c s is a reference vapour concentration (we may use the values for c e given in Table 1 as values for c s ). In the thin-film limit, the curvature κ * of the liquid-gas interface is given by The expression (15) for the evaporation rate E is then replaced by In the outer region, we deduce from (130) that the Kelvin effect may be neglected at leading order provided that σ 1. This is true provided that the contact-set radius R is much larger than a critical value R 1 . We report the values of R 1 for hexane, isopropanol, and HFE-7100 in Table 2, and we see that the assumption R R 1 is essentially always satisfied in practice.
In (assuming the slope of the drop in the inner region to be of order-unity size), we deduce from (130) that the Kelvin effect may be neglected at leading order provided that σ Pe k −3/2 . This is true provided that the contact-set radius R is much smaller than a critical value R 2 . We report the values of R 2 for hexane, isopropanol, and HFE-7100 in Table 2, and we see that the assumption R R 2 is satisfied for all but very large drops. Thus, for drops with a contact-set radius R such that R 1 R R 2 (and subject to the caveat on the slope mentioned above), we reach the conclusion that it is reasonable to neglect the Kelvin effect compared to kinetic effects at leading order. We note that this conclusion agrees with that made by [20] for a closely related model.