Fluctuation-driven dynamics of liquid nano-threads with external hydrodynamic perturbations

Abstract

Instability and rupture dynamics of a liquid nano-thread, subjected to external hydrodynamic perturbations, are captured by a stochastic lubrication equation (SLE) incorporating thermal fluctuations via Gaussian white noise. Linear instability analysis of the SLE is conducted to derive the spectra and distribution functions of thermal capillary waves influenced by external perturbations and thermal fluctuations. The SLE is also solved numerically using a second-order finite difference method with a correlated noise model. Both theoretical and numerical solutions, validated through molecular dynamics, indicate that surface tension forces due to specific external perturbations overcome the random effects of thermal fluctuations, determining both the thermal capillary waves and the evolution of perturbation growth. The results also show two distinct regimes: (i) the hydrodynamic regime, where external perturbations dominate, leading to uniform ruptures, and (ii) the thermal-fluctuation regime, where external perturbations are surpassed by thermal fluctuations, resulting in non-uniform ruptures. The transition between these regimes, modelled by a criterion developed from linear instability theory, exhibits a strong dependence on the amplitudes and wavenumbers of the external perturbations.

1. Introduction

The interfacial dynamics of liquid nano-threads plays a crucial role in modern fluid-based techniques, including in-fibre particle production (Kaufman et al. 2012), fabrication of structures in micro-/nano-electromechanical systems (Li et al. 2015), and nano-printing (Zhang et al. 2016). Experimentally observing the dynamics at the nanoscale is often challenging, highlighting the significance of modelling and simulation in unravelling the underlying physics.

Classical models describing the macroscale dynamics of liquid threads typically consist of two stages (Eggers & Villermaux 2008): (i) the linear dynamics of instability, and (ii) the nonlinear dynamics leading to rupture. Theoretical foundations for linear instability were laid by two pioneers: Plateau (1873) deduced the critical wavelength ($\lambda _{crit}$) below which all interface disturbances decay, and Lord Rayleigh (1878) identified the fastest-growing mode ($\lambda _{max}$) by applying the normal mode expansion to the axisymmetric Navier–Stokes equations. Concerning the nonlinear dynamics, various scaling laws have been developed to describe the final pinch-off stage in three typical scenarios: the inertial regime (Chen & Steen 1997; Day, Hinch & Lister 1998), the viscous regime (Papageorgiou 1995), and the viscous-inertial regime (Eggers 1993). Experimental confirmations of these regimes (Castrejón-Pita et al. 2015; Lagarde, Josserand & Protière 2018) further indicate that transitions between them are notably intricate. Building upon these classical models for the Rayleigh–Plateau (RP) instability and rupture, recent studies have employed specific actuations to introduce external perturbations to manipulate the interfacial dynamics of liquid threads/jets at the macroscale, facilitating the generation of uniform droplets (Yang et al. 2019; Zhao et al. 2021b; Mu et al. 2023).

However, at the nanoscale, classical theories prove inadequate due to the emergence of new physical mechanisms (Kavokine, Netz & Bocquet 2021). One significant factor is thermal fluctuations caused by random molecular motions. Their influence on the interfacial dynamics of liquid threads was first emphasised numerically (Koplik & Banavar 1993) and experimentally (Shi, Brenner & Nagel 1994). Subsequently, Moseler & Landman (2000) conducted molecular dynamics (MD) simulations of nano-jets to reveal a distinctive ‘double-cone’ shape near the rupture point, contrasting the long neck observed at the macroscale. The discrepancy attributed to thermal fluctuations was elucidated using a stochastic lubrication equation (SLE) derived by applying the slender-body approximation to the governing equations of fluctuating hydrodynamics (Landau, Lifshitz & Pitaevskij 1987). This approach captures the influence of thermal fluctuations by incorporating stochastic flux terms. Furthermore, thermal fluctuations have proven to play a vital role in other nanoscale interfacial hydrodynamics, such as fluid mixing (Kadau et al. 2007), droplet coalescence (Perumanath et al. 2019), bounded film flows (Zhang, Sprittles & Lockerby 2019; Zhao, Zhang & Si 2023), and moving contact lines (Liu et al. 2023).

For the linear instability of liquid nano-threads, it was initially demonstrated that the classical instability criterion remains applicable at the nanoscale (Min & Wong 2006; Tiwari et al. 2008). Subsequent research by Gopan & Sathian (2014) indicated that thermal fluctuations affect the dynamics only during the last stage of breakup. However, Mo et al. (2016) contested this, asserting that the growth rates of thermal capillary waves deviate considerably from classical theories. In a study by Zhao, Sprittles & Lockerby (2019), a framework was developed for modelling the linear instability of interfaces in the presence of thermal fluctuations, named the SLE-RP. At the nanoscale, the SLE-RP shows that the criterion of the classical RP instability can be violated, and $\lambda _{max}$ predicted by classical theories is significantly modified (notably becoming time-dependent), recently supported by numerical simulations for the governing equations of fluctuating hydrodynamics (Barker, Bell & Garcia 2023).

Concerning the nonlinear dynamics, Eggers (2002) first derived a similarity solution from the SLE to describe the nonlinear dynamics of nano-thread rupture. This solution successfully replicated the double-cone profile observed by Moseler & Landman (2000), and presented a power law governing the progression of the minimum thread radius to rupture: $h_{min} \sim (t_r-t)^{0.418}$, where $t_r$ represents the rupture time. This power law, distinct from the macroscale counterparts, was subsequently verified experimentally (Hennequin et al. 2006; Petit et al. 2012) and numerically (Arienti et al. 2011; Zhao et al. 2020b). However, a recent study by Zhao, Lockerby & Sprittles (2020a) challenged this power law, demonstrating its validity only under the condition of ultra-low surface tension.

