2 Mathematical modelling

2.1 Stochastic lubrication equation

Consider a cylindrical polar coordinate system (figure 1), with $z$ -axis along the centre line of an axisymmetric initially cylindrical liquid volume of length $L$ and radius $r_{0}$ . In the lubrication approximation, the dynamics is described by the radius of the cylinder $r=h(z,t)$ and the axial velocity $u(z,t)$ , so that the SLE (Moseler & Landman Reference Moseler and Landman2000) is given by

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}u=-uu^{\prime }-\frac{p^{\prime }}{\unicode[STIX]{x1D70C}}+3\unicode[STIX]{x1D708}\frac{(h^{2}u^{\prime })^{\prime }}{h^{2}}+A\frac{(h{\mathcal{N}})^{\prime }}{h^{2}}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}h=-uh^{\prime }-u^{\prime }h/2, & \displaystyle\end{eqnarray}$$

with the full Laplace pressure retained

(2.3) $$\begin{eqnarray}p=\unicode[STIX]{x1D6FE}[h^{-1}(1+(h^{\prime })^{2})^{-1/2}-h^{\prime \prime }(1+(h^{\prime })^{2})^{-3/2}],\end{eqnarray}$$

and primes denoting differentiation with respect to $z$ . In these equations, $\unicode[STIX]{x1D708}$ is the liquid’s kinematic viscosity, $\unicode[STIX]{x1D70C}$ is its density, $\unicode[STIX]{x1D6FE}$ is the surface-tension coefficient and ${\mathcal{N}}$ denotes a standard Gaussian (white) random variable, capturing the thermal fluctuation effects, where $A=\sqrt{6k_{B}T\unicode[STIX]{x1D708}/(\unicode[STIX]{x1D70C}\unicode[STIX]{x03C0})}$ is a dimensional parameter, $k_{B}$ is the Boltzmann constant and $T$ is the liquid temperature. If there are no fluctuations (i.e. $A=0$ ) the classical lubrication equations (LE) (Eggers & Dupont Reference Eggers and Dupont1994) are recovered.

Figure 1. Schematic of the cylinder instability.

2.2 Stability analysis

For the linear stability analysis, we take $h(z,t)=r_{0}[1+\hat{r}(z,t)]$ , with $\hat{r}(z,t)\ll 1$ a dimensionless perturbation. Substituting this into the SLE and linearising gives

(2.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}\hat{r}}{\unicode[STIX]{x2202}t^{2}}+\frac{\unicode[STIX]{x1D6FE}}{2\unicode[STIX]{x1D70C}}\left(\frac{\hat{r}^{\prime \prime }}{r_{0}}+r_{0}\hat{r}^{\prime \prime \prime \prime }\right)-3\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}\hat{r}^{\prime \prime }}{\unicode[STIX]{x2202}t}=-\frac{A}{2r_{0}}{\mathcal{N}}^{\prime \prime }.\end{eqnarray}$$

Note, this linearisation also implies an extent to the strength of fluctuations, which we will discuss later. A finite Fourier transform is applied to (2.4) to obtain

(2.5) $$\begin{eqnarray}\frac{\text{d}^{2}R}{\text{d}t^{2}}+\unicode[STIX]{x1D6FC}\frac{\text{d}R}{\text{d}t}+\unicode[STIX]{x1D6FD}R=\frac{A}{2r_{0}}k^{2}N,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}=3\unicode[STIX]{x1D710}k^{2}$ and $\unicode[STIX]{x1D6FD}=(\unicode[STIX]{x1D6FE}/2\unicode[STIX]{x1D70C})(r_{0}k^{4}-(k^{2}/r_{0}))$ and the transformed variables are defined as follows:

(2.6a,b ) $$\begin{eqnarray}R(k,t)=\int _{0}^{L}\hat{r}(z,t)\text{e}^{-\text{i}kz}\,\text{d}z\quad \text{and}\quad N(k,t)=\int _{0}^{L}{\mathcal{N}}(z,t)\text{e}^{-\text{i}kz}\,\text{d}z.\end{eqnarray}$$

The solution of (2.5) is linearly decomposed into two parts:

