Dynamics of liquid nano-threads: fluctuation-driven instability and rupture

The instability and rupture of nanoscale liquid threads is shown to strongly depend on thermal fluctuations. These fluctuations are naturally occurring within molecular dynamics (MD) simulations and can be incorporated via fluctuating hydrodynamics into a stochastic lubrication equation (SLE). A simple and robust numerical scheme is developed for the SLE that is validated against MD for both the initial (linear) instability and the nonlinear rupture process. Particular attention is paid to the rupture process and its statistics, where the ‘double-cone’ profile reported by Moseler & Landmann [Science, 2000, 289(5482): 1165-1169] is observed, as well as other distinct profile forms depending on the flow conditions. Comparison to the Eggers’ similarity solution [Physical Review Letters, 2002, 89(8): 084502], a power law of the minimum thread radius against time to rupture, shows agreement only at low surface tension; indicating that surface tension cannot generally be neglected when considering rupture dynamics.

universal regime, where viscous, capillary and inertial forces are all comparable [14]. These classical theories are supported by numerous macroscopic experiments [6] and there is recent evidence that the transitions between these regimes is rather complex [9,15,16].
Motivated by the emergence of micro/nano-fluidic technology (e.g. lab-on-a-chip devices [17], 3D printing [18], and nano-device fabrication [19]), the validity of these classical theories at the microscale and nanoscale has recently been brought into question. In a pioneering article, Moseler and Landman [1] used molecular dynamics (MD) simulations of nanoscale jets to discover a rupture profile not described by macroscopic theory: a 'double-cone profile' was observed at rupture, a phenomenon attributed to thermal fluctuations that are negligible at usual engineering scales. The importance of thermal fluctuations at these scales has been confirmed by physical experiments using specially prepared low-surface-tension liquid-liquid combinations that enhance fluctuations [20,21] and further MD simulations [22,23].
The double-cone profile was also predicted by solutions to stochastic lubrication equations (SLE) [1]. The SLE are derived by applying a lubrication approximation to the Landau-Lifshitz Navier-Stokes equations (LLNS), which are the Navier-Stokes equations plus stochastic terms used to model the thermal fluctuations [24]. The SLE have two very significant advantages over MD: (i) they are far less computationally intensive to solve numerically and (ii) they are amenable to theoretical analysis that provide more insight into the process. In relation to the latter, in Zhao et al [25], a framework was developed for modelling (linear) interface instability in the presence of thermal fluctuations (named the SLE-RP). At the nanoscale, the SLE-RP shows that Plateau stability can be violated and k max predicted by classical theories is significantly modified (notably becoming time dependent). Interestingly, for rupture dynamics, Eggers [2] has derived a 'nanoscale' similarity solution that incorporates thermal fluctuations and is able to reproduce the double-cone profile observed by Moseler & Landman [1]. However, the SLE has never been focussed on this case in order to study the accuracy of this solution and the validity of the assumptions made (e.g. negligible surface tension).
Numerical solutions to the SLE offer much broader applicability than the analytic results described above, and can be obtained at a small fraction of the computational cost of MD. A set of closely related stochastic equations are well studied for thin-film flows [26][27][28][29][30] and thermally activated vapor-bubble nucleation [31]. But, surprisingly, there are no detailed numerical SLE studies in the literature for liquid thread rupture. In previous work (e.g. [1,23]) only qualitative comparisons between MD and selected SLE realisations were presented. However, the SLE are stochastic, and many independent solutions are needed in order to (i) understand the statistics of the rupture process, and (ii) verify that these statistics are well described by the SLE (in comparison to MD). Filling this gap in knowledge is the primary contribution of the work presented here.
The article is laid out as follows. In Section II the SLE are introduced ( §II B), a simple, yet robust, scheme for their numerical solution is proposed ( §II C), and its convergence characteristics demonstrated ( §II D). In Section III, numerical SLE solutions are verified against known analytical results and validated against MD calculations (introduced in §II A); first for initial (linear) instability growth ( §III A), and second for nonlinear growth of disturbances to the point of rupture ( §III B). In Section IV we use the SLE solver to (i) provide a deeper understanding of the impact of fluctuations on rupture dynamics and (ii) reach cases that will be computationally intractable for MD.

