Hostname: page-component-7c8c6479df-7qhmt Total loading time: 0 Render date: 2024-03-29T00:29:39.126Z Has data issue: false hasContentIssue false

Revisiting the Rayleigh–Plateau instability for the nanoscale

Published online by Cambridge University Press:  07 January 2019

Chengxi Zhao
Affiliation:
School of Engineering, University of Warwick, Coventry CV4 7AL, UK
James E. Sprittles*
Affiliation:
Mathematics Institute, University of Warwick, Coventry CV4 7AL, UK
Duncan A. Lockerby*
Affiliation:
School of Engineering, University of Warwick, Coventry CV4 7AL, UK
*
Email addresses for correspondence: J.E.Sprittles@warwick.ac.uk, D.Lockerby@warwick.ac.uk
Email addresses for correspondence: J.E.Sprittles@warwick.ac.uk, D.Lockerby@warwick.ac.uk

Abstract

The theoretical framework developed by Rayleigh and Plateau in the 19th century has been remarkably accurate in describing macroscale experiments of liquid cylinder instability. Here we re-evaluate and revise the Rayleigh–Plateau instability for the nanoscale, where molecular dynamics experiments demonstrate its inadequacy. A new framework based on the stochastic lubrication equation is developed that captures nanoscale flow features and highlights the critical role of thermal fluctuations at small scales. Remarkably, the model indicates that classically stable (i.e. ‘fat’) liquid cylinders can be broken at the nanoscale, and this is confirmed by molecular dynamics.

Type
JFM Rapids
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© 2019 Cambridge University Press

1 Introduction

Understanding the surface-tension-driven break-up of liquid cylinders and their subsequent disintegration into droplets is crucial to a range of technologies such as ink-jet printing (Basaran, Gao & Bhat Reference Basaran, Gao and Bhat2013), fibre manufacture (Shabahang et al. Reference Shabahang, Kaufman, Deng and Abouraddy2011) and drug delivery (Mitragotri Reference Mitragotri2006). The classical theoretical results are due to Plateau (Reference Plateau1873), who showed that a cylinder of radius $r_{0}$ is unstable to sufficiently long wavelength disturbances $\unicode[STIX]{x1D706}>\unicode[STIX]{x1D706}_{crit}=2\unicode[STIX]{x03C0}r_{0}$ and Rayleigh (Reference Rayleigh1878), who conducted a linear stability analysis of the inviscid dynamics to find a fastest growing mode with wavelength $\unicode[STIX]{x1D706}_{max}=9.01r_{0}$ (or wavenumber $k_{max}=0.697/r_{0}$ ). The Rayleigh–Plateau (RP) theoretical framework has been subject to numerous generalisations (e.g. Rayleigh (Reference Rayleigh1892), Tomotika (Reference Tomotika1935)) and accurately describes macroscopic experiments (Eggers & Villermaux Reference Eggers and Villermaux2008).

With an increasing interest in micro and nano fluid dynamics, driven by emerging technologies in these domains (Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004), there has been cause to reassess the validity of classical analyses at smaller scales and question the absolute truth of long-held intuition. For liquid flows, molecular dynamics (MD) allows one to perform nanoscale ‘experiments’ and these have been shown to qualitatively reproduce RP-like instability (Koplik & Banavar Reference Koplik and Banavar1993). However, MD studies of long cylinders (Kawano Reference Kawano1998; Gopan & Sathian Reference Gopan and Sathian2014) have been unable to conclusively (i.e. quantitatively) validate classical theory.

It is well recognised that although $k_{max}$ can give a first approximation of the drop size formed, the actual break-up is a nonlinear asymmetric ‘universal’ solution (Eggers Reference Eggers1993) that resembles a drop connected to a thin cylinder. Interestingly, at the nanoscale, MD showed instead the emergence of a new double-cone break-up profile (Moseler & Landman Reference Moseler and Landman2000), generated by the presence of thermal fluctuations of molecules whose effects are negligible for standard liquids at the macroscale. This behaviour was described using a ‘stochastic lubrication equation’ (SLE) derived by applying the lubrication approximation to the equations of fluctuating hydrodynamics (Landau, Lifshits & Pitaevskii Reference Landau, Lifshits and Pitaevskii1980). The importance of the fluctuations has been confirmed by experiments (Hennequin et al. Reference Hennequin, Aarts, van der Wiel, Wegdam, Eggers, Lekkerkerker and Bonn2006; Petit et al. Reference Petit, Rivière, Kellay and Delville2012), analytical models (Eggers Reference Eggers2002) and subsequent MD (Choi, Kim & Kim Reference Choi, Kim and Kim2006; Kang & Landman Reference Kang and Landman2007).