(2.7) $$\begin{eqnarray}R=R_{LE}+R_{fluc}.\end{eqnarray}$$

The first part is the solution to the homogenous form of (2.5) (i.e. with $A=0$ ) with some stationary initial disturbance (i.e. $R=R_{i}$ and $\text{d}R/\text{d}t=0$ at $t=0$ ). The solution to the homogeneous ordinary differential equation is straightforward to obtain:

(2.8) $$\begin{eqnarray}R_{LE}=R_{i}\text{e}^{-at/2\unicode[STIX]{x1D70F}}\left[\cosh (ct/2\unicode[STIX]{x1D70F})+\frac{a\sinh (ct/2\unicode[STIX]{x1D70F})}{c}\right]\!,\end{eqnarray}$$

where

(2.9a,b ) $$\begin{eqnarray}a=3(kr_{0})^{2}\left(\frac{\ell _{\unicode[STIX]{x1D708}}}{r_{0}}\right)^{1/2}\quad \text{and}\quad c=2\left[\frac{9}{4}\frac{\ell _{\unicode[STIX]{x1D708}}}{r_{0}}(kr_{0})^{4}+\frac{(kr_{0})^{2}}{2}(1-(kr_{0})^{2})\right]^{1/2},\end{eqnarray}$$

and characteristic flow scales for time $\unicode[STIX]{x1D70F}=(\unicode[STIX]{x1D70C}r_{0}^{3}/\unicode[STIX]{x1D6FE})^{1/2}$ and length $\ell _{\unicode[STIX]{x1D708}}=\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}^{2}/\unicode[STIX]{x1D6FE}$ have been introduced. This is a solution to the standard lubrication equations (there is no fluctuating component), and is thus denoted $R_{LE}$ in (2.7). Note, for the rest of this paper, the initial modal disturbance $R_{i}$ is treated as a complex random variable with zero mean.

The second component of the solution arises from solving the full form of (2.5) with zero initial disturbance; this part of the solution is solely due to fluctuations, and is thus denoted $R_{fluc}$ . This is obtained by the homogeneous equation’s impulse response,

(2.10) $$\begin{eqnarray}H(k,t)=2\unicode[STIX]{x1D70F}\text{e}^{-at/2\unicode[STIX]{x1D70F}}\sinh (ct/2\unicode[STIX]{x1D70F})/c,\end{eqnarray}$$

which due to the linear, time-invariant nature of the system, allows us to write

(2.11) $$\begin{eqnarray}R_{fluc}=\frac{Ak^{2}}{2r_{0}}\int _{0}^{t}N(k,t-{\mathcal{T}})H(k,{\mathcal{T}})\,\text{d}{\mathcal{T}}.\end{eqnarray}$$

The modal amplitude $R\,(=R_{LE}+R_{fluc})$ is a complex random 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, in order to obtain information on how disturbances associated with each wavenumber develop in time, allowing the identification of unstable and fastest growing modes, the root mean square (r.m.s.) of $|R|$ is sought (equivalent to the standard deviation of $|R|$ ):

(2.12) $$\begin{eqnarray}|R|_{rms}=\sqrt{\overline{|R_{LE}+R_{fluc}|^{2}}}=\sqrt{\overline{|R_{LE}|^{2}}+\overline{|R_{fluc}|^{2}}},\end{eqnarray}$$

where, from (2.8),

(2.13) $$\begin{eqnarray}\overline{|R_{LE}|^{2}}=\overline{|R_{i}|^{2}}\text{e}^{-at/\unicode[STIX]{x1D70F}}\left[\cosh (ct/2\unicode[STIX]{x1D70F})+\frac{a\sinh (ct/2\unicode[STIX]{x1D70F})}{c}\right]^{2},\end{eqnarray}$$

and since $N$ is uncorrelated Gaussian white noise, and the variance of the norm of the white noise $\overline{|N|^{2}}=L$ , equations (2.10) and (2.11) combine to give:

