Molecular simulation of thin liquid films: Thermal fluctuations and instability

The instability of a thin liquid film on a solid surface is studied both by molecular dynamics simulations (MD) and a stochastic thin-film equation (STF), which models thermal fluctuations with white noise. A linear stability analysis of the STF allows us to derive a power spectrum for the surface fluctuations, which is quantitatively validated against the spectrum observed in MD. Thermal fluctuations are shown to be critical to the dynamics of nanoscale films. Compared to the classical instability mechanism, which is driven by disjoining pressure, fluctuations (a) can massively amplify the instability, (b) cause the fluctuation wavelength that is dominant to evolve in time (a single fastest-growing mode does not exist), and (c) decrease the critical wavelength so that classically stable films can be ruptured.


I. INTRODUCTION
The spontaneous rupture and dewetting of thin liquid films deposited on substrates is of broad engineering interest due to a myriad of applications in coatings [1,2], lubricants [3], chemical sensors [4], and microfluidics [5]. Dewetting experiments of polymeric [6][7][8][9][10] and metal films [11] have revealed three different instability mechanisms [12]: spinodal dewetting due to undulation growth, thermal nucleation of random holes, and heterogeneous nucleation around defects. While the rupture of a film can be influenced by all three modes simultaneously [13], the breakup of relatively thin films, where the disjoining pressure is strong, is dominated by the spinodal dewetting mode [6][7][8][9]. In this regime, thermally excited capillary waves can be amplified by the disjoining pressure, but in competition with the restoring force of surface tension, such that only disturbances above a critical wavelength grow. The presence of unstable modes can be used as an indication that spinodal dewetting is the primary dewetting mechanism for a particular liquid film [6,10,11,13].
Thin film flows, though complicated by their free surface, can usually be described relatively simply by the thin-film equation (TF), a lubrication approximation to the Navier-Stokes equations (NS) [14]. Williams and Davis [15] derived the TF equation by treating the disjoining pressure as a body force on the film. Later, Oron [16] and Becker [17] obtained numerical solutions to this equation and found good agreement with experiments for the evolution of dewetting patterns. growing mode and a critical wavelength (perturbations with a wavelength larger than this value grow) [15]. The critical wavelength is often used to predict whether a film of a particular size is stable, and the fastest growing mode is used to estimate the rupture time and number of holes or drops after film rupture [6,13].
As the thickness of thin liquid films reaches the nanoscale, thermal fluctuations may play an important role in the instability process. In fact, their importance in many kinds of freesurface flows at small scales has been generally recognized. For example, Moseler and Landman [18] found a double-cone rupture profile for nanojets in molecular dynamics (MD) simulations, in contrast to the long-thread profile predicted by the classical lubrication equation simplified from NS equations. It was shown that the deficiency of this lubrication equation for describing nanojet dynamics can be remedied by adding a noise term of appropriate strength to the equation [18]. Later, based on this equation, Eggers found thermal noise leads to a self-similar profile that describes the most probable breakup mode and dominates the pinch-off of nanojets [19]. Recently, Zhao et al. [20] derived a spectrum for growing surface waves of a nanojet due to the Rayleigh-Plateau instability, modified by the inclusion of thermal fluctuations, and found good agreement with MD simulations. Experimentally, Aarts et al. [21] were able to observe capillary waves driven by thermal fluctuations by using mixtures that create ultralow interfacial tension. Davidovitch et al. [22] studied the spreading of drops on solids under the influence of thermal noise and showed that it increased contact line speeds.
In terms of the thin-film flows studied here, Grün et al. [23] derived a stochastic variant of the thin film equation (STF) and their numerical solutions demonstrated that noise can accelerate the rupture of thin films. In later work, Mecke and Rauscher [24] derived a spectrum for the capillary waves in thin-film flows and showed that noise can shift the spectrum from exhibiting an exponential to a power law decay with increasing wave numbers for large wave vectors. This finding was later investigated in experiments on the dewetting of thin polymer films [25]. However, due to experimental limitations, a direct comparison between the theoretical and the experimental spectrum could not be made. Nesic et al. [26] conducted detailed numerical solutions of the STF equation and found that thermal fluctuations reduce the number of formed droplets and lead to larger variability in their size and spatial distribution. Diez et al. [27] considered the noise to be spatially correlated and found that, to produce close agreement between their theoretical spectrum and experimental data for metal films, different correlation lengths for the nonwhite noise had to be employed.
So far, experimental studies on the effects of thermal fluctuations on thin film flows are limited due to the technical difficulties associated with the spatiotemporal scale. As such, MD simulations are a natural and convenient tool to investigate the thin-film problem. Whilst there have been previous MD studies of thin-film rupture [13,28], they have concentrated on the nonlinear stage of instability with the aid of classical theory (no fluctuations).
In this paper, we conduct MD simulations of the rupture of liquid films on substrates and focus on the influence of thermal fluctuations on the initiation of instability. MD results are compared to a power spectrum derived from the STF equation in order to establish the applicability and usefulness of the STF at the nanoscale. While MD simulations are molecularly accurate, their computational intensity is restrictive for simulating most physical systems. This motivates the study of the STF, which also incorporates thermal fluctuations, as an efficient simulation tool for studying the dewetting of liquid films at the nanoscale. This paper is organized as follows. In Sec. II, the MD model of a nanoscale liquid film on a solid surface is introduced. In Sec. III, the stochastic governing equation is presented and the power spectrum for the surface modes of the liquid film is derived. The spectra from MD simulations, along with comparison with the analytical solutions, are shown in Sec. IV. We conclude our findings in Sec. V.