The influence of thermal fluctuations on liquid nano-threads leads to their breakup into droplets with sizes spanning a broad range (Gopan & Sathian 2014; Xue et al. 2018), presenting challenges for potential nanoscale applications. Despite extensive investigations into manipulating interfacial dynamics at the macroscale (Yang et al. 2019; Zhao et al. 2021b; Mu et al. 2023), the exploration of such external perturbations at the nanoscale has been relatively limited. Fowlkes et al. (2012) investigated dewetting of striped liquid films with prescribed perturbations using MD simulations. The external perturbations with wavelength $\lambda > \lambda _{crit}$ were found to determine the RP instability of the liquid films and lead to uniform droplets, while perturbations with wavelength $\lambda < \lambda _{crit}$ proved ineffective in controlling droplet sizes. But the effects of the thermal fluctuations were not quantitatively analysed in their work. Shah et al. (2019) explored the instability of ultra-thin films driven by both thermal fluctuations and drainage due to the curvature of the initial interface profile (similar to an external perturbation). The competition between these two effects yields two regimes of the instability: the dimple-dominated regime and the fluctuation-dominated regime. Notably, only the tangential curvature (for calculating surface tension) was considered in the planar liquid film, while both tangential and circumferential curvatures are crucial in the interfacial dynamics of liquid threads. Surface tension forces due to the latter serve as the dominant driving force for the instability, distinguishing it significantly from planar liquid films. Hence the physics governing the interaction between external perturbations and thermal fluctuations, and their collective impact on the interfacial dynamics of nano-threads, remains uncertain.

In this study, the SLE and MD are employed to investigate the effects of external perturbations on the interfacial dynamics of liquid nano-threads. The paper is organised as follows. In § 2, we introduce the models of fluctuating hydrodynamics of liquid nano-threads, where the SLE is first presented in § 2.1, followed by analytical solutions of the linearised SLE (§ 2.2) and schemes for the numerical solutions of the SLE (§ 2.3). Subsequently, details of MD simulations are introduced in § 3. Results in § 4 show the influence of the external perturbations and thermal fluctuations on (i) the thermal capillary wavelengths (§ 4.1), and (ii) the evolution of perturbation growth (§ 4.2). Two instability regimes are also defined with boundaries of regime conversions, presented in § 4.3.

2. Fluctuating hydrodynamics modelling

In this section, the SLE, as a simple mathematical model, is introduced (§ 2.1) to pursue theoretical solutions (§ 2.2) and numerical solutions (§ 2.3) of the dynamics of liquid nano-threads.

2.1. Stochastic lubrication equation

We consider an axisymmetric liquid nano-thread in a cylindrical coordinate system, with its axis aligned along the $z$ direction (figure 1), with $h_0$ the initial radius of the thread. The SLE was derived by Moseler & Landman (2000) via applying a lubrication approximation to the axisymmetric fluctuating hydrodynamic equations, allowing the dynamics of the interface to be described by the thread radius $h(z,t)$ and the axial velocity $u(z,t)$.

Figure 1. Schematic of an axisymmetric liquid nano-thread with capillary waves.

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:

(2.1a–e)\begin{equation} h=\frac{h^*}{h_0}, \quad t=\frac{t^*}{\sqrt{\left.\rho h_0^3 \right/ \gamma}}, \quad u=\frac{u^*}{\sqrt{\gamma / \rho h_0}}, \quad p=\frac{p^*}{\gamma / h_0}, \quad N=\frac{N^*}{\sqrt[4]{\left.\gamma \right/ \rho h_0^{5}} }, \end{equation}

where $h^*$, $t^*$, $u^*$ and $p^*$, respectively, indicate dimensional thread radius, time, velocity and pressure. The variables without asterisks represent corresponding dimensionless ones (note that the dimensional material parameters are not given asterisks). Here, $\rho$ is the density, and $\gamma$ is the surface tension. To model thermal fluctuations, a stochastic term $N(z,t)$ is introduced, standing for a Gaussian white noise that obeys the fluctuation–dissipation theorem. Its mean and autocovariance are respectively calculated as

$$\begin{align} \begin{cases} \left\langle N(z,t) \right\rangle & = 0, & (2.2a) \\[12pt] \left\langle N(z,t)\,N(\acute{z}, \acute{t}) \right\rangle &= \delta(z-\acute{z})\,\delta(t-\acute{t}), & (2.2b)\end{cases} \end{align}$$

where $\delta$ represents a unit impulse function, and $\langle {\,\cdot\,} \rangle$ denotes ensemble averages. Here, $\acute {t}$ and $\acute {z}$ could be infinitesimally close to the original values in time or space. The dimensionless SLE can be written as

$$\begin{align} \begin{cases}\frac{\partial u}{\partial t} & =-u u^{\prime}-p^{\prime}+3\,\textit{Oh}\,\frac{\left(h^{2} u^{\prime}\right)^{\prime}}{h^{2}}+\sqrt{\frac{6\,\textit{Oh}}{\rm \pi}}\,\textit{Th}\,\frac{(h N)^{\prime}}{h^{2}}, & (2.3a) \\[12pt] \frac{\partial h}{\partial t} & =-u h^{\prime}-\frac{u^{\prime} h}{2}, & (2.3b) \end{cases} \end{align}$$

where the primes denote partial derivatives with respect to $z$. In (), one dimensionless quantity is the Ohnesorge number $\textit{Oh}$, which relates the viscous forces to inertial and surface tension forces, i.e. $\textit{Oh} = \mu / \sqrt {\rho \gamma h_0}$. Here, $\mu$ is the dynamic viscosity. Another dimensionless quantity is the thermal-fluctuation number $\textit{Th} = l_{T} / h_0$, representing the relative intensity of interfacial fluctuations. Here, $l_{T}= \sqrt {k_{B} T / \gamma }$ is the characteristic thermal-fluctuation length, with $k_{B}$ being the Boltzmann constant, and $T$ the temperature. When $\textit{Th}=0$, the SLE reduces to the deterministic lubrication equation (LE) proposed by Eggers & Dupont (1994). Additionally, the full Laplace pressure in () is determined from the principal curvatures

(2.4)\begin{equation} p=\frac{1}{h \sqrt{1+\left(h^{\prime}\right)^{2}}}-\frac{h^{\prime \prime}}{\left[1+\left(h^{\prime}\right)^{2}\right]^{{3}/{2}}}. \end{equation}

2.2. Linear instability analysis

For the linear instability, we take $h(z,t) = 1 + \tilde {h}(z,t)$ and assume that $\tilde {h}(z,t) \ll 1$. Substituting this into () and ignoring higher-order terms gives the linearised SLE

(2.5)\begin{equation} \frac{\partial^{2} \tilde{h}}{\partial t^{2}}-3\,\textit{Oh}\,\frac{\partial\tilde{h}^{\prime\prime}}{\partial t}+\frac{1}{2}\left(\tilde{h}^{\prime\prime}+\tilde{h}^{\prime\prime\prime\prime}\right) =-\sqrt{\frac{3\,\textit{Oh}}{2{\rm \pi}}}\,\textit{Th}\,N^{\prime\prime}. \end{equation}