II. MODELS AND NUMERICS
In this section we introduce the stochastic lubrication equations and present a numerical method for their solution. Given the stochastic and nonlinear nature of the equations, we dedicate some time to demonstrating the scheme's robustness to increasing numerical resolution. First, though, we give details of the Molecular Dynamics simulations presented in this paper, which are used in Section III for validation of the numerical SLE solutions.

A. Molecular Dynamics
The MD simulations of this work are performed in LAMMPS [32] on nanoscale threads of water. The simulation box extends 10r 0 , 10r 0 , and L in the x, y and z directions, respectively, where r 0 is the initial thread radius and L is the thread length. The liquid thread is placed in the centre of the domain, and there are periodic boundary conditions imposed in all three directions.
The initial configuration is cut from a liquid bulk, created from an equilibrium NVT simulation with a Nosé-Hoover thermostat set at a specific temperature. The same ensemble and thermostat is used for the main simulations.
The Green Kubo method [33,34] is applied to calculate the dynamic viscosity by integration of the time-autocorrelation function of the off-diagonal elements of the pressure tensor where V bulk is the volume of the bulk fluid, k B is the Boltzmann constant, and T is temperature. The pressure tensor components are obtained using the definition of [35] and the angular brackets indicate the expectation.
The surface tension is calculated from the profiles of the components of the pressure tensor in a simple liquid-vapour system, using the mechanical definition [36]: where L z is the length of the MD domain, and subscripts n and t denote normal and tangential components, respectively. These methods are well known and have been validated by comparing MD results (µ and γ) with experimental data over a range of temperatures (see [37,38]). In the present work, liquid water is chosen because of its wide applications and its ability to create a large range of material properties [39] due to its stable liquid phase over a wide temperature range. The detailed properties (e.g. temperature, surface tension, viscosity) will be listed in the relevant sections. For the instability validation cases in Sec. III A (requiring long cylinders) a coarse-grained water molecule model, known as mW [40], is adopted (for computational efficiency); whereas for the breakup validation cases in Sec. III B, the TIP4P/2005 water model [41] is used (thus achieving a more accurate result). Selected MD realisations of both models are shown in Fig. 1.

