Requirements for the Formation and Shape of Microscopic Precursors

Apr 15, 2016 - While the mass transport mechanisms and dynamics of molecular precursors, ultrathin films that precede spreading droplets, are well und...
0 downloads 11 Views 3MB Size
Download PDF
Article pubs.acs.org/Langmuir

Requirements for the Formation and Shape of Microscopic Precursors in Droplet Spreading Rolf E. Isele-Holder* and Ahmed E. Ismail*,† Aachener Verfahrenstechnik: Molecular Simulations and Transformations and AICES Graduate School, RWTH Aachen University, Schinkelstraße 2, 52062 Aachen, Germany ABSTRACT: While the mass transport mechanisms and dynamics of molecular precursors, ultrathin films that precede spreading droplets, are well understood, the requirements for their formation and the reasons for the occurrence of different precursor shapes remain unclear. In this work, we study these requirements using molecular dynamics simulations of spreading droplets and extensive free energy computations. For a simple simulation model, we demonstrate that with growing droplet−substrate attraction, spreading passes succesively through regimes with no precursor, a monolayer precursor, and a continuously growing precursor. We show that the onset of layer formation and the changes in the precursor shape correlate with the free energy of layer formation. On the basis of these findings, we show that a positive spreading coefficient is sufficient but not necessary for precursor formation. Microscopic films (see ref 4 for a comprehensive review), in contrast, are characterized by a thickness of only a few molecules, making the continuum view invalid, as confirmed by the inability of continuum approaches to model microscopic precursors.13,14 Continuum approaches cannot successfully capture the apparently universally valid diffusive spreading dynamics of these microscopic films.4 While the mass transport mechanisms and spreading dynamics of microscopic precursors are well understood from numerous experimental, theoretical, and simulation studies,15−19 two fundamental questions remain. First, what are the requirements for precursor formation? Second, if a molecular precursor evolves, what controls the shape of the observed precursor? The purpose of this study is to address these questions. Because the onset of precursor formation has often been observed to coincide with the complete wetting transition,19−23 where the spreading coefficient S becomes positive, it is often argued that S > 0 is a requirement for the formation of microscopic precursors. This view, however, is challenged by observations of molecular precursors in the partial wetting regime.24,25 Moreover, there is no theoretical justification for why the spreading coefficient S, a macroscopic quantity, should control the onset of a microscopic film. In fact, the only theoretical model that comments on the onset of formation of microscopic precursors is that of Burlatsky et al.15,26 This model predicts that precursors form if a microscopic spreading coefficient s, which is equal to the negative of the energy

1. INTRODUCTION Droplets of even the simplest liquids can exhibit a wealth of phenomena close to the contact line.1,2 One such phenomenon is the occurrence of thin precursors3,4 that can appear in various shapes, such as single layers of constant thickness; continuously growing precursors; thick, pancakelike films; and multiple distinct layers, also known as terraced precursors.5 Precursor formation has a number of applications, including the spreading of cells on substrates,6 the mobilization of water droplets on hydrophobic surfaces with thermocapillary actuation,7 droplet transport on nanowires,8 and use in nanotechnology, where precursors can either serve a desired purpose or be an undesired contaminant. Because spreading dynamics are controlled by the microscopic region at the droplet contact line,2 precursor formation is relevant not only for specific technical applications but also in a wider context of spreading, particularly for scenarios in which the macroscopic spreading behavior is controlled by the presence (or absence) of a precursor film. Models of precursor films are typically subdivided into mesoscopic and microscopic films. The latter are also sometimes termed diffusive or molecular precursors.2 Mesoscopic films are sufficiently thick, typically between 30 Å and 1 μm, such that the continuum view is applicable.1 Using the concept of the disjoining pressure,9,10 it is possible to describe the dynamics and shape of these mesoscopic wetting films (cf. ref 1 for a detailed review). Exploiting the fact that the spreading coefficient S is the integral of the disjoining pressure, one can show that a positive spreading coefficient is a requirement for mesoscopic precursor formation,1,2,4 although spreading can occur in some cases for a slightly negative value of the spreading coefficient.11,12 © XXXX American Chemical Society

Received: March 1, 2016 Revised: April 14, 2016

A

DOI: 10.1021/acs.langmuir.6b00807 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir required to transport a fluid particle from the bulk to the tip of the film, becomes positive. Thus, while the first question mentioned above has been previously addressed, with the macroscopic spreading coefficient S identified as a crucial factor, the exact role of S, as well as the actual requirements for precursor formation, remains unclear, and the second question mentioned above has been hitherto completely untouched in the scientific literature. In this study, we use molecular dynamics (MD) simulations to provide answers for the two open questions. Section 2 describes the setup of the spreading simulations and the extensive free energy computations required to understand the onset of precursor formation. Results in section 3 show that the onset of precursor formation does not have to correlate with the wetting transition but that both the onset and the precursor shape correlate with the free energy of layer formation. On the basis of these observations, we show why the onset of precursor formation can correlate with the wetting transition in section 4. We summarize our findings in section 5.

