Molecular dynamics simulation of nanosized water droplet spreading on chemically heterogeneous surfaces via complex surface trapping

with Boltzmann simulations: A wetting state is dynamically trapped ABSTRACT The wetting kinetics of water droplets on chemically heterogeneous surfaces is important in several industrial technologies, such as biomedicine and microfluidics. Surfaces with different wettabilities can be designed to control the spread of droplets. In this study, nanosized water droplet spreading on chemically heterogeneous surfaces was investigated using molecular dynamics simulations. Chemically heterogeneous surfaces with different wetting patterns were investigated, and the equivalent spreading radius and dynamic contact angle during the spreading process were analyzed. Results showed that droplet spreading is mainly dependent on the area fractions of hydrophobic and hydrophilic regions and the shape of the wetting pattern has a minor influence on the spreading process. The dynamic contact angle can be well predicted by molecular kinetics theory. The static contact angle data remarkably deviate from Cassie’s equation, while they agree better with the modified Cassie’s equation as a function of the hydrophobic length fraction, indicating that the wetting pattern has a substantial influence in the vicinity of the contact line.


I. INTRODUCTION
When a droplet is deposited on a solid surface, it is influenced by surface tension and viscous force, and spontaneously spreads toward equilibrium. The wetting kinetics of droplet spreading on a solid surface is important to numerous practical applications, such as evaporation, 1 biomedicine, 2 microfluidics, 3-7 and self-cleaning. 8,9 Droplet spreading on a flat and homogeneous substrate follows a power law, R b ∼ t n , where R b is the spreading radius, t is the spreading time, and n is the spreading index depending on the spreading mechanism. Tanner's law 10 assumes that viscous dissipation is the dominant resistance of wetting and gives n = 1/10, whereas molecular kinetic theory (MKT) 11,12 considers the nonhydrodynamic energy dissipation at the contact line as the dominant resistance and gives n = 1/7. In addition, droplet spreading on a smooth surface is faster than the prediction of Tanner's law at the initial stage. [13][14][15][16] The index is found to be 1/2 for complete wetting.
As the droplet reaches a static condition, the equilibrium contact angle θ can be described by Young's equation: 17 cos where θ Y is Young's contact angle, γ is the surface tension coefficient, and the subscripts S, L, and G respectively refer to the solid, liquid, and gas phases. However, heterogeneities widely exist in real life, and no strictly flat and homogeneous surface exists. A surface can be physically and chemically heterogeneous. Physically heterogeneous surfaces refer to the surfaces with physical structures, e.g., the lotus leaf is waterrepellent due to the papillose epidermal cells on its surface, which form papillae or microasperities; they are artificial surfaces decorated with a micropillar array to obtain specific wettability in the applications of biomedicine, 2 microfluidics, [3][4][5][6][7] and self-cleaning. 8,9 Chemically heterogeneous surfaces refer to the surfaces with heterogeneous wetting patterns, e.g., contaminated surfaces which exhibit ARTICLE scitation.org/journal/adv different wetting behaviors; they are artificial surfaces with stripes of hydrophilic and hydrophobic regions realizing anisotropic wetting. Spreading of droplets on physically and chemically heterogeneous surfaces can be rather different from that on smooth surfaces. Spreading of droplets on physically heterogeneous surfaces has been extensively investigated. Nanostructured surfaces are typical physical heterogeneous surfaces. The droplet may behave at the Cassie 18 or Wenzel 19 state on the nanostructured surfaces. In the Wenzel state, the upper part of the droplet collapses, while the bottom part penetrates into the gaps between nanostructures, forming a fringe film; 20 in the Cassie state, the liquid stays above the nanostructures with air trapped in the gaps. The spreading radius can be well described by the exponential correlation, with the spreading index ranging from 0.108 to 1/3. [20][21][22][23] The spreading index depends on several parameters such as viscosity, surface tension, and solid fraction 24 and varies with the sizes of the droplet and the nanostructure on the surface as well as the geometry of the nanostructure.
Compared to physically heterogeneous surfaces, droplet behavior on chemically heterogeneous surfaces is still insufficiently understood. When a droplet is placed on a chemically heterogeneous surface, the three-phase contact line of the droplet is distorted by the regions of different wettability, with the microscopic contact angle changing along it. 25,26 Several studies have also pointed out that the wetting characteristics, such as the static contact angle, are not highly dependent on the length scale of the alternating lyophilic and lyophobic regions for droplets of several micrometers if they cover multiple regions. [27][28][29] Meanwhile, other results have shown that some changes can occur with the change of length scale. [30][31][32][33] When sliding on a surface with alternating lyophilic and lyophobic regions, the wetting behavior of droplet relies on the size of regions even if the droplet is sufficiently large. As for the nanoscale, the droplet cannot be regarded as a continuous liquid. Whether the size and shape of the wetting pattern influence the spreading process and whether the power law can be applied remain unclear. As a droplet reaches the equilibrium state on chemically heterogeneous surfaces, Cassie's equation is widely used to describe the static contact angle 18 of the droplet. Adao et al. showed that the static contact angle of a nanodroplet only depends on the area fraction of the hydrophobic region and is predicted well by Cassie's equation. 34 Several investigations pointed out that Cassie's equation holds for droplets sufficiently large compared to the length scale of heterogeneity. [27][28][29]35 Extrand 36 and Gao and McCarthy 37 claimed that the static contact angle of droplets relies on the condition of vicinity of the triple contact line instead of the whole wetted area. Several studies demonstrated that the length scale of chemical heterogeneity, deposition position, and wetting pattern shape mainly influence the static condition of the droplet on chemically heterogeneous surfaces. [38][39][40][41] However, more effort is still needed for better understanding the droplet behavior on chemically heterogeneous surfaces. The validity of the power law for droplet spreading on chemically heterogeneous surfaces needs to be examined, and the influence of surface properties, especially near the triple contact line, on the spreading of droplet needs to be understood.
The molecular dynamics (MD) method can be applied to investigate the kinetics of nanodroplet spreading, which is otherwise difficult to handle in experiments. In the present study, the spreading of water nanodroplets on chemically heterogeneous surfaces with wetting patterns of alternating hydrophilic and hydrophobic regions was investigated through MD simulations. The spreading radius and dynamic contact angle were compared to previous experimental and theoretical studies and were analyzed to understand the influence of surface properties on nanodroplet spreading, and the influence of the length scale of chemical heterogeneity, wetting pattern, and fraction of hydrophobic surface is discussed.