A Fourier transform within $[0, L]$ is then applied for (2.5) to give a second-order ordinary differential equation

(2.6)\begin{equation} \frac{\mathrm{d}^2 H}{\mathrm{d} t^2}+3\,\textit{Oh}\,k^2\,\frac{\mathrm{d} H}{\mathrm{d} t}+\frac{k^4-k^2}{2}\, H=\sqrt{\frac{3\,\textit{Oh}}{2{\rm \pi}}}\,\textit{Th}\, k^2 \xi, \end{equation}

where

(2.7a,b)\begin{equation} H(k,t) = \int_{0}^{L}\tilde{h}(z,t)\ \mathrm{e}^{-\mathrm{i} kz} \,\mathrm{d} z, \quad \xi(k,t) = \int_{0}^{L}N(z,t)\ \mathrm{e}^{-\mathrm{i} kz} \,\mathrm{d} z. \end{equation}

Here, $k$ is the wavenumber. The transformed variable $H$ represents the spectrum of the thermal capillary waves of the interfaces. Its final solution is expressed as follows (see Appendix A for derivation):

(2.8)\begin{equation} \left|H\right|_{rms} = \sqrt{\left\langle \left|H_{LE}\right|^2 \right\rangle + \left\langle \left|H_{fluc}\right|^2 \right\rangle}, \end{equation}

where

$$\begin{align} \begin{cases} \left\langle \left|H_{LE}\right|^2 \right\rangle = \left|H_{i}\right|^2 {\rm e}^{-at} \left[ \cosh\left( \frac{bt}{2} \right) + \frac{a}{b}\sinh\left( \frac{bt}{2} \right) \right]^2, & (2.9a) \\[12pt] \left\langle \left|H_{fluc}\right|^2 \right\rangle = \frac{3L\,\textit{Oh}}{\rm \pi}\,\textit{Th}^2\,k^4\,\frac{a^2-b^2-a^2 \cosh (b t)-a b \sinh (b t)+b^2\,{\rm e}^{a t}}{a b^2\left(a^2-b^2\right) {\rm e}^{a t}}. & (2.9b) \end{cases} \end{align}$$

Here, the subscript ‘rms’ represents root mean square, $a=3\,\textit{Oh}\,k^2$, and $b=$ $\sqrt {(9\,\textit{Oh}^2 -2)k^4+2k^2}$. The solution is linearly decomposed into two terms: the hydrodynamic component $H_{LE}$, and the thermal-fluctuation component $H_{fluc}$, where $H_{LE}$ is the solution of the homogeneous form of (2.6), representing the classical (deterministic) RP instability, and $H_{fluc}$ arises from solving the full form of (2.6) with zero initial disturbances, representing the fluctuation-drive instability. In (2.9a), $H_{i}$ is the initial spectrum, determined by the capillary waves at $t=0$. With initial perturbation waves, $\tilde {h}(z, 0)=A_0\sin (k_0 z)$, we have

(2.10)\begin{equation} \left| H_{i} \right| = \left|\int_{0}^{L}A_0\sin(k_0z)\ \mathrm{e}^{-\mathrm{i} kz}\,\mathrm{d} z\right| = 2A_0k_0\left| \frac{\sin(kL/2)}{k^2-k_0^2} \right|. \end{equation}

2.3. Numerical scheme

In this work, we solve the SLE numerically with periodic boundary conditions using the MacCormack method (MacCormack 2003), a simple second-order explicit finite difference scheme in both time and space. At each time level, the solution is represented by two arrays, $\{h_i\}^M_{i=1}$ and $\{u_i\}^M_{i=1}$, where $M$ denotes the number of mesh points. The time-derivative terms are approximated as $(h_i^{t+1}-h_i^t) / \Delta t$ and $(u_i^{t+1}-u_i^t) / \Delta t$. The numerical method follows a two-step process, beginning with a predictor step

(2.11)\begin{equation} \begin{pmatrix} \bar{u}^{t+1}_i \\[2pt] \bar{h}^{t+1}_i \end{pmatrix} = \begin{pmatrix} u^t_i \\[2pt] h^t_i \end{pmatrix} + \boldsymbol{D}\left(u^t_i, h^t_i\right) \Delta t, \end{equation}

and a corrector step

(2.12) \begin{equation} \begin{pmatrix} u^{t+1}_i \\[2pt] h^{t+1}_i \end{pmatrix} = \begin{pmatrix} u^t_i \\[2pt] h^t_i \end{pmatrix} + \frac{\Delta t}{2} \left[ \boldsymbol{D}\left(\vphantom{\bar{u}^{t+1}_i, \bar{h}^{t+1}_i}u^t_i, h^t_i\vphantom{\bar{u}^{t+1}_i, \bar{h}^{t+1}_i}\right)+\bar{\boldsymbol{D}}\left( \bar{u}^{t+1}_i, \bar{h}^{t+1}_i \right) \right], \end{equation}

where $\bar {u}_i^{i+1}$ and $\bar {h}_i^{t+1}$ denote the ‘provisional’ values at time level $t+1$, and $\boldsymbol {D}$ encompasses all the partial spatial derivative terms on the right-hand side (expressions for $\boldsymbol {D}$ are listed in Appendix C).

For the stochastic term $N(z,t)$, the autocovariance of uncorrelated fluctuations can be approximated numerically by a two-dimensional rectangular (boxcar) function, non-zero over a time step ($\Delta t$) and grid spacing ($\Delta z$), expressed as $N(z,t) \approx N_i^t = \mathcal {N} / \sqrt {\Delta t\,\Delta z}$. Here, $\mathcal {N}$ signifies computer-generated random numbers, following a normal distribution with zero mean and unit variance. However, this model has been shown to be prone to numerical instability for the SLE (Zhao et al. 2020a), problems that are exacerbated as $\Delta z$ and $\Delta t$ become smaller, and the amplitude of noise becomes larger.

To establish a robust numerical scheme, we adopt the methodology introduced by Zhao et al. (2022), combining a spatially and temporally correlated noise model. This integration enforces correlations in the noise beneath the spatial correlation length $l_{c}$ and the temporal correlation length $t_{c}$. The uncorrelated behaviour is then approximated by taking the limit of these lengths approaching zero, ensuring that their numerical resolution remains accurate throughout the limiting process. The stochastic term $N(z,t)$ is expanded using separation of variables in the $Q$-Wiener process $W(z,t)$ (Grün, Mecke & Rauscher 2006; Diez, González & Fernández 2016), as follows:

