Joint effect of advection, diffusion, and capillary attraction on the spatial structure of particle depositions from evaporating droplets

TL;DR

This paper presents a simplified simulation model to analyze how advection, diffusion, and capillary forces influence the formation of particle chains and deposits in evaporating droplets, aligning with experimental observations.

Contribution

It introduces a computational model that captures the combined effects of advection, diffusion, and capillary attraction on particle deposition patterns during droplet evaporation.

Findings

Chains form due to capillary attraction near the fixing radius.

Close packing occurs if evaporation time exceeds diffusion ordering time.

Chains develop at the final evaporation stage when the fixing radius approaches the droplet center.

Abstract

A simplified model is developed, which allows us to perform computer simulations of the particles transport in an evaporating droplet with a contact line pinned to a hydrophilic substrate. The model accounts for advection in the droplet, diffusion and particle attraction by capillary forces. On the basis of the simulations, we analyze the physical mechanisms of forming of individual chains of particles inside the annular sediment. The parameters chosen correspond to the experiments of Park and Moon [Langmuir 22, 3506 (2006)], where an annular deposition and snakelike chains of colloid particles have been identified. The annular sediment is formed by advection and diffusion transport. We find that the close packing of the particles in the sediment is possible if the evaporation time exceeds the characteristic time of diffusion-based ordering. We show that the chains are formed by the end of the evaporation process due to capillary attraction of particles in the region bounded by a fixing radius, where the local droplet height is comparable to the particle size. At the beginning of the evaporation, the annular deposition is shown to expand faster than the fixing radius moves. However, by the end of the process, the fixing radius rapidly outreaches the expanding inner front of the ring. The snakelike chains are formed at this final stage when the fixing radius moves toward the symmetry axis.

Equations

$$h(r,t)=\theta (t)\frac{R^2-r^2}{2R},$$

$$R_f(t)=\sqrt{R^2-\frac{4r_{p}R}{\theta (t)}}.$$

$$\theta (t)={\theta }_0\left(1-\frac{t}{t_{\max }}\right),$$

$$F\approx \frac{2\pi \sigma Q_1{Q}_2}{L},$$

$$F\approx \frac{2\pi \sigma Q^2}{L}.$$

$$\delta l_{\mathrm{cap}}=\frac{F_{23}-F_{12}}{m_p}\frac{\delta t^2}{2}.$$

$${\bar{v}}_r\approx -\frac{1}{rh}\underset{0}{\overset{r}{\int }}\left(\frac{J}{{\rho }_l}+\frac{\partial h}{\partial t}\right)rdr,$$

$$J(r)\approx \frac{2}{\pi }\frac{D_v(c_s-c_{\infty })}{\sqrt{R^2-r^2}}$$

$$m=\frac{\pi {\rho }_l{R}^3{\theta }_0}{4}=2\pi \underset{0}{\overset{R}{\int }}\underset{0}{\overset{t_{\max }}{\int }}Jrdrdt=4D_v(c_s-c_{\infty })R{t}_{\max },$$

$${\bar{v}}_r=\frac{R}{4\tilde{r}(t_{\max }-t)}\left[\frac{1}{\sqrt{1-{\tilde{r}}^2}}-\left(1-{\tilde{r}}^2\right)\right],$$

$$t_{\max }\approx \frac{{\rho }_l{h}_{0}R}{\pi D_v{c}_s\left(1-H\right)\left(0.27\theta +1.3\right)}\approx 0.3\text{ s,}$$

$$w=R\sqrt{\frac{\varphi }{4p}}{\left[1-{\left(1-\frac{t}{t_{\max }}\right)}^{3/4}\right]}^{2/3},$$

Full text

I Introduction

The processes of heat and mass transfer occurring in droplets and films are of interest in various applications such as fuel droplet evaporation and combustion in engines Al Qubeissi et al. (2018); Hänichen et al. (2019), interaction of droplets with surfaces of varying wettabilities in ink-jet printing Sprittles and Shikhmurzaev (2012), obtaining stable ultrathin film surfaces on the basis of polar liquids Gordeeva and Lyushnin (2017), removing nanoparticles from a solid surface Mahdi et al. (2015), and many other applications.