parameters were selected so that interactions are softer at shorter separations compared to an LJ potential, as shown in Figure 1. As a

Figure 1. LJ and modified Buckingham potentials (left) for the liquid− liquid interactions. The Buckingham potential is less repulsive at shorter particle distances. Equilibrated slab system of the fluid (right) at T = 0.6ε/kB. The condensed phase is liquid, and the vapor phase is almost empty.

2. METHODS We simulate the spreading of a nonvolatile liquid on a crystalline substrate. The solid−liquid interaction energy is used as a control parameter to sample different precursor regimes. Results from the spreading simulations are compared to the spreading coefficient obtained from the dry-surface simulation method (DSSM)27 and the free energy for layer formation obtained from the two-phase thermodynamic (2PT) method.28−30 2.1. Simulation Model. Our work requires a model that is both nonvolatile so that effects caused by evaporation and condensation can be neglected and simple enough computationally so that the employed free energy computation methods provide highly accurate results. Low volatility has been achieved in previous MD studies by using chain molecules18,20,21,25 and metallic liquids.31 The three-body interactions in the EAM potential for metals and the difficulties in properly sampling the phase space of the chain molecules, however, make these models too demanding for the study at hand. The simplest computational model is the single-center Lennard-Jones (LJ) fluid. This model, however, is volatile above its melting temperature32 and thus unsuitable. We overcome these shortcoming of the LJ, chain, and metallic models by introducing a new potential for nonvolatile, simple liquids. We emphasize that the purpose of this study is not to examine properties of the model defined by this new potential. The motivation for defining this potential is to have a model that is both nonvolatile and sufficiently simple such that the spreading coefficient and the free energy of layer formation can be accessed. The units in our system are defined by a reference energy ε, mass m, distance σ, and derived time unit τ = (mσ2/ε)1/2. The interaction between particles i and j is described by a modified Buckingham potential33 U=

⎡ rij ⎢ 6 e λ(1 − σij ) − 6 1 − λ ⎢⎣ λij ij εij

⎛ σij ⎞6 ⎤ D ⎜⎜ ⎟⎟ ⎥ + ij 12 ⎥ ⎝ rij ⎠ ⎦ rij

result, the melting point of the liquid is reduced and simulations can be performed at lower temperatures, which results in reduced fluid volatility. Effects due to evaporation and condensation through the gas phase can therefore be neglected. A lack of volatility in three dimensions can still allow for volatility along the substrate in two dimensions. This aspect is briefly discussed in section 3.1. 2.2. Spreading Simulations and Analysis. Because simulations of cylindrical droplets show the same behavior observed for spherical droplets while being computationally less demanding for a given droplet radius,21 we used cylindrical droplets in our simulations, with the axis of the cylinder pointing in the y-direction. The substrate’s surface normal points in the z-direction, so that spreading occurs in the x-direction. Simulations were performed with periodic boundary conditions in all directions. The substrate was a fcc crystal with a unit length of 1.5874σ and the (111) vector normal to the surface pointing in the z-direction. The substrate was composed of eight layers, with the lowest layer held rigid in all simulations. For the spreading simulations, the substrate extended 600.747σ and 33.6739σ in the xand y-directions, respectively, and contained 142772 particles. The liquid droplet was composed of 180000 beads and was created using PACKMOL.34 The substrate and droplet were equilibrated separately for 50τ and 1000τ at a reduced temperature of T = 0.6ε/kB using a Nosé−Hoover thermostat35 with a damping factor τT = 0.5τ. The equilibrated droplet was then positioned 2.5σ above the equilibrated substrate. The first five substrate layers above the rigid layer were coupled to a Langevin thermostat36 with a damping factor τD = 0.3τ. No thermostat was used for the rest of the simulation because it provides the greatest similarity to experiment.18 During the equilibration and production, the particle positions and velocities were updated using a two-stage rRESPA integrator,37 with interactions within a distance of 3.0σ shifted to zero starting at 2.5σ and computed every 0.005τ. Long-range dispersion interactions were computed using a PPPM 38,39 solver with ik differentiation, interpolation order P = 5, grid spacing h ≈ 0.75σ, Ewald parameter β = 1.0/σ, and real-space cutoff rc = 3.0σ and were evaluated every 0.02τ. The repulsive terms were truncated at rc. The spreading simulations are analyzed using our method presented in ref 40. Liquid layers are defined from the natural layering in fluids that occurs close to solid surfaces. The position of the droplet’s liquid−vapor interface is determined for each layer from a twodimensional histogram of the liquid density. A circular arc is fit to this interface, with the interface close to the solid surface being neglected. The horizontal distance from the interface position of layer i to the circular fit is the measured precursor width li. In addition, we measure