(2.13)\begin{equation} N(z,t) = \frac{\partial W(z,t)}{\partial t} = \sum^{\infty}_{q =-\infty} \chi_q\,\dot{c}_q(t)\, g_q(z) . \end{equation}

Here, $\chi _q$ represents the eigenvalues of the correlation function $F_{c}$:

(2.14)\begin{equation} \chi_q = \int_{0}^{L} F_{c}(z) \exp({-2{\rm \pi}\mathrm{i} q z/L})\,\mathrm{d} z, \end{equation}

where $q$ represents an integer sequence. The expressions for $F_{c}$ and $g_q$, and other details regarding this model, can be found in Appendix D. The coefficient $\dot {c}_q(t)$, representing a temporally correlated noise process, is modelled using a straightforward linear interpolation between uncorrelated random noise at the endpoints of the temporal correlation interval, as proposed by Zhao et al. (2020a). Therefore, the final discretised expression of the noise term becomes

(2.15)\begin{equation} N^t_i= \frac{1}{\sqrt{t_{c}}}\sum^{({M+1})/{2}}_{q=-({M+1})/{2}} \chi_q\, \mathcal{N}_q\, g_q(z) . \end{equation}

In this work, we set the dimensionless grid size $\Delta z= 0.03$ and time step $\Delta t = 5 \times 10^{-6}$ (the dimensional parameters corresponding to MD simulations are $\Delta z^* = 0.11 \ \mathrm {nm}$ and $\Delta t^* = 0.13 \ \mathrm {fs}$, respectively). The correlation lengths are $l_{c} = 5\,\Delta z=0.15$ and $t_{c} = 10\,\Delta t = 5 \times 10^{-5}$. A discussion of the influence of $l_{c}$ and $t_{c}$ is presented in Appendix E. The boundary conditions are set as periodic.

3. Molecular dynamics

The MD simulations of this work are performed using the open source package LAMMPS (Thompson et al. 2022). The simulation box ($40\ \mathrm {nm} \times 40\ \mathrm {nm} \times 432 \ \mathrm {nm}$ in the $x$, $y$ and $z$ directions, respectively) has periodic boundary conditions imposed in all directions. A nano-thread of water, initially with radius $h_0=3.6\ \mathrm {nm}$, is positioned at the centre of the simulation domain. The thread length is equal to the length of the simulation box in the $z$ direction, i.e. the dimensionless thread length $L=120$. External perturbations are introduced through a sinusoidal function $h(z,t) = h_0+A_0\sin ( k_0 z )$, where $A_0$ denotes the initial amplitude, and $k_0$ represents the (angular) wavenumber (figure 2). Initial configurations with various external perturbations are generated from equilibrium simulations of a liquid bulk in the canonical (NVT) ensemble, employing the Nosé–Hoover thermostat at $T = 300 \ \mathrm {K}$. The interactions between water molecules are described by a coarse-grained force field, the mW potential (Molinero & Moore 2009). The final density of the bulk at equilibrium is $997\ \mathrm {kg}\ \mathrm {m}^{-3}$. Considering the ultra-low density of vapours predicted by the mW model, the nano-threads are simulated in a vacuum environment with time step $2.5 \ \mathrm {fs}$. The canonical ensemble with the Nosé–Hoover thermostat is employed again to keep $T = 300 \ \mathrm {K}$.

Figure 2. Initial configuration for MD simulations: a liquid nano-thread with an external perturbation.

To determine $\textit{Oh}$ and $\textit{Th}$, viscosity and surface tension of the liquid are required to be extracted from MD. Here, we employ the Green–Kubo relation (Green 1954; Kubo 1957) to calculate the dynamic viscosity

(3.1)\begin{equation} \mu=\frac{V}{k_{B} T} \int_0^{\infty} \left\langle \,p_{ij}(t)\,p_{ij}(0) \right\rangle {\rm d}t,\quad i\neq j, \end{equation}

where $V$ is the volume of a liquid bulk, $p_{ij}$ represents the off-diagonal elements of the pressure tensor, and $\langle \,p_{ij}(t)\,p_{ij}(0) \rangle$ is the autocorrelation function of $p_{ij}$. In addition, a liquid layer lying in the $x$–$y$ plane is used to estimate the surface tension. Resorting to pressures on the two free surfaces, we have the expression (Kirkwood & Buff 1949)

(3.2)\begin{equation} \gamma=\int_{-\infty}^{\infty}\left[\,p_{n}(z)-p_{t}(z)\right]{\rm d}z, \end{equation}

where the normal and tangential pressure components are defined as $p_{n} = p_{zz}$ and $p_{t} = ( p_{xx}+p_{yy} ) / 2$, respectively. Finally, we have $\mu = 3.14\times 10^{-4} \ \mathrm {Pa}\ \mathrm {s}$ and $\gamma = 6.5\times 10^{-2} \ \mathrm {N}\ \mathrm {m}^{-1}$ at $T = 300\ \mathrm {K}$, leading to $\textit{Oh} = 0.65$ and $\textit{Th} = 0.07$.

4. Results and discussions

In this section, the analytical and numerical solutions of the SLE are validated by MD results, showing the influence of the thermal fluctuations and external perturbations on the thermal capillary waves (§ 4.1) and evolution of perturbation growth (§ 4.2). Two instability regimes are also defined with boundaries of regime conversions presented in § 4.3.

4.1. Thermal capillary waves

To examine the theoretical solution in § 2.2, we conduct MD simulations on long threads with various wavenumbers $k_0$ and amplitudes $A_0$: case 1 ($A_0=0$), case 2 ($A_0=0.1$, $k_0={\rm \pi} /6$) and case 3 ($A_0=0.2$, $k_0=2{\rm \pi} /5$). For each case, 30 independent MD simulations (realisations) are performed to gather statistics.

Figures 3(a,c,e) illustrate the MD snapshots of cases 1–3. Specifically, case 1 represents a situation with no external perturbations, while cases 2 and 3 involve different external perturbations. In case 1 (figure 3a), perturbations arising from thermal fluctuations grow over time, generating significant capillary waves, and eventually lead to the final rupture at $t_4$. Since this fluctuation-driven instability is naturally stochastic, the liquid threads break up into non-uniform droplets. In contrast, the external perturbation in case 2 grows, despite being disturbed by the thermal fluctuations, and ultimately leads to a uniform rupture similar to the macroscale cases actively controlled (figure 3c). The external perturbation with a larger wavenumber (compared to $k_0$ in case 2) decays rapidly and is then overwhelmed by fluctuation-driven perturbations (figure 3e), leading to an irregular breakup pattern similar to that in case 1.