II. MOLECULAR DYNAMICS SIMULATIONS OF A THIN LIQUID FILM
Molecular dynamics simulations are used to simulate the rupture of a liquid film on a substrate. These simulations are performed in LAMMPS [29]. The domain contains three phases with the liquid bounded by the vapor above and the solid below, as shown in Fig. 1(a). The liquid of the film is argon, simulated with the standard Lennard-Jones (LJ) 12-6 potential: where ll denotes liquid-liquid interactions and i j represents pairwise particles. For argon, the energy parameter ε ll , the length parameter σ ll , and atomic mass are 1.67 × 10 −21 J, 0.34 nm, and 6.63 × 10 −26 kg, respectively [30]. The temperature of this system is kept at T = 85 K or T * = 0.7ε ll /k B (* henceforth denotes LJ units and k B is the Boltzmann constant). At this temperature, the mass density of liquid argon is 1.40 × 10 3 kg/m 3 and number density n * = 0.83/σ 3 ll [31]. The number density of the vapor phase is 1/400n * [31]. As such, in the continuum model of Sec. III, the vapor is assumed to be dynamically passive and has no effect on the rupture of the liquid film. The substrate is platinum with a face centered cubic (fcc) structure and its 100 surface in contact with the liquid. The platinum mass density is 21.45 × 10 3 kg/m 3 with an atomic mass of 3.24 × 10 −25 kg [32]. The solid substrate is assumed to be rigid, which saves considerable computational cost. The liquid-solid interactions are also modeled by the same 12-6 LJ potential with ε ls = 0.7ε ll and σ ls = 0.8σ ll , which creates partial wetting of argon on the substrate. The cutoff distance, beyond which the intermolecular interactions are omitted, is chosen as r c * = 5.5σ ll . In order to compare with the predictions of a continuum model, the transport properties of liquid argon are calculated by simulating the thermal motion of a thick layer of liquid argon. It is found that the surface tension γ = 1.52 × 10 −2 N/m (obtained from the virial expression [31]) and the dynamic viscosity μ = 2.44 × 10 −4 kg/(ms) (obtained from the Green-Kubo relation [33]).
The initial dimensions of the liquid film (L x , L y , h) in Fig. 1(a) are chosen so that L x L y , with L x = 313.90 nm and L y = 3.13 nm, making the 3D MD simulations quasi-2D, which allows us to consider large aspect ratio films and compare to 2D theories. Three cases with different film thickness h 0 = (a) 1.18 nm, (b) 1.57 nm, and (c) 1.96 nm are considered. Thus L x is also much larger than h, enabling the system to be described by lubrication theory. The lateral size of the substrate is the same as that of the liquid film and has a thickness h s = 0.78 nm (composed of five layers of platinum atoms). The vapor above the liquid film has a thickness of 15.70 nm.
We initialize each MD simulation as follows. First, the liquid film and vapor are equilibrated separately in periodic boxes at T * = 0.7ε ll /k B . The liquid film is then deposited on the substrate and the vapor on top of the film. After assembly, the positions and velocities of the liquid and vapor atoms are updated with a Nosé-Hoover thermostat. Periodic boundary conditions are applied in the x and y directions of the system, while vapor particles are reflected specularly in the z direction at the top boundary of the system (see Fig. 1).
To obtain the position of the liquid-vapor interface, we first locate all liquid molecules by calculating the number density within one sphere radius of each molecule. The threshold density for defining a liquid molecule is chosen as 0.5n * as the interface molecules have half of their spherical volume in vapor and half in liquid. Liquid molecules are then subdivided into columns and the molecules with the highest z coordinate in each column are found; this defines the liquid-vapor interface. After determining the 2D interface, we average over the y direction (since these simulations are expected to be quasi-2D). Thus the final interface is 1D and allows comparison with the STF equation.
After defining the interface h(x, t ), the undulations δh are obtained by subtracting the initial thickness h 0 from h. Then we perform a discrete Fourier transform of the data and average over a number of realizations [70, 50, and 30 times for cases (a)-(c)] to obtain the spectra presented in Sec. IV.