II. SIMULATION METHOD
A. Molecular dynamics model The interactions between water molecules were described by the extended simple point charge (SPC/E) model, 42 in which the interactions between two water molecules are calculated as follows: The first term in the right side refers to the Lennard-Jones potential for oxygen-oxygen interactions with ε OO = 0.6502 kJ/mol and σ OO = 3.166 Å, and the second term in the right side refers to the Coulomb potential. r OO is the distance between oxygen atoms of two water molecules; rij is the distance between atoms i and j; qi and qj are the charges of atoms i and j, respectively; and ε 0 is the vacuum permittivity. The SHAKE algorithm 43 was used to keep the water molecules rigid, and the long-range Coulombic force was calculated using the particle-particle-particle-mesh technique. 44 The embedded atom model (EAM) potential 45 was used to describe the interaction between gold atoms, and the standard 12-6 Lennard-Jones potential was used for the water-Au interaction, in which the interaction constant ε O-Au is modified to model a hydrophilic or hydrophobic substrate and the distance parameter σ O-Au is fixed. A cutoff distance of 0.95 nm was applied for Au-water and waterwater non-bonded interactions. All simulations were performed using the LAMMPS package. The time step in all simulations was set to 1 fs and the temperature was set to 298 K. The simulation continues for each condition until the droplet reaches a stable condition, and the total simulation time is 6-12 ns depending on the wetting pattern.
To simulate droplet spreading on a solid surface, periodic boundary conditions were applied in all directions of the simulation box with a region size of 244.8 × 244.8 × 170.0 Å 3 . Water droplets with 5000-15 000 water molecules and diameters of 6.8-10 nm and Au substrates with 5 layers of lattices (lattice constant of 4.08 nm) were applied in the simulation. The droplets and Au substrates were simulated to equilibrium before placing the droplets on the surface of the substrates. The surface of the substrates is the {100} surface of Au. The deposition position of the droplet is chosen to be the corner of the island patch on the surfaces for all the simulations in this study in order to avoid the influence of deposition.