(1)

where rij is the distance between i and j and λij controls the shape of the curve. If Dij = 0, the original Buckingham potential, where εij and σij are the depth and location of the potential well, respectively, is recovered. The last term avoids the unphysical maximum for small rij in the original Buckingham potential. For the substrate, we use εss = 5ε, λss = 12, σss = 1.122σ, and D = 0.05εσ12. For the liquid, we use εll = ε, λll = 9, σll = 1.1σ, and Dss = 0.03εσ12. For the solid−liquid interaction, we use λsl = 12, σsl = 1.1σ, and Dsl = 0.03εσ12, while εsl is varied in the simulations. We use ms = 2m and ml = m for the mass of the solid and liquid particles, respectively. The substrate−particle interaction was chosen so that a stable solid phase is formed at the simulated temperature, while the liquid−liquid B

DOI: 10.1021/acs.langmuir.6b00807 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir the effective surface density of the liquid to characterize the shape of the precursor. See ref 40 for further details. 2.3. Computation of the Spreading Coefficient. The spreading coefficient S describes the wettability of a liquid on a substrate. S is defined via S = γs − γsl − γl

molecular simulations, it has not been previously computed. Here we approximate the free energy of forming precursors by computing the change in free energy required to deposit a layer of liquid from a bulk reservoir on a bare solid substrate. In this approach, the deposited layer is a substitute for the precursor and the bulk reservoir is representative of the droplet. Free energy changes are computed using the two-phase thermodynamic method28−30 using the implementation of Lin et al.43 This method splits the densities of states, which can be computed from the velocity autocorrelation function, into solid and gas contributions (hence the name “two-phase”). For each contribution, thermodynamic quantities can be accessed directly using either a harmonic oscillator for the solid or the ideal gas approximation for the gas. Its accuracy has been demonstrated for various test cases and was found to be especially striking for single-center Lennard-Jones particles.28 Given that the model studied here resembles single-center Lennard-Jones particles, we expect the method to provide accurate results. The free energy of the liquid state is computed from bulk systems of fluid; 1024 particles are randomly positioned in a cubic box of length 12σ. The system is then equilibrated in the NPT ensemble at zero pressure and temperature T = 0.6ε/kB for 250τ using Nosé−Hoover barostats and thermostats with τT = 0.5τ and τp = 5τ. Afterward, the simulation was run for 500τ in the NVT ensemble. A Verlet integrator with a time step of 0.005τ was used. Every 50τ, the particle velocities for computing the velocity autocorrelation functions for the 2PT method were stored with a frequency 0.001τ for a period of 10τ. The free energy of the fluid-laden substrate was computed from a simulation cell with an eight-layer solid slab consisting of 3736 atoms covered with liquid in a 23.3σ × 22.5σ × 46σ box. The simulation setup is depicted in Figure 3. Simulations were performed for different

(2)

as a balance of the solid, solid−liquid, and liquid surface energies (γs, γsl, and γl, respectively). While γl can be computed easily from MD simulations, γs and γsl cannot be readily accessed. S < 0 in the partial wetting regime, where the equilibrium shape of a droplet is a spherical or cylindrical cap, and S > 0 in the complete wetting regime, where the droplet spreads out to form a thin film. S = 0 marks the transition between the wetting regimes. γl is computed from slab geometry simulations of the liquids with the surface of the slab being perpendicular to the z-direction with 16000 fluid particles in a 23.3σ × 22.5σ × 60.0σ box similar to the setup in Figure 1. Simulations were run with a Nosé−Hoover thermostat with damping factor τd = 0.5τ. After equilibration for 500τ, surface tensions were computed from a subsequent run over 20000τ using the approach of Kirkwood and Buff.41,42 The remaining simulation settings were similar to those used in the spreading simulations. The difference between γs and γsl was determined using the DSSM,27 in which the potential between the substrate and the liquid is gradually switched off. This is realized in our simulations by scaling εsl and Dsl. Simulations were executed with 16000 fluid particles on a solid substrate with 3736 particles in a 23.3σ × 22.5σ × 86.0σ box, with an eight-layer solid slab, a liquid slab, and a vapor slab, as shown in the left image of Figure 2. The path was sampled for values of εsl