One of the most important and actively discussed problems is connected to studying structures of colloidal particles, which emerge on the surface of an evaporating sessile droplet and remain on the substrate after drying Larson (2014); Thiele (2014); Erbil (2015); Sefiane (2014); Brugarolas et al. (2013); Brutin and Starov (2018); Tarafdar et al. (2018); Giorgiutti-Dauphiné and Pauchard (2018). One of the examples is the effect of evaporative contact line deposition, the so-called coffee-ring effect Deegan et al. (1997, 2000); Mampallil and Eral (2018). While a droplet is drying on the substrate, capillary flows carry the colloid particles toward the three-phase boundary. In this case, formation of an annular deposition is observed if the contact line was pinned throughout the entire process. According to the experimental results, coffee-rings may have an inner structure Marín et al. (2011), and colloid particles can also merge into chains and other geometric shapes Park and Moon (2006) (Fig. 1). However, the existing models which describe this phenomenon Deegan et al. (2000); Fischer (2002); Vodolazskaya and Tarasevich (2011); Ozawa et al. (2005) do not allow predicting the inner structure and location of the particles relatively to each other.

Effects of this kind are studied experimentally in Park and Moon (2006); Marín et al. (2011); Zhang et al. (2013); Yunker et al. (2011, 2013); Mayarani et al. (2017); Li et al. (2008); Andac et al. (2019); Govor et al. (2004a, b); Govor (2009). Marín et al. primarily explained the mechanism of these effects by convective transfer and diffusion Marín et al. (2011). Park and Moon obtained circular depositions and stains as a result of evaporating picoliter droplets based on water ink and different solvent mixtures Park and Moon (2006). Such depositions are characterized by a well-ordered hexagonal structure of particles. In case of the annular deposition, the images show the presence of independent chains of particles within the ring.

The role of the Marangoni flow in case of a heated substrate is discussed in Zhang et al. (2013). Ternary mixtures of solvent, polymer and semiconductor nanocrystal were investigated in Miller et al. (2014). Yunker et al. studied the influence of particle shape on the resulting structure of deposition Yunker et al. (2011, 2013). Depositions of microgel soft particles are characterized by a gradual order-to-disorder transition Mayarani et al. (2017), which is not abrupt like the hard colloid depositions Mayarani et al. (2017); Marín et al. (2011). The transition from order to disorder is also observed in the structures obtained by a dip-coating technique Kaplan et al. (2015); Mehraeen et al. (2015); Colosqui et al. (2013). Li et al. used solvents with dendrimers and obtained annular structures of a periodic thickness Li et al. (2008). Callegari et al. studied the influence of active matter (moving bacteria) on the dynamics of the coffee-rings growth Andac et al. (2019). A series of papers is dedicated to the formation of nanoparticle chains along the contact line of micron-size droplets Govor et al. (2004a, b); Govor (2009).

The network deposition patterns inside the annular sediment have been obtained experimentally for a drying sessile droplet in Li et al. (2017), where a tailored substrate with a circular hydrophilic domain is used. The evolution of a dry patch is found to be divided into three stages: rupture initiation, dry patch expansion, and drying of the residual liquid. The authors argue that in their case the capillary attraction is important, while the Marangoni flows are considerably large only in the first phase of evaporation, when the contact angle is relatively large.

The capillary attraction of the particles in liquid films was studied in Kralchevsky and Nagayama (1994); Denkov et al. (1992, 1993); Park and Furst (2011); Lotito and Zambelli (2017); Tang and Cheng (2018). Based on such interaction between the particles, a lithographic method of self-assembly of nanostructures was developed Liddle et al. (2004).

There are several discrete models, which explicitly describe the dynamics of each particular particle. These models are usually based on Monte Carlo or molecular dynamics methods Liao et al. (2000); Gundabala (2003); Lebedev-Stepanov and Vlasov (2013); Mehraeen et al. (2015); Shi and Webb (2016); Akkus et al. (2019); Tang and Cheng (2018). Some Monte Carlo based models are lattice models Zhang et al. (2016); Kim et al. (2011); Jung et al. (2014), and some are not Petsi et al. (2010); Andac et al. (2019); Kolegov (2019). The flow velocity is often calculated analytically for simple particular cases Lebedev-Stepanov and Vlasov (2013); Kim et al. (2011); Petsi et al. (2010); Jung et al. (2014); Crivoi and Duan (2014); Kolegov (2019). The obtained velocity field is taken into account for calculations of the particles dynamics. Callegari et al. Andac et al. (2019) use a constant velocity value. Mobile bacteria are modeled by the Monte Carlo method Andac et al. (2019). Lebedev-Stepanov and Vlasov took account of the adhesion of particles and roughness of the substrate Lebedev-Stepanov and Vlasov (2013). The formation of branched nanoparticle aggregates in the drying sessile droplet is modeled in Zhang et al. (2016), where the two-dimensional kinetic Monte Carlo Ising-like approach is used. The results pertain to those experiments, where nanoparticles are considered, and the coffee ring effect is absent. The branched aggregates begin to form when the thickness of the ultrathin film reaches several angstroms and local ruptures appear. On such a scale, the process can be accompanied by the formation of other structures, such as networks, worm-like structures and so on Stannard (2011). The models Kim et al. (2011); Crivoi and Duan (2014) describe the formation of annular depositions in the process of droplet evaporation using a cubic lattice. Petsi et al. compared the Monte Carlo and lattice Boltzmann equation methods for a semi-cylindrical geometry of the droplet and also provided a numerical description of the process of particles deposition on the substrate as a result of the transfer by the compensatory fluid flow and their Brownian motion Petsi et al. (2010). The computational experiment was carried out for different wetting angles and three-phase boundary modes. Jung et al. used a model with a circular lattice to describe the coffee-ring effect mathematically Jung et al. (2014). A computational method which allows one to calculate the capillary attraction forces between particles in a thin film is presented in Zhou et al. (2006).