B. Wetting patterns
The wetting patterns on chemically heterogeneous surfaces are shown in Fig. 1 hydrophilic region, the constant is ε O-Au,1 = 5.260 kJ/mol with a corresponding static contact angle of approximately 25 ○ , whereas for the hydrophobic region, the constant is ε O-Au,2 = 0.482 kJ/mol with a corresponding static contact angle of approximately 140 ○ . Although a contact angle of 140 ○ is hard to obtain in experiment, it is still used in the simulation in order to have a large difference between the two regions of a wetting pattern, so the effect of wetting pattern on droplet spreading can be prominently shown. The geometric parameters for 15 surfaces with different patterns are listed in Table I, where the Type column shows the pattern types as listed in Fig. 1. w, l, and d are the length parameters defined in Fig. 1, and a 0 is the lattice constant. f is the area fraction of the hydrophobic region on the surface and is obtained by dividing the hydrophobic atom number by the total atom number on the surface.

C. Equivalent spreading radius and dynamic contact angle
The wetting kinetics of water droplet spreading can be characterized by the equivalent spreading radius R b and dynamic contact angle θ d . The following steps are applied to calculate the equivalent spreading radius R b and dynamic contact angle θ d . First, the liquid droplet is divided into several horizontal layers. The layer thickness is chosen such as to maximize the layer number, and the water molecule number in each layer is sufficient to provide a uniform density. Then, the density of water molecules for each layer is calculated as a function of spatial coordinates, and a cutoff density is chosen to identify the region occupied by liquid molecules. The equivalent spreading radius R b is the equivalent radius of the layer that is 5 Å from the surface, and the area of the wetting region is equal to the circle area of the equivalent spreading radius. The layer thickness and cutoff density need to be chosen so that their changes will not influence R b and θ d . In this study, the layer thickness is set to be 0.4 Å, and the cutoff density is half of the density of the bulk.
To obtain the dynamic contact angle of the droplet, the equivalent radius of each layer at different heights of the droplet is plotted in Fig. 2 and a circular fitting is then applied to obtain the curve of the droplet. Note that when fitting the curve, the equivalent radius below a height of 5 Å is not used to avoid the abnormal change caused by the solid-liquid interaction. Although the droplet deforms greatly in the vicinity of the triple contact line, the droplet shape above 5 Å from the surface agrees well with a spherical cap. Extending the spherical cap to the surface, the angle between the cap and the surface is obtained as the dynamic contact angle θ d of the droplet.

