Nanoscale thin-film flows with thermal fluctuations and slip

The combined effects of thermal fluctuations and liquid-solid slip on nanoscale thin-film flows are investigated using stochastic lubrication equations (SLEs). The previous no-slip SLE for films on plates is extended to consider slip effects and a new SLE for films on fibers is derived, using a long-wave approximation to fluctuating hydrodynamics. Analytically derived capillary spectra, which evolve in time, are found from the new SLEs and compared to molecular dynamics simulations. It is shown that thermal fluctuations lead to the generation and growth of surface waves, and slip accelerates this growth. SLEs developed here provide useful tools to study nanoscale film dewetting, nanofiber coating, and liquid transport using nanofibers.


I. INTRODUCTION
The emergence of microfluidic and nanofluidic technologies has focused attention on the nature of liquid flows at the nanoscale. Nanoscale flows can often behave in a manner that is qualitatively different from those at the macroscale. One important difference is the breakdown of the deterministic description at the nanoscale and the increasingly dominant role of thermal fluctutaions of liquid molecules. For example, at the nanoscale, the position of a contact line in equilibirum is not fixed, as deterministic models predict, but instead fluctuates with a Gaussian probability distribution [1]; the breakup of liquid nanojet in molecular dynamics simulations (MD) shows a double-cone rupture profile, in contrast to the longthread profile predicted by deterministic equations [2], and deterministic models fail to accurately predict the rupture time of a dewetting polymer nanofilm [3,4].
Stochastic modeling is thus essential to capture nanoscale flow physics, and broadly there are two options: a model based on Langevin theory or one based on fluctuating hydrodynamics (FH). Langevin theory is based on fluctuation-dissipation theorem and its applications to stochastic modeling include Brownian motion, single-file water transport in carbon nanotubes (CNT) [5], fluctuating contact lines [1,6], fluctuating contact angles [7], and more. The equations of FH, proposed by Landau and Liftshiz [8], are stochastic versions of the Navier-Stokes (NS) equations, with thermal fluctuations modelled by an additional random stress tensor. The applications of FH include the study of fluctuations in a dilute gas [9], fluctuations in diffusively mixing fluids [10], and the rupture of nanoscale liquid jets and films [2,4].
Thin-film flows are characterized by disparities of length scale in different dimensions, i.e., the ratio of film width to film length is very small: = h/λ 1. This allows the adoption of a long-wave approximation to derive lubrication equations (LE) from the full governing equations, reducing the dimensionality of the problem. In terms of macroscopic thin-film flows, which are well studied, deterministic LEs have been developed from NS: for liquid jets (which we here refer to as the "jet LE") to study breakup of liquid threads [11]; for planar liquid films on plates (the "planar-film LE") [12], which has important applications in dewetting of polymer films [3,13]; and for an annular liquid film on a fiber (the "annular-film LE") [14][15][16][17], which is used to predict surface morphologies of falling liquid films down a fiber [14,15] and fiber coating [16,17]. LEs are much simpler to solve than the NS and can be easily extended to consider many interesting effects, such as electric field forces [18,19], Marangoni stresses [20,21], and evaporation [22].
For thin-film and jet flows at the nanoscale, stochastic lubrication equations (SLEs) are obtained from FH, also using the long-wave approximation. For example, Moseler and Landman [2] proposed an SLE for nanojets (the "jet SLE") that is able to predict the double-cone rupture profile observed in MD, unlike the jet LE. Later, based on the jet SLE, Eggers [23] confirmed that thermal noise leads to a symmetric self-similar rupture profile. Grün et al. [4] and Davidovitch et al. [24] derived a SLE for planar liquid films on substrates (the "planar-film SLE") able to reconcile the discrepancy in dewetting time between experiments [3] and solutions to the planar-film LE. The rupture of thin films has subsequently been widely investigated by numerical solutions to the planarfilm SLE [4,[25][26][27][28], showing thermal fluctuations indeed speed up the instability growth. Linear stability analyses of the jet SLE and planar-film SLE lead to capillary spectra of surface waves [29][30][31][32][33], which have been investigated in MD simulations [31][32][33] and experiments [30]. The analytical solutions to SLEs in previous work [32,33] show thermal fluctuations can massively amplify the growth of waves and cause the dominant wavelength to evolve in time, in contrast to a constant value predicted by LEs. Both the jet SLE and planar-film SLE can be viewed as extensions of the jet LE and planar-film LE with an additional, appropriately scaled, noise term. However, an extension of the annular-film LE to a stochastic version (an "annular-film SLE") is performed for the first time in this article.
Another aspect that makes nanoscale flows differ from macroscale flows is the breakdown of the no-slip condition and the greater importance of slip effects on flows. For example, the fast water transport inside CNTs, which is far beyond the predictions of no-slip hydrodynamics, can be attributed to slip at the water-carbon interface [34]. This flow-rate enhancement makes CNTs promising membrane materials for ultrafiltration devices [35,36]. Despite lots of evidence of slip in MD simulations and experiments, measured by nonequilibrium simulations or shear-driven methods, physical explanations for the origin of slip were not initially clear. To resolve this issue, Bocquet and Barrat [37] proposed a Green-Kubo-type expression for slip length, to identify slip as a property derivable within thermal equilibrium, like viscosity, providing a way to measure slip length in equilibrium simulations. Notably, the slip length measured by nonequilibrium simulations, like in Couette flows, find values that are constant at low-shear rates but grow rapidly at high-shear rates [38]. The constant slip at low-shear rates is thus where nonequilibrium and equilibrium measurements should coincide [39].
It is clear that at the nanoscale both slip and thermal fluctuations play an important role. However, until now, these effects have never been explored simultaneously for thin-film flows, due to the lack of proper tools: The SLE proposed in Refs. [4,24] does not model slip. In this paper, the planar-film SLE is extended to consider slip effects and a new annularfilm SLE is derived, also with slip modelled. These newly developed SLEs are validated against MD simulations and used to study the effects of thermal fluctuations and slip on capillary wave evolution of thin liquid films. This paper is organized as follows. In Sec. II, we briefly present the equations of FH. In Sec. III and Sec. IV, we derive the planar-film SLE and annular-film SLE, with slip effects, from FH, respectively. In Sec. V, the static spectra and evolving spectra of capillary waves are derived. In Sec. VI, the MD model of nanoscale liquid films, on both plates and fibers, is introduced. Section VII, shows the comparison between MD and SLE and includes discussion on the effects of thermal fluctuations and slip. We conclude our findings and outline future directions of research in Sec. VIII.

II. FLUCTUATING HYDRODYNAMICS
The equations of FH for incompressible constant-density flow are given by [8]: FIG. 1. Sketch of a liquid film on a plate, where h = h(x, t ) is the film thickness, λ is the characteristic length, u is the x component of liquid velocity, and is the liquid-solid slip length. The film has depth L y in the y direction (into the page).
Here u, ρ, t, p, and μ are the velocity field, density, time, pressure, and viscosity, respectively. The random stress tensor τ captures thermal fluctuations, modelled by white noise with zero mean and covariance where x = (x, y, z) and · · · is the ensemble average.

III. PLANAR-FILM SLE WITH SLIP EFFECTS
For a two-dimensional (2D) film on a substrate ( Fig. 1) with thermal fluctuations but without slip, Grün et al. [4] derived the planar-film SLE. To do this, a long-wave approximation to the equations of FH was adopted, i.e., ε = h 0 /λ 1. We extend this analysis to consider the effect of slip at the liquid-solid interface, using a similar long-wave approximation to the equations of FH, but with a special treatment for the slip boundary condition.
In terms of a 2D film on a plate, Eq. (1) in Cartesian coordinates is ∂u ∂x and Eq. (2) is Here u and w are the x and z components of velocity, and ψ is a 2D random stress tensor with zero mean and covariance given by The factor 1/L y appears because the planar films we consider here are quasi-2D (L y L x ), allowing all variables of 053105-2 interest to be averaged over the y direction; in particular, ψ = 1 L y L y 0 τ dy. For boundary conditions, at z = h, we have the dynamic condition: where σ is the hydrodynamic stress tensor, γ is the surface tension, φ(h) is the disjoining pressure, and n is the normal vector at the surface, The kinematic condition at z = h is ∂h ∂t At z = 0, Navier's slip boundary condition is where is the slip length. The substrate is impermeable so that To get a lubrication equation from governing equations Eqs. (4)-(12), we need to establish the leading-order terms of Eqs. (4)-(12) by their asymptotic expansion in , for which we use the scaling relations suggested by Grün et al. [4], before then returning to the dimensional form of the equations for ease of interpretation: Here nondimensional variables are upper case and u 0 is the characteristic velocity. Using these scaling relations, Eq. (5) and Eq. (6) become Here the Reynolds number Re = ρu 0 h 0 /μ is assumed to be order one or smaller which is usually valid for nanoscale flows. The leading-order form of the dimensional momentum equations Eq. (5) and Eq. (6) are thus On the other hand, the leading-order form of the incompressible condition, Eq. (4), remains unchanged.
In the normal direction to the surface, the scaled dynamic condition, Eq. (8), is (18) and the scaled dynamic condition in the tangential direction to the surface is Therefore, the leading order of the dimensional dynamic boundary conditions in the normal and tangential direction are seen to be The leading order of the kinematic condition is unchanged. For the slip boundary Eq. (11), its scaled form is Thus, the leading-order slip boundary depends on the ratio /h 0 . In the present work, we consider the slip is at the order of h 0 so that the leading order of slip boundary condition is unchanged from Eq. (11).
To obtain the lubrication equation, Eq. (16) is integrated in the z direction, and using Eq. (21) gives Integrating Eq. (23) from 0 to z then gives 053105-3 and u| z=0 is determined as Thus the expression for velocity u is Integrating Eq. (4) in z from 0 to h and using the boundary condition Eq. (10) and (12) give By substituting Eq. (26) into (27) one can obtain where the mobility M(h) = 1 3 h 3 + h 2 . For the above double integral, integration by parts leads to Before simplifying the noise terms in Eq. (29), the covariance of ψ zx | z=0 has to be determined, which we will show is given by Here η = μ/ is the so-called friction factor. Equation (30) is inferred from Bocquet and Barrat's Green-Kubo type expression of friction factor [37], which is where F f is the friction force at the wall (whose area is L x L y ) along the x direction. As which suggests that the friction factor is related to the stochastic shear stress at the boundary. Notably, this expression for friction factor appears analogous to the Green-Kubo expression for bulk shear viscosity, namely, μ = , from which one can obtain the covariance of bulk shear stress [40]. Thus the covariance of ψ zx | z=0 , namely, Eq. (30), is obtained from Eq. (32). Now to simplify the noise terms in Eq. (29), we use the method provided in Appendix A for the simplification of the stochastic integral. The two noise terms in Eq. (29) are respectively simplified to The √ in Eq. (34) comes from Eq. (30). The ξ 1 and ξ 2 have zero mean and covariance Note that the noise term in the bulk is assumed to be uncorrelated to the value at the surface, Eq. (37), so that the two noise terms can be combined together: Here ξ 3 has the same covariance as ξ 1 and ξ 2 . Thus we derive the planar-film SLE with slip effects as is the driving pressure, and the noise ξ (the subscript "3" is omitted) has zero mean and covariance ξ ( . Notably, to derive the planar-film SLE including slip effects, it is essential to include a random friction force at the liquid-solid interface as defined by Eq. (30). This finding was also noticed in the study of Brownian motion of solid particles within fluids when slip effects were prominent [41]. We also note that there are several other ways [4,24,27,42] to simplify the stochastic integral in Eq. (29). The fluctuation-dissipation theorem is satisfied by the SLE: The variance of the noise (the square of the prefactor to the noise term) is equal to the mobility, which appears in the diffusion term (the first term of the right-hand side of the SLE). For efficient numerical solutions of stochastic partial differential equations such as this, the reader is referred to the numerical integration schemes discussed in Refs. [43,44].
The generalization of this 2D planar-film SLE to a 3D version is quite straightforward, and one can obtain Here Q(h) = 1 3 h 3 + h 2 , p = −γ h + φ, and the noise ς has zero mean and covariance,

IV. ANNULAR-FILM SLE WITH SLIP EFFECTS
Consider now axisymmetric flow of an annular liquid film on a fiber, see Fig. 2. Using the same procedure as for the planar thin-film case, at leading order in = h 0 /λ 1, the incompressible condition in cylindrical coordinates is Here w and u are the axial and radial components of velocity, respectively. The governing momentum equations in cylindrical coordinates (with only leading-order terms) become Here ζ rz is the rz component of the random stress tensor ζ , which has zero mean and covariance, as we will show now. The covariance of τ in cylindrical coordinates is Due to the axisymmetry of the flow the variables of interest are averaged over the azimuthal direction, e.g., ζ = 1 2π 2π 0 τ dθ . Thus the covariance for the azimuthally averaged fluctuating stress tensor is At the free surface r = h, the dynamic condition in the normal direction of the surface is Notably, the term γ ∂ 2 h ∂z 2 is not at leading order, but conventionally in this field it is still included in the pressure in an attempt to extend the validity of the model [11,15]. In the tangential direction, The kinematic condition at r = h is ∂h ∂t At the substrate surface r = a, the slip boundary is and the impermeability condition is Integrating Eq. (43) twice and using two boundary conditions Eq. (49) and (51) lead to an expression for velocity where the mobility G(h) is A simplification of the double integral in Eq. (55) leads to From Bocquet and Barrat's expression of friction factor (in cylindrical coordinates), the covariance of ζ rz | r=a is found  to be Now we are left with the task of finding an equivalent stochastic partial differential equation that does not contain integrals. Using the method in Appendix A, the two noise terms in Eq. (57) are simplified respectively to Here the noise β 1 and β 2 have covariance, Equation (63) results from the assumption that the bulk noise is uncorrelated with the noise at the boundary, allowing the two noise terms to be combined: Here β 3 has the same covariance as β 1 and β 2 . Thus the stochastic lubrication equation for an annular film on a fiber with slip effects (the annular-film SLE) is derived as Here p = −γ ( ∂ 2 h ∂z 2 − 1 h ) + φ, and the noise β (the subscript "3" is omitted) has zero mean and covariance ). Notably, the derived annular-film SLE can be viewed as an extension of the existing annular-film LE [15]; an extension that constitutes the addition of an appropriately scaled noise term.

V. STATIC AND EVOLVING CAPILLARY WAVE SPECTRA
The study of fluctuating capillary waves has long been the interest of researchers, as they are essential to modern theories of surface physics [45,46] and they contain information relating to the properties of liquid-solid systems [47], providing a noninvasive method of measuring transport properties. Here, the newly developed SLEs are used to investigate the effects of thermal fluctuations and the slip length on the temporal growth of capillary waves of thin liquid films (planar and annular) in terms of their spectra.

A. Static spectra
The static spectrum of capillary waves on a planar film can be determined by the equipartition theorem; this forms the basis of classical capillary wave theory (CWT) [45,46,48]. For a planar film influenced by the interfacial potential , whose derivative is disjoining pressure φ = ∂ ∂h , the free energy f 1 (the subscript "1" denotes that variables are for planar films henceforth) due to the change of surface area (in 2D configuration) is For small deformations (∂h/∂x 1), Let us define the Fourier transform of δh = h(x) − h 0 as δh = δhe −iqx dx. From Parseval's theorem, we can express Eq. (67) in terms of Fourier modes: As each summand appears quadratically, it has the same energy 1 2 k B T , from equipartition theorem, so that Thus, the well-known CWT for a planar film is derived [45,46,48]: It is easy to extend consideration to an axissymmetric annular film. The free energy for an annular film, f 2 (the subscript "2" denotes that variables are for annular films henceforth), is so that in Fourier space Applying the equipartition theorem, we obtain the static spectrum of an annular film, (73)

B. Time-dependent spectra
The static spectrum from CWT does not describe how capillary waves develop in time. For example, in the absence of an interfacial potential, CWT cannot reveal how an initially smooth surface of a planar film with | δh| 2 = 0 evolves into a rough surface with | δh| To resolve this problem, we perform a linear stability analysis of the newly developed planar-film SLE and annular-film SLE to obtain time-dependent spectra (for details, see Ref. [33]). For a planar film, with Eq. (40), one can obtain Here | δh(q, 0)| 2 is the initial condition of the film surface and the dispersion relation is For wave numbers with ω 1 (q) > 0 and in the long-time limit t → ∞, Eq. (74) simplifies to Eq. (70). Thus the static spectrum from CWT is a special case of the time-dependent spectrum derived here. We also note that the right-hand side of Eq. (74) consists of two terms. The first term is from the deterministic part of the SLE (namely, the LE) and the second term comes from the stochastic part of the SLE [33].
For an annular film, the time-dependent spectrum derived from the annular-film SLE is where the dispersion relation is In the following sections we use S i = | δh(q, t )| 2 i (i = 1, 2) for notational simplicity.

VI. MOLECULAR DYNAMICS SIMULATION
To validate the time-dependent spectra obtained from the planar-film SLE and the annular-film SLE, MD simulations are performed of nanoscale liquid films on planar and cylindrical substrates, using the open-source code LAMMPS [49]. Previously, we have validated the no-slip planar-film SLE [4] using MD simulations as well [33]. All simulation domains contain three phases with the liquid bounded by the vapor The liquid of the film is argon, simulated with the standard Lennard-Jones (LJ) 12-6 potential: 6 ], where ll denotes liquid-liquid interactions and i j represents pairwise particles. The energy parameter ε ll is 1.67 × 10 −21 J and the length parameter σ ll is 0.34 nm. The temperature of this system is kept at T = 85 K or T * = 0.7ε ll /k B (* henceforth denotes LJ units). At this temperature, the number density of liquid argon is n * l = 0.83/σ 3 ll . The number density of the vapor phase is 1/400n * l . For a planar substrate shown in Figs. 3(a) and 3(b), the solid is platinum with a face centered cubic (fcc) structure and its 100 surface in contact with the liquid. The platinum number density is n * s = 2.60/σ 3 ll . For a cylindrical substrate, the one in Figs. 3(c) and 3(d) is generated by cutting out a cylinder from a large cuboid of platinum. All solid substrates are assumed to be rigid, which saves considerable computational cost. The liquid-solid interactions are also modelled by the same 12-6 LJ potential with σ ls = 0.8σ ll for the length parameter. For planar films, three different values of the energy parameter are used, in order to generate varying slip lengths: Case (P1) ε ls = 0.65ε ll , Case (P2) ε ls = 0.35ε ll , and Case (P3) ε ls = 0.2ε ll . For an annular film on a fiber: Case (A) ε ls = 0.7ε ll . The cut-off distance, beyond which the intermolecular interactions are omitted, is chosen as r c * = 5.5σ ll . The surface tension is γ = 1.52 × 10 −2 N/m (obtained from the virial expression) and dynamic viscosity μ = 2.87 × 10 −4 kg/(ms) (obtained by the Green-Kubo method).
The initial dimensions of a planar liquid film [see Fig. 3(a)] are L x = 313.90 nm, L y = 3.14 nm and h 0 = 3.14 nm so that L x L y , making the 3D MD simulations quasi-2D, which allows us to consider large aspect ratio films and compare to 2D theories. The substrate has a thickness of h s = 0.78 nm (composed of five layers of platinum atoms). The initial size of the annular film [see Fig. 3(c)] has film length L x = 229.70 nm, outer radius h 0 = 5.74 nm. As the fiber is cut from a box of platinum with fcc structure, its size is defined The liquid-vapor interface is defined as the usual equimolar surface. After the interface is extracted from MD simulations (see Ref. [33] for details), the interfacial profile is averaged over the y direction for a planar film and over θ direction for an annular film. Then Fourier transforms of the interfacial profile are performed and averaged over a number of independent realizations (65 realizations for a planar film and 10 realizations for an annular film) to get the time-dependent spectra of capillary waves S shown in Fig. 4, for planar films, and in Fig. 6 for annular films.

