{'fraction_non_alphanumeric': 0.04025889967637541, 'fraction_numerical': 0.01581445523193096, 'mean_word_length': 4.273751280009102, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evaporation and fluid dynamics of a sessile drop of capillary size

Theoretical description and numerical simulation of an evaporating sessile drop are developed. We jointly take into account the hydrodynamics of an evaporating sessile drop, effects of the thermal conduction in the drop and the diffusion of vapor in air. A shape of the rotationally symmetric drop is determined within the quasistationary approximation. Nonstationary effects in the diffusion of the vapor are also taken into account. Simulation results agree well with the data of evaporation rate measurements for the toluene drop. Marangoni forces associated with the temperature dependence of the surface tension, generate fluid convection in the sessile drop. Our results demonstrate several dynamical stages of the convection characterized by different number of vortices in the drop. During the early stage the street of vortices arises near a surface of the drop and induces a non-monotonic spatial distribution of the temperature over the drop surface. The initial number of near-surface vortices in the drop is controlled by the Marangoni cell size which is similar to that given by Pearson for flat fluid layers. This number quickly decreases with time, resulting in three bulk vortices in the intermediate stage. The vortices finally transform into the single convection vortex in the drop, existing during about 1/2 of the evaporation time.


I. INTRODUCTION
Evaporation of a liquid drop in an ambient gas is a well-known yet not completely solved problem of the classical physics (see, for example, [1,2]). Past decade was marked by significant advance in experiment and progress in understanding several key aspects of evaporation process, in particular, the vapor diffusion from the sessile drop surface and the hydrodynamic effects within the evaporating drops [3,4,5,6,7,8,9]. It was found, in particular, that the evaporating flux density is inhomogeneous along the surface and diverges on approach to the pinned contact line [3,4]. The resulting inhomogeneous mass flow modifies the temperature distribution over the drop surface and, hence, the Marangoni forces associated with the temperature dependent surface tension. The convection inside a droplet [8,9,10,11,12,13,14,15,16,17,18,19] appears quite different from the classical Marangoni convection in the systems with a simple flat geometry [20,21]. Thermal conductivity of the substrate can also influence the formation of flows within a liquid drop since it is the magnitude of the conductivity which determines the sign of the tangential component of the temperature gradient at the surface close to the contact line and, therefore, the direction of the convection [7]. The observation of the distinct stages of the evaporation process [22,23,24] has revealed that the longest and dominating regime of the evaporation process is the constant contact area mode, where the contact line is pinned. As further the contact line gets depinned, the different regime, the constant contact angle mode switches on. Finally, the drying mode follows, in which the height, the contact area and the contact angle rapidly decrease with time.
Understanding the details of the drop dynamics and evaporation is crucial to many applications involving this process, including preparing ultra-clean surfaces [25,26,27,28], protein crystallography [29,30], the studies of DNA stretching behavior and DNA mapping methods [31,32,33], developing methods for jet ink printing [34,35,36], and in many other fields (see, for example, [37]). The process of evaporation of the drop with the colloidal suspension in it is of interest for methods of fabrication of various structures on the substrate. One of the examples is the effect of evaporative contact line deposition, the so-called coffee-ring effect [3,4,6,38,39,40]. Other important example is the self-assembly of long-range-ordered nanocrystal superlattice monolayer [41,42,43].
While the fundamentals of quasistationary evaporation behavior are established, to some extent [8,9], the full dynamical description of the liquid drop evaporation is yet not available. In this work we undertake a step towards such a description offering a quantitative approach which enables us to account for all the relevant components of the evaporation process, namely, the fluid dynamics, the vapor diffusion, and the spatial temperature distribution in an evaporating sessile drop. The calculation are carried out according to the following scheme: (i) We derive the equation yielding the shape of the sessile drop and find how the evaporation rate is influenced by the deviations of the drop shape from the ideal spherical cap caused by the gravitational force. (ii) We solve the diffusion equation describing the vapor kinetics and calculate the local evaporation flux from the surface of the drop as well as the resulting evolution of the drop evaporation rate with time. At this stage we take into account the nonstationary effects in a vapor diffusion and find that the dynamical corrections to the evaporation rate do not vanish exponentially, but decay as ∝ 1/ √ t. (iii) We solve the thermal conduction equation and find the spatial temperature distribution in the drop. (iv) We solve the Navier-Stokes equations and derive the velocity field corresponding to the convection within the drop. (v) Finally we arrive at a full description of the time evolution of the temperature and the fluid convection in the drop and identify characteristic stages of the thermocapillary convection. To this end we take into account inertial terms in the Navier-Stokes equations, including the time derivatives.
To ensure self consistency of the derivation of the physical characteristics of the system we include the temperature variation over the drop surface (since the surface tension depends on temperature) into the boundary conditions for the fluid dynamics in the drop. Furthermore, we take into account that the velocity field influences the thermal conduction, since the velocities enter the thermal conduction equation. And, finally, we take into account that the local evaporation flux is related to the heat transfer and, hence, to the temperature gradient on the drop surface.
While the developed approach is quite general and thus applies to a wide variety of the evaporating problems, we focus here on the experimental situation corresponding to evaporation of the toluene sessile drop of the capillary size. The typical values of the characteristic parameters at T 0 = 295K are taken from [44] and are presented in the lower part of Table I. The initial values which we use in our calculations are given in the upper part of Table I and correspond to the evaporating toluene drop containing the colloidal solution of gold nanoparticles. These nanoparticles stick to the drop surface making islands there and form a monolayer nanocrystal on the substrate after toluene dries out [41,42,43]. We have identified three major distinct dynamic stages of the Marangoni convection during the evaporation of the toluene drop. During the initial stage, vortices appear near the surface of the drop, during the second stage they grow, coalesce and rapidly migrate into the drop bulk, and, finally, a single vortex survives governing all the fluid dynamics of the thermocapillary effect in the system. We derive the temperature profiles, evaporation rates, and the velocity distributions for each stage of the convection and uncover the important role of the convective heat transfer and the effect of inertial terms in the Navier-Stokes equations in the evaporation dynamics.
The paper is organized as follows. In Sec. II we present basic equations and main analytical results along with the quantitative description of main dynamical processes. In particular, Sec. II A deals with the shape of a sessile drop. Various aspects of the vapor diffusion and the drop evaporation rates are discussed in Sec. II B-II D. The characteristic Marangoni numbers for the drop are found in Sec. II E. Sec. II F makes use of a self-consistent approach for calculating the velocity field and the spatial temperature. In Sec. III we present our main results and discussion, including the experimental and theoretical data for the evaporation rates and contact angles, the dynamics of vortex structure and the temperature profile. Sec. IV contains discussion. The details of the developed numerical approach are described in Appendix A.

II. BASIC EQUATIONS AND FORMULAS
A. The sessile drop shape During the evaporation process a drop loses its mass, hence with time its volume decreases and the shape changes. In this section we present the quasistationary method for calculating the shape of a sessile drop based on the hydrostatic Laplace equation. The approach holds as long as viscous forces, which generally enter the boundary condition for the pressure, are small. The ratio of viscous to capillary forces is characterized by the dimensionless number C a = ηv/σ, where v is the characteristic value of the velocity. In our case C a ≈ 2 · 10 −4 ≪ 1.
The Laplace equation states that the pressure difference ∆p, taken at the different sides of the surface of the liquid in an arbitrary point equals σk, where k = 1/R 1 + 1/R 2 is the mean curvature of the surface and σ is the surface tension [45]. Taking into account the gravity, one finds that at the surface of the drop where the shape of the drop surface is defined by the relation z = f (x, y) and a = 2σ/(ρg) is the capillary constant. The curvature of the surface of the drop is, in its turn, also expressed in terms of the function f (x, y). Indeed, k = −Tr(G −1 Q), where G and Q are matrices of the first and the second quadratic forms of the surface [46]. Plugging in the surface equation z = f (x, y) we find Therefore, We consider a sessile drop with the axial symmetry, where z = f (r), r = x 2 + y 2 . Then Eq. (3) assumes the form where prime denotes the derivative with respect to r. It is convenient to parametrize and separate the variables as [47,48]: where s is the arc length on the drop surface taken from the apex. The curvature radii are R 1 and R 2 , where R 1 is the curvature radius in the meridian rz-plane, and R 2 is the second curvature radius. Both curvature radii can be expressed in terms of the angle φ between the normal vector to the drop surface and the symmetry axis. Namely, Also, ∂r/∂s = cos φ, ∂z/∂s = − sin φ and r(0) = z(0) = φ(0) = 0. Introducing the curvature radius R 0 at the apex, we can rewrite the Laplace equation as and therefore Hence, for the column vector y = (r(s), φ(s), z(s)) T , where T is the transposition operation, we have the set of first-order differential equations with Cauchy boundary conditions: f (s, y) = cos φ, where the last relationship ensures that φ(s)/r(s) → 1/R 0 for s → 0. The Cauchy problem (9)-(11) can be solved using the fourth order Runge-Kutta method or the Adams-Bashforth method [47]. Therefore, given the values of R 0 and s max , one finds the drop shape in the form of y(s) = (r(s), φ(s), z(s)) T . The quantities that are measured in the experiment are the mass of the drop and the radius r 0 of the substrate to which the drop is pinned. Given the values of R 0 and r 0 , one can easily find s max , and, therefore, m. Inverting the function m(R 0 , r) numerically, one can obtain the value of R 0 and a drop shape for a given mass m and contact line radius r 0 . We note that Eqs. (9)-(11) are exact for a sessile drop with axial symmetry. The evaporation rate for a given drop shape can be obtained as where J(r(s)) is the mass evaporated per second from unit area of the surface, i.e. the local evaporation flux.

B. Evaporation of the drop
Consider a sessile droplet resting on a flat substrate. The vapor concentration above the droplet is time dependent and inhomogeneous during the evaporation process. Since the diffusion of the vapor from a near-surface layer is slower than the evaporation [45], the vapor concentration at the surface of the drop is assumed to equal the saturation value u s . Far above the drop, the toluene vapor concentration is negligible. The dynamics of the vapor concentration in the surrounding atmosphere is described by the diffusion equation We carry out our calculations (see Sec. III and Appendix A) taking into account both time dependence of vapor concentration and the deviations of the sessile drop shape from a spherical cap. It is instructive to begin discussing the problem with the more simple case allowing for an analytical description. To this end we notice that if evaporation can be viewed as an adiabatic process in a sense that the vapor concentration adjusts fast enough to the change of the drop size (and shape) i.e. on the time scales much less than the droplet evaporation time t f , then the diffusion equation (14) can be replaced by the Laplace's equation ∆u = 0, and the evaporation kinetics can be considered as a steady state process. Indeed the time t p required for the Brownian particle to pass the characteristic length r 0 is t p = r 2 0 /D ≪ t f under the experimental conditions t p ≈ 0.2 s and t f ≈ 500 s. At the same time it interesting as well to consider particular dependencies of the local evaporation rates on time. Our simulation results based on Eq. (14) for the drop with the fixed surface can be fitted with the power-law time dependence which is almost exact with the accuracy within 1% for t > 0.5 s and |r 0 −r| > 0.01 cm (see Fig. 1). Here J(r, t) = |D∇u| is the local evaporation flux on the drop surface, and constant A = 0.966 is the only fitting parameter. For the case of the spherical cap the asymptotic value of the evaporation flux density J(r, ∞) can be related to J(ξ, θ) from Eq. (20). Eq. (15) and Fig. 1 confirm the above estimation that the vapor concentration becomes stationary under the condition t ≫ r 2 0 /D. At the same time, for t = t f the second term in (15) still exceeds 1% of the first term. The time dependence J(r, t) − J(r, ∞) ∝ 1/ √ t in Eq. (15) is known for an isotropic diffusion when the vapor concentration is kept saturated on the fixed surface of the sphere [2]. This is also valid for a diffusion from a flat plane [45]. Our results demonstrate that such dependence also takes place for an inhomogeneous diffusion from a fixed surface of the sessile drop.
We further derived the numerical results which take into account the time dependent drop profile (see Fig. 5). In particular, they show that the corrections for the evaporation rate due to the nonstationary effects may be up to 5% of the resulting value. Deegan et al. [4] reported an analytical solution for a stationary spatial distribution of the vapor concentration for a drop with the shape of a spherical cap. The problem was shown to be mathematically equivalent to that solved by Lebedev [49], who obtained the electrostatic potential of a charged conductor having the shape defined by two intersecting spheres. Modifying the equations derived in [4,49], we arrive at the following formula for the vapor concentration u(α, β) above the surface (25) of a spherical drop: Here θ is the drop contact angle, −π + θ ≤ β ≤ π − θ, u ∞ = 0, α and β are the toroidal coordinates: At The evaporation flux at the surface is and, therefore, After substituting ξ = r/r 0 and (18) in Eq. (20), one obtains J(r). Fig. 2 shows log J(r) for drops with different values of the contact angle θ. The total drop mass is and the total evaporation rate is While it were possible to simplify the expression for the drop evaporation rate adapting the sessile drop shape to that of the spherical cap of the same mass (see Sec. II D), and then integrate Eqs. (20), (22), we will employ a more general method and develop a description of the evaporation rate which takes into account altogether spatial and time variations of the vapor diffusion and the sessile drop shape.

C. Estimation for the evaporation time
In this section we derive estimates for the evaporation time of the drop, which in spite of their relative simplicity, offer, nevertheless, a fairly good description of the experimental results. Let us consider a spherical evaporating drop of radius r 0 = 0.2 cm, so that the vapor density is saturated at the drop surface and vanishes far away from the drop. The diffusion of the toluene vapor controls the evaporation process. The outward flow of the vapor through the concentric sphere around the drop of the radius R is J = −D · 4πR 2 du/dR and does not depend on R. Here u is the vapor density on the surface of the sphere of radius R, and D is the diffusion constant for the toluene vapor.
Since the vapor concentration far away from the drop is negligible, we obtain u = A/R, A = r 0 u s and J = 4πDr 0 u s . On the other hand, the flow of molecules from the fluid drop is J = −ρd(4πr 3 0 /3)/dt. Therefore, r 0 dr 0 /dt = −Du s /ρ, and the evaporation time of the drop can be written as Eq. (23) is a well-known result of the classical Maxwell's theory of the evaporation (see, for example, [2]). The diffusion constant for toluene vapor can be estimated as Here v = 8RT /(πµ) = 26035 cm/s is the average velocity for the thermal motion of the vapor molecules, λ is the mean free path which may be estimated as follows. We replace the toluene molecules by the impermeable balls of the diameter d ∼ 1 nm (we found the estimate for the diameter by adding together the lengths of the chemical bonds and taking into account the molecule geometry). Therefore, the mean free path is λ = 1/(π √ 2nd 2 ) = kT /(π √ 2P 0 d 2 ) ≈ 2.7 · 10 −5 cm. Here n is the density of the toluene molecules in the toluene vapor. If we substitute the obtained values into the Eq. (24), we get the estimate D ∼ 0.23 cm 2 /s. This qualitative estimate is about in 1.5 times larger than the result obtained from comparison of our numerical results and the experimental data in Sec III.
Substituting the value of D to the Eq. (23), we get the evaporation time of the spherical drop: t ≈ 594 s.

D. Approximating a sessile drop surface with spherical caps
Spherical caps offer a very simple model allowing for carrying out analytical calculations to the end and are thus widely used for the modeling of the droplet. The equation for the surface is Figure 3 shows surface of the sessile drop and three approximating spherical caps with the parameters listed in Table II D. The radius of the substrate is r 0 = 0.2 cm for all the drops. The first spherical surface has the contact angle of the sessile drop. The second spherical surface has the mass of the sessile drop. The third spherical drop has the height of the sessile drop.
As seen from the Fig. 3 and Table II D, the characteristics m, h and θ of all three approximating spherical caps are quite close to those of the sessile drop. This is not the case for local parameters such as the local curvature of the sessile drop. As follows from Eq. (1), for the local curvature to be approximately constant along the drop surface, the condition B o = ρghr 0 /(2σ sin θ) ≪ 1 has to be satisfied. Here B o is dimensionless number which is analogous to the Bond number.  Table II D).    Table II D. time, it is not surprising that the total evaporation rate, integrated over the surface of the drop, does not depend much on whether the drop is sessile or spherical. With the values from Table I, one may numerically integrate (22) and obtain the evaporation rates for the three spherical drops: dm 1 /dt = 20.52 µg/s, dm 2 /dt = 19.76 µg/s, dm 3 /dt = 19.47 µg/s. The difference is within a few per cent.