III. STOCHASTIC THIN-FILM EQUATION
The hydrodynamics of a macroscopic thin liquid film on a substrate can be captured by the thin-film equation [14]: a well-known approximation to the NS equations. This approximation is based on a long-wave assumption, namely, h/λ 1, where h/λ is the ratio of the film thickness h to film characteristic length λ [14]. At the nanoscale, thermal fluctuations lead to a stochastic stress tensor added to the NS equations, which is described by Landau-Lifshitz Navier-Stokes (LLNS) [34] equations. Applying the same lubrication approximation to LLNS, the one-dimensional stochastic thinfilm equation (STF) can be derived (for details, see [23]): Here, φ is the disjoining pressure and ψ is the random stress tensor component, which has the form of uncorrelated white noise, both in time and space. The average of this noise ψ satisfies ψ = 0 and its covariance is given by where μ is dynamic viscosity, δ is the Dirac delta function, and r = (x, z). Notably, the noise amplitude differs from the form usually presented, with a factor 1/L y reflecting that the noise is an averaged value in the y direction and our MD simulations are quasi-2D [22]. The noise term in Eq. (2) is a stochastic integral, which makes the STF equation complex to solve. Grün et al. [23] simplified the problem considerably by proving Eq. (2) is equivalent to a stochastic partial differential equation where ψ is white noise with zero mean and covariance,

A. Disjoining pressure
There are many forms of disjoining pressure reported in literature and different choices of disjoining pressure cause different flow behavior in a thin film [35,36]. The disjoining pressure for the system considered here is the classical φ = A/(6π h 3 ), obtained by integrating the attractive part of 12-6 LJ potential over a semi-infinitely extended substrate (the repulsive part makes a negligible contribution to the linear stability growth), where A is the Hamaker constant equal to the difference between the Hamaker constant of the liquid film itself and the liquid-solid interactions, i.e., A = A ll − A ls [28,37]. A simplistic explanation for this form is that a liquid film will rupture itself due to the attraction of its two surfaces, but the substrate will hinder this by attracting the bottom surface of the film. Finally, due to the finite thickness of the substrate in MD, the disjoining pressure takes the form This is because the top surface of a liquid film, which is located at h from a semi-infinite substrate, has the pressure φ 1 = A ls /(6π h 3 ) and if located at h + h s away from the substrate, the pressure is 3 ]. Thus the contribution from a substrate with thickness h s should be φ 1 − φ 2 . Added to the disjoining pressure of the liquid film itself, we arrive at Eq. (6). The Hamaker constant A ll is 4.5 × 10 −20 J, from the expression 4π 2 ε ll σ 6 ll ρ l 2 , and A ls is 2.61 × 10 −20 J, from the expression 4π 2 ε ls σ 6 ls ρ l ρ s (ρ l and ρ s are number density of liquid and solid in SI units) [38]. Then, in the conventional theory, which contains disjoining pressure but no thermal fluctuations, the critical wavelength and fastest growing wavelength are λ c = −4π 2 γ /(∂φ/∂h) and λ max = −8π 2 γ /(∂φ/∂h), respectively [15]. For the thickest film [case (c)], which has the largest critical wavelength and fastest growing wavelength, λ c is evaluated to be 48.33 nm and λ max = 68.35 nm. Thus the chosen length of the film L x is long enough to contain multiple waves of the fastest growing wavelength for all the cases, and the chosen width of the film L y is small enough to suppress any wave growth in the y direction for all the cases, ensuring the simulations remain quasi-2D.