A. Spectra of planar films
Figures 4(a)-4(c) show spectra of planar films with two different slip lengths. Since the film thickness is larger than the cut-off distance, there are no interactions between the film surface and solid surface. As such, the disjoining pressure φ vanishes, which leads to ω 1 > 0 for all wave numbers. This means that, according to classic theory, the surface should remain smooth for an initially smooth surface (| δh(q, 0)| 2 = 0) shown in Fig. 1(a). However, as illustrated by Fig. 1(b), the free surface becomes rough and the spectra of surface wave evolve with time; see MD results (triangles) in Figs. 4(a) and 4(b). Notably, this growth is toward an upper limit, which is the static spectrum L x L y k B T γ q 2 (orange dash-dot). Meanwhile, this growth is driven purely by thermal fluctuations-the deterministic term in Eq. (74) makes no contribution to the spectrum.
Notably, the no-slip planar-film SLE [4] cannot accurately predict the MD results, as shown in the inset of Fig. 4(b). However, using the measured slip length ( = 0.68 nm for P1, = 3.16 nm for P2, and = 8.77 nm for P3) and effective thickness of fluid domain h 0 = 2.90 nm (see Appendix C for measurements of slip length and the definition of effective thickness, arising from the fact that the hydrodynamic boundary does not align with the solid surface), analytical spectra from the newly derived planar-film SLE with slip (solid lines) agree well with MD results [see Fig. 4(a) and Fig. 4(b)]. Notably, the agreement gets worse when slip length becomes larger than the film thickness [see Fig. 4(c)]. In the derivation of the SLE, the slip length is assumed to be of the same order as the film thickness-so we can infer that the observed worsening agreement is a result of this assumption becoming progressively less valid. By comparing Fig. 4(a) and Fig. 4(b), slip is shown to speed up the evolution of capillary spectra. However, from Eq. (74), the spectra with different slip lengths will approach the same static spectrum, as t → ∞.
The dominant wave number q d decreases with time and is found to be from Eq. (74) by requiring ∂S 1 ∂q = 0. This fascinating, entirely nonclassical, result is confirmed by the MD results in Fig. 5. Figure 6 shows evolving spectra for the capillary waves of an annular film. For wave number qh 0 > 1, the MD spectra (triangles) at different times collapse to the static spectra L x 2πh 0 k B T γ (q 2 −1/h 2 0 ) . However, for qh 0 < 1, the Laplace pressure from the circumferential curvature results in a negative dispersion relation (see the inset of Fig. 6) so that the amplitude will grow unboundedly until the film ruptures, where beads are formed as shown in Fig. 1(d). Using the measured slip length = 0 with the hydrodynamic boundary position at a = 2.60 nm and the surface position at h 0 = 5.74 nm (see Appendix C), the analytical spectrum from the annular-film SLE (solid lines) is in reasonable agreement with MD results shown in Fig. 6. The difference between MD and the annular-film SLE in terms of spectral amplitude increases with time, and this is because the annular-film SLE overpredicts the dispersion relation, as shown in the inset. By the comparison of the dispersion relation to the full governing equations based on Stokes flow [15], it is known that the accuracy of the deterministic annular-film LE depends on the ratio a/h 0 with improving agreement for larger a/h 0 . Since the dispersion relation of the annular-film SLE is the same as the one of the annular-film LE, it is expected that the annular-film SLE has better accuracy for larger a/h 0 as well (currently, a/h 0 = 0.45).