Figure 3. Solid substrates (pink) covered with thin layers of liquid (cyan) in simulations to compute the free energy of layer formation with the 2PT method: (left) low particle density and (right) maximal particle density used in this study.

Figure 2. Simulation setup for the computation of Δγ with fluid particles (cyan) on a solid substrate (pink). Box edges are depicted as black lines. Switching off the interactions between the solid and the liquid particles is equivalent to pulling the solid and liquid particles apart and thereby removing the solid−liquid interface and generating a solid and a liquid surface.

values of the density of liquid particles ρsurf, where ρsurf = nsurf/Asurf is the ratio of the number of liquid particles to the surface area. The number of liquid particles nsurf was varied between 0 and 1200. Simulations are performed with εsl = 1.3ε. A Nosé−Hoover thermostat with damping factor τT = 0.5τ and a Verlet integrator with a time step of 0.005τ was used to equilibrate the system in the NVT ensemble for 250τ. Afterward, simulations were run for 1000τ, with particle velocities stored for 10τ every 50τ with a frequency of 0.001τ to compute the velocity autocorrelation functions. Velocity autocorrelation functions and the derived thermodynamic quantities were computed separately for each set of particle velocities. Mean values and statistical uncertainties are computed from the results obtained for each separate velocity autocorrelation function. The internal energy is taken directly from the simulations, while the entropy is obtained from the 2PT method. The specific volume of the liquid on the covered substrate, which is required as an input for the 2PT method, was assumed to be equal to that of the pure liquid. Planck’s constant in reduced units was chosen as h = 0.185 σ(mε)1/2, in agreement with ref 28. The change in free energy with varying substrate energy compared to the reference energy of εsl = 1.3ε was computed using the thermodynamic integration scheme already used for computing the spreading coefficient. Differences compared to the description above are that we used a Langevin thermostat instead of the Nosé−Hoover thermostat to enhance the sampling, that the equilibration was done for 125τ and the production for 250τ, that a Verlet integrator with a

between 0.06ε and 2.0ε with spacing 0.02ε. For εsl ≥ ε, Dsl = 0.03εσ12, in agreement with spreading simulations. For ε < εsl, Dsl scales linearly to zero with the same factor as εsl. For εsl = 0, the contribution to the change in free energy is taken as zero.27 At each εsl, the simulations are equilibrated for 500τ. Production runs are performed for 2000τ, during which the change in potential energy with respect to the change in εsl and Dsl is evaluated every 2.5τ. The temperature was controlled using a Nosé−Hoover thermostat with τT = 0.5τ. Other settings were as described above. Using this method, the solid−liquid interface is removed and separate solid and liquid interfaces are created, as illustrated in Figure 2. The result from this method is thus

Δγ = γl + γs − γsl

(3)

The spreading coefficient can then be determined as

S = Δγ − 2γl

(4)

2.4. Computation of the Free Energy for Layer Formation. Given that processes occur spontaneously if the associated change in free energy is negative, it is logical to assume that precursors will form if the change in the free energy of forming a precursor is negative: the formation of layers on the substrate should be energetically favorable. Even though this change in free energy is in principle accessible from C

DOI: 10.1021/acs.langmuir.6b00807 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir time step of 0.005τ was used, and that the derivative of the potential was evaluated every 0.625τ.

The solid−liquid interaction energy for the onset of the multilayer precursor is 1.6ε < εsl,mult < 1.8ε, as can be seen from panels c and d of Figure 5. Because of the noise in the measured li, and because precursor growth is very slow close to the onset of single-layer precursor formation, the solid−liquid interaction energy εsl,single at which a single layer starts to separate from the droplet cannot easily be determined from the raw data of li in panels a and b of Figure 5. The strong noise in the data stems from the uncertainties of the circular fit to the droplet surface. By eliminating of the cause for the noise in the data, an accurate value for εsl,single is determined. Uncertainties in the fit of the circular arc affect the layer width li for all layers in the same way. A large fitted droplet radius leads to small values for the width of all layers, li, and vice versa. The noise in the data can therefore be reduced by subtracting a reference layer that is not in the immediate vicinity of the substrate. The uncertainties originating from the fit of the circular arc to the droplet shape are thereby subtracted to the same extent from all reported li. As can be seen from the the modified data in panels e and f of Figure 5, a single layer separates from the droplet at εsl = 1.3ε, whereas no layer separates at εsl = 1.25ε. The onset for a single-layer precursor is thus at 1.25ε < εsl,single < 1.3ε. Additional information is required to fully characterize the observed precursors. Figure 5 reveals that precursors separate, but for the multilayer plot, the question of whether the observed precursor is continuously growing or whether multiple distinct layers, the so-called terraced precursor, evolve is unanswered. The effective surface density of the droplet in Figure 6 reveals that the precursor at interaction energies εsl > εsl,mult is continuously growing and not terraced. The profile of the effective droplet height in Figure 6 also allows us to evaluate the importance of two-dimensional (2D) volatility of the liquid along the substrate surface. A liquid with high 2D volatility would show a very slow and smooth decay of the effective droplet height as a function of the spreading direction. For our system, however, the effective droplet height decays steeply to zero at the tip of the precursor, demonstrating that our model has low 2D volatility.