(2.14) $$\begin{eqnarray}\displaystyle \overline{|R_{fluc}|^{2}} & = & \displaystyle \left(\frac{Ak^{2}}{2r_{0}}\right)^{2}\overline{\left|\int _{0}^{t}N(k,t-{\mathcal{T}})H(k,{\mathcal{T}})\,\text{d}{\mathcal{T}}\right|^{2}},\nonumber\\ \displaystyle & = & \displaystyle \left(\frac{Ak^{2}}{2r_{0}}\right)^{2}\left(\int _{0}^{t}\overline{|N(k,t-{\mathcal{T}})|^{2}}H(k,{\mathcal{T}})^{2}\,\text{d}{\mathcal{T}}\right),\nonumber\\ \displaystyle & = & \displaystyle \left(\frac{Ak^{2}}{2r_{0}}\right)^{2}L\int _{0}^{t}H^{2}\,\text{d}{\mathcal{T}},\nonumber\\ \displaystyle & = & \displaystyle \ell _{fluc}^{2}b\frac{(a^{2}-c^{2})-a^{2}\cosh (ct/\unicode[STIX]{x1D70F})-ac\sinh (ct/\unicode[STIX]{x1D70F})+c^{2}\text{e}^{at/\unicode[STIX]{x1D70F}}}{ac^{2}(a^{2}-c^{2})\text{e}^{at/\unicode[STIX]{x1D70F}}},\end{eqnarray}$$

where non-dimensional $b=(3/\unicode[STIX]{x03C0})(L/r_{0})(kr_{0})^{4}(\ell _{\unicode[STIX]{x1D708}}/r_{0})^{1/2}$ , and the thermal capillary length $\ell _{fluc}=(k_{B}T/\unicode[STIX]{x1D6FE})^{1/2}$ gives the characteristic length scale of the fluctuations.

2.3 Convergence to the classical model

From (2.12), fluctuations can be seen to be negligible when the thermal capillary length is very much shorter than the initial modal amplitude; i.e. $R\rightarrow R_{LE}$ as $\ell _{fluc}^{2}/\overline{|R_{i}|^{2}}\rightarrow 0$ . We refer to this classical limit as the LE-RP model (as distinct from the SLE-RP). In order to compare with standard stability analysis, we want to find the fastest growing mode at a long time after any initial disturbance. As such we take $t\rightarrow \infty$ to give:

(2.15) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \overline{|R_{LE}|^{2}}\rightarrow \overline{|R_{i}|^{2}}\text{e}^{(c-a)t/\unicode[STIX]{x1D70F}}\left(\frac{c+a}{2c}\right)^{2},\\ \displaystyle \overline{|R_{fluc}|^{2}}\rightarrow \frac{\ell _{fluc}^{2}\,b}{2ac^{2}(a^{2}-c^{2})}[2c^{2}-\text{e}^{(c-a)t/\unicode[STIX]{x1D70F}}(a^{2}+ac)],\end{array}\right\} & & \displaystyle\end{eqnarray}$$

with $c-a\geqslant 0$ for $kr_{0}\leqslant 1$ (which is the case as $t\rightarrow \infty$ ). Therefore, a functional form of (2.12) is:

(2.16) $$\begin{eqnarray}{\mathcal{R}}(k,t)={\mathcal{F}}_{1}(k)\text{e}^{{\mathcal{G}}(k)t}+{\mathcal{F}}_{2}(k),\end{eqnarray}$$

where

(2.17) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\mathcal{G}}(k)=(c-a)/\unicode[STIX]{x1D70F},\\ \displaystyle {\mathcal{F}}_{1}(k)=\overline{|R_{i}|^{2}}\left(\frac{c+a}{2c}\right)^{2}+\frac{\ell _{fluc}^{2}\,b}{2c^{2}(c-a)},\\ \displaystyle {\mathcal{F}}_{2}(k)=\frac{\ell _{fluc}^{2}\,b}{a(a^{2}-c^{2})}.\end{array}\right\}\end{eqnarray}$$

In order to find the maximum of ${\mathcal{R}}(k,t)$ , which defines the fastest growing mode $k=k_{max}$ , we calculate the derivative of (2.16) with respect to $k$ to obtain