III. RESULTS AND DISCUSSION
A. Effect of droplet size Spreading of water droplets with different sizes on 15 surfaces was evaluated. The numbers of molecules are 5000, 9000, and 15 000 for droplets with diameters R 0 of 68 Å, 84 Å, and 98 Å, respectively.  Table I.
The result shows that the dimensionless equivalent spreading radius and dynamic contact angle change when the number of molecules increases from 5000 to 9000, whereas they are almost identical for droplets with 9000 and 15 000 molecules. This finding indicates that the droplet size has no remarkable influence on the spreading process on No. 2 surface when the droplet molecules are more than 9000, which is consistent with other surfaces in this study. This result also agrees with several other investigations that wetting characteristics are not highly dependent on the length scale of alternating lyophilic and lyophobic regions for droplets of several micrometers if they cover multiple regions. [27][28][29] In this study, droplets with 9000 molecules are large enough and are used in the following simulations. Figure 4 shows the equivalent spreading radius changes with time on No. 2, 5, 6, and 7 surfaces. Note that the area fraction of the hydrophobic region sequentially increases in the order of No. 2, 5, 7, and 6 surfaces. The equivalent spreading radius and dynamic contact angle are dependent on the area fraction of the hydrophobic region, and the spreading index is almost the same before the droplet reaches a stable condition. The equivalent spreading radius decreases as the area fraction of hydrophobic region increases, and the dynamic contact angle increases as the area fraction of the hydrophobic region increases. The spreading process becomes stable earlier on surfaces with a higher area fraction of the hydrophobic region. The spreading index is around 1/2 at the early stage and decreases to around 1/7 at the late stage. The spreading index at the early stage is identical to that from the experimental study of Stapelbroek et al. 46 for the early stage of droplet spreading on either physically or chemically heterogeneous surfaces, and the spreading index at the late stage agrees with that predicted by MKT, which is valid for the late stage of droplet spreading. [20][21][22][23]47 The results indicate that the spreading index at the early or late stages is independent of the chemical properties of surface, which agrees with previous experimental and theoretical study. 40 The area fractions for samples 2, 5, 7, and 6 are 1/9, 4/9, 5/9, and 8/9, respectively. In Fig. 4(b), the solid lines correspond to MD simulation and the dashed lines correspond to MKT calculation. tend to penetrate into the hydrophilic region. When the area fraction is 1/9 or 4/9, the hydrophobic regions are individual islands separated by the hydrophilic regions, and droplet spreading is mainly dominated by the hydrophilic regions. When the area fraction increases to 5/9, the difference occurs in the late stage of spreading. The droplet on the surface of type B stops spreading earlier than that on the surface of types C and F, with a smaller equivalent spreading radius and lager contact angle. There is a network of hydrophilic regions with isolated hydrophobic regions on surface Nos. 9 and 15, which can promote the spreading process, while a wetting pattern on surface No. 7 in the contrary can hinder the spreading process. This leads to the consequence that the droplet on surface No. 7 stops its spreading process earlier than that on the surface Nos. 9 and 15, with a smaller equivalent spreading radius and larger contact angle.

C. Effect of wetting pattern
When the area fraction increases to 8/9, the difference mainly occurs at a time of 100 ps. The droplet on the surfaces No. 6 and No. 10 spreads faster than that on surface No. 14 at that time. The energy dissipation is proportional to the surface tension, the perimeter of the moving contact line, and the moving velocity of contact line 13. The energy dissipation should be equal to the energy input during droplet spreading. Therefore, with an identical energy input, a longer perimeter of the moving contact line, with a larger dissipation force at the contact line, leads to a lower moving velocity of the contact line. In the early stage, the length of the contact line is mainly influenced by the curvature of the liquid bridge close to the substrate which is identical for all the three substrates, so the spreading processes are the same. In the next stage, the length of the contact line increases and is significantly affected by the local wetting pattern. The contact line is almost circular on surface No. 6 and is elliptical on surface No. 10 with close lengths, while it is more distorted and longer on surface No. 14 due to the triangular shape of the isolated hydrophilic regions. The longer perimeter of the contact line on surface No. 14 leads to a lower moving velocity and thus, a smaller equivalent spreading radius and a larger contact angle. In the late stage, the droplet well spreads on the substrate with the contact line across a larger area, and the shape of the contact line is mainly affected by the overall properties of the substrate. Since the three surfaces have the same area fraction and network of the hydrophilic region, the contact lines on these surfaces are similar, leading to similar changes with time in the equivalent radius or dynamic contact angle during the late stage of spreading.
From the analysis, it is seen that the shape of the network has a relatively important influence on the spreading process when the hydrophobic regions are connect to a network on the surface of the substrate.

D. Effect of pattern period
In this study, the pattern period d was chosen to be smaller than the droplet diameter so that the droplet can cover several pattern periods during the spreading process. Figure 7 shows the equivalent spreading radius and dynamic contact angle of the droplet as functions of spreading time for droplets on surface Nos. 1-4 with the same area fraction of 1/9 but different pattern periods. The equivalent spreading radius and dynamic contact angle show almost the same trend regardless of the difference in the pattern period, indicating that the pattern period has a small contribution to the spreading process when the droplet covers several periods. The results agree with previous studies in which the wetting characteristics are slightly dependent on the absolute size of the pattern period when a droplet covers multiple pattern periods. [27][28][29]