B. Spectra of annular films
As for the dominant wave number for a planar film, the dominant wave number for an annular film also decreases with time but approaches a constant value which is q dc h 0 = √ 2/2 from the annular-film LE. It is interesting to know at what timescale the q d reaches the constant value q dc . By requiring we obtain 2e 2ω 2 (q d )t ω 2 (q d )t[2 − 1/(q d h 0 ) 2 ] + 1 − e 2ω 2 (q d )t = 0. When q d h 0 is close to √ 2/2, this expression is simplified to Thus, the timescale for the dominant wave number to reach q dc h 0 = √ 2/2 is . However, it is not feasible to extract the dominant wave number from current MD simulations to confirm Eq. (79), as the number of data points is too sparse, as shown in Fig. 6. One would need to simulate a much longer film to have data dense enough to verify Eq. (79), and the computational cost of this is currently prohibitively high.

VIII. CONCLUSION
To investigate the combined effects of thermal fluctuations and slip on thin-film flows at the nanoscale, we derive stochastic lubrication equations for thin films on both plates and fibers using the long-wave approximation to fluctuating hydrodynamics, which are validated by molecular simulations in terms of capillary wave spectra. To complete the derivation, it is essential to include a random stress at the liquid-solid interface and its covariance is related to the slip length. The derivation process is general, as we demonstrated for films on different substrate geometries, and can be extended to study other thin-film flows at the nanoscale, like free liquid films [50].
The analysis of SLEs here is linear, since our focus is on the wave growth at early stages. In the future, it will be interesting, using the developed SLEs, to investigate how thermal fluctuations and slip influence the dynamics during the nonlinear stages of film growth, such as on rupture profile, rupture time and droplet distribution. The SLEs developed here also provide useful tools to study nanoscale film dewetting, nanofiber coating and liquid transport using nanofibers [51] where thermal fluctuations and slip are expected to play an important role.
Despite the successes, there are a number of ways in which the models derived here can be improved to predict fluctuating capillary waves. The SLEs are derived using a long-wave approximation, which means the dynamics of capillary waves with short wavelengths are inaccurate. To derive the SLEs, the assumption of small slip length is adopted, which limits their application to relatively small slip systems. A more general model, beyond the current lubrication framework, is needed to predict the growth of capillary waves in these more general cases and presented in the future work. (A10)