The objective of the present work is to analyze physical mechanisms of the formation of standalone particle chains located within an annular monolayer deposition Park and Moon (2006) (Fig. 1), where the particle size is of the order of micrometers. Our hypothesis is that this effect is caused by the particles’ capillary attraction, which occurs at different times in those areas of the droplet where the thickness of the liquid layer reduces to the size of particles. To verify this mechanism, we employ a simplified mathematical model, which accounts for joint consideration of advection, diffusion, and capillarity and allows us to perform computer simulations of the transport of particles.

II Methods

II.1 Problem formulation

Let us consider an evaporating sessile droplet of colloid solution (Fig. 2). The radius of the spherical particles of $\mathrm{SiO_{2}}$ $r_{p}\approx$ 0.35 $\mu$m Park and Moon (2006).

We estimate the sedimentation time for a single particle by dividing the droplet height $h_{0}=R\theta/2\approx 5$ $\mu$m by the Stokes velocity $U_{0}=2r_{p}^{2}\Delta\rho g/9\eta\approx 0.5$ $\mu$m/s Russel et al. (1991), where $\Delta\rho=1.65$ g/cm3 is the density difference between the particles and the water, $g$ is acceleration due to gravity, and $\eta$ is the fluid viscosity. Hence, the sedimentation time is of the order of $t_{s}\approx 10$ s. The advection further increases the sedimentation time, which exceeds the evaporation time $t_{max}\leq 10$ s. For this reason, we do not take account of sedimentation.

Although the particles of silicon dioxide may acquire a charge during their chemical synthesis Behrens and Grier (2001), we assume that the screening is significant due to a large concentration of ions in the aqueous solution. For this reason, we disregard the electrostatic interactions Linse (2005); Liang et al. (2007).

The shape of a sessile droplet profile is a spherical cap due to the predominance of capillary forces over gravity because of the small volume of liquid ($V\approx 72$ pl Park and Moon (2006)). Such geometry facilitates working with a cylindrical coordinate system ($r,\varphi,z$). The direction of axis $r$ is parallel to the substrate. The symmetry center of the droplet corresponds to $r=0$. Axis $z$ is perpendicular to the substrate. The surface of the solid base corresponds to $z=0$. The three-phase boundary is fixed throughout the entire process at $r=R$. The contact line radius is $R\approx$ 60 $\mu$m. The substrate is made from a hydrophilic material, and the droplet is flattened. Its contact angle $\theta\leq$ 10∘ Park and Moon (2006). In such cases, one or two layers of particles are formed Park and Moon (2006). In this work, we only consider the case with a monolayer formation from particles.

We take into account radial flow of the liquid, directed from the symmetry axis toward the droplet periphery. This compensatory flow is known as capillary flow and is a result of intense evaporation near the contact line (Fig. 2). We do not take account of circulating flows of the liquid, because the liquid layer is thin (see Sec. II.2.3).

II.2 Description of the model

II.2.1 Capillary attraction

It is convenient to use both polar ($r,\varphi$) and Cartesian coordinates ($x,y$). Since thickness of the liquid layer $h$ in the flattened droplet is much smaller than $R$, the particle dynamics can be described in a horizontal plane, $xy$. Despite the fact that we study the case with a fixed three-phase boundary, the $xy$ plane also has a movable boundary, $R_{f}$. Let the fixing radius $R_{f}$ define a boundary, where the particle size and the local droplet height are comparable ($d_{p}\approx h$), provided that $R_{f}<R$. Here $d_{p}=2r_{p}$. In a thin droplet, the particles cannot reach the contact line, because the local droplet height is very small in the vicinity of the contact line. The particles will not move farther than $R_{f}$, since beyond this boundary $h<d_{p}$. The surface tension forces will restrain them. Let us write the approximate expression for the shape of the droplet surface Popov (2005)