E. Static contact angle
The static contact angle is the contact angle of the droplet at equilibrium. It usually takes 1200 ps or more to reach the equilibrium state, in which the contact angle and equivalent spreading radius no longer change. Cassie's equation is widely used to describe the static contact angle of the droplet on chemically heterogeneous surfaces, 18 where θ 0 1 and θ 0 2 are the static contact angles of the droplet on two areas with different chemical properties on a surface and f 1 and f 2 are the fractions of the two areas on the surface. The validity of Cassie's equation has been questioned and investigated for many years. To consider more detailed molecular interactions including the additivity of molecular polarizabilities, dipoles, and charges, the modified Cassie's equation is proposed as 48 As seen in Fig. 8(a), the static contact angle calculated by MD simulations deviates from the predictions of Cassie's equation and modified Cassie's equation, which may be attributed to the distorted shape of the droplet bottom on the substrates. Figure 9 shows the visualization performed using visual molecular dynamics (VMD), 50 which shows the shape of droplets on substrates 5 and 7. As shown, the water droplet tends to spread in the hydrophilic area but is hindered in the hydrophobic area. As a result, the droplet deviates from a spherical cap and the triple contact line is distorted, which is different from the usual droplets. The regions where the triple contact line stays could be important to the static contact angle. The triple contact line alternately passes the hydrophobic and hydrophilic regions of the substrate. The ratio of the total length of the hydrophobic sections to the whole length of the triple line is hereby defined as the hydrophobic length fraction f l . The hydrophobic length fraction f l measures the surface properties near the triple contact line more properly than the area fraction of the hydrophobic region f on the substrate. Figure 8(b) shows the static contact angle as a function of the hydrophobic length fraction f l . Compared to Fig. 8(a), the difference in the results between the MD simulation and Cassie's equation or the modified Cassie's equation becomes smaller, indicating that the regions that the triple line passes are more important to the static contact angle than the whole substrate. Besides, in both Figs. 8(a) and 8(b) the modified Cassie's equation predicts the results of MD simulations better than Cassie's equation because it has considered the molecular interactions with different regions on the substrate.
The results agree with the experimental investigation of Extrand 36 and Gao and McCarthy 37 and the theory proposed by Xu and Wang. 49 They claimed that the static contact angle of the droplet does not solely rely on the average properties of the entire substrate. The region in the vicinity of the contact line is crucial for the static contact angle. Therefore, the interfacial energy in the vicinity of the distorted contact line, instead of the average interfacial energy of the surface, should be adopted to explain the behavior of the static contact angle. Further investigations are still needed to better understand the influence of the contact line on droplet spreading.

F. Dynamic contact angle
The dynamic contact angle of the droplet changes during the entire spreading process as shown in Fig. 4(b) and can be predicted by MKT. According to this theory, the macroscopic behavior of the contact line depends on the overall statistics of the individual molecular displacements that occur within the three-phase zone. 51 The spreading velocity v can be described as follows: where κ 0 is the equilibrium frequency; λ is the average distance between the centers of adsorption sites; γ LV is the surface tension of the liquid; θ 0 and θ d are the static and dynamic contact angles, respectively; n is the number of adsorption sites per unit area; k B is Boltzmann's constant; and T is the temperature. For a small droplet, gravity can be neglected and the shape of the droplet can be assumed to be a spherical cap. 47,52 The geometry of the spherical cap provides the following relation: where r is the base radius of the droplet, V is the droplet volume, and θ d is the dynamic contact angle. Assuming no evaporation occurs, the volume is constant and hence, Equations (5) and (7) can be linked to a partial differential equation with two parameters: a = 2κ 0 λ and b = γ LV /2nk B T. The parameter a can be regarded as the moving speed of liquid molecules at the contact line, and the parameter b represents the liquid and solid properties. The equations can be solved numerically to obtain the relation between the dynamic contact angle and spreading time.
In this study, the dynamic contact angle data of all samples are used to fit the two parameters. Figure 4(b) shows the fitting curves using MKT. Given that the site density n and other parameters in b are independent of the interactions between the solid and the liquid, b should be the same in all the simulation cases. The differential equations were solved using a fourth-order Runge-Kutta algorithm. The difference between the simulation data and theoretical curve was minimized to obtain the best fit. In this study, the best fit gave b = 1.74, and Table II shows the best fit of a. The parameter a varies highly nonlinearly with the area fraction of substrates, suggesting the pattern type (shown in the Type column) has an important influence on the spreading process.