Figure 3. (a,c,e) The MD snapshots, and (b,d,f) the comparisons between theoretical (solid lines) and numerical (dashed lines with circles) spectra in cases 1–3. (a,b) Results of case 1 at six instants: 0.00, 3.27, 10.76, 18.24, 21.98 and 29.93. (c,d) Results of case 2 at six instants: 0.00, 1.40, 5.14, 8.88, 15.90 and 22.45. (e,f) Results of case 3 at six instants: 0.00, 3.27, 10.76, 18.24, 23.85 and 32.73. The liquid threads illustrated here are all truncated at $z = L/2$.

The phenomena presented above can be further explained quantitatively by the spectra in figures 3(b,d,f). Following the approach used by Zhao et al. (2021a), the profile $h(z,t)$ in each MD realisation is extracted from axially distributed annular bins based on a threshold value of particle number density. A discrete Fourier transform is applied to $h(z,t)$ to get the spectra. We then ensemble average the spectra at each instant over the realisations, and take the square root to produce the numerical spectra of cases 1–3 in figure 3 (dashed lines with circles). Good agreement with the theoretical model for the spectra can be found for all the cases, confirming the validity of (2.8).

The spectra of case 1 (figure 3b) display a modal distribution with a certain bandwidth at each time instant, explaining why the liquid thread exhibits a non-uniform breakup. In cases 2 and 3, the external perturbations lead to initial spikes in the spectra, representing the initial conditions of the hydrodynamic component $H_{LE}$, modelled by $H_{i}$ of (2.10). Note that (2.10) can cause ‘spectral leakage’ (Proakis & Manolakis 1996), which leads to noise in the initial spectra. To avoid this problem and compare with the results from the discrete Fourier transform, further processing on (2.10) is required (see Appendix B for details). The spike in case 2 increases rapidly and indicates that the hydrodynamic component dominates the entire dynamics, resulting in the formation of uniform droplets after the rupture. However, the spike in case 3 decays drastically and is overwhelmed by the fluctuation modes, denoting that thermal fluctuations re-dominate the instability. The difference between cases 2 and 3 can be explained by the classical RP theory (Plateau 1873; Lord Rayleigh 1878), where perturbations of short wavelength $\lambda < \lambda _{crit}$ would dissipate.

The dominant wavenumbers $k_{max}$ of the instability are also extracted from the peaks of spectra, illustrated in figure 4(a). For case 1, $k_{max}$ decreases monotonically to a constant, which has been pointed out by Zhao et al. (2019). Since the external perturbation dominates the instability of case 2, its $k_{max}$ maintains an invariant value, i.e. $k_{max} = k_0={\rm \pi} /6$. In case 3, $k_{max}$ remains equal to $k_0=2{\rm \pi} /5$ at the early stage. When the peak of $k_0$ is surpassed by the instability modes due to the thermal fluctuations, $k_{max}$ return to the trajectory observed in case 1.

Figure 4. The MD (circles) and theoretical (solid lines) results of temporal evolutions of (a) dominant wavenumber $k_{max}$ and (b) surface roughness $\varTheta$, in cases 1 (black), 2 (blue) and 3 (red).

Moreover, we define the evolution of surface roughness $\varTheta (t)$ via integrating the square of $\tilde {h}(z,t)$ over the entire spatial domain to measure the development of the thermal capillary waves. According to Parseval's theory, $\varTheta (t)$ can be expressed as

(4.1)\begin{equation} \varTheta^2 = \frac{1}{L} \left\langle \int_{0}^{L}\tilde{h}^2 \,\mathrm{d} z \right\rangle = \frac{1}{{\rm \pi} L}\int_0^{\infty}\left|H\right|_{rms}^2 {\rm d}k . \end{equation}

Since $|H|_{rms}$ consists of two components, the surface roughness is also divided into

(4.2)\begin{equation} \varTheta^2 = \varTheta_{LE}^2 + \varTheta_{fluc}^2 = \frac{1}{{\rm \pi} L}\left( \int_0^{\infty}\left\langle \left|H_{LE}\right|^2 \right\rangle{\rm d}k +\int_0^{\infty}\left\langle \left|H_{fluc}\right|^2 \right\rangle{\rm d} k \right). \end{equation}

Figure 4(b) illustrates the evolution of the roughness versus time in cases 1–3. The roughness $\varTheta$ increases constantly with time in both cases 1 and 2. Case 2 exhibits a higher initial growth rate due to external perturbations. In case 3, the surface roughness initially decreases due to the dissipation of initial hydrodynamic perturbation, and then increases driven by thermal fluctuations.

The spectra above present the interfacial dynamics only in the frequency domain, i.e. $|H|_{rms}(k,t)$. To gain a better understanding of dynamics in the spatial domain, we propose a distribution function of the perturbation amplitudes, $P(\hat {h})$, where $\hat {h}$ is introduced to represent a possible value of the random perturbation amplitudes $\tilde {h}$. The perturbations at the linear stage can be divided into two independent components: $\tilde {h} = \tilde {h}_{LE} + \tilde {h}_{fluc}$, where $\tilde {h}_{LE}$ represents waves generated by the classical RP instability, and $\tilde {h}_{fluc}$ accounts for waves from thermal fluctuations. So $P(\hat {h})$ can be modelled by the convolution of the probability distributions of each component (Rice 2007), expressed as

(4.3)\begin{equation} P= P_{LE} \otimes P_{fluc} . \end{equation}

Here, ‘$\otimes$’ denotes convolution.

To get the expression for $P_{LE}$, we introduce the cumulative distribution function

(4.4)\begin{equation} F_{\tilde{h}} ( \hat{h} ) = \int_{-\infty}^{\hat{h}} P_{LE} \,\mathrm{d} \tilde{h}. \end{equation}

Based on the classical RP instability, $\hat {h} = A(t)\sin ( k_0 \hat {z} )$, where $A$ grows or decays exponentially from the initial value $A_0$. So $\hat {h}$ and $\hat {z}$ have a one-to-one functional relationship and are piecewise monotonic. It is easy to get an inverse function, $\hat {z} = \mathrm {arcsin}(\hat {h}/A)/k_0$. According to the approach of the distribution function transformation (Papoulis & Pillai 2002), $F_{\tilde {h}}(\hat {h})$ is equal to the cumulative distribution function of $\hat {z}$, i.e. $F_z(\hat {z})$. When $z$ follows a uniform distribution, we have