[TABLE]

where radial coordinate is $r=\sqrt{x^{2}+y^{2}}$. From (1) we derive the dependence of the fixing radius on time

[TABLE]

The droplet height and the contact angle reduce during the evaporation process, and, hence, the fixing radius $R_{f}(t)$ is changed as well (Fig. 2). It should be noted that $R_{f}$ is changed insignificantly for the majority of time. Only at the end of the process, this radius starts shrinking. At $t\to t_{\mathrm{max}}$ the value of $R_{f}\to 0$. Contact angle $\theta$ decreases linearly over time, $\theta\to 0$ at $t\to t_{\mathrm{max}}$ Hamamoto et al. (2011), where $t_{\mathrm{max}}$ is the time of full evaporation. This dependence can be expressed as follows:

[TABLE]

where $\theta_{0}$ is the initial value of the contact angle Popov (2005). The exact value of $\theta_{0}$ is not determined in the experiment Park and Moon (2006). It is only known that it does not exceed 10∘ for the case of the hydrophilic substrate. For the calculations, we will use $\theta_{0}=$ 10∘.

In the region of boundary $R_{f}$, the particles are subject to the capillary attraction Denkov et al. (1992); Kralchevsky and Nagayama (1994); Lotito and Zambelli (2017)

[TABLE]

where $L$ is the distance between particles, and $Q_{k}$ is the capillary charge of the $k$th particle, and $\sigma$ is the surface tension coefficient of the liquid. Expression (4) resembles Coulomb’s law. Here, $Q_{k}$ characterizes the ability of a particle to deform the free surface of the liquid, $Q_{k}=r_{p}\sin\psi_{k}$. The radius of any particle is assumed to be constant value $r_{p}$. The meniscus slope angle $\psi_{k}$ may vary depending on the location of a particle. To carry out assessment, let us consider the average of $\psi$ (Fig. 3). In this case, $Q_{1}=Q_{2}=Q$, where $Q=r_{p}\sin\psi$. Then, expression (4) can be written as follows:

[TABLE]

Let us study a three-particle system as an example (Fig. 3) and assume that angle $\psi=$ 1∘.

We assume that the left and the right particles are attached to the substrate. Let us determine the displacement of the central particle over a time step of $\delta t=$ 10*-4* s. The choice of the time step value is discussed in Sec. II.3. Let the distance between the first and second particles be $L_{12}=11r_{p}$, and the distance between the second and third particles be $L_{23}=12r_{p}$. Now, let us determine the values of $F_{12}$ and $F_{23}$ by substituting $L_{12}$ and $L_{23}$ in (5). The second particle displacement caused by capillary forces is

[TABLE]

Mass of a particle $m_{p}=V_{p}\rho_{p}\approx$ 0.5 pg, where $V_{p}$ is the volume of a spherical particle, and $\rho_{p}$ is the density of silicon dioxide. We obtain the value of relation $\delta l_{\mathrm{cap}}/L_{12}\approx-1$. This means that the second particle shifts tightly to the first particle, since the distance between them is smaller and the capillary effect is stronger. It will occur more rapidly at a larger angle $\psi$.

II.2.2 Advection

Based on the mass conservation law Deegan et al. (2000) and taking into consideration the lubrication approximation $\sqrt{1+(\partial h/\partial r)^{2}}\approx 1$ we express the height-averaged velocity of radial flow of the liquid $\bar{v}_{r}(r,t)$,

[TABLE]

where density of the liquid $\rho_{l}\approx 10^{3}$ kg/m3, and $J=J(r,t)$ is the local evaporation flux with a dimension of kg/(m2 s).

If we consider the evaporation as only a diffusion of vapor into air and assume the process is quasistationary, then, by solving the Laplace equation for the distribution of vapor concentration $\Delta c=0$, we find $J$. In this approximation, $J$ only depends on $r$. The boundary condition at the liquid-vapor interface is $J(r)=D_{v}\left.(\mathbf{n}\cdot\nabla c)\right|_{z=h}$, where $\mathbf{n}$ is the normal vector to the droplet surface, and $D_{v}$ is the coefficient of vapor diffusion. Analytical solution

[TABLE]

had been obtained earlier for $\theta\to 0$ Deegan et al. (2000); Popov (2005); Hu and Larson (2005). Here $c_{s}$ is the concentration of saturated vapor, and $c_{\infty}$ is the concentration of vapor in the surrounding air far from the droplet surface.