Remarkably, despite failures to match classical RP theory to MD and a long recognition that ‘fluctuations substantially distort the shape of the cylinder’ (Koplik & Banavar Reference Koplik and Banavar1993), no one has considered how the RP instability mechanism is modified in the presence of thermal fluctuations. This is despite recognition of their importance to enhancing instabilities in thin film flows (Grün, Mecke & Rauscher Reference Grün, Mecke and Rauscher2006; Diez, González & Fernández Reference Diez, González and Fernández2016). Note, the focus of this work is on how thermal fluctuations affect stability, which is a related but distinct question to how thermal fluctuations affect break-up, as studied in Moseler & Landman (Reference Moseler and Landman2000).

In the next section, we develop an SLE-RP framework for incorporating fluctuations into the stability analysis of liquid cylinders at the nanoscale. The MD simulation method is introduced in § 3 to validate the SLE-RP. Results in § 4 show that the SLE-RP allows us to identify deviations from the classical picture, most notably a loss of the clear-cut Plateau stability criterion and very significant modifications of $k_{max}$ .

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.

5 Conclusions

The SLE has been shown to be a remarkably robust tool for extending classical RP arguments into the nanoscale, predicting the loss of a clear-cut Plateau stability boundary and the modification of $k_{max}$ due to thermal fluctuations, at a fraction of the computational cost of the benchmark MD. More generally, the underlying formulation based on fluctuating hydrodynamics has been confirmed as a useful framework for investigating nanoscale flows; whilst the deterministic nature of classical hydrodynamics is lost, nanoflows can apparently often be well described with (stochastic) partial differential equations without resorting to particle methods.

Promisingly, the effects of thermal fluctuations on capillary flows have previously been observed in experiments (Hennequin et al. Reference Hennequin, Aarts, van der Wiel, Wegdam, Eggers, Lekkerkerker and Bonn2006) using ultra-low surface-tension liquid–liquid systems that massively enhance $\ell _{fluc}$ to bring thermal fluctuations into the optically observable range. It is our hope that such systems can be used to experimentally verify our new predictions, most notably the violation of Plateau stability. Potential extensions of the framework are numerous. For example, to consider phase change phenomena, which are known to affect the break-up of nanocylinders (Kang & Landman Reference Kang and Landman2007), and to explore related flow configurations, such as liquid jets (Moseler & Landman Reference Moseler and Landman2000) where an intriguing open problem is to understand the interplay between fluctuations within the jet and perturbations imposed at the orifice in an attempt to control break-up.

Acknowledgements

Useful discussions with Dr L. Gibelli and Dr J. C. Padrino are gratefully acknowledged. This work was supported by the Leverhulme Trust (Research Project Grant) and the EPSRC (grants EP/N016602/1, EP/P020887/1 & EP/P031684/1).

References