(4.5)\begin{equation} P_{LE}( \hat{h} ) = \frac{\mathrm{d} F_{z}\left( \hat{z} \right)}{\mathrm{d} \hat{h}} = \frac{\mathrm{d} F_z}{\mathrm{d} \hat{z}} \left| \frac{\mathrm{d} \hat{z}}{\mathrm{d} \hat{h}} \right| = \frac{1}{{\rm \pi}\sqrt{A^2 - \hat{h}^2}}, \end{equation}

where $\hat {h} \in (-A, A)$. Note that the final expression for $P_{LE}$ is independent of the wavenumbers ($k_0$) of the initial perturbations.

Additionally, it is challenging to pursue a theoretical model of $P_{fluc}$ mainly due to the complexity of (2.5). So we extract numerically predicted $P_{fluc}$ from the MD simulations of case 1, where $P_{LE}$ can be neglected. We collect all the values of $\tilde {h}(z)$ from 30 realisations, then plot the numerical distributions of $\tilde {h}$ by the histograms in figures 5(a)–5(c). Promisingly, $\tilde {h}_{fluc}$ is observed to follow a Gaussian distribution with mean zero. Moreover, the standard deviation of $\tilde {h}_{fluc}$ is equal to $\varTheta _{fluc}$ in (4.1). Therefore, we have a ‘semi-theoretical’ model

(4.6)\begin{equation} P_{fluc}(\hat{h}) = \frac{1}{\sqrt{2{\rm \pi}}\,\varTheta_{fluc}}\exp\left( -\frac{\hat{h}^2}{2\varTheta_{fluc}^2} \right). \end{equation}

From a theoretical perspective, $\hat {h} \in (-\infty, \infty )$ in (4.6). However, $\hat {h}$ typically falls within $[-3 \varTheta _{fluc}, 3 \varTheta _{fluc}]$, accounting for a $99.7\,\%$ confidence interval. Here, $| 3 \varTheta _{fluc} | < 1$, ensuring that the perturbation amplitude is always smaller than the thread radius. The model is validated by the good agreement between its predictions and MD results in figures 5(a)–5(c). A more rigorous derivation of $P_{fluc}$ involves pursuing the Fokker–Planck equation of (2.5), which will be a subject of our future research. Combining (4.3), (4.5) and (4.6) gives us the final theoretical expression for $P(\hat {h})$.

Figure 5. Temporal evolutions of $P$ from MD simulations (histograms) and the theory (red solid lines) in cases 1–3. (a–c) Results of case 1 at three instants: 3.27, 10.76 and 18.24. (d–f) Results of case 2 at three instants: 1.40, 5.14 and 8.88. (g–i) Results of case 3 at three instants: 3.27, 10.76 and 18.24.

Figures 5(d)–5(i) compare the numerical results from MD simulations with the theoretical distributions predicted by (4.3) for cases 2 and 3. In case 2, the distribution function maintains a bimodal curve, signifying that the interface can largely preserve the sinusoidal feature. Similar to the trend in case 1, the two spikes also propagate outwards as the thermal capillary waves develop. In case 3, the initial spikes dissipate and ultimately merge into a Gaussian curve, aligning with the observations in figure 4.

4.2. Evolution of perturbation growth

Besides the distribution of wavelengths (wavenumbers) investigated in § 4.1, the growth of the perturbations, particularly in the cases with uniform rupture, is explored in this subsection. MD simulations are performed on long threads with initial amplitude $A_0=0.2$ and various wavenumbers: $k_0={\rm \pi} /15$ for case 4, $k_0={\rm \pi} /12$ for case 5, and $k_0=2{\rm \pi} /15$ for case 6. Additionally, numerical simulations for the SLE are also conducted to compare with the MD results and provide deeper insights into the evolution of perturbation growth.

Figure 6 displays two selected realisations at three instants from both the MD and SLE simulations of case 6. Though the SLE predictions deviate slightly from the ‘double-cone’ profile documented by Moseler & Landman (2000), they agree well with the MD results qualitatively. To study the evolution of the perturbations, we focus on the temporal evolution of the minimum (over $z$) thread radius, i.e. $h_{min}(t)$. To get statistics, we conduct multiple independent realisations: 30 for MD, and 100 for the SLE.

Figure 6. Evolution of interface profiles extracted from the selected SLE (top) and MD (bottom) simulations. The minimum thread radius $h_{min}$ and the rupture point are marked in (b,c), respectively.

According to the classical theory (Eggers & Villermaux 2008), the growth of the perturbations can be divided into the linear and nonlinear stages. Figure 7(a) illustrates the ensemble-averaged perturbation growth ($1-h_{min}$) at the linear stage, extracted from both the numerical solutions of the SLE and the MD results. For cases 4–6, good agreement is observed at all instants for both mean values and standard deviations (from the thermal fluctuations), further validating the numerical solutions of the SLE. Interestingly, the perturbation is found to grow exponentially, approximately following the relation $1-h_{min} \sim \mathrm {e}^{\omega t}$. Despite the presence of the thermal fluctuations, the growth rate $\omega$ is close to the analytical results (dashed lines in figure 7a), predicted by the dispersion relation of the LE (Eggers & Dupont 1994):

(4.7)\begin{equation} \omega = \frac{k}{2}\sqrt{9\,\textit{Oh}^2\,k^2 + 2(1-k^2)} - \frac{3}{2}\,\textit{Oh}\,k^2. \end{equation}

These observations further explain the occurrence of uniform breakup. The surface tension forces, induced by the external perturbations with specific wavenumbers, overcome the random effects due to thermal fluctuations, determining the final form of the thermal capillary waves. Moreover, the influence of $\textit{Oh}$ and $\textit{Th}$ is investigated using the SLE solver with $k_0$ and $A_0$ from case 6. We set $\textit{Th}=0.07$ in figure 7(b), and $\textit{Oh}=0.65$ in figure 7(c). Figure 7(b) shows that the growth rates of the perturbations, which decline with increasing $\textit{Oh}$, agree well with the predictions of (4.7), further confirming the dominant roles of the surface tension forces induced by the external perturbations. However, when $\textit{Th}$ increases, the growth rate deviates from the predictions of (4.7), indicating that thermal fluctuations regain a significant role. Notably, each realisation evolves over different time periods, so the ensemble average can account only for the shortest time across all realisations. When thermal fluctuations become crucial, the variance in evolution time is larger (i.e. the minimum of rupture time is smaller), hence the trajectory for the case with $\textit{Th}=0.24$ is quite short.