E. The Marangoni convection
A fluid convection within a drop caused by the temperature dependent surface tension which persists in the drop during the all stages of an evaporation process is called the Marangoni convection, and was first observed in its classical form by Bénard in a process of formation of the characteristic hexagonal convection patterns in flat fluid films. A theoretical description of the Marangoni convection was developed by Pearson [21]. The Marangoni number M a (see Eq. (26) below) characterizes a relative importance of surface tension forces caused by the temperature variation and viscous forces, leading to the Benard-Marangoni instability if M a exceeds a critical value M c . For a flat fluid film, M c ≈ 83 [21]. To estimate the Marangoni number for a drop we exploit a formal similarity between a drop and a flat layer. We extract the temperature difference between the substrate plane and the apex of the drop from our simulation results as ∆T ≈ 1 K. The drop height h, the thermal conductivity k, the kinematic viscosity ν, the specific heat c V , and the rate of the change of the surface tension with the temperature, −dσ/dT , are taken from the Table I. Then the value of the Marangoni number at the beginning of the drop evaporation process is According to experimental data, a turbulent regime arises for considerably larger values of the Marangoni number. For example, in the water drop the turbulence takes place for M a > 22000 [50]. As the Reynolds number for toluene drop, R e ≈ 62, is also comparatively small, the flow in our problem is of a laminar character. A competition between the buoyancy-induced flow and the thermocapillary convection is characterized by the dimensionless number B = ρgh 2 β/(7.14σ ′ ) [21]. In our case B ≈ 0.02 ≪ 1, and the buoyancy-driven convection is negligibly small.
For a flat liquid layer the height of the Marangoni cell coincides with the layer thickness h, whereas the transverse size of the cell is determined by the characteristic instability wavelength λ. If the Marangoni number well exceeds the critical value, then λ = 2πh/α and α = M a /8 [21].
Our simulations show that near the surface a "street" of vortices forms at the initial stage of the convection (see Sec. III). The size of the vortices (i.e. the distance over which they extend into the bulk of the drop), h d , constitutes a noticeable fraction of the height h of the drop. With the assumption that the instability wavelength λ and the Marangoni number can be described similarly to the case of flat liquid layer, one can take the value h d for the layer thickness and obtain λ = 2π √ hh d /α. Taking as an example the number of the near-surface vortices N = 9 and s max = 0.256 cm, we find λ = s max /N ≈ 0.028 cm and This although crude estimate turns out to be in a good agreement with our simulations for the Marangoni convection (see Figs. 6,7,8 and their discussions).