B. Spectra of surface waves
To compare the surface waves in our MD simulations with those of the STF equations, it is convenient to investigate them in Fourier space. Thus it is natural to derive the spectrum of the surface waves for the liquid film from the STF equations. We start by rewriting Eq. (4) and Eq. (5) in terms of uncorrelated Gaussian white noise N: where the mean of N is zero and autocovariance We linearize Eq. (7) using h = h 0 + δh, with the assumption that δh represents small deviations from the initial film thickness h 0 and the noise amplitude is also assumed to be small, as explained in [24]. By expanding Eq. (7) to first order and combining with Eq. (6), the linearized STF equation is Taking the Fourier transform of Eq. (9) using leads to Here ω(q) is the dispersion relation of the deterministic TF equation: The solution of Eq. (12) can be represented as the linear superposition of two contributions where δh flu is the contribution purely caused by thermal fluctuations and δh det is the solution to the deterministic part of Eq. (12), i.e., ∂ δh ∂t = ω(q) δh, obtained as below: where the initial disturbance is δh(q, 0); here this is the Fourier transform of the liquid surface found in MD simulations at t = 0. To find the contribution of the fluctuating component to the spectrum, we determine the impulse response of the linear system ∂ δh Performing a Laplace transform of Eq. (16) using g(q, s) = ∞ 0 δh(q, t )e −ts dt with zero initial disturbance δh(q, 0) = 0 gives , (17) so that from the inverse Laplace transform, the impulse response is simply As δh is both a random and complex variable, the root mean square (rms) of its norm is sought, which, from Eq. (14), is given by (as the average of δh flu is zero), where from Eq. (15) and from Eq.
Here we have used | N (q, t )| 2 = L x , due to finite length of the discrete Fourier transform used in MD simulations. Thus we derive the spectrum of surface waves of bounded film flow as A similar expression for the surface wave spectrum can be found in Refs. [24,27], but the derivation provided here is potentially more intuitive. We note that, from capillary wave theory, the average magnitude of each mode of waves can be determined from equipartition theorem [21,39]. Here is the interface energy related with disjoining pressure by φ = −∂ /∂h. For large wave numbers ω(q) < 0 and, in the long-time limit L y , which is consistent with capillary wave theory.