ARTICLE
scitation.org/journal/adv  Figure 4(b) shows the dynamic contact angle curves obtained by MD simulations and MKT. MKT well predicts the dynamic contact angle despite some small deviations. The assumption that the droplet is a spherical cap is probably the main cause of deviations.
In Table II, the parameter a slightly varies as f increases. The molecular interaction between the liquid (water) and solid (Au) atoms is weaker in the hydrophobic region than that in the hydrophilic region. Increasing the hydrophobic area fraction f leads to an increase in the hydrophobic region area and thus, a decrease in the solid-liquid interaction. As a result, the liquid molecules are less bound to the solid sites and it is easier for them to jump with a high jump frequency between sites on the solid surface. 47 In this study, the average distance between the centers of solid sites λ can be assumed to be constant due to the same atomic distribution in the substrate. Therefore, according to the definition a = 2κ 0 λ, the increase in the jump frequency κ 0 of liquid molecules on a solid surface leads to the increase of parameter a. This well explains the increase of a with the increase of f, whereas even with the same hydrophobic area fraction f (e.g., surface Nos. 5 and 11, or Nos. 7 and 9), different wetting patterns (shown in the Type column) also influence the molecular interaction at the contact line and cause the difference in the molecular interaction and thus, the differences in the jump frequency κ 0 and parameter a.
Note that a has a remarkable increase as f changes from 5/9 to 8/9. This shows that the moving speed of liquid molecules increases quickly as the surface of the substrate approaches to a pure hydrophobic one. In other words, for a purely hydrophobic surface, adding a small fraction of hydrophilic area will cause a remarkable decrease in the moving speed of liquid molecules on the surface. To qualitatively explain this phenomenon, further investigations are need on the moving mechanism of the contact line on a chemically heterogeneous surface.

IV. CONCLUSIONS
In this study, the spreading of water droplets on chemically heterogeneous substrates is investigated through MD simulations. The spreading process of droplets is analyzed to obtain the equivalent spreading radius and dynamic contact angle, and the following conclusions are obtained in this study.
1. The diameter of droplet does not affect the spreading process on chemically heterogeneous surfaces when the droplet is large enough to cover multiple (>3-4) periods of patterns on the surfaces.
2. The spreading index at the early or late stage is independent of the chemical distribution of surfaces in this study, which agrees with the previous experimental study of Stapelbroek et al. 46 and MKT. The shape of the wetting pattern and the hydrophobic area fraction have a relatively minor effect on the spreading index in the power law. 3. The advancing contact line of the droplet is hindered by the hydrophobic region and tends to penetrate into the hydrophilic region. The distribution of hydrophilic and hydrophobic regions and whether these regions are connected strongly influence the shape of the droplet. 4. The static contact angle as a function of area fraction of the hydrophobic region remarkably deviates from the original and modified Cassie's equation, while it agrees better with the modified Cassie's equation as a function of hydrophobic length fraction, indicating the regions near the triple contact line are more important to the static contact angle than the whole substrate. 5. The dynamic contact angle can be predicted by MKT. The increasing trend of a is attributed to the decreasing solid-liquid interaction. As the liquid molecules are less bound to the solid sites, they jump more frequently and thus, their moving speed on the surface increases. Besides, the pattern also influences the movement of molecules on the surface. However, the specific explanation for the change of the dynamic parameter a with area fraction and wetting pattern is still missing. Further studies are needed to understand the molecule jumping mechanism on different wetting patterns.