F. Hydrodynamic equations
The basic equations inside the drop are the Navier-Stokes equations, the continuity equation for the incompressible fluid, and the equation for thermal conduction Here ∆ = ∂ 2 /∂r 2 + ∂/r∂r + ∂ 2 /∂z 2 , ν = η/ρ is kinematic viscosity, κ = k/(ρc p ) is thermal diffusivity. Other terms in the thermal conduction equation ν/(2c p ) (∂v i /∂x k + ∂v k /∂x i ) 2 , where sum over i and k is assumed, are negligibly small.
Applying the "rot " operation to the both sides of equation (28), one excludes the pressure p and obtains Therefore, where the vorticity γ is introduced by We define the stream function ψ, such that One has ∂ 2 ψ/(∂r∂z) = v r + r(∂v r /∂r) = −r(∂v z /∂z), therefore v r /r + ∂v r /∂r + ∂v z /∂z = 0, i.e. the velocities obtained with Eqs.(34) will automatically satisfy the requirement (29). Applying the Laplace operator to the stream function, one obtains ∆ψ = rγ − 2v z , which can be transformed to a more convenient form∆ψ = rγ with the modified operator∆ that differs from the Laplace operator by the sign of the term ∂/(r∂r):∆ The method of the numerical solution is as follows: at each step, inside the drop 1. We solve Eq. (32) with the proper boundary conditions and find γ(r, z); 2. We solve Eq. (35) and find ψ(r, z). Then we find velocities with Eqs.(34); 3. We solve Eq. (30) and obtain the new boundary conditions for Eq. (32).
Details of the numerical approach are given in Appendix A. The proper boundary conditions for quantities γ and ψ, satisfying Eqs. (32) and (35) respectively, are derived in Appendix B. They take the form γ = 0 for r = 0; γ = ∂v r /∂z for z = 0; γ = dσ/(ηds) + 2v τ dφ/ds on the surface of the drop; ψ = 0 at all boundaries: at the surface of the drop, at the axis of symmetry of the drop (r = 0), and at the bottom of the drop (z = 0).
Here dσ/ds = −σ ′ ∂T /∂s is the derivative of surface tension along the surface of the drop, where the distribution of temperatures, which is found with Eq. (30) is taken into account, and is the experimental dependence of surface tension σ on temperature T . The boundary conditions for the Eq. (30) are ∂T /∂r = 0 for r = 0; T = T 0 for z = 0; on the surface of the drop. Here Q 0 (r) = LJ(r) is the rate of heat loss per unit area of the upper free surface, J(r) is the local evaporation flux, T 0 is the temperature of the substrate, and n is a normal vector to the surface of the drop. The relation Q 0 (r) = LJ(r) implies that the heat flow from ambient air towards the drop surface is negligible. This is the case if the temperature difference between the drop surface and the air far from the drop is less than LDu s /k [2]. This is well satisfied in the problem in question.