(2.18) $$\begin{eqnarray}\frac{1}{t\text{e}^{{\mathcal{G}}t}}\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}k}=\frac{1}{t}\frac{\text{d}{\mathcal{F}}_{1}}{\text{d}k}+\frac{\text{d}{\mathcal{G}}}{\text{d}k}{\mathcal{F}}_{1}+\frac{1}{t\text{e}^{{\mathcal{G}}t}}\frac{\text{d}{\mathcal{F}}_{2}}{\text{d}k},\end{eqnarray}$$

and set $\unicode[STIX]{x2202}{\mathcal{R}}/\unicode[STIX]{x2202}k=0$ . As $1/t$ and $\text{e}^{-{\mathcal{G}}t}$ vanish as $t\rightarrow \infty$ , and ${\mathcal{F}}_{1}(k_{max})\neq 0$ , the equation for determining $k_{max}$ is simply

(2.19) $$\begin{eqnarray}\left.\frac{\text{d}{\mathcal{G}}}{\text{d}k}\right|_{k=k_{max}}=\frac{1}{\unicode[STIX]{x1D70F}}\left.\frac{\text{d}(c-a)}{\text{d}k}\right|_{k=k_{max}}=0,\end{eqnarray}$$

which is the same as found by Eggers & Dupont (Reference Eggers and Dupont1994) who neglected fluctuations entirely. However, as break-up occurs in a finite time, both terms in (2.12) could play a role in determining $k_{max}$ at any instant, with the second term becoming more important as $\ell _{fluc}^{2}/\overline{|R_{i}|^{2}}$ increases (all else being constant).

3 Simulations

To test our hypothesis, MD simulations (in figure 2) are performed in LAMMPS (Plimpton Reference Plimpton1995) on long cylinders $L/r_{0}=160$ of three different radii: Cylinder 1 ( $r_{0}=5.76$  nm, $2.1\times 10^{6}$ particles), Cylinder 2 ( $r_{0}=2.88$  nm, $2.8\times 10^{5}$ particles) and Cylinder 3 ( $r_{0}=$ 1.44 nm, $4.6\times 10^{4}$ particles). The simulation box ( $57$ nm $\times$ $57$ nm $\times \,L$ in the $x$ , $y$ and $z$ directions, respectively) has periodic boundary conditions imposed in all directions and is filled with Lennard-Jones (LJ) fluid whose atomic interactions are modelled using the standard 12-6 pair potential with the cutoff distance $5.5\unicode[STIX]{x1D70E}$ ; interaction potentials for argon have been employed (Lee, Barker & Pound Reference Lee, Barker and Pound1974). The initial configuration is created from equilibrium simulations of separate liquid-only and vapour-only systems (using a canonical ensemble (NVT) with a Nosé–Hoover thermostat) with respective densities of $1398~\text{kg}~\text{m}^{-3}$ and $3.22~\text{kg}~\text{m}^{-3}$ , which correspond to the saturated liquid and vapour densities at a temperature of 84.09 K (Fernandez, Vrabec & Hasse Reference Fernandez, Vrabec and Hasse2004). The same ensemble and thermostat is used for the full simulations.

The kinematic viscosity is found using the Green–Kubo method (Green Reference Green1954; Kubo Reference Kubo1957) to be $\unicode[STIX]{x1D708}=1.76\times 10^{-7}~\text{m}^{2}~\text{s}^{-1}$ . This technique calculates dynamic viscosity by integration of the time-autocorrelation function of the off-diagonal elements of the pressure tensor $P_{ij}$ so that

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D707}=\frac{V_{bulk}}{k_{B}T}\int _{0}^{\infty }\left\langle P_{ij}(t)\cdot P_{ij}(0)\right\rangle \,\text{d}t\quad (i\neq j),\end{eqnarray}$$

where the pressure tensor is obtained using the definition of Kirkwood & Buff (Reference Kirkwood and Buff1949) and the brackets indicate the expectation value. $V_{bulk}$ represents the volume of the bulk fluid under consideration. The surface tension is calculated from the profiles of the components of the pressure tensor in a simple liquid–vapour system (Trokhymchuk & Alejandre Reference Trokhymchuk and Alejandre1999):

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}=\frac{1}{2}\int _{0}^{L_{z}}[P_{n}(z)-P_{t}(z)]\,\text{d}z,\end{eqnarray}$$