B. Stochastic lubrication equations (SLE)
The SLE were derived by Moseler and Landman [1], who applied a lubrication approximation to the axisymmetric LLNS (also known as the stochastic Navier-Stokes equations) allowing the dynamics of the interface to be described by the thread radius h(z, t) and a velocity u(z, t) (see Fig. 2). Modified stochastic lubrication equations, where a dependence of the evaporation-condensation flux on the presence of a surrounding gas is incorporated into the continuity equation, have yielded results agreeing well with MD [23]. However, in this work we restrict our attention to the orginal SLE. To identify the governing dimensionless parameters, we use the following variables as scales of length, time, velocity and pressure, based on (but not confined to) a balance of inertial and surface-tension forces: whereh,t,ũ andp represent the dimensional interface height, time, velocity and pressure, respectively (note that the dimensional material parameters are not given tildes). The dimensionless SLE are written as follows: with the full Laplace pressure retained: where primes denote differentiation with respect to z, and N is a standard Gaussian random variable -a model for thermal fluctuations. The non-dimensional quantity Oh = µ/ √ ργr 0 is the Ohnesorge number, which relates the viscous forces to inertial and surface-tension forces (where µ is the liquid dynamic viscosity, γ is surface tension, and ρ is liquid density).
The Ohnesorge number is all that is needed to characterise the dynamics of free macroscopic threads, but here we obtain an additional dimensionless quantity: the thermal-fluctuation number, Th = l T /r 0 , to express the relative intensity of interface fluctuations, where l T = k B T /γ is the characteristic thermal fluctuation length and r 0 is the initial thread radius.
One of the purposes of the simulations presented in Sec. IV is to explore the combined influence of Oh and Th on different aspects of thread dynamics.
C. Numerical scheme for the SLE By construction, after enforcing the fluctuation-dissipation balance, the covariance of the stochastic term in Eq. (5) is where the presence of a Dirac delta function ensures infinitely small temporal/spatial correlation functions (i.e. the noise term is temporally and spatially uncorrelated). To represent N numerically, we introduce computer-generated random numbers, N t i , that are normally distributed with zero mean and unit variance. The delta function in Eq. (7) can be approximated by a 2D rectangular (boxcar) function (in t and z) that is non-zero over a time step ( t) and grid spacing ( z). The amplitude of the rectangular function, 1/( t z), is such that the integral properties of the delta function are preserved, i.e. ∞ −∞ ∞ −∞ δ(z, t) dz dt = 1 [43]. The complete noise term is thus discretized by Equation (8) provides robust and accurate numerical performance when used in conjunction with linear equations, e.g. one-dimensional LLNS [44] or for the linearised SLE. However, the full SLE are nonlinear (including the stochastic driving force: (hN ) /h 2 ), which creates stability issues that exacerbate as z and t become smaller and the amplitude of noise becomes larger (see Eq. (8)). Consequently, for some cases, it is impossible to achieve a spatially and temporally resolved result (i.e. one that converges as z → 0 and t → 0).
As a straightforward solution to this problem, we propose a numerical method where, beneath a certain scale, the noise becomes spatially and temporally correlated; thus remaining finite as z → 0 and t → 0. MD results show that this 'correlation scale' is much smaller than any scale of interest to the current paper, but we cannot be sure this will always be the case.
While this solution is largely pragmatic in nature, it actually reflects the physics better than uncorrelated noise. Motivated by these MD results, we introduce a correlation time scale, T c , and correlation length scale, L c , into our SLE simulations. These correlation scales must be larger than the time step and grid spacing, respectively. Beneath the correlation scales, a simple linear interpolation between uncorrelated random noise at the end points of the correlation interval is applied (as illustrated in Fig. 3 b for temporal noise). The uncorrelated noise (which is interpolated between) has a mean of zero and a variance 1/(T c L c ). See Appendix B for a derivation.
Note, a spatially correlated, but temporally uncorrelated, noise model has been developed by Grün [26] for the stochastic thin-film equation. The approach we present above has the advantage of simplicity, compares well to the results of Grün's model and has the capability to produce spatio-temporally correlated noise as well.
In order to solve the full nonlinear SLE, we use the MacCormack method [45], a simple second-order explicit finite difference scheme in both time and space. The solution at each time level is defined by two arrays, Here, n is the number of mesh points. The time-derivative terms are approximated by (h t+1 The numerical method proceeds in two steps: a predictor step: and a corrector step: where u t+1 i and h t+1 i are "provisional" values at time level t+1, and F represents all the partial spatial derivative terms on the right-hand side (expressions for F are listed in Appendix C).

D. Time-step and grid-size convergence
In order to test the integrity of the SLE numerical approach introduced above, we consider the simulation of a short thread (L = 10, Oh = 1.07 and Th = 0.11) with an increasingly fine time step and grid spacing (note, for this case, a model using uncorrelated noise would not converge).

III. NUMERICAL VERIFICATION AND VALIDATION
Having demonstrated the convergence, the numerical solutions are verified and validated by the analytical models and MD simulations at both the linear stage for instability ( § III A) and the nonlinear stage for rupture ( § III B).

A. Linear instability and thermal capillary waves
In this section, the analytical SLE-RP [25] has been employed as a benchmark for the numerical solutions of the SLE introduced above (Sec. II C). The SLE-RP can be written in a dimensionless form as follows [see Appendix D for derivation]: Here, k is the dimensionless wave number, L is the dimensionless length of the thread, and R i is the initial model disturbance (a complex random variable with zero mean).
We also perform MD simulations for comparison, where we adopt a coarse-grain water molecule model, mW [40], to limit the computational resources required for such long threads; in all cases, L = 100. The inital radius r 0 and temperature T is selected to obtain specific Oh and Th (shown in Table I   The results in Fig. 6 a show a modal distribution (spectrum) that varies with time. For small wavenumbers (k < 1), the spectrum becomes sharper with time, while the spectrum at high wavenumbers (k >∼ 2) is static over these time scales; i.e., it quickly reaches its asymptotic limit. This limit can be obtained from Eq. (11), for k > 1, by taking t → ∞: This is consistent with the theory for thermal capillary waves in thin-film flows [46,47],  is adopted, with liquid properties as listed in Table II. Here, L = 12.
Our first comparison, Fig. 7, is for the time evolution of the minimum (over z) thread    surface tension (note we vary a dimensional quantity here, in order to connect most transparently with the assumption in Eggers' work). Notably, when surface tension is stronger, the breakup is faster than the analytical prediction, which suggests that the destabilising effect of surface tension can also contribute to the thread dynamics near to rupture. This limit of applicability might also explain the deviation between Eggers' similarity solution and MD results in previous studies [49,50].
Although the agreement for h min with Eggers' similarity solution is good for low γ, we were unable to find any agreement between either MD or SLE with the associated universal profile. The reason is currently unclear and should be the subject of future investigation.
The ensemble-averaged profiles plotted in Figure 9 show good overall agreement between the MD and SLE for three time instances leading to rupture. The limitation of bin sizes in the MD data prevents a more detailed comparison of the profile shape. In particular, it is not clear whether the finer features seen in the SLE (namely, the V-notch -or 'widow's peak' -near to the minimum) is physical because these local features reach the molecular scale and cannot be reliably extracted from the MD.

IV. EXPLOITING SLE: RUPTURE BEYOND MD
Having established its predictive capability, in this section we use the nonlinear SLE solver to further explore the impact of fluctuations on rupture dynamics, over a broader range of conditions than have been studied previously and for cases that are too computationally demanding to consider with MD. We start, in §IV A, by exploring the shape of the thread at rupture whilst in §IV.B, we focus on the time evolution of the point of the thread's minimum thickness.

A. Rupture profiles
Moseler and Landmann [1] were the first to demonstrate, using MD, that thermal fluctuations could lead to a symmetric double-cone rupture profile, and that SLE solutions were also able to capture this (whereas deterministic equations cannot). We reproduce this result for Thread 5 (see Table II for parameters) in Figure 10. in a far broader range of conditions than accessible to MD; as is illustrated in Figure 11.  estimated (the red dashed lines). Notably, it is the most probable profile which Eggers [2] computes from the Fokker-Planck equation for the SLE. In the cases considered here there is little difference between the mean and the estimated most-probable profile.
The bottom right-hand corner profile in Fig. 12  it is, seemingly, a non-trivial interplay of effects as we would expect when inertia, viscosity, surface tension and fluctuations all play a role.
The next most important observation is that, for low Oh, the impact of noise results in an asymmetric mean thread profile at rupture (see top right-hand corner image of Fig. 12).
The double-cone profile observed in Moseler and Landman [1] is not observed here. Instead, we see a quite distinct rupture shape (a drop and funnel), on average, which looks more like typical rupture profiles seen macroscopically when a drop breaks off from a thread. This behaviour is not surprising, as the chance of two points pinching off at precisely the same instance becomes slim when we have fluctuations (and indeed this kind of 'perfect' pinch off is difficult to reproduce experimentally macroscopically as well). We stress that these flow conditions are not accessible by our MD simulations at present; the SLE calculations are essential to provide this insight.

B. Evolution of minimum thread radius
In the classical picture there is the potential for multiple transitions between distinct 'dynamic regimes' (defined by Oh) leading to rupture [9,15,16]. The three main regimes, described in Section I, are the viscous regime (V-regime), the inertial regime (I-regime), and the universal viscous-inertial regime (VI-regime). These regimes are characterised by a power-law (linear for the first and third) evolution of minimum thickness with time to rupture, at rates given by various analytical results [9,15].
On top of this already complex situation, thermal fluctuations can introduce yet another regime (here referred to as the F-regime), which generates non-linear (power-law) evolution of minimum thread radius. For moderate Th (>∼ 0.1) and non-negligible Oh, fluctuations appear to dominate the entire thread evolution (see, e.g., the non-linear evolution in Fig. 7 a).
However, at lower Th, we can observe transitions from the classical behaviour to one that is fluctuation dominated as the rupture process progresses.
In Fig. 13 we compare the SLE with the classical model (LE) and various analytical results for fixed Th = 0.02 at Oh = 0.01, 1 and 100. Experimentally, this corresponds to using fluids of a range of different viscosities (with fixed surface tensions) for the same breakup configuration. As was the case for the rupture profiles presented above, at low Oh ( Fig.13 a there is seemingly no impact of fluctuations on the time evolution of h min (i.e. V-regime come from Refs. [14] and [12], respectively. The power law solution in the F-regime for (c) is obtained from Ref. [2].
there is no discernible difference between the SLE and LE). For larger Oh (see Fig. 13 c), however, a clear transition between macroscopic and a fluctuation-dominated regime can be observed. At early times the evolution is described by a linear time dependence, derived by Papageourgiou [12,13] for the V-regime; at later times (in the F-regime) the evolution matches the power law proposed by Eggers [2] and greatly accelerates the breakup process.
Of course, of the two numerical methods presented in the figure, only the SLE can capture both. When Oh = 1 (Fig. 13 b), the dynamics become more complicated. The LE predicts a transition from the V-regime to the VI-regime in the final stages, which has been proved experimentally [15] and numerically [9] at macroscopic scales. However, this transition does not occur in the presence of thermal fluctuations, according to the SLE. Instead there exists a similar transition from the V-regime to a new regime as was the case for the large Oh case, but the power-law exponent does not match that found by Eggers (which, as explored in Section 3, is possibly due to the assumption in his analytic treatment that surface tension is unimportant in the final stages leading to rupture).  Oh (or Th). However, we would need to get many decades of h min to determine the precise crossover point, which is not available from our current simulations.

V. DISCUSSION AND FUTURE DIRECTIONS
In this work, a numerical solver of the SLE has been developed with a new simple scheme proposed for the noise term. Based on validation from MD for both instability and the rupture of liquid nano-threads, this solver is demonstrated to be a powerful tool for studying the interface dynamics of nano-threads; and operating over a thousand times faster than MD.
Furthermore, it allows us to operate in the regions of parameter space where analytic models are outside their limits of applicability and MD is impractical either due to (i) exorbitant computational cost or (ii) limits in the molecular properties available from known potentials.
Whilst this article provides new understanding of interface dynamics, it opens up several new avenues of enquiry, that we discuss here.

A. Correlation Scales
In this article the use of a correlation scale is motivated from two angles. First, the computational SLE scheme is seen to be unable to converge unless these are introduced, with huge spikes on the free surface observed that seemingly prematurely rupture the thread and/or destroy the numerical accuracy. Second, MD suggests that correlation scales exist, and as one may expect these are typically on the molecular scale. These issues motivate a number of different questions. From a modelling viewpoint, the incorporation of molecular correlation scales within a continuum model should be treated with caution, and thus one may interpret a continuum limit as when the correlation length goes to zero. However, in some cases the correlations seem to have a profound effect on thin film dynamics simulations and experiments [28], so that these issues are far from trivial. In terms of numerical and then take the correlation scale to zero (as one may expect for the continuum limit)?
Such questions are related to the development of robust and efficient numerical schemes for SPDE problems. Here, we focus on simplicity, both with the use of an explicit time-stepping scheme and in the linear interpolation of noise. More complex schemes exist, where the noise is represented in terms of appropriate basis functions [26,29], but although these are more mathematically rigorous approaches, we found them to give the same dynamic behaviours at increased computational cost. Clearly there is scope for more work in this direction, particularly as one considers the possibility of developing 2D schemes, as even though the SLE is much cheaper than MD, it is still more computationally burdensome than deterministic methods (as a minimum, due to the requirement for ensembles)

B. Similarity Solutions
By considering a liquid with sufficiently small surface tension, we were able to recover the power law predicted by Eggers' similarity solution [2], as previously also identified in specially designed experiments with colloid-polymer mixtures [20]. However, no agreement was obtained between the SLE and similarity solution for the universal profiles predicted in [2] and surface tension was seen to influence the power law even at physical values, as seen in Figures 8 and 13. Therefore, it remains an open problem to derive similarity solutions for the breakup that incorporate surface tension, building on the new framework, considering the most probable breakup in a stochastic process, as developed by Eggers for this class of flows.

C. Transition prediction
Liquid thread rupture is a multiscale phenomenon with complicated transitions between different regimes. Numerical solutions in Fig. 13 and 14 show a transition from the V-regime to fluctuation-dominated regimes. However, it remains unclear whether there exists a scaling law between the crossover (transition) point and Oh. To answer this question, we need to develop more accurate and efficient schemes (e.g. higher-order schemes and implicit time marching methods) to capture many decades of dynamics, which could be the subject of future work.

D. Experimental analysis
In nanoscale breakup phenomena, one not only has small spatial scales, but also small temporal ones for the problems of interest. This is in contrast to thin film dynamics, where typical time scales are macroscopic when highly viscous films are considered [52,53] and therefore experimental analysis that temporally resolves features becomes possible. For the problems considered in this article, it seems most likely that experimental verification would come first from the ultra-low surface tension liquids developed in [20,21] that make Th moderate even at the microscale where one can perform imaging. There are many potential directions for experimental analysis to take, but a starting point would be to more carefully consider the rupture profiles and scaling of h min to see how this compares to our predictions.

Appendix A: Molecular models
In this section, we introduce the two molecular models used in this work: (i) mW [40] and (ii) TIP4P/2005 [41].
In mW model, the hydrogen-bonded structure of water is mimicked through the introduction of a non-bond angular dependent term that encourages tetrahedral configurations.
The model contains two terms: (i) φ ij depending on the distances between pairs of atoms (represented by r ij and s ik ); and (ii) φ ijk depending on the angles formed by triplets of atoms (represented by θ ijk ). The full expression is given by where is the depth of the potential well, σ represents the particle diameter, A, B, p, q, χ and κ respectively give the form and scale to the potential, and θ 0 represents the tetrahedral angles. All the parameters are presented in Table III.
where ε 0 represents the vacuum permittivity, q i and q j are the atomic charges. We list all the parameters for different atoms in Table IV.
Appendix B: Derivation for the noise variance In this section, a mathematical procedure is proposed to derive the variance of our numerical correlated noise model in the main text. To simplify the problem further, we only use temporal fluctuations in this example. The variance of the fluctuation is where A f is the noise amplitude. The mean of the noise, N , over a period T N is where the variance of N can be obtained by The interpolated noise proposed in the main article (illustrated in Fig. 3 b) can be written as where X i are normal distributed random numbers with zero mean, Π is the hat function, and τ i is defined as The term N can be calculated by When T N T c , this is approximately: so the variance of N is then: Because Var(N ) should equal the theoretical result in Eq.(B1), we can obtain that: A similar procedure can be applied to the spatial noise, and can be used to explain why we use 1/(T c L c ) to replace the 1/( t z) in Sec. II B.
Appendix C: F for the MacCormack method Two differential operators, f and b are introduced to represent forward and backward difference respectively: F is discretized by the forward difference for the predictor step, written as while the backward method is applied for F: (C3)  R is both a random and complex variable with zero mean. Note, R LE is also random as it develops from a random initial condition, but is uncorrelated with R fluc (both have zero mean). So the root mean square (rms) of R is sought, which from Eq.(D3) is given by where |R LE | 2 can easily be obtained from Eq.(D4).