IV. RESULTS AND DISCUSSION
In this section, the surface undulation spectrum obtained from MD simulations is compared with the analytical spectrum derived above. Figure 1 shows an MD simulation of a flat liquid film [case (a)], in which perturbations spontaneously grow [ Fig. 1(b)] and subsequently rupture the film [ Fig. 1(c)]. Figure 2 shows the rapidly growing amplitude of certain low wave number disturbances in MD simulations, which suggests that the rupture of these liquid films is mainly due to the spinodal instability. Notably, the analytical spectra (solid red lines in Fig. 2) compare very well with MD simulations. In Fig. 2(a), the inset shows the deterministic spectrum for h 0 = 1.18 nm. As can be seen, the amplitude of waves in the deterministic spectrum is far below that of MD results, showing that conventional models cannot be relied on at such scales. On the other hand, the stochastic spectrum agrees well with MD, suggesting thermal fluctuations substantially amplify the underlying instability and promote the rupture process. The deterministic spectrum in the inset of Fig. 2(a) indicates that the dominant wave number q max (wave number with peak amplitude) is constant over time, while both MD simulations and the stochastic spectra show clearly that this evolves. The dominant mode at different times is extracted from MD simulations and plotted in Fig. 3 (triangles) along with values obtained from the stochastic spectrum (solid red lines) and the deterministic spectrum (dashed black lines). The results demonstrate that thermal fluctuations induce a dominant wave number much higher than classical predictions and that this gradually decreases to the classical result over time, but in these cases, not before rupture of the film. It is also interesting to note that, while the stochastic spectra are strictly only valid in the early linear stages of the rupture process, the dominant wave number matches well with MD results even close to rupture.
Apart from wave amplitude and the dominant wave number, thermal fluctuations can also affect the critical wave number below which waves grow. From the classic theory, the critical wave number (wavelength) for the bounded film studied here is given by q c = − 1 γ ∂φ ∂h (λ c = 2π q c ) and it is 0.35 × 10 9 (17.95 nm), 0.2 × 10 9 (31.41 nm), and 0.13 × 10 9 (48.33 nm) for cases (a)-(c), respectively. We perform here MD simulations with a film length considerably smaller (i.e., conventionally stable): 13 nm, 24 nm, and 36 nm for cases (a)-(c), respectively. Interestingly the result in Fig. 4 shows that a spontaneous rupture still occurs so that the critical wave number has been significantly altered. To identify wave numbers which grow in time we consider the critical wave number q c to be defined by dS dt | q=q c = 0 and find that where | δh(q, 0)| 2 is the initial condition of the film and k B T γ is the square of thermal length which premultiplies the new term due to thermal fluctuations. The expression indicates that, as observed in Fig. 4, thermal noise will increase the critical wave number, making liquid films of a certain thickness more susceptible to rupture by shorter-wave undulations. The parameters | δh(q, 0)| = 0.0977 × 10 −18 m 2 , 0.1328 × 10 −18 m 2 , and 0.1625 × 10 −18 m 2 are directly obtained from MD simulations for cases (a)-(c), respectively, to find that q c = 5.80 × 10 9 (λ c = 1.08 nm), 5.79 × 10 9 (λ c = 1.08 nm), and 5.79 × 10 9 (λ c = 1.08 nm). Therefore, it is seen that at the nanoscale the critical wave number can become independent of film height as thermal fluctuations [second term in Eq. (24)] overwhelm the conventional instability mechanism of disjoining pressure [first term in Eq. (24)]. Given the consistency between the MD results and the theoretical analysis, we can have some confidence that thermal fluctuations have the potential to rupture conventionally stable film geometries.

V. CONCLUSIONS
Thermal fluctuations play an important role in different types of free-surface flow at the nanoscale. In this article, we investigate their effects on the instability of a thin liquid film on a substrate using MD simulations as an experimental probe and, analytically, solving stochastic thin film equation (STF), to provide a deeper insight into the underlying physics. While thermal fluctuations are intrinsically captured in MD simulations, the STF equation models these fluctuations using an appropriately scaled white noise. To facilitate the comparison between MD simulations and the STF equation, we derive a stochastic spectrum of surface waves and show close agreement between the analytical result and MD. By comparison with a deterministic (fluctuation-free) result, we conclude that thermal fluctuations are critical to the nature of the instability of thin-film flows: they significantly intensify the amplitude of undulations, render the dominant wave number time dependent, and decrease the critical wavelength. Thus our work indicates that the consideration of thermal fluctuations is essential when investigating the behavior of liquid films at the nanoscale. A potentially related phenomenon occurs in the study of water transport in carbon nanotubes (CNT), where recent experiments show continuous flowing water in a channel tends to break up and form a consecutive void-water structure [40]. The current work reveals the effects of thermal fluctuations on the initial stages of instability growth; however, there are also abundant interesting flow dynamics in the later stages of the process where thermal fluctuations could play a role. For example, future research could consider the coarsening dynamics after the film has ruptured into droplets [41,42].