Let us substitute (1), (3) and (7) in (6) and then integrate. Considering the fact that the mass of the droplet $m$ is equal to the mass of the entire evaporated liquid Popov (2005),

[TABLE]

we find

[TABLE]

where $\tilde{r}=r/R$ Popov (2005); Hu and Larson (2005). Displacement of the particle by the radial fluid flow over one time step is $\delta l_{\mathrm{adv}}=\bar{v}_{r}\delta t$.

Equation (8) is a rough approximation for an actual fluid velocity, so the suggested model is only of a phenomenological nature. The process can be considered as quasistationary when $R^{2}/D_{v}\ll t_{\mathrm{max}}$ Popov (2005). In order to further refine the suggested model, one has to consider $J$, which depends on time $t$, as well. For example, dependencies $J(r,t)$ are suggested in Murisic and Kondic (2008); Barash et al. (2009); Cazabat and Guéna (2010); Semenov et al. (2017a).

II.2.3 Marangoni flows

It can be shown that the capillary flow dominates the Marangoni effect in the fixing radius region. As seen in Eq.(8), the velocity of capillary flow substantially increases both on approaching the contact line and during the final phase of the evaporation process. Since for large contact angles the fixing radius (2) is close to the contact line, the capillary flow dominates the Marangoni flow there. The velocities of Marangoni flow are proportional to $\theta^{2}$, and decrease rapidly during the final phase of the evaporation process. Hence, at small contact angles, the capillary flow dominates the circulatory one as well.

It should be noted that the thermocapillary instability is known to occur only when the Marangoni number $\mathrm{Ma}=-\sigma^{\prime}\Delta Th_{0}/(\eta\kappa)$ exceeds a certain critical value. Here $\Delta T$ is the temperature difference between the apex and the substrate, $\sigma^{\prime}=d\sigma/dT$, and $\kappa$ is the thermal diffusivity. The threshold value for a flat fluid film is about $80$ Pearson (1958). The onset of the Marangoni convection has been also confirmed in MacDonald and Ward (2012) for the particular case of droplets with contact angle $\theta=90^{\circ}$. Using the known estimates of the evaporation flux Hu and Larson (2002), we obtain the value $\mathrm{Ma}\approx 0.4$ for $\theta=10^{\circ}$ and $\mathrm{Ma}\approx 0.03$ for $\theta=3^{\circ}$. Although the onset is not known for droplets with small contact angles, it is possible that for such small values of $\mathrm{Ma}$ the thermocapillary flow does not occur at all.

For this reason, we do not take into account Marangoni flows in the simulations.

II.2.4 Diffusion

Using the generalized Stokes-Einstein equation Yiantsios and Higgins (2006); Russel et al. (1991), we estimate the diffusion coefficient of the particles $D\approx 2\cdot 10^{-13}$ m2/s. The diffusion displacement distance is $\delta l_{\mathrm{dif}}=\sqrt{2D\delta t}\approx$ 6 nm. We obtain the relation: $\delta l_{\mathrm{cap}}/\delta l_{\mathrm{dif}}\approx 615$. Hence, capillary attraction of the particles prevails over diffusion, which can be neglected on achievement of fixing radius $R_{f}$ by the particle. The characteristic time of diffusion-based ordering of particles is $t_{d}=d_{p}^{2}/D\approx$ 2.45 s. The evaporation time of the droplet depends on the temperature, humidity, air pressure, concentration of colloidal particles, and other parameters. The precise value of evaporation time $t_{\mathrm{max}}$ is not given in Park and Moon (2006). For a pure liquid, there is an approximation

[TABLE]

where the coefficient of diffusion is $D_{v}\approx 2.4\cdot 10^{-5}$ m2/s, concentration of saturated vapor $c_{s}\approx 17.3\cdot 10^{-3}$ kg/m3 Larson (2014). The value of the contact wetting angle is $\theta=\pi/18$ and height of the flattened droplet $h_{0}\approx R\tan(0.5\theta)\approx$ 5 $\mu$m. In the experiment Park and Moon (2006), relative humidity of the surrounding air is $H\approx 0.4$. A number of experiments Okuzono et al. (2010); Kim and Weon (2018) demonstrate that the evaporation rate can decrease with increasing concentration of colloidal particles. For this reason, we consider $t_{\mathrm{max}}$ as a parameter. In calculations, we will use $t_{\mathrm{max}}=$ 1 s $<t_{d}$ and $t_{\mathrm{max}}=$ 10 s $>t_{d}$.