3. RESULTS After positioning the liquid droplet above the substrate and starting the simulations, the long-range dispersion interactions attract the liquid to the substrate and the droplet starts to spread, as shown in Figure 4.

Figure 4. Spreading of the liquid droplet (cyan) on a substrate (pink) with εsl = 1.6ε at times 0τ, 4200τ, and 8200τ (from top to bottom, respectively).

3.1. Precursor Regimes and Transitions. The first step is to define the observed spreading regimes and the transition points. As shown in Figure 5a, the layers li do not separate from the main part of the droplet for weak solid−liquid attractions. No precursor forms in this regime. At an increased interaction energy εsl = 1.6ε (cf. Figure 5c), only the first layer l1 separates from the droplet; a single-layer molecular precursor forms at this intermediate interaction energy. Further increasing the interaction energy leads to multiple layers li separating from the droplet and growing over time (cf. Figure 5d).

Figure 5. Evolution of the first four layers of the atomistic fluid for different solid−liquid interactions εsl. (a) No precursor forms for low substrate energies. (b) The onset of a monolayer precursor cannot be detected from the raw data of li. (c) For intermediate energies at εsl = 1.6ε, a single layer separates. (d) Multiple layers separate for substrate energies εsl ≥ 1.8ε. (e and f) Evolution of the first four layers minus a reference layer (to reduce the noise in the data) reveals the onset of single-layer formation at 1.25ε < εsl,single < 1.3ε. D

DOI: 10.1021/acs.langmuir.6b00807 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

Figure 6. Effective droplet thicknesses for different εsl. Black points denote the measured droplet thicknesses. Red lines are circular arcs fitted to the droplet shape. Dashed lines indicate the positions of the substrate. Top panels show data for the entire droplet with equal scaling in x- and zdirections. Bottom panels show selected regions with alternate scaling. The snapshots correspond to the final snapshot from each simulation. No precursor forms at the lowest substrate energy; a single layer forms for εsl = 1.6ε, and a continuously growing precursor develops for εsl = 2.0ε.

3.2. Precursor Requirements. The spreading coefficient is shown in Figure 7 as a function of the solid−liquid interaction

Figure 8. Change in free energy of depositing layers of the atomistic fluid from a bulk reservoir on a bare substrate over surface density ρsurf and different values of εsl. The symbols are larger than the statistical uncertainties (2σ environment). Figure 7. Spreading coefficient as a function of the substrate energy. The region in which the onset of precursors occurs is shaded. The width of the red line for the spreading coefficient corresponds to the statistical uncertainty (2σ environment) from the simulations. Note that the onset of precursor formation does not correspond to the nonwetting−wetting transition at which S = 0.

for εsl = 1.8ε at ρsurf ≈ 1.5/σ2. The free energy thus predicts that formation of thicker layers instead of a single-layer precursor is advantageous, which is in agreement with the transition from a single-layer precursor to a continuously growing precursor in the spreading simulations.

energy εsl. The onset of precursor formation 1.25ε < εsl,single < 1.3ε does not coincide with the transition from the partial to the complete wetting regime at εsl ≈ 1.4ε. Note that the results for spreading coefficient S are corrupted because of crystallization of the droplet material in simulations for determining the spreading coefficient at εsl > 1.34ε. Results for S at εsl ≤ 1.34ε are, however, correct. The observation that the wetting transition and the onset of precursor formation do not coincide is thus not affected by the crystallization. Results for the free energy of layer formation from a bulk liquid on a bare substrate ΔA are depicted in Figure 8. ΔA shows multiple characteristic features that relate to the observations from the spreading simulations. (i) For interaction energies εsl ≤ 1.2ε, ΔA > 0, and thus, no precursors should form at low interaction energies, in agreement with the spreading simulations. (ii) For low and intermediate substrate energies, free energy ΔA shows a minimum that drops below zero for the first time at εsl ≈ 1.3ε, in agreement with the onset of a single-layer precursor at this εsl. (iii) This minimum is located at a surface density ρsurf ≈ 1.0/σ2, equivalent to a substrate covered by a fluid monolayer. (iv) While for εsl ≈ 1.6ε the minimum at ρsurf ≈ 1.0/σ2 is a global minimum in the direction of the surface density, ΔA drops below the minimum