A. Evaporation rates
The evaporation rate of sessile drops is controlled mainly by the vapor diffusion in the surrounding atmosphere [4].
Here we determined the diffusion coefficient of toluene vapor in air by measuring the toluene evaporation rate. The main conditions and parameters of the experiment are as follows. The sessile drop of 10µl of toluene is lying on a substrate and evaporates to the atmosphere during ≈ 500 s. The contact line of the drop is pinned to the edge of the substrate. The weight of the drop was measured every 10 seconds. As substrates, we used silicon wafers with a 100 nm thick amorphous silicon nitride layer.
The results of the measurements of drop evaporation rates are presented in Fig.5. Open circles are experimental values for the pure toluene evaporation and open triangles are values for the gold nanoparticles colloid. The solid line is the simulation result obtained in terms of Eq. (13) within the full numerical scheme. As seen in Fig. 5, it agrees well with the experimental data. Comparing the experimental data and the simulation results permits us to find the diffusion coefficient of the vapor, D, the only parameter controlling evaporation, as D = 0.1449 cm 2 /s. It is worth noting that the rough estimate mentioned in Sec. II C is larger by the factor of about 1.5.
One can divide the evaporation process into three characteristic parts. Our simulation describes the main part, i.e. the quasi-stationary diffusion-limited evaporation, which lasts for about 400 seconds. While there is a good agreement of the simulation and experimental data during the diffusion-limited regime of evaporation, this is not the case for the drying stage, i.e. the final 50-100 seconds. At the drying stage, the experimental evaporation rate drops rapidly but not abruptly. The main physical reason for such a behavior is the depinning of the contact line. The interface area then shrinks over time. The two distinct parts in the drying stage of evaporation can be identified. The first part, that occurs in the t=400−510 time window can be described as dm/dt∝(t 0 −t) α that fits well the experimental data with t 0 =550(2), α=0.41 (2) for the evaporation of the pure toluene and t 0 =522(3), α=0.28 (8) for the colloid evaporation. The last part of the drying stage is a pronounced exponential decay of the evaporation rate dm/dt∝ exp(−t/t d ), with the decay time t d =20.5(3) for the colloid evaporation. The exponential behavior during this stage qualitatively follows from the fact that change of the drop surface area A controls here the evaporation rate evolution: dA/dt ∝ dm/dt ∝ A. Inset in Fig. 5 represents the experimental mass variation during the drying stage of the evaporation process. drop. For a toluene drop, this regime quickly arises and evolves up to t ≈ 0.3 s. This is quite a short time period as compared with the total evaporation time ≈ 550 s, but it admits an experimental study. The vortices grow, the number of vortices decreases, and eventually they evolve into the Marangoni convection cells in the bulk of the drop. This can be seen in Figs. 6,7, where the vortex structure contains six pairs of near-surface vortices and a corner vortex, and the temperature displays just six humps at the surface. As was shown in Sec. II E, this number of vortices in the drop is controlled by the Marangoni cell size which is similar to that given by Pearson for flat fluid layers. The existence of near-surface vortices and the associated humps in the profile of the surface temperature, become more pronounced with the decrease of the viscosity of the liquid. There are no near-surface vortices when the viscosity increases more than in four times as compared with the toluene drop.
As seen in Fig. 7, the extrema of the surface temperature correspond to the change of sign of the tangential component of velocities at the surface. The reason for this behavior is that the fluid flow moves from the higher to the lower temperature regions of the surface, because, according to (36), the surface tension decreases with increasing the temperature. The flows result in a redistribution of the temperature due to the convective heat transfer. The number of the surface temperature humps and the number of the near-surface vortices decrease during their evolution.
The initial conditions chosen when generating Fig. 6 were: the vanishing velocities, constant room temperature, and the vanishing vapor concentration. For a description of a real experiment they have to be slightly modified to include the weak stochastic distribution of the surface temperature and velocities. We have carried out such calculations using a small-scale stochastic initial conditions. This has changed the particular behavior of the fluid dynamics only on the initial stage of the process, where a large number of small-scale surface vortices arise. Then the surface vortex structure quickly evolves into exactly the same one as we obtained for the basic initial conditions. This result demonstrates the generic character of the near-surface vortex regime in Fig. 6 at the early stage of the formation of the Marangoni convection.
An initial velocity field within the sessile drop can also be strongly disturbed right after the drop has fallen down on the substrate, or for some other reason. We model such a situation by choosing random initial conditions for the bulk velocity field, which are in agreement with the continuity equation (29). We find that strong disturbances of the bulk velocities can noticeably modify the initial stage of the drop dynamics, but they do not modify its main stage. For the initial random velocities of the order of 5 cm/s (which well exceeds the typical velocity in the vortex 1 cm/s), the difference between the dynamics of the disturbed and resting drops disappears after t ≈ 0.5 s. This verifies the stability of the large-scale drop dynamics with respect to disturbances of the initial temperature and velocity fields.  During the enlargement of the near-surface vortices, their number decreases, and the convection involves the bulk of the drop. As a result, for t ≈ 0.45 s, three bulk vortices control the velocity and temperature fields in the drop, as seen in left panels in Figs. 8,9. During the coexistence of three vortices, the corner vortex starts growing at the expense of the other two vortices, and eventually at t ≈ 2.0 s it occupies the whole drop volume. A spatial dependence of the temperature along the drop surface is nonmonotonic, if the drop contains more than one vortex (see Fig. 9). Right panels in Figs. 8,9 demonstrate how in the single-vortex regime effects of Marangoni forces drive liquid along the surface to the apex, where the fluid penetrates along the symmetry axis in the depth of the drop.  Fig. 10, the single-cell pulsating state breaks into two-center (and later three-center) pulsating structure. Eventually at t ≈ 300 s a quasistationary state with three vortices arises.
The numerical calculations of the fluid dynamics were tested with several different mesh sizes. The respective results are qualitatively identical and show reliable convergence of the quantitative characteristics. For example, the single-vortex regime was found to arise at 3.48 s, 2.48 s, 2.2 s, 2.06 s, 2.01 s for 100×100, 150×150, 200×200, 250×250, 300×300 mesh elements covering a half of the drop cross-section.  If the thermal conductivity of a substrate is large compared to that of the liquid, then the temperature can be maintained practically constant at the substrate-fluid interface. This is the case, in particular, for the silicon nitride substrates used in experiments [41,42,43]. The silicon nitride is a material with high thermal conductivity, approximately in three orders larger than for the toluene. For this reason, the boundary condition for the temperature distribution at the substrate can be reduced to the constant temperature. Heat transfer between the substrate and the drop plays an important role in establishing the temperature profile in the drop. This also excludes a possibility for the reversal of the Marangoni convection [7], taking place for substrates with relatively small thermal conductivity.
The characteristic scale of the temperature variation on the drop surface can be easily estimated if one disregards the effect of velocities on the thermal conduction. According to the evaporation rate data, the heat loss per unit area of the drop near apex is Q 0 = LJ(r = 0) ≈ 0.027 W/cm 2 during the early stage of the evaporation process. Therefore, as follows from the boundary condition (37), we obtain |∂T /∂n| = 22.88 K/cm. The calculations show that the temperature dependence is almost linear on z at r = 0 in disregarding the velocity term in thermal conduction. This permits to obtain the temperature difference between the substrate and the apex of the drop as δT ≈ |∂T /∂n| r=0 h 0 = 3.0087K. This estimate is in a good agreement with the temperature profile along the surface shown in Fig. 11, where the effect of the velocity field on the temperature in Eq. (30) is disregarded. The temperature variation in Fig. 11 has larger amplitude and takes the monotonic form, as compared with the respective curve in Fig. 7. This demonstrates that the fluid flow noticeably modifies the temperature variations in the drop.
Similar consideration for the evaporating water drop studied in [8] results in Q 0 = 0.093 W/cm 2 , |∂T /∂n| = 15.49 K/cm, δT ≈ 0.56 K. This estimation for the temperature difference exceeds in about 24 times the respective numerical results obtained in [8] (see Figs. 1 and 2 there). We consider this discrepancy as an internal contradiction of [8] leading to a significant underestimation of the velocity values in the water drop. Our numerical calculations confirm that δT ≈ 0.6 K for the water drop of [8]. At the same time, our calculations exactly reproduce the velocity field shown in Figs. 4 and 5 in [8], if one takes a given surface temperature distribution of [8]. One should also note that the role of convective term in thermal conduction equation is substantially less for the water drop of [8] due to a small size of the drop.