Brownian motion of the particles is simulated using a Monte Carlo method. We denote a random angle between a particle displacement and the $x$ axis as $\alpha$, $\alpha\in[-\pi;\pi)$. Then the displacement of a particle is described by the vector $(\delta x,\delta y)^{T}=(\delta l_{\mathrm{dif}}\cos\alpha,\delta l_{\mathrm{dif}}\sin\alpha)^{T}$.

II.3 The algorithm

The pseudocode of the algorithm is given in the Appendix. We use the Mersenne Twister generator in order to generate pseudorandom numbers Matsumoto and Nishimura (1998). Efficient implementations for the current computer architectures are contained in modern software libraries for random number generation Barash and Shchur (2011); Guskova et al. (2016); Barash and Shchur (2014).

The particles can be in different states during the calculation. This approach resembles a finite state machine (Fig. 4). The Supplemental Material video [77] shows results of the simulations and displays the states by different colors. The particles subject to advection and diffusion are marked in green. Yellow particles are only subject to diffusion. Black particles are affected by capillarity. Red particles are the ones that stopped their motion and are attached to the substrate. Conventionally, below we use this color notation of the states, because it is identical to the color notation in the Supplemental Material video.

Initially, according to the color designation, the particles are green by default. A particle becomes black on reaching fixing radius $|R_{f}-r|\leq r_{p}$. If the particle offsets beyond boundary $R_{f}$ toward $R$ when $r>R_{f}+r_{p}$, it becomes red.

The particles within the forming ring are more affected by diffusion than by fluid flow. This can be seen, e.g., in the supplemental video in Yunker et al. (2011). We assume that a green particle becomes yellow if its attempt to move was failed due to a collision with a red or black particle. In this context, collision is an overlap of the particles, i.e., in the two-dimensional model representation, an overlap of the circles of radius $r_{p}$. The particles are of a solid material and are not deformed. The regions of space occupied by the particles cannot cross even partially.

If the yellow particle has moved away by a small distance from $R_{f}$ due to the Brownian motion toward the symmetry axis (e.g., $r<R_{f}-5d_{p}$), it becomes green again. In other words, the liquid flow starts moving it toward the periphery again.

The capillarity of particles on the basis of assessments using formula (5) in Sec. II.2.1 will be simulated as follows. Over one time step, a black particle is attracted and moves closely to the nearest black or red particle in the neighborhood (within the surrounding region of a radius $R_{n}$). For a calculation, we take $R_{n}=20r_{p}$ Park and Furst (2011). The value of time step $\delta t=$ 10*-4* s was chosen on the basis of a series of computation experiments to satisfy the Einstein relation for the mean square displacement $\langle\delta L_{\mathrm{dif}}^{2}\rangle=2Dt_{\mathrm{max}}$, where $\delta L_{\mathrm{dif}}$ is the total displacement during the period $t_{\mathrm{max}}$, and the averaging is performed over all particles.

We do not consider collisions of the particles. Since the Monte Carlo method updates the particle motion at each time step, its stochastic character excludes a possibility to describe collisions of particles given such a small time step $\delta t$. The diffusion displacement $\delta l_{\mathrm{dif}}$ during a time step is about 1% of the particle size. Particle displacement by advection $\delta l_{\mathrm{adv}}$ usually does not exceed 1.5% of the particle size. This results in a negligibly small portion of colliding particles during one step.

III Results and discussion

We chose the number of particles $N_{p}=9000$ such that the width of the formed annular deposition corresponded to the value from the experiment Park and Moon (2006). The time $t_{\mathrm{max}}$ was considered as a parameter of the simulations. Figure 5 represents the obtained structures of depositions at the drying time of $t_{\mathrm{max}}=$ 1 s (left) and 10 s (right). Figure 6 represents the radial profile of the particle density distribution $\rho$, which is one of the key physical quantities studied. Here, $\rho$ is the number of particles per unit area ($\rho$ is normalized to $N_{p}/R^{2}$). The profile is demonstrated in two columns of panels for droplets with two substantially different characteristic values of $t_{\mathrm{max}}$, respectively. Five rows of panels in each of the two columns show the radial profile calculated for five combinations of the main effects considered, namely, calculated by taking into account solely the advection [Figs. 6(a) and 6(b)], or the effects of diffusion and capillary attraction [Figs. 6(c) and 6(d)], or the advection and capillary attraction [Figs. 6(e) and 6(f)], or the advection and diffusion [Figs. 6(g) and 6(h)], or the advection, diffusion and capillary attraction jointly [Fig. 6(i) and 6(j)]. In parallel with these ten panels, the panels of Fig. 5 show the corresponding final spatial particle structures formed in each case. Such a representation allows one to see the relative role of each of the effects in forming the spatial particle structure on the substrate at four characteristic times chosen, that correspond to four curves in each of the panels in Fig. 6. The calculation error does not exceed the size of the marker. Each calculation was repeated ten times. The lines connecting the markers are shown for convenience only and do not have a physical meaning.