4. DISCUSSION From the excellent agreement of the spreading simulations and the free energy computations, we conclude that the free energy of layer formation controls the onset as well as the shape of precursors, as one would naturally expect. The important question regarding the role of the spreading coefficient S, however, remains unexplained by these results. Even though the wetting transition and the onset of precursor formation do not coincide for the model studied here, analyzing our results resolves the role of the spreading coefficient for precursor formation. In particular, we clarify why the wetting transition and the onset of precursor formation can coincide. The key to understanding the role of the spreading coefficient is to understand the relationship between ΔA and S. As shown in Figure 9, increasing ρsurf corresponds to increasing the thickness of the liquid slab deposited on the substrate. In the limit ρsurf → ∞, the deposited layer is so thick that the substrate−liquid interface and the liquid−vapor interface do not interfere. The change in free energy of creating the system from a bulk liquid and a bare substrate thus corresponds to the removal of a solid surface and the creation of solid−liquid and vapor−liquid interfaces: −γs + γsl + γl = −S. Thus, for ρsurf → ∞, ΔA → −S. E

DOI: 10.1021/acs.langmuir.6b00807 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

onset of a molecular monolayer precursor in the spreading simulations matches the free energy of forming a single molecular layer becoming favorable. The onset of the multilayer precursor corresponds to the free energy of forming multiple molecular layers falling below the energy of forming a single layer. The observations show that ΔA controls the existence and the shape of microscopic precursors. Finally, it is shown that the free energy of layer formation for layers of infinite thickness is equivalent to the negative spreading coefficient. This relationship between ΔA and S reveals why the onset of microscopic precursor formation can, but does not have to, correlate with the wetting−nonwetting transition, that a precursor has to exist in the complete wetting regime, and that precursors can exist in the partial wetting regime only if ΔA as a function of the surface density of the liquid ρsurf passes through a minimum. A positive spreading coefficient S is thus sufficient but not necessary for microscopic precursor formation.

Figure 9. Schematic explanation of the relationship between ΔA and S: pink for solid and cyan for fluid. The free energy changes as the number of particles on the surface increases. Starting from the bare substrate with the surface free energy γs over a substrate covered with a thin layer of fluid that has an effective surface energy γc to a substrate covered with a thick layer. For the thick layer, the upper and lower interfaces of the liquid do not interfere and thus correspond to a liquid and a solid−liquid interface with surface energies γl and γsl, respectively, and a bulklike phase between the interfaces. The change in free energy per unit area when moving from the left system to the right system is thus γl + γsl − γs = −S.

■

ΔA goes to −S for large values of ρsurf. Additionally, ΔA = 0 for ρsurf = 0; the limits of ΔA over ρsurf are thus fixed. It immediately follows that the change in the nonwetting−wetting transition has to match the onset of precursors if ΔA does not exhibit a minimum as one moves toward increasing ρsurf values at a fixed εsl. In this case, ΔA cannot drop below zero in the partial wetting regime and has to drop below zero in the complete wetting regime; the wetting transition and the onset of precursor formation thus have to match in this case. This finding also implies that a precursor in the partial wetting regime can exist only if ΔA has a minimum in the direction of ρsurf that drops below zero, and that a precursor has to exist in the complete wetting regime. Finally, our observations can be linked to findings for mesoscopic precursor films that can be described using the concept of disjoining pressure. Thick microscopic films with ρsurf → ∞ correspond to mesoscopic films. Using again that ΔA → −S for ρsurf → ∞, and that ΔA < 0 is required for precursor formation, our analysis predicts that very thick microscopic films, i.e., mesoscopic films, can exist only in the complete wetting regime where S > 0, which is in agreement with predictions using the concept of disjoining pressure for mesoscopic films. Lastly, our observations can be directly connected to mesoscopic precursor films via the disjoining pressure. Because ρsurf is directly related to the thickness of the fluid layer, the derivative dΔA/dρsurf is an effective microscopic analogue of the macroscopic disjoining pressure. Our computations are thus a natural extension of the theory of disjoining pressure for mesoscopic films to microscopic precursors.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Present Address †

A.E.I.: Department of Chemical and Biomedical Engineering, West Virginia University, P.O. Box 6102, Morgantown, WV 26506-6102. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Francisco Fontenele and André Bardow for helpful discussions and Brooks D. Rabideau and Tod A. Pascal for support with the 2PT method. Financial support for REI from the Deutsche Forschungsgemeinschaft (German Research Foundation) through Grant GSC 111 and access to computer resources provided by JARA-HPC are gratefully acknowledged.