D. Contact angles
The evaporation process with the pinned contact line is accompanied with an increase of the oblateness of the drop shape. Figure 12 displays our results for time dependence of the contact angle of sessile drop. Our results agree with the time-dependent contact angle determined experimentally under the identical initial drop parameters for the colloidal solution of gold nanoparticles in the toluene drop [42]. In [42] highly ordered nanoparticle islands were found to form on the drop surface. For this reason their orientation angle, which was determined experimentally close to the contact line with small angle X-ray scattering, coincides with the drop contact angle.  We have developed the approach for studying the evaporation and fluid dynamics of a sessile drop of a capillary size and applied it for the description of the toluene drop evaporation. In [41,42,43] the evaporating toluene drop containing colloidal solution of gold nanoparticles was used for realizing the self-assembly of nanocrystal monolayer. During the evaporation process, the nanoparticle islands were found to form on the drop surface [42]. For understanding the island formation as well as nanoparticle segregation, a study of the fluid and evaporation dynamics is necessary. In particular, the convective flows described in this paper are important for driving the nanoparticles to the contact line and the drop surface.
The numerical simulations were carried out and the system of the diffusion equation for the vapor, the thermal conduction equation for the temperature and the Navier-Stokes equations for the fluid flow in the drop is solved. The shape of the drop is controlled by the quasistationary Laplace equation, which includes the effect of the gravitational forces. Their role in forming a profile of a nonspherical sessile drop is characterized by dimensionless number B o = ρghr 0 /(2σ sin θ), which is analogous to the Bond number. Since in our case B o ≈ 0.4, effect of gravitational forces, in general, should be taken into account. We have found that deviations of the sessile drop shape from the spherical cap are noticeable in local drop characteristics such as the local curvature and the evaporation flux density. At the same time the integrated over the drop surface characteristics, like the rate of the mass loss, are well described by the spherical cap approximation.
The experimental and simulation results for the time dependent drop evaporation rate agree well during the main longest stage of the evaporation process. The calculated evolution of the contact angle agrees with the results of [42]. In studying nonstationary diffusion equation, we found that the time-dependent corrections to the stationary local evaporation rates of the sessile drop do not vanish exponentially, but decay much slowly ∝ 1/ √ t with time and, in general, can be noticeable. This kind of behavior was established earlier only for homogeneous diffusion from a surface.
Our solution demonstrates the presence of several time stages in the evolution of the Marangoni convection in the drop. The main quasistationary fluid flow, containing only one vortex in the drop, is formed during the longest period of the evaporation process. Disregarding the time derivatives and the dynamical behavior of the quantities would result directly in a stationary state with one vortex in the sessile drop for a given contact angle. This kind of single-vortex solutions was obtained in previous theoretical studies of the sessile drop evaporation [8,9]. Taking into account an explicit time dependence of all quantities, permitted us to identify a sequence of early dynamical stages of the Marangoni convection in the drop, containing a street of near-surface vortices which then transforms to the state with three bulk vortices. The stability of the large-scale drop dynamics with respect to perturbations of the initial conditions is identified. Also, for a flattened evaporating drop we found under the conditions in question a three-vortex state instead of the single-vortex one. We establish an important role of inertial terms in the Navier-Stokes equations and convective heat transfer terms in the thermal conduction equation.
Our further plans in developing the present approach include investigations of effects of surface surfactants and nanoparticles, effects of a nonstationary dynamics of the drop profile and a detailed study of the influence of substrate properties on the drop dynamics.
Authors thank H.M. Jaeger for fruitful and encouraging discussions. This research was supported by the U.S. Department of Energy Office of Science under the contract DE-AC02-06CH11357 and by the Program of Russian Academy of Sciences. L.Yu. Barash thanks the Dynasty Foundation for the partial support of his research.