where $L_{z}$ is the length of the MD domain, and subscripts ‘ $n$ ’ and ‘ $t$ ’ denote normal and tangential components, respectively. We have found that to obtain results that are close to experimental data (Trokhymchuk & Alejandre Reference Trokhymchuk and Alejandre1999), a cutoff distance significantly larger than the popular choice of 2.5 $\unicode[STIX]{x1D70E}$ is required; a choice of 5.5 $\unicode[STIX]{x1D70E}$ (used throughout) has proved to be sufficiently accurate. Finally, we have $\unicode[STIX]{x1D6FE}=1.42\times 10^{-2}~\text{N}~\text{m}^{-1}$ and $\ell _{fluc}=0.2859$  nm.

Figure 2. Molecular dynamics simulation of the Rayleigh–Plateau instability showing a liquid cylinder (Cylinder 1) breaking into droplets.

4 Results and discussion

To gather statistics, multiple independent MD simulations (realisations) (Cylinder 1: $30$ , Cylinder 2: $45$ and Cylinder 3: $100$ ) are performed. For each realisation, a discrete Fourier transform of the interface position (which is extracted from axially distributed annular bins based on a threshold density) is performed and then an ensemble average at each time allows us to produce the results shown in figure 3 (dashed lines) and figure 4 (red dots). Using the initial conditions from MD to extract $\overline{|R_{i}|^{2}}$ , remarkably good agreement with the SLE-RP is obtained, giving us confidence that our approach captures the essential physics.

Figure 3. The r.m.s. of dimensionless modal amplitude versus dimensionless wavenumber; a comparison of ensemble-averaged MD data (dashed lines) and (2.12) (solid lines) at various time instances. (a) Results of Cylinder 1; the inset shows selected MD realisations; (b) Results of Cylinder 2, the three timesteps are: $0.23$ (red), $0.47$ (blue), $0.70$ (black)/ns; (c) Results of Cylinder 3, the three timesteps are: $0.062$ (red), $0.125$ (blue), $0.100$ (black)/ns.

The MD results in figure 3 illustrate that there exists a modal distribution which varies with time, becoming sharper as it evolves. Extracting $k_{max}$ from these data yields the plots in figure 4, which confirm that $k_{max}$ tends to the Eggers & Dupont result as $t\rightarrow \infty$ . However, figure 4 also shows that $k_{max}$ at the average break-up time (which ultimately determines drop size) is consistently overpredicted by Rayleigh’s inviscid result, as seen in previous MD, and underpredicted by the Eggers & Dupont model (valid across all values of viscosity) – particularly for the smallest radius (Cylinder 3) where $k_{max}=0.52/r_{0}$ in the MD and $k_{max}=0.35/r_{0}$ from Eggers & Dupont. The modal analysis based purely on the LE-RP also underpredicts the MD data and fails to exhibit the dominant short wavelength modes observed at early times. In contrast, the SLE-RP curves in both figures 3 and 4 give excellent agreement with the MD and underline the critical role of thermal fluctuations in the instability mechanism at the nanoscale. Note, however, that the agreement between the SLE-RP and the MD is less good for the smallest cylinder. This can potentially be explained by the decreasing validity with decreasing scale of the weak-fluctuation assumption in our SLE-RP analysis. Crudely, $l_{fluc}/r_{0}\ll 1$ for the linear analysis to be valid. Here, $l_{fluc}/r_{0}=0.05$ , 0.1 and 0.2, for the three cylinders, listed largest to smallest.

Figure 4. Evolution in time $t$ of the wavenumber with greatest amplitude ( $k_{max}$ ) for (a) Cylinder 1 (b) Cylinder 2 and (c) Cylinder 3. Red dots and solid lines are the maximum predicted by MD (interpolation) and the SLE-RP respectively. Average break-up time, $\bar{t}_{b}$ , is obtained from the MD data.