■

REFERENCES

(1) De Gennes, P.-G. Wetting: statics and dynamics. Rev. Mod. Phys. 1985, 57, 827−863. (2) Bonn, D.; Eggers, J.; Indekeu, J.; Meunier, J.; Rolley, E. Wetting and spreading. Rev. Mod. Phys. 2009, 81, 739−805. (3) Heslot, F.; Fraysse, N.; Cazabat, A.-M. Molecular layering in the spreading of wetting liquid drops. Nature 1989, 338, 640−642. (4) Popescu, M. N.; Oshanin, G.; Dietrich, S.; Cazabat, A.-M. Precursor films in wetting phenomena. J. Phys.: Condens. Matter 2012, 24, 243102. (5) Cazabat, A.-M.; Fraysse, N.; Heslot, F.; Levinson, P.; Marsh, J.; Tiberg, F.; Valignat, M. P. Pancakes. Adv. Colloid Interface Sci. 1994, 48, 1−17. (6) Douezan, S.; Guevorkian, K.; Naouar, R.; Dufour, S.; Cuvelier, D.; Brochard-Wyart, F. Spreading dynamics and wetting transition of cellular aggregates. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 7315− 7320. (7) Zhao, Y.; Liu, F.; Chen, C.-H. Thermocapillary actuation of binary drops on solid surfaces. Appl. Phys. Lett. 2011, 99, 104101. (8) Huang, J. Y.; Lo, Y.-C.; Niu, J. J.; Kushima, A.; Qian, X.; Zhong, L.; Mao, S. X.; Li, J. Nanowire liquid pumps. Nat. Nanotechnol. 2013, 8, 277−281. (9) Derjaguin, B. V. Theory of the capillary condensation and other capillary phenomena taking into account the disjoining effect of longchain molecular liquid films. Zh. Fiz. Khim. 1940, 14, 137.

5. SUMMARY We report molecular dynamics simulations of the spreading of a simple, nonvolatile fluid on a crystalline substrate. With increasing solid−liquid interaction energy, the droplet passes through three precursor regimes: no precursor, a single-layer precursor, and a continuously growing precursor. The spreading coefficient S and the free energy of layer formation ΔA are computed from extensive free energy simulations. For the examined model, the nonwetting−wetting transition does not match the onset of precursor formation. In contrast, the free energy of layer formation correlates accurately with the observations from the spreading simulations. The F