Based on the simulation results, time $t_{\mathrm{max}}$ influences the forming structure only in the cases when both diffusion and advection are taken into account [Figs. 5(g)–5(j)]. The fluid flow carries the particles toward the periphery. A ring with unordered particles is formed near the contact line over time $t_{\mathrm{max}}=$ 1 s. Thus, if $t_{\mathrm{max}}<t_{d}$, the particles do not have enough time to pack densely due to the Brownian motion. We have calculated the packing fraction $p\approx 0.79$. The coffee rings formed during $t_{\mathrm{max}}=$ 10 s, are characterized by almost hexagonal package of particles. In this case, the ratio of sediment area in the ring to ring area is approximately 0.84. This ratio is about 0.91 in hexagonal packaging. The difference is due to visible defects and local subdomains in the obtained structure. In particular, the defects are induced by the annular geometry of the region. Over $t_{\mathrm{max}}>t_{d}$, the particles have enough time to densely pack due to diffusive mixing. Thus, as a result of transfer of the particles by the advective flow alone, the total density of the particles near the periphery is not so high [$\rho\approx$ 0.56 in Figs. 6(a) and 6(b)], while under the additional influence of diffusion $\rho\approx$ 0.81 for $t_{\mathrm{max}}=$ 1 s [Figs. 6(g) and 6(i)] and $\rho\approx$ 0.86 for $t_{\mathrm{max}}=$ 10 s [Figs. 6(h) and 6(j)]. The influence capillarity produces is that, instead of a uniform value of $\rho$ in the central part of the deposition [Fig. 6(g) and 6(h)], the final density distribution of the particles grows linearly from the symmetry axis ($\tilde{r}=0$) toward the ring boundary ($\tilde{r}\approx 0.7$) in Figs. 6(i) and 6(j). We find that $\rho(\tilde{r}=1)=0$, since the particles do not move beyond boundary $R_{f}(0)$ toward $R$ (Fig. 6).

In the case of studying only advection or its combination with capillarity, treelike structures of particles are observed within the ring [Figs. 5(a), 5(b), 5(e), and 5(f)]. Such dendritic shapes are not characteristic to diffusing particles. In the absence of the Brownian motion, the particles do not mix within the ring. A particle driven by the fluid flow stops moving when it encounters an obstacle. Then, this particle can represent an obstacle for another particle. This is the essence of the dendritic shape formation mechanism.