Figure 7. Evolution of the perturbation growth at the linear stage. (a) The ensemble-averaged SLE predictions (solid lines) and MD results (circles) for cases 4 (black), 5 (red) and 6 (blue). Here, $\textit{Oh} = 0.65$ and $\textit{Th} = 0.07$. The dashed lines represent the growth rates predicted by the dispersion relation (4.7) for the three specific wavenumbers in cases 4–6: $\omega = 0.108$ for $k_0 = {\rm \pi}/15$ (black), $\omega = 0.124$ for $k_0 = {\rm \pi}/12$ (red), and $\omega = 0.148$ for $k_0 = 2{\rm \pi} /15$ (blue). The error bars and shadows in the inset represent one standard deviation (either side of the mean) for MD and the SLE, respectively. (b) The SLE predictions for three values of $\textit{Oh}$: $\omega = 0.244$ for $\textit{Oh}=0.1$ (black), $\omega = 0.168$ for $\textit{Oh}=0.5$ (red), and $\omega = 0.053$ for $\textit{Oh}=2.5$ (blue). Here, $\textit{Th}=0.07$, $k_0 = 2{\rm \pi} /15$ and $A_0 = 0.2$. (c) The SLE predictions for three values of $\textit{Th}$: 0.08 (black), 0.16 (red) and 0.24 (blue). We have $\omega = 0.148$ for $\textit{Oh}=0.165$ (black). Here, $k_0 = 2{\rm \pi} /15$ and $A_0 = 0.2$.

To investigate the nonlinear evolution near rupture, $h_{min}$ extracted from simulations is plotted against time to rupture, $t_r-t$, shown in figure 8(a). The nonlinear dynamics of cases 4–6 is found to be nearly identical as approaching the rupture point, indicating that external perturbations do not affect the rupture dynamics despite their significant impacts on the evolution at the linear stage. Additionally, the inset of figure 8(a) suggests that a power law might govern the progression of the minimum thread radius to rupture: $h_{min} \sim (t_r-t)^\alpha$. However, the power law does not satisfy either the thermal-fluctuation-dominated power law, $\alpha =0.418$ (Eggers 2002), or the surface-tension-dominated one, $\alpha =1$ (Eggers 1993). Instead, it lies between the two, indicating that both fluctuations and surface tension forces contribute to the dynamics during the rupture stage. Additionally, figure 8(b) shows the ensemble-averaged rupture profiles of cases 4–6. The overall interface shape varies due to the influence of external perturbations with different wavelengths, whereas profiles near the rupture overlap, further supporting the conclusion in figure 8(a) that external perturbations do not impact the rupture dynamics.

Figure 8. (a) Minimum thread radius ($h_{min}$) against time to rupture ($t_r - t$) for cases 4 (black), 5 (red) and 6 (blue): comparison between ensemble-averaged MD results (circles) and SLE calculations (solid lines). The inset illustrates $h_{min}(t_r-t)$ on a logarithmic scale. The two dashed lines represent two power laws of similarity solutions for the rupture dynamics (Eggers 1993, 2002). (b) Ensemble-averaged rupture profiles of cases 4 (black), 5 (red) and 6 (blue): comparison between ensemble-averaged MD results (circles) and SLE calculations (solid lines).

4.3. Regime transition

Based on the results and discussions in §§ 4.1 and 4.2, the final interface profiles of the liquid nano-threads are determined by both thermal fluctuations and external perturbations. Figure 9 illustrates the influence of $k_0$ and $A_0$ of the external perturbations on the interface profiles extracted from the MD simulations.

Figure 9. Interface profiles after rupture with various wavenumbers $k_0$ and amplitudes $A_0$ of external perturbations. (a–c) The initial amplitudes of the external perturbations are fixed, i.e. $A_0=0.15$. (d–f) The wavenumbers of external perturbations are fixed, i.e. $k_0={\rm \pi} /5$. (a) Results with $k_0={\rm \pi} /30$ at four instants: 0.09, 8.42, 16.83 and 28.52. (b) Results with $k_0={\rm \pi} /6$ at four instants: 0.09, 4.68, 8.88 and 17.77. (c) Results with $k_0=3{\rm \pi} /10$ at four instants: 0.09, 6.55, 13.09 and 22.91. (d) Results with $A_0= 0.025$ at four instants: 0.09, 8.88, 17.77 and 30.40. (e) Results with $A_0= 0.1$ at four instants: 0.09, 6.08, 12.16 and 20.11. (f) Results with $A_0= 0.175$ at four instants: 0.09, 4.21, 8.42 and 16.83.

In figures 9(a)–9(c), the initial amplitude of the external perturbations is fixed ($A_0=0.15$) with various wavenumbers. The uniform breakup appears only in the case with $k_0= {\rm \pi}/6$. According to (4.7), the dimensionless growth rate in case 5 is $0.146$, much larger than those with $k_0= {\rm \pi}/30$ ($\omega =0.064$) and $k_0= 3 {\rm \pi}/10$ ($\omega =0.028$). Starting from the same amplitude, the external perturbation with larger growth rate is better able to overwhelm the effects of the thermal fluctuations, leading to the results in figures 9(a)–9(c).

Figures 9(d)–9(f) show the impact of different initial amplitudes with the same wavenumber ($k_0={\rm \pi} /5$), where the external perturbations have the same growth rate. The maximum $A_0$ is found to enhance the hydrodynamic component of the instability, generating uniform droplets after the rupture, shown in figure 9(f), while the minimum $A_0$ is overwhelmed by the thermal fluctuations, leading to the non-uniform breakup in figure 9(d). Interestingly, the rupture in figure 9(e) is ‘quasi-uniform’ with only one droplet coalescence. Note that this result is extracted from one selected realisation. Uniform breakup can also be found in other realisations of the case with $A_0= 0.1$, indicating a transition regime from the non-uniform breakup to the uniform breakup.

According to the simulation results in the preceding sections, two principal instability regimes can be summarised, providing a framework to describe different breakup patterns: (i) the ‘hydrodynamic regime’, characterised by the generation of uniform droplets, and (ii) the ‘thermal-fluctuation regime’, associated with non-uniform breakup. To distinguish the regimes, a parameter $\phi$ is introduced to quantify the relative intensity of the hydrodynamic component due to the external perturbations and thermal-fluctuation component, written as