Explicit method for temperatures inside the drop
For solving Eq. (30) the following explicit finite difference method is used: Here T n ij is the temperature at the n-th time step, i, j are the coordinates on the regular mesh, h r ,h z are steps along the respective axis, h t is the time step, For r = 0 we have ∂T /∂r = 0, therefore for small r one has T = α + O(r 2 ), i.e. ∂ 2 T /∂r 2 = ∂T /(r∂r). Hence, instead of (A2) we have the following formula: For z = 0 we have T i0 = T 0 .

Boundary interpolation for temperatures inside the drop
Let G = ∂T /∂n = −LJ(r)/k, and D(i, j) is the mesh point which is close to the surface. The point D is inside the drop and at least one of its nearest neighbors has to be outside the drop. In linear approximation T = a + br + cz in a vicinity of the point D(i, j). We denote the temperatures at points B(i − 1, j) and C(i, j − 1) as u B and u C correspondingly. Then The solution of the set of equations is The coefficients a, b, c allow to calculate the temperature at the point D(i, j) using the formula T i,j = a + bh r i + ch z j. Also, the coefficients a, b, c allow to find the temperature at the surface of the drop near the point D.

Alternating direction implicit method for vapor density
In order to cover quite a large region for the vapor density calculation with Eq. (14), this is convenient to use irregular mesh. We use the following mesh: r i = ih r for i < n r , r i = h r (n r + 1.04 i−nr+79 − 1.04 79 ) for i ≥ n r , z j = jh z for j < n z , z j = h z (n z + 1.04 j−nz+79 − 1.04 79 ) for j ≥ n z . Here n r = n z = 700, n r h r = R, n z h z = h. This is the mesh with sufficiently small steps near the drop surface and with exponentially increasing steps, which are chosen in accordance with the asymptotic decay of the vapor density far away from the drop.
We denote the distances between the mesh point (i, j) and its nearest neighbors as a = r i − r i−1 , b = r i+1 − r i , c = z j+1 − z j and d = z j − z j−1 . Then the finite difference representations for the second derivatives are where u n ij is the vapor density at mesh point (i, j) for the n-th time step. We apply the alternating direction method to Eq. (14) with the above notations. In the first part of the method one takes r-derivative implicitly. Then the finite difference representation of Eq. (14) is For given vapor density at time step n it is convenient to rewrite this expression as Here For each j the tridiagonal matrix algorithm is used to solve the set of equations (A13) for the vapor density at the time step n + 1/2. In the second part of the method one takes z-derivative implicitly and represents Eq. (14) as For given vapor density at time step n + 1/2 this expression takes the form where the coefficients are given in (A14)-(A17). For each i the tridiagonal matrix algorithm is used to solve the set of equations (A19) for the vapor density at time step n + 1.