If the system is only affected by diffusion and capillarity, the depositions also have dendritic structures [Fig. 5(c) and 5(d)]. The difference is that such structure occupies the entire surface area, which was previously in contact with the liquid layer. Particle density distribution is almost uniform ($\rho\approx 0.35$) in this case [Fig. 6(c) and 6(d)]. Besides, the “tree branches” turn out sparse, since they are not formed due to advection. The particles are uniformly distributed in the entire volume and fluctuate chaotically at the beginning of the process. As the fixing radius decreases, the particles occurring at the fixing boundary ($\left|R_{f}-r\right|\leq r_{p}$), are drawn to each other and stop their movement. They are in the state of rest for several reasons. Firstly, the local droplet height is smaller than the particle size, therefore, the particles lie on the substrate. Secondly, the capillary forces press them to the substrate, because a thin film wets these particles. At the beginning of the process, the motion of the boundary $R_{f}$ is slow, and separate strips of particles are formed. These strips are located near the periphery perpendicular to the radial direction. There is a small dip in the particle density distribution [Fig. 6(c) and 6(d)] right before the deposit, e.g., for $t=0.6t_{\mathrm{max}}$. The dip is located exactly at the fixing radius area, where the local droplet height is comparable to the particle size at each time, as can be seen in the Supplemental Material video 111See Supplemental Material at [http://link.aps.org/

supplemental/10.1103/PhysRevE.100.033304](http://link.aps.org/supplemental/10.1103/PhysRevE.100.033304) for the simulation videos. The dip can be explained by the capillary attraction acting in this small region, while the tree-like structure of particles is still not finalized and the density is still increasing.

It should be noted that in cases when advection is taken into account but capillary attraction is not, uniform distribution of the particles is observed in the central region [Fig. 5(a), 5(b), 5(g), and 5(h)]. These are the particles not yet carried by the flow toward the annular deposition at the droplet periphery. When capillarity is added to the model, we see how these particles form chains [Fig. 5(e), 5(f), 5(i), and 5(j)]. Capillarity does not affect the structure inside the forming ring except drawing the particles, which were in the proximity $R_{f}$, toward each other. Most of the time, the change of $R_{f}(t)$ occurs slower than the growth of ring width $w(t)$. By the end of the process, the situation becomes opposite, $R_{f}(t)<R_{f}(0)-w(t)$ at $t>t_{\mathrm{cri}}$. Then the chains start forming from the particles as they appear at fixing boundary $R_{f}$ (Fig. 7).

According to the numerical results, critical time $t_{\mathrm{cri}}\approx(0.74\pm 0.06)t_{\mathrm{max}}$. It is possible to obtain the time $t_{\mathrm{cri}}$ from the following equality $R_{f}(t_{\mathrm{cri}})=R_{f}(0)-w(t_{\mathrm{cri}})$. An analytical formula for $w(t)$ was obtained in Deegan (2000),

[TABLE]

where $\phi$ is the volume fraction of colloidal particles. Let us consider the critical time for a particular example, when three effects affecting mass transfer are taken into account ($t_{\mathrm{max}}=1$ s). As shown in Fig. 8, the functions $R_{f}(t)$ and $R_{f}(0)-w(t)$ overlap at $t_{\mathrm{cri}}\approx 0.74t_{\mathrm{max}}$: $R_{f}(t_{\mathrm{cri}})=R_{f}(0)-w(t_{\mathrm{cri}})=r_{\mathrm{cri}}\approx 0.7R$. It is seen that the value $t_{\mathrm{cri}}$ corresponding to the intersection of the $R_{f}(t)$ and the numerically obtained $R_{f}(0)-w(t)$ is only slightly smaller than $0.74t_{\mathrm{max}}$ (Fig. 8).

The result of modeling the dynamics of the inner front of the annular sediment is quantitatively different from the analytical prediction of (9) at the final stage (Fig. 8). The value of the function $R_{f}(0)-w(t)$ at $t>t_{\mathrm{cri}}$ remains constant. This means that the width of the ring does not increase further in time, as in the case of using the formula (9), because Eq. (9) does not take into account the existence of the boundary $R_{f}$.

Branching in the polar direction is mostly suppressed in the region where radial motion of the boundary $R_{f}$ is faster than the characteristic velocity of advection and diffusion (60 $\mu$m/s), while the particle density is low ($\rho<0.2$) (see Figs. 5, 6, and 8).

Considering jointly advection, diffusion, and capillarity for $t_{\mathrm{max}}=10$ s allowed obtaining a numerical result [Fig. 5(j)], which is qualitatively consistent with the experimental data in Fig. 1. The structure of the deposition shows close-packed particles in the annular part and snakelike chains of particles in the central region. The difference is that these chains are short in the experiment. The numerical results show long chains of particles. It is most likely connected to the used approximation in the simulation of the particles capillarity.

IV Conclusions

The formation of both close-packed annular deposition at the periphery of the droplet and the snake-like chains of particles in the central region has been described based on the suggested model. A major part of the colloidal particles is transferred toward the periphery due to advection. The fluid flow toward the contact line originates from the nonuniform evaporation rate along the free surface of the droplet. This results in the formation of the ring of colloid particles. The hexagonal structure within the ring is found to be induced by the diffusion-based ordering of the particles when the evaporation time $t_{\mathrm{max}}$ exceeds the characteristic time $t_{d}$. Chains of particles are formed at the final stage in the flattened droplet as a result of the particles capillarity, when the fixing boundary outreaches the inner boundary of the expanding ring, $R_{f}(t)<R_{f}(0)-w(t)$.

The presented model jointly takes into account the effects of advection, diffusion, and capillarity on the particle motion. We consider only the case when the concentration of colloidal particles and the droplet size allow a formation of a monolayer of particles. Our simplified model and the corresponding numerical simulations allowed us to study some of the effects discovered in the experiment Park and Moon (2006). Since the results obtained are in a qualitative agreement with the observations, we believe that the model in this case accounts for the major physical mechanisms of the deposition microstructure formation.

Acknowledgements.

This work is supported by Grant No. 18-71-10061 from the Russian Science Foundation. The authors would like to thank Yuri Yu. Tarasevich, Irina V. Vodolazskaya, Renat K. Akhunzhanov, and Yury A. Budkov for useful discussions and remarks.

V Appendix

The pseudocode of the algorithm is given in the following Algorithm 1.