Intriguingly, the early time behaviour in figure 4 indicates that $k_{max}r_{0}$ can be greater than unity, violating the classical stability criterion of Plateau. Therefore, it seems possible that ‘fat’ cylinders, whose length is below the classical critical stability ( $L<L_{crit}=2\unicode[STIX]{x03C0}r_{0}$ ), may be unstable in the presence of fluctuations. To test the hypothesis, we consider (2.12) at the critical point, i.e. when $kr_{0}=1$ to obtain

(4.1) $$\begin{eqnarray}\overline{|R|^{2}}=\overline{|R_{i}|^{2}}+\frac{L}{4\unicode[STIX]{x03C0}r_{0}}\ell _{fluc}^{2}\left(\frac{\ell _{\unicode[STIX]{x1D708}}}{r_{0}}\right)^{1/2}\frac{6at/\unicode[STIX]{x1D70F}+3\text{e}^{-2at/\unicode[STIX]{x1D70F}}(4\text{e}^{at/\unicode[STIX]{x1D70F}}-1)-9}{a^{3}}.\end{eqnarray}$$

Notably the contribution from LE equations (the first term on the right-hand side of (4.1)) is a constant; so, according to the classical model, the initial disturbance neither decays or grows. Hence, it is critically stable. However, the second term (due purely to fluctuations) grows in proportion to $t$ as $t\rightarrow \infty$ , giving a potential mechanism for break-up. This suggests that cylinders of the critical length, and perhaps shorter, are likely to be unstable at the nanoscale. To verify these conclusions, we perform a further series of MD experiments for cylinders of two radii ( $r_{0}=1.44,2.88$  nm) slightly shorter than the critical length $L_{crit}$ so that all classically unstable (long) wavelengths are suppressed by the domain size. This has been performed using two different flow configurations, one in which periodic boundary conditions are applied and the other in which the liquid cylinder is bounded at each end by a solid wall.

The four cases in figure 5 (cylinders of two radii for each end condition) all have length $L=6r_{0}<L_{crit}$ , and thus satisfy Plateau’s stability criterion. However, all break-up in a finite time, supporting our conclusion from SLE-RP that fat cylinders can indeed become unstable at the nanoscale. Notably, the break-up shapes resemble the double-cone profiles first observed by Moseler & Landman (Reference Moseler and Landman2000).

Having established the possibility of violating the Plateau criterion at the nanoscale, in figure 6 we show the average break-up time of such cylinders using 50 independent MD simulations for each data point (the standard deviation is indicated). There are two intuitive observations: (i) for the smaller radius cylinder, the break-up (which is partly or wholly due to fluctuations) occurs significantly faster than for the larger radius cylinder; (ii) as the aspect ratio of the cylinder decreases (i.e. becomes fatter) and crosses the classical stability limit ( $L_{crit}/L>1$ ), the average break-up time increases dramatically, and so does the variance in the measured times. This is because, at lower aspect ratios, the stabilising effect of surface tension becomes stronger and one has to wait longer (on average) for the ‘perfect storm’ of fluctuations to arrive that will overcome these and rupture the cylinder. This could explain why previous MD (Min & Wong Reference Min and Wong2006) appears to support the classical criterion: to violate the Plateau stability one must either be close to $L_{crit}/L=1$ or potentially wait a long time. Notably, while this is a ‘long time’ in MD simulations, from the macroscopic perspective the time scales on which classical stability is lost are miniscule.

Figure 5. Selected realisations for the break-up of classically stable cylinders (i.e. those satisfying the Plateau stability criterion). The two simulations on the left satisfy periodic boundary conditions, those on the right are bounded by a wall (in blue). Non-dimensional time $t^{\ast }=t/\unicode[STIX]{x1D70F}$ .

Figure 6. The non-dimensional break-up time ( $t_{b}^{\ast }=t_{b}/\unicode[STIX]{x1D70F}$ ) of short nanocylinders near the classical stability boundary, obtained from MD. $L_{crit}/L=2\unicode[STIX]{x03C0}r_{0}/L=kr_{0}$ . Error bars indicate one standard deviation of $t_{b}^{\ast }$ either side of the mean.