The boundary interpolation for vapor density
Consider a mesh point D(i, j) close to the surface. The point D is outside the drop and at least one of its nearest neighbors is inside the drop. In linear approximation u(r, z) = a + br + cz in a vicinity of the point D(i, j). Consider mesh points B(i + 1, j), C(i, j + 1) and a point A of the drop surface near D. Then one has The solution of the set of equations is In the first part of the alternating direction method the calculations of rows proceed towards smaller values of j. For this reason one should consider here u(C) and u(A) = u s as given quantities, whereas u(B) and u(D) are unknown. It is convenient under these conditions to represent (A21)-(A24) as and obtain explicit expressions for a 0 , a 1 , b 0 , b 1 , c 0 , c 1 . This results in linear relation between u(D) and u(B) which can be transformed to the form This completes the set of equations (A13) for the tridiagonal matrix algorithm. Here To carry out the boundary interpolation in the second part of the alternating direction method, similar expressions can be derived to relate u(D) and u(C).
The obtained values b and c also allow to find the local evaporation rate at the surface point A: The total drop evaporation rate is found from local values (A27) with numerical integration of (13).

Alternating direction implicit method for stream function
The equation to be numerically solved is (compare with (35)) 1 ω The finite difference representation for the second derivatives on a regular mesh arê Taking r-derivative implicitly in the first part of the alternating direction method, we have for Eq. (A28) For a given ψ at time step n it is convenient to rewrite this expression as where For each j the tridiagonal matrix algorithm is used to solve the set of equations (A31) and obtain ψ at time step n + 1/2. In the second part of the method one takes z-derivative implicitly and represents Eq. (A28) as For given ψ at time step n + 1/2 this expression takes the form c ′′ j ψ n+1 i,j−1 + (d ′′ j + 2/(ωh t ))ψ n+1 ij + e ′′ j ψ n+1 i,j+1 = c ′ i ψ n+1/2 i−1,j + (d ′ i + 2/(ωh t ))ψ n+1/2 ij where the coefficients are defined in (A32), (A33). For each i the tridiagonal matrix algorithm is used to solve the set of equations (A35) for the stream function ψ at time step n + 1.

Boundary interpolation for stream function
Consider a mesh point D(i, j) close to the surface. The point D is inside the drop and at least one of its nearest neighbors is outside of the drop. In linear approximation ψ(r, z) = a + br + cz in a vicinity of the point D(i, j). For mesh points B(i − 1, j), C(i, j − 1) and a point A of the drop surface near D one gets a + br A + cz A = 0, a + br B + cz B = ψ(B), a + br C + cz C = ψ(C).
The coefficients a, b, c allow to calculate γ at the point D(i, j) using the formula γ(D) = a + bh r i + ch z j.
Therefore, ψ = const on the surface of the drop, and, hence, ψ = const at all boundaries. Since only the derivatives of the stream function ψ enter expressions for physical quantities, the particular value of the constant does not have any observable consequences. For this reason, one can put ψ = 0 throughout the boundary.