APPENDIX B: MD INITIALIZATION
For a planar film, a liquid film with thickness h 0 = 3.14 nm and vapor are equilibrated separately in periodic boxes at T = 85 K. Then the film is deposited above the substrate and with a vapor on top of the film. Because there exists a gap (a depletion of liquid particles) between the solid and liquid, arising from the repulsive force in the LJ potential, it is necessary to deposit the liquid above the substrate by some distance. The thickness of the gap is found to be about 0.2 nm after which the liquid-solid system reaches equilibrium, so that we choose a deposit distance d = 0.2 nm. This makes the position of the film surface at h 0 + d = 3.34 nm, initially. For an annular film, liquid molecules and vapor molecules are equilibrated separately, in cuboid periodic boxes at T = 85 K. Then an annular film is cut out from the cuboid box with an outer radius at h 0 = 5.74 nm and an inner radius above the fiber radius, with an interval 0.2 nm. Then the fiber is put into the annular film and vapor placed to surround the film. Notably, in this case, the position of film surface is still at h 0 = 5.74 nm, initially.

APPENDIX C: MEASUREMENTS OF SLIP LENGTH
Since we have chosen varying liquid-solid interactions for three cases (P1, P2, and A), it is necessary to know what slip length those cases have. Slip length is thus measured by nonequilibrium simulations of pressure-driven flow past a substrate surface as shown by the MD snapshots in the top-left corner of Fig. 7 (for a planar film) and Fig. 8 (for an annular film). The pressure gradient is created by applying a body force g to the fluid. The generated velocity distribution is u(z) = ρg 2μ (z − z 1 )(2z 2 − z 1 − z) + u s for a planar film. Here z 1 and z 2 are positions of the hydrodynamic boundary (HB) and free surface (FS) for a planar film, respectively, and u s is the slip velocity at the HB. For an annular film, the axisymmetric velocity profile is u(r) = − ρg 4μ [r 2 − r 2 1 − 2r 2 2 log(r/r 2 )] + u s , where r 1 and r 2 are positions of the HB and FS for this system.
The precise location of two boundary positions for each system is not trivial since there is an interfacial zone between the two different phases (solid-liquid and liquid-vapor) as demonstrated by the density distribution( the orange line in Fig. 7 and Fig. 8). For the HB, research has shown it is located inside the liquid, between the first-peak density and the second-peak density, rather than being located at the solid surface [37,52], by comparing the analytical solution and MD measurements of the correlations of momentum density (an offset which matters when the interfacial layer has comparable width to the film). In line with this finding, we choose the position of the HB at the first valley of the density distribution: z * 1 = 1.3σ for a planar film and r * 1 = 7.65σ for an annular film (see Fig. 7 and Fig. 8). The position of FS is determined in the standard way by the location of the equimolar surface where density is 0.5n * l , with z * 2 = 9.8σ for a planar film and r * 2 = 16.55σ for an annular film (see Fig. 7 and Fig. 8).
After locating the boundary, the slip velocity is obtained by fitting velocity profiles of MD data (symbols) with analytical expressions of velocity (solid black lines) as shown in Fig. 7. The slip length is the distance between the HB and the position where the the linear extrapolation of the velocity profile vanishes. Figure 7 is, in particular, for the measurement of the slip length of case P2 where the slip length is measured to be * = 9.3σ (3.16 nm). Two different values of driving force g * = 0.01 and g * = 0.006 are used to demonstrate that the measured slip length is a constant, independent of the driving force (g * 0.01). However, as the inset shows, the slip length does become force (shear)-dependent for g * 0.01, which is beyond the current consideration [38]. The driving forces in the free-surface flows studied for capillary waves are small enough that the assumption of a constant slip length holds. In the same way, = 0.68 nm is measured for case P1 and = 8.77 nm for case P3. For an annular film, as shown in Fig. 8, the slip length is = 0 nm (no-slip) for case A.
We note that as the HB does not align with the edge of the solid, the effective thickness of the fluid domain simulated for capillary waves is different to its initial thickness. For a planar film, as the position of the initial free-surface is at 3.34 nm (see Appendix B) and the HB is at z 1 = 0.44 nm, the effective thickness of a planar film is 2.9 nm. For an annular film, this means a = r 1 = 2.6 nm and outer radius h 0 is 5.74 nm.