Basaran, O. A., Gao, H. & Bhat, P. P. 2013 Nonstandard inkjets. Annu. Rev. Fluid Mech. 45, 85113.Google Scholar
Choi, Y. S., Kim, S. J. & Kim, M. U. 2006 Molecular dynamics of unstable motions and capillary instability in liquid nanojets. Phys. Rev. E 73, 016309.Google Scholar
Diez, J. A., González, A. G. & Fernández, R. 2016 Metallic-thin-film instability with spatially correlated thermal noise. Phys. Rev. E 93 (1), 013120.Google Scholar
Eggers, J. 1993 Universal pinching of 3D axisymmetric free-surface flow. Phys. Rev. Lett. 71 (21), 34583460.Google Scholar
Eggers, J. 2002 Dynamics of liquid nanojets. Phys. Rev. Lett. 89, 084502.Google Scholar
Eggers, J. & Dupont, T. F. 1994 Drop formation in a one-dimensional approximation of the Navier–Stokes equation. J. Fluid Mech. 262, 205221.Google Scholar
Eggers, J. & Villermaux, E. 2008 Physics of liquid jets. Rep. Prog. Phys. 71, 179.Google Scholar
Fernandez, G. A., Vrabec, J. & Hasse, H. 2004 A molecular simulation study of shear and bulk viscosity and thermal conductivity of simple real fluids. Fluid Phase Equilib. 221, 157163.Google Scholar
Gopan, N. & Sathian, S. P. 2014 Rayleigh instability at small length scales. Phys. Rev. E 90, 033001.Google Scholar
Green, M. S. 1954 Markoff random processes and the statistical mechanics of time-dependent phenomena. II. Irreversible processes in fluids. J. Chem. Phys. 22, 398413.Google Scholar
Grün, G., Mecke, K. & Rauscher, M. 2006 Thin-film flow influenced by thermal noise. J. Stat. Phys. 122 (6), 12611291.Google Scholar
Hennequin, Y., Aarts, D., van der Wiel, J. H., Wegdam, G., Eggers, J., Lekkerkerker, H. N. W. & Bonn, D. 2006 Drop formation by thermal fluctuations at an ultralow surface tension. Phys. Rev. Lett. 97, 244502.Google Scholar
Kang, W. & Landman, U. 2007 Universality crossover of the pinch-off shape profiles of collapsing liquid nanobridges in vacuum and gaseous environments. Phys. Rev. Lett. 98, 064504.Google Scholar
Kawano, S. 1998 Molecular dynamics of rupture phenomena in a liquid thread. Phys. Rev. E 58, 44684472.Google Scholar
Kirkwood, J. G. & Buff, F. P. 1949 The statistical mechanical theory of surface tension. J. Chem. Phys. 17, 338343.Google Scholar
Koplik, J. & Banavar, J. R. 1993 Molecular dynamics of interface rupture. Phys. Fluids A 5, 521536.Google Scholar
Kubo, R. 1957 Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems. J. Phys. Soc. Japan 12, 570586.Google Scholar
Landau, L. D., Lifshits, E. M. & Pitaevskii, L. P. 1980 Course of Theoretical Physics, 2nd edn. vol. 9, chap. 9, pp. 8691. Pergamon Press.Google Scholar
Lee, J. K., Barker, J. A. & Pound, G. M. 1974 Surface structure and surface tension: perturbation theory and Monte Carlo calculation. J. Chem. Phys. 60, 19761980.Google Scholar
Min, D. & Wong, H. 2006 Rayleigh’s instability of Lennard-Jones liquid nanothreads simulated by molecular dynamics. Phys. Fluids 18, 024103.Google Scholar
Mitragotri, S. 2006 Current status and future prospects of needle-free liquid jet injectors. Nat. Rev. Drug Discov. 5 (7), 543548.Google Scholar
Moseler, M. & Landman, U. 2000 Formation, stability, and breakup of nanojets. Science 289, 11651169.Google Scholar
Petit, J., Rivière, D., Kellay, H. & Delville, J.-P. 2012 Break-up dynamics of fluctuating liquid threads. Proc. Natl Acad. Sci. USA 109 (45), 1832718331.Google Scholar
Plateau, J. A. F. 1873 Statique Expérimentale et Théorique des Liquides Soumis Aux Seules Forces Moléculaires, vol. 2. Gauthier-Villars.Google Scholar
Plimpton, S. 1995 Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117, 119.Google Scholar
Rayleigh, Lord 1878 On the instability of jets. Proc. Lond. Math. Soc. 1, 413.Google Scholar
Rayleigh, Lord 1892 XVI. On the instability of a cylinder of viscous liquid under capillary force. Phil. Mag. J. Sci. 34, 145154.Google Scholar
Shabahang, S., Kaufman, J. J., Deng, D. S. & Abouraddy, A. F. 2011 Observation of the Plateau–Rayleigh capillary instability in multi-material optical fibers. Appl. Phys. Lett. 99 (16), 161909.Google Scholar
Stone, H. A., Stroock, A. D. & Ajdari, A. 2004 Engineering flows in small devices: microfluidics toward a lab-on-a-chip. Annu. Rev. Fluid Mech. 36, 381411.Google Scholar
Tomotika, S. 1935 On the instability of a cylindrical thread of a viscous liquid surrounded by another viscous fluid. Proc. R. Soc. Lond. A 150, 322337.Google Scholar
Trokhymchuk, A. & Alejandre, J. 1999 Computer simulations of liquid/vapor interface in Lennard-Jones fluids: some questions and answers. J. Chem. Phys. 111, 85108523.Google Scholar
Figure 0

Figure 1. Schematic of the cylinder instability.

Figure 1

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

Figure 2

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.

Figure 3

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.

Figure 4

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 5

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.