(4.8)\begin{equation} \phi(t) = \left.{\displaystyle \int_{0}^{\infty}\sqrt{\left\langle \left|H_{LE}\right|^2 \right\rangle}\ {\rm d} k} \right/ {\displaystyle \int_{0}^{\infty}\sqrt{\left\langle \left|H_{fluc}\right|^2 \right\rangle}\ {\rm d}k}. \end{equation}

Note that $\phi$ is time-dependent. We set $\phi (t_r)=1$ as the boundary separating the hydrodynamic and thermal-fluctuation regimes. When $\phi (t_r)>1$, the external perturbations dominate the instability, while the thermal fluctuations exert more significant influence when $\phi (t_r)<1$. For the fixed values of $\textit{Oh}$ and $\textit{Th}$, contours of $\phi$ are generated as a regime map based on $k_0$ and $A_0$, illustrated in figure 10. Here, the distribution of $t_r$ in the regime map is fitted using a third-order polynomial based on the numerical results from the SLE.

Figure 10. Regime maps at (a) $\textit{Oh}=0.65$ and $\textit{Th}=0.07$, (b) $\textit{Oh}=0.10$ and $\textit{Th}=0.07$, (c) $\textit{Oh}=0.65$ and $\textit{Th}=0.04$. The regime maps are depicted using contours of $\phi$ and symbols representing the numerical results obtained from (a) MD and (b,c) the SLE. Circles, triangles and crosses denote the hydrodynamic, thermal-fluctuation and transition regimes, respectively.

Figure 10(a) presents the regime map for the MD results ($\textit{Oh} = 0.65$ and $\textit{Th} = 0.07$). Besides the cases presented in figure 9, more MD simulations with different values of $k_0$ and $A_0$ are performed to support the criterion of the regime map. Promisingly, the regime boundary (black solid line) from (4.8) generally matches the MD results represented by symbols (circles for the hydrodynamic regime, and triangles for the thermal-fluctuation regime), except for the four circles at the bottom. The crosses suggest the transition regime, which emerges near the boundary, i.e. $\phi (t_r) \approx 1$. This is consistent with the results of the case in figure 9(e). Moreover, the bottom of the boundary indicates the optimum wavenumber ($k_0= 0.49$) for the hydrodynamic regime, closely matching the dominant mode predicted by the classical RP theory of (4.7).

Figures 10(b) and 10(c) depict the influence of $\textit{Oh}$ and $\textit{Th}$ on the boundaries in the regime maps. Since MD is not available for arbitrary values of $\textit{Oh}$ and $\textit{Th}$, numerical solutions of the SLE are employed to confirm the regime map across a broad range of $\textit{Th}$ and $\textit{Oh}$, specifically $\textit{Oh} = 0.10$ and $\textit{Th} = 0.07$ in figure 10(b), and $\textit{Oh} = 0.65$ and $\textit{Th} = 0.04$ in figure 10(c). Comparison between figures 10(a) and 10(b) reveals that reducing $\textit{Oh}$ results in a rightward shift of the regime boundary. This trend is further presented in figure 11(a) and can be explained by the classical RP theory, where the dominant wavenumber of the instability increases as $\textit{Oh}$ decreases. Notably, the bottom point of the boundary in figure 11(a) also exhibits a slight upward movement. The main reason is that $\textit{Oh}$ not only affects the hydrodynamic component but also modifies the intensity of thermal fluctuations, as shown in (2.9b). Examining figures 10(a,c) and 11(b), the regime boundary is observed to move downwards as $\textit{Th}$ decreases, indicating that it becomes easier to enter the hydrodynamic regime with weaker thermal fluctuations.

Figure 11. Regime boundaries at (a) different $\textit{Oh}$ values ($\textit{Oh}_1=0.10$, $\textit{Oh}_2=0.50$, $\textit{Oh}_3=2.50$ and $\textit{Oh}_4=12.50$) with $\textit{Th} = 0.07$, and (b) different $\textit{Th}$ values ($\textit{Th}_1=0.03$, $\textit{Th}_2=0.05$, $\textit{Th}_3=0.07$ and $\textit{Th}_4=0.09$) with $\textit{Oh} = 0.65$.

5. Conclusions

In this paper, the SLE and MD are utilised to explore the influence of external perturbations and thermal fluctuations on the dynamics of liquid nano-threads.

Linear instability analysis is performed to derive a theoretical model for the spectra of thermal capillary waves, influenced by both thermal fluctuations and external perturbations. This model, validated by MD simulations, reveals the instability mode of a spike from a specific external perturbation and a continuous curve due to thermal fluctuations, corresponding to the uniform and non-uniform ruptures, respectively. An analytical model is then established for the two typical distributions of thermal capillary waves: bimodal distribution for uniform waves, and Gaussian distribution for stochastic ones. Besides the formulation of thermal capillary waves, the evolution of perturbation growth, particularly in cases with uniform rupture, is also investigated. The results of uniform rupture show that the perturbation grows exponentially at the linear stage, approximately following the classical linear theory proposed by Eggers & Dupont (1994), indicating the dominant roles of surface tension forces arising from the external perturbation with specific wavenumbers. However, the nonlinear evolution near rupture, determined jointly by surface tension forces and thermal fluctuations, is observed not to be affected by the external perturbations. Finally, two distinct regimes are defined to characterise the instability: (i) the hydrodynamic regime, marked by uniform droplets controlled by external perturbations, and (ii) the thermal-fluctuation regime, exhibiting a stochastic breakup pattern. A criterion is proposed to draw a regime map based on the perturbation amplitude ($A_0$) and wavenumber ($k_0$). The boundaries of these regimes, validated by MD and SLE simulations, are obtained, including a transition area observed.

While this paper provides new understanding of interfacial dynamics, it opens up several new avenues of enquiry. One avenue involves deriving the Fokker–Planck equation of the SLE, a deterministic equation describing the probability density function of $\tilde {h}$. The utilisation of the Fokker–Planck equation holds the promise not only to fortify the mathematical underpinnings of the distribution function in § 4.1 but also to provide additional theoretical insights into the nonlinear dynamics of liquid nano-threads. Another avenue is extending this study to a more practical fluid configuration, i.e. a liquid nano-jet. Despite the performed MD simulations for nano-jets (Moseler & Landman 2000; Choi, Kim & Kim 2006; Kang, Landman & Glezer 2008), the introduction of external perturbations, widely employed at the macroscale (Yang et al. 2019; Mu et al. 2023), remains unexplored for actively controlling the breakup of nano-jets.