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. Reference Kaufman, Tao, Shabahang, Banaei, Deng, Liang, Johnson, Fink and Abouraddy2012), fabrication of structures in micro-/nano-electromechanical systems (Li et al. Reference Li, Zhao, Mao, Wu and Xu2015), and nano-printing (Zhang et al. Reference Zhang, He, Li, Xu and Li2016). 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 Reference Eggers and Villermaux2008): (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 (Reference Plateau1873) deduced the critical wavelength (
$\lambda _{crit}$) below which all interface disturbances decay, and Lord Rayleigh (Reference Rayleigh1878) 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 Reference Chen and Steen1997; Day, Hinch & Lister Reference Day, Hinch and Lister1998), the viscous regime (Papageorgiou Reference Papageorgiou1995), and the viscous-inertial regime (Eggers Reference Eggers1993). Experimental confirmations of these regimes (Castrejón-Pita et al. Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, Lister and Basaran2015; Lagarde, Josserand & Protière Reference Lagarde, Josserand and Protière2018) 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. Reference Yang, Qiao, Mu, Zhu, Xu and Si2019; Zhao et al. Reference Zhao, Wan, Chen, Chao and Xu2021b; Mu et al. Reference Mu, Qiao, Ding and Si2023).
However, at the nanoscale, classical theories prove inadequate due to the emergence of new physical mechanisms (Kavokine, Netz & Bocquet Reference Kavokine, Netz and Bocquet2021). 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 Reference Koplik and Banavar1993) and experimentally (Shi, Brenner & Nagel Reference Shi, Brenner and Nagel1994). Subsequently, Moseler & Landman (Reference Moseler and Landman2000) 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 Reference Landau, Lifshitz and Pitaevskij1987). 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. Reference Kadau, Rosenblatt, Barber, Germann, Huang, Carlès and Alder2007), droplet coalescence (Perumanath et al. Reference Perumanath, Borg, Chubynsky, Sprittles and Reese2019), bounded film flows (Zhang, Sprittles & Lockerby Reference Zhang, Sprittles and Lockerby2019; Zhao, Zhang & Si Reference Zhao, Zhang and Si2023), and moving contact lines (Liu et al. Reference Liu, Zhao, Lockerby and Sprittles2023).
For the linear instability of liquid nano-threads, it was initially demonstrated that the classical instability criterion remains applicable at the nanoscale (Min & Wong Reference Min and Wong2006; Tiwari et al. Reference Tiwari, Reddy, Mukhopadhyay and Abraham2008). Subsequent research by Gopan & Sathian (Reference Gopan and Sathian2014) indicated that thermal fluctuations affect the dynamics only during the last stage of breakup. However, Mo et al. (Reference Mo, Qin, Zhao and Yang2016) contested this, asserting that the growth rates of thermal capillary waves deviate considerably from classical theories. In a study by Zhao, Sprittles & Lockerby (Reference Zhao, Sprittles and Lockerby2019), 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 Reference Barker, Bell and Garcia2023).
Concerning the nonlinear dynamics, Eggers (Reference Eggers2002) 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 (Reference Moseler and Landman2000), 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. Reference Hennequin, Aarts, van der Wiel, Wegdam, Eggers, Lekkerkerker and Bonn2006; Petit et al. Reference Petit, Rivière, Kellay and Delville2012) and numerically (Arienti et al. Reference Arienti, Pan, Li and Karniadakis2011; Zhao et al. Reference Zhao, Zhou, Zhang, Chen, Liu and Wang2020b). However, a recent study by Zhao, Lockerby & Sprittles (Reference Zhao, Lockerby and Sprittles2020a) 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 Reference Gopan and Sathian2014; Xue et al. Reference Xue, Sbragaglia, Biferale and Toschi2018), presenting challenges for potential nanoscale applications. Despite extensive investigations into manipulating interfacial dynamics at the macroscale (Yang et al. Reference Yang, Qiao, Mu, Zhu, Xu and Si2019; Zhao et al. Reference Zhao, Wan, Chen, Chao and Xu2021b; Mu et al. Reference Mu, Qiao, Ding and Si2023), the exploration of such external perturbations at the nanoscale has been relatively limited. Fowlkes et al. (Reference Fowlkes, Horton, Fuentes-Cabrera and Rack2012) 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. (Reference Shah, van Steijn, Kleijn and Kreutzer2019) 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 (Reference Moseler and Landman2000) 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:
\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 (Reference Eggers and Dupont1994). Additionally, the full Laplace pressure in () is determined from the principal curvatures
\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
\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
\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
\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):
\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
\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 Reference MacCormack2003), 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
\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
\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. Reference Zhao, Lockerby and Sprittles2020a), 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. (Reference Zhao, Liu, Lockerby and Sprittles2022), 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 Reference Grün, Mecke and Rauscher2006; Diez, González & Fernández Reference Diez, González and Fernández2016), as follows:
\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}$:
\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. (Reference Zhao, Lockerby and Sprittles2020a). Therefore, the final discretised expression of the noise term becomes
\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. Reference Thompson2022). 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 Reference Molinero and Moore2009). 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 Reference Green1954; Kubo Reference Kubo1957) to calculate the dynamic viscosity
\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 Reference Kirkwood and Buff1949)
\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. (Reference Zhao, Zhao, Si and Chen2021a), 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 Reference Proakis and Manolakis1996), 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 Reference Plateau1873; Lord Rayleigh Reference Rayleigh1878), 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. (Reference Zhao, Sprittles and Lockerby2019). 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
\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
\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 Reference Rice2007), expressed as
\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
\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 Reference Papoulis and Pillai2002), 
$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
\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
\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 (Reference Moseler and Landman2000), 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 Reference Eggers and Villermaux2008), 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 Reference Eggers and Dupont1994):
\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 Reference Eggers2002), or the surface-tension-dominated one, 
$\alpha =1$ (Eggers Reference Eggers1993). 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 Reference Eggers1993, Reference Eggers2002). (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
\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 (Reference Eggers and Dupont1994), 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 Reference Moseler and Landman2000; Choi, Kim & Kim Reference Choi, Kim and Kim2006; Kang, Landman & Glezer Reference Kang, Landman and Glezer2008), the introduction of external perturbations, widely employed at the macroscale (Yang et al. Reference Yang, Qiao, Mu, Zhu, Xu and Si2019; Mu et al. Reference Mu, Qiao, Ding and Si2023), remains unexplored for actively controlling the breakup of nano-jets.
Acknowledgements
Useful discussions with Dr K. Mu and Dr J. Liu are gratefully acknowledged.
Funding
This work was supported by the National Natural Science Foundation of China (grant nos. 12202437, 12027801, 12388101), the Youth Innovation Promotion Association CAS (grant no. 2018491), the Fundamental Research Funds for the Central Universities, the Opening fund of State Key Laboratory of Nonlinear Mechanics, the China Postdoctoral Science Foundation (grant no. 2022M723044) and the International Postdoctoral Fellowship (grant no. YJ20210177).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of the spectrum function
Equation (2.6) is a linear equation with time-invariant coefficients, so it satisfies the superposition principle of solutions, which supports us to decompose the full solution into two components:
\begin{equation} H=H_{LE}+H_{fluc}. \end{equation}
Considering the initial interface shape, and assuming the initial velocity of the interface to be zero, the initial conditions are 
$H|_{t=0}=H_{i}$ and 
${\mathrm {d} H/\mathrm {d} t} |_{t=0}=0$, where 
$H_{i}$ represents the spectrum of the initial interface profile 
$h(z,0)$. The first term on the right-hand side can be obtained by solving the homogeneous form of (2.6) with these initial conditions:
\begin{equation} H_{LE}=H_{i}\exp\left( \frac{-a}{2}\,t \right) \left[ \cosh\left( \frac{b}{2}t \right) + \frac{a}{b}\sinh\left( \frac{b}{2}\,t \right) \right], \end{equation}where
\begin{equation} a=3\,\textit{Oh}\,k^2, \quad b=\sqrt{\left(9\,\textit{Oh}^2 -2\right)k^4+2k^2}. \end{equation}The second term on the right-hand side of (A1) could be calculated from the convolution of the excitation function (i.e. the inhomogeneous term) and the impulse response of (2.6):
\begin{equation} H_{fluc}=\sqrt{\frac{3\,\textit{Oh}}{2{\rm \pi}}}\,\textit{Th}\,k^2 \int_{0}^{t} \xi(k, t-\tau)\,G(k, \tau) \,\mathrm{d}\tau. \end{equation}
The impulse response 
$G(k, t)$ could be obtained by solving the equation
\begin{equation} \frac{\mathrm{d}^2 G}{\mathrm{d} t^2}+ 3\,\textit{Oh}\,k^2\,\frac{\mathrm{d} G}{\mathrm{d} t}+ \frac{k^4-k^2}{2}\,G=\delta(t), \end{equation}so we have
\begin{equation} G(k, t) = \frac{2}{b}\exp\left( \frac{-a}{2}\,t \right)\sinh\left( \frac{b}{2}\,t \right). \end{equation} As 
$H$ is a complex random variable with a zero mean, we should analyse it statistically, i.e. seek its root mean square from (A1):
\begin{equation} \left|H\right|_{rms} = \sqrt{\left\langle \left| H_{LE}+H_{fluc} \right|^2 \right\rangle} = \sqrt{\left\langle \left|H_{LE}\right|^2 \right\rangle + \left\langle \left|H_{fluc}\right|^2 \right\rangle}, \end{equation}
where the cross-term is erased since 
$H_{LE}$ and 
$H_{fluc}$ are orthogonal. Then applying the same operation to (A2) readily yields
\begin{equation} \left\langle \left|H_{LE}\right|^2 \right\rangle =\left|H_{i}\right|^2\mathrm{e}^{-at} \left[ \cosh\left( \frac{b}{2}\,t \right) + \frac{a}{b}\sinh\left( \frac{b}{2}\,t \right) \right]^2. \end{equation} Given that 
$\xi (k,t)$ is an uncorrelated Gaussian white noise, we derive 
$\langle | \xi (k,t) |^2 \rangle = L$, which leads to
\begin{align} \left\langle \left|H_{fluc}\right|^2 \right\rangle &= \frac{3\,\textit{Oh}}{2{\rm \pi}}\,\textit{Th}^2\,k^4 \left\langle \left|\int_{0}^{t} \xi(k, t-\tau)\,G(k, \tau) \,\mathrm{d}\tau\right|^2 \right\rangle\nonumber\\ &= \frac{3\,\textit{Oh}}{2{\rm \pi}}\,\textit{Th}^2\,k^4 \int_{0}^{t} \left\langle \left|\xi(k, t-\tau)\right|^2 \right\rangle G(k, \tau)^2 \,\mathrm{d}\tau\nonumber\\ &= \frac{3\,\textit{Oh}}{2{\rm \pi}}\,\textit{Th}^2\,k^4 L \int_{0}^{t} G(k, \tau)^2 \,\mathrm{d}\tau \nonumber\\ &= \frac{3L\,\textit{Oh}}{\rm \pi}\,\textit{Th}^2\,k^4\,\frac{a^2-b^2-a^2 \cosh (b t)-a c \sinh (b t)+b^2\ \mathrm{e}^{a t}}{a b^2\left(a^2-b^2\right) \mathrm{e}^{a t}}. \end{align}
Organising all the above results, we have the spectrum function 
$|H|_{rms}$ described in (2.8) and (2.9).
Appendix B. Spectral leakage
The Fourier transform over a finite range 
$[0, L]$ introduces spectral leakage, leading to the prediction of numerous irrelevant modes (sidelobes) besides 
$k_0$ as depicted in figure 12 (Proakis & Manolakis Reference Proakis and Manolakis1996). However, in the discrete Fourier transform (DFT), a finite signal is extended periodically, resulting in a discrete spectrum (Proakis & Manolakis Reference Proakis and Manolakis1996) where the sidelobes cannot be captured. To align with the outcomes obtained from the DFT, we eliminate these sidelobes and retain only the main lobe. The peak value of this main lobe is 
$| H_{i} |_{max} = | H_{i} |_{k=k_0} = A_0L/2$, and its bandwidth is 
$4{\rm \pi} /L$.
Figure 12. Spectra of a sinusoidal wave 
$A_0\sin ( k_0 z )$ with 
$z \in [0, L]$, where 
$A_0=0.06$, 
$k_0=3{\rm \pi} /20$ and 
$L=120$. The results are obtained from (2.10) (blue solid line) and a numerical discrete Fourier transform (DFT) (red dashed line with circles).
Appendix C. Details of the MacCormack method
Two differential operators, 
$\Delta _{f}$ and 
$\Delta _{b}$, are introduced to represent the forward and backward differences, respectively:
\begin{equation} \Delta_{f} f = \frac{f_{i+1} - f_i}{z_{i+1} - z_i}, \quad \Delta_{b} f = \frac{f_i - f_{i-1}}{z_i - z_{i-1}}. \end{equation}
Then 
$\boldsymbol {D}$ is discretised by the forward difference for the predictor step, written as
\begin{equation} \boldsymbol{D}\left(u^t_i,h^t_i\right) = \begin{pmatrix} D_1 \\ D_2 \end{pmatrix}, \end{equation}where
\begin{equation} \begin{cases} \displaystyle D_1 =-u^t_i\,\Delta_{f} u^t_i -\Delta_{f} p^t_i + \frac{3\,\textit{Oh}}{\left(h^t_i\right)^2}\, \frac{\left(h^t_{i+1}\right)^2 \Delta_{f} u_i^t - \left(h^t_{i}\right)^2 \Delta_{b} u_i^t}{z_{i+1}-z_i}\\[12pt] \quad\qquad{}+ \sqrt{\dfrac{6}{\rm \pi}}\,\dfrac{\textit{Th}\,\sqrt{\textit{Oh}}} {\left(h^t_i\right)^2}\,\Delta_{f} \left(h^t_i N^t_i\right), \\[15pt] D_2 =-u_i^t\,\Delta_{f} h^t_i - \dfrac{1}{2} h^t_i\,\Delta_{f} u^t_i. \end{cases} \end{equation}
The backward difference is applied for 
$\bar {\boldsymbol {D}}$:
\begin{equation} \bar{\boldsymbol{D}}\left(\bar{u}^{t+1}_i, \bar{h}^{t+1}_i\right)= \begin{pmatrix} \bar{D}_1 \\ \bar{D}_2 \end{pmatrix}, \end{equation}where
\begin{equation} \begin{cases} \bar{D}_1 =-\bar{u}^{t+1}_i\,\Delta_{b} \bar{u}^{t+1}_i -\Delta_{b} \bar{p}^{t+1}_i + \dfrac{3\,\textit{Oh}}{\left(\bar{h}^{t+1}_i\right)^2}\,\dfrac{\left(\bar{h}^{t+1}_{i}\right)^2 \Delta_{f} \bar{u}_i^{t+1} - \left(\bar{h}^{t+1}_{i-1}\right)^2 \Delta_{b} \bar{u}_i^{t+1}}{z_{i}-z_{i-1}} \\[12pt] \quad\qquad{}+\sqrt{\dfrac{6}{\rm \pi}}\,\dfrac{\textit{Th}\,\sqrt{ \textit{Oh} }}{\left(\bar{h}^{t+1}_i\right)^2}\, \Delta_{b}\left(\bar{h}^{t+1}_i N^t_i\right), \\[15pt] \bar{D}_2 =-\bar{u}_i^{t+1}\,\Delta_{b} \bar{h}^{t+1}_{i} - \dfrac{1}{2}\bar{h}^{t+1}_{i}\,\Delta_{b} \bar{u}^{t+1}_i. \end{cases} \end{equation}Appendix D. Spatially correlated noise model
In this appendix, we introduce the spatially correlated noise model, first proposed by Grün et al. (Reference Grün, Mecke and Rauscher2006) for nanoscale bounded films, where an exponential correlation function is employed:
\begin{equation} F_{c}(z,l_{c})= \begin{cases} \displaystyle \frac{1}{X} \exp\left(-\frac{1}{2} \left[ \frac{L}{l_{c}} \sin\left( \frac{{\rm \pi} z}{L}\right)\right]^2 \right) & \text{for } l_{c} > 0, \\[12pt] \delta(z) & \text{for } l_{c} = 0. \end{cases} \end{equation}
Here, 
$l_{c}$ is the spatial correlation length, 
$L$ is the domain length, and 
$X$ is such that 
$\int _0^L F_{c}(z,l_{c}) \,\mathrm {d} z = 1$. Diez et al. (Reference Diez, González and Fernández2016) calculated the integral and found that 
$\chi _q$ could be expressed by the Bessel function
\begin{equation} \chi_q = \left.\mathrm{I}_q(\alpha) \right/ \mathrm{I}_0(\alpha), \end{equation}where
\begin{equation} \alpha = \left( \frac{L}{2 l_{c}} \right)^2, \quad k = \frac{2{\rm \pi} q}{L}. \end{equation}
Figure 13(a) shows the eigenvalue spectra for several values of 
$l_{c}$. Note that for 
$l_{c} \rightarrow 0$ (i.e. 
$\alpha \rightarrow \infty$), we have 
$\chi _q \rightarrow 1$ for all 
$q$, leading to the limiting case of the white (uncorrelated) noise.
Figure 13. (a) Linear spectra of eigenvalues for several values of 
$l_{c}$ from (D2). Here, the wavenumber is 
$k=2 {\rm \pi}q/L$. (b) Spatially correlated noise with different 
$l_{c}$.
The term 
$g_q$ corresponds to the set of orthonormal eigenfunctions according to
\begin{equation} g_{q}(z)=\begin{cases} \displaystyle \sqrt{\frac{2}{L}} \cos\left( \frac{2{\rm \pi} qz}{L} \right) & \text{for } q > 0, \\[12pt] \displaystyle \sqrt{\frac{1}{L}} & \text{for } q = 0, \\[12pt] \displaystyle \sqrt{\frac{2}{L}} \sin\left( \frac{2{\rm \pi} qz}{L} \right) & \text{for } q < 0. \end{cases} \end{equation}Combining with the temporal correlated model proposed by Zhao et al. (Reference Zhao, Lockerby and Sprittles2020a), the final discretised expression of the noise term is
\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}
Samples of 
$N^t_i$ at one time instant are illustrated in figure 13(b) with different spatial correlation lengths. Note that a larger 
$l_{c}$ leads to smooth large-wavelength and small-amplitude noise.
Appendix E. Influence of the correlation lengths
In this appendix, we investigate the influence of the correlation length on the dynamics at both linear and nonlinear stages.
For the linear instability, we conduct the SLE simulations by using different correlation lengths for a long thread with parameters from case 1 (
$L=120$, 
$\textit{Oh}=0.65$, 
$\textit{Th}=0.07$ and 
$A_0=0$). Using an approach similar to that employed in § 4.1, 50 independent realisations are performed to gain statistics. A DFT of the interface position is then applied to get the ensemble-averaged spectra. Figure 14 illustrates the influence of both the spatial correlation length 
$l_c$ and the temporal one 
$t_c$. Comparing figures 14(a) and 14(b), the spatial correlation length is not found to have a significant impact on spectra when 
$l_c \leq 20 \Delta z$. However, when 
$l_c = 80 \Delta z$, there is a notable reduction in the spectrum at high wavenumbers compared to the theoretical results, suggesting that a larger correlation length suppresses capillary waves driven by thermal fluctuations. Additionally, figures 14(d)–14(f) indicate that for the time step (
$\Delta t = 5 \times 10^{-6}$) used in this paper, 
$t_c$ within the range of 
$1000 \Delta t$ have no significant impact on the instability results.
Figure 14. Comparisons between theoretical (solid lines) and numerical (dashed lines) spectra from the SLE at three instants: 
$t_1 = 6$ (black), 
$t_2 = 11$ (blue) and 
$t_3 = 18$ (red). (a–c) Influence of 
$l_c$ with 
$t_c=10 \Delta t$. (d–f) Influence of 
$t_c$ with 
$l_c=5 \Delta z$.
Furthermore, we examine the impacts of correlation lengths on the interface profiles, particularly at the nonlinear stage. Given that figure 14 demonstrates the minimal effects of 
$t_c$, only the influence of 
$l_c$ is explored here. To reduce computational costs, we consider the simulations of a short thread (
$L=12$, 
$\textit{Oh}=0.65$, 
$\textit{Th}=0.07$ and 
$A_0=0.4$) with various spatial correlation lengths. Multiple SLE simulations (100 for each case) are then performed to get ensemble-averaged profiles. The results in figure 15 indicate that when 
$l_c \leq 20 \Delta z$, the overall interface profiles are not significantly affected, aligning with the findings at the linear stage in figure 14. However, the local interface morphology near 
$h_{min}$ is found to be affected by the spatial correlation lengths, indicating that 
$l_c$ in the numerical model represents the smallest spatial scale of thermal fluctuations.
Figure 15. Ensemble-averaged interface profiles at two instants (
$t_1 = 1.5$ and 
$t_2 = 4.5$) influenced by 
$l_c$.
Therefore, we can conclude that variations in correlation length within a certain range do not significantly affect the computational results. They only influence the local behaviours of fluctuating hydrodynamics below the correlation length. In this study, we choose two relatively small correlation lengths, 
$l_{c} = 5 \Delta z=0.15$ and 
$t_{c} = 10 \Delta t = 5 \times 10^{-5}$, which approximately correspond to the molecular scale and a time scale of one femtosecond in MD simulations, respectively. These parameters essentially preserves the true physical characteristics of thermal fluctuations in physical space.