DOI: 10.1021/acs.langmuir.6b00807 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir (10) Derjaguin, B. V.; Churaev, N. V. The current state of the theory of long-range surface forces. Colloids Surf. 1989, 41, 223−237. (11) Indekeu, J. O.; Ragil, K.; Bonn, D.; Broseta, D.; Meunier, J. Wetting of alkanes on water from a Cahn-type theory: effects of longrange forces. J. Stat. Phys. 1999, 95, 1009−1043. (12) Indekeu, J. O.; Ragil, K.; Broseta, D.; Bonn, D.; Meunier, J. New Approaches to Problems in Liquid State Theory; Springer: Berlin, 1999; pp 337−344. (13) De Gennes, P.-G.; Cazabat, A.-M. Étalement d’une goutte stratifiée incompressible. C. R. Acad. Sci. 1990, 310, 1601−1606. (14) Abraham, D. B.; Collet, P.; De Coninck, J.; Dunlop, F. Langevin dynamics of spreading and wetting. Phys. Rev. Lett. 1990, 65, 195−198. (15) Burlatsky, S. F.; Oshanin, G.; Cazabat, A.-M.; Moreau, M. Microscopic model of upward creep of an ultrathin wetting film. Phys. Rev. Lett. 1996, 76, 86−89. (16) Abraham, D. B.; Cuerno, R.; Moro, E. Microscopic model for thin film spreading. Phys. Rev. Lett. 2002, 88, 206101. (17) Bekink, S.; Karaborni, S.; Verbist, G.; Esselink, K. Simulating the Spreading of a Drop in the Terraced Wetting Regime. Phys. Rev. Lett. 1996, 76, 3766−3769. (18) De Coninck, J.; D’Ortona, U.; Koplik, J.; Banavar, J. R. Terraced Spreading of Chain Molecules via Molecular Dynamics. Phys. Rev. Lett. 1995, 74, 928−931. (19) Voué, M.; Valignat, M. P.; Oshanin, G.; Cazabat, A.-M.; De Coninck, J. Dynamics of Spreading of Liquid Microdroplets on Substrates of Increasing Surface Energies. Langmuir 1998, 14, 5951− 5958. (20) Heine, D. R.; Grest, G. S.; Webb, E. B., III Spreading dynamics of polymer nanodroplets. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 68, 061603. (21) Heine, D. R.; Grest, G. S.; Webb, E. B., III Spreading dynamics of polymer nanodroplets in cylindrical geometries. Phys. Rev. E 2004, 70, 011606. (22) Wu, C.; Qian, T.; Sheng, P. Droplet spreading driven by van der Waals force: a molecular dynamics study. J. Phys.: Condens. Matter 2010, 22, 325101. (23) Voué, M.; De Coninck, J. Spreading and wetting at the microscopic scale: recent developments and perspectives. Acta Mater. 2000, 48, 4405−4417. (24) Tiberg, F.; Cazabat, A.-M. Spreading of Thin Films of Ordered Nonionic Surfactants. Origin of the Stepped Shape of the Spreading Precursor. Langmuir 1994, 10, 2301−2306. (25) Heine, D. R.; Grest, G. S.; Webb, E. B., III Surface Wetting of Liquid Nanodroplets: Droplet-Size Effects. Phys. Rev. Lett. 2005, 95, 107801. (26) Burlatsky, S. F.; Oshanin, G.; Cazabat, A.-M.; Moreau, M.; Reinhardt, W. P. Spreading of a thin wetting film: Microscopic approach. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1996, 54, 3832−3845. (27) Leroy, F.; Müller-Plathe, F. Dry-Surface Simulation Method for the Determination of the Work of Adhesion of Solid-Liquid Interfaces. Langmuir 2015, 31, 8335−8345. (28) Lin, S.-T.; Blanco, M.; Goddard, W. A., III The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: Validation for the phase diagram of LennardJones fluids. J. Chem. Phys. 2003, 119, 11792. (29) Lin, S.-T.; Maiti, P. K.; Goddard, W. A., III Two-phase thermodynamic model for efficient and accurate absolute entropy of water from molecular dynamics simulations. J. Phys. Chem. B 2010, 114, 8191−8. (30) Pascal, T. A.; Lin, S.-T.; Goddard, W. A., III Thermodynamics of liquids: standard molar entropies and heat capacities of common solvents from 2PT molecular dynamics. Phys. Chem. Chem. Phys. 2011, 13, 169−181. (31) Webb, E. B., III; Grest, G. S.; Heine, D. R. Precursor Film Controlled Wetting of Pb on Cu. Phys. Rev. Lett. 2003, 91, 236102. (32) Lotfi, A.; Vrabec, J.; Fischer, J. Vapour liquid equilibria of the Lennard-Jones fluid from the NpT plus test particle method. Mol. Phys. 1992, 76, 1319−1333.

(33) Borodin, O.; Smith, G. D. Development of many-body polarizable force fields for Li-battery components: 1. Ether, alkane, and carbonate-based solvents. J. Phys. Chem. B 2006, 110, 6279−92. (34) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. Packmol: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157−2164. (35) Shinoda, W.; Shiga, M.; Mikami, M. Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 134103. (36) Schneider, T.; Stoll, E. Molecular-dynamics study of a threedimensional one-component model for distortive phase transitions. Phys. Rev. B: Condens. Matter Mater. Phys. 1978, 17, 1302−1322. (37) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible multiple time scale molecular dynamics. J. Chem. Phys. 1992, 97, 1990−2001. (38) Isele-Holder, R. E.; Mitchell, W.; Ismail, A. E. Development and application of a particle-particle particle-mesh Ewald method for dispersion interactions. J. Chem. Phys. 2012, 137, 174107. (39) Isele-Holder, R. E.; Mitchell, W.; Hammond, J. R.; Kohlmeyer, A.; Ismail, A. E. Reconsidering Dispersion Potentials: Reduced Cutoffs in Mesh-Based Ewald Solvers Can Be Faster Than Truncation. J. Chem. Theory Comput. 2013, 9, 5412−5420. (40) Isele-Holder, R. E.; Ismail, A. E. Classification of precursors in nanoscale droplets. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2016, 93, 043319. (41) Tolman, R. C. Consideration of the Gibbs Theory of Surface Tension. J. Chem. Phys. 1948, 16, 758−774. (42) Kirkwood, J. G.; Buff, F. P. The Statistical Mechanical Theory of Surface Tension. J. Chem. Phys. 1949, 17, 338−343. (43) Lin, T.-S.; Pascal, T. A.; Blanco, M.; Goddard, W. A., III California Institute of Technology, Pasadena, California. Two-Phase Thermodynamics (2PT) Program for Entropy and Free Energies of Complex Systems, Materials and Process Simulation Center, 2010.

G

DOI: 10.1021/acs.langmuir.6b00807 Langmuir XXXX, XXX, XXX−XXX