2. Modified stochastic lubrication equation

Figure 1.

Schematic of the liquid thread instability. Here, $z$ is the axial coordinate, $r$ is the radial coordinate. The thread is infinitely long but only a portion with length $L$ is taken for study. The blue solid line depicts the profile of the liquid thread. The parameter $h$ is the radius of the thread as a function of $z$ and $t$ and $k$ is the wavenumber of the perturbation wave. The black dashed line depicts the profile of the cylindrical liquid thread at $t=0$ . The value of $h_0$ is the radius of this initial cylindrical thread and $C$ is a condensation parameter which will be defined below. The red dashed line depicts the time-dependent basic state that we will explain and use in § 4.

We consider an axisymmetric liquid thread in a cylindrical coordinate system, as depicted in figure 1. The initially cylindrical liquid thread has an initial radius $h_0$ . The thread is infinitely long but we take only a finite length $L$ for study. The SLE (Moseler & Landman Reference Moseler and Landman2000) that describes the evolution of thread profile and axial flow velocity reads

(2.1) \begin{equation} \frac {\partial h}{\partial t} = - uh^{\prime} - \frac {u^{\prime}h}{2}, \end{equation}

(2.2) \begin{equation} \frac {\partial u}{\partial t} = -u u^{\prime} - \frac {p^{\prime}}{\rho } + 3\nu \frac {(h^2u^{\prime})^{\prime}}{h^2} + \sqrt {\frac {6k_BT\nu }{\rho \pi }}\frac {(h\mathcal{N})^{\prime}}{h^2}, \end{equation}

where $h(z,t)$ is the thread radius, $u(z,t)$ is the axial flow velocity, both being functions of the axial coordinate $z$ and time $t$ , $p=\gamma \kappa$ is the pressure, $\kappa =h^{-1}(1+(h^{\prime})^2)^{-1/2}-h^{\prime\prime}(1+(h^{\prime})^2)^{-3/2}$ is the sum of the principal curvatures, $\rho$ is the fluid density, $\nu$ is the kinematic viscosity, $\gamma$ is the surface tension, $T$ is the temperature, $k_B$ is the Boltzmann constant, $\mathcal{N}$ is white noise with correlation function $\langle \mathcal{N}(z,t) \mathcal{N}(z^{\prime},t^{\prime})\rangle = \delta (z-z^{\prime})\delta (t-t^{\prime})$ and $\langle \boldsymbol{\cdot }\rangle$ denotes the ensemble average. The prime notation indicates derivative against $z$ , except when the prime is on $t$ or $z$ where it simply means another time or location variable.

Then we consider the evaporation–condensation process. There are two predominant evaporation–condensation models: the diffusion-limited model (Wilson & D’Ambrosio Reference Wilson and D’Ambrosio2023) and the one-sided model (Burelbach, Bankoff & Davis Reference Burelbach, Bankoff and Davis1988; Ajaev Reference Ajaev2005). The diffusion-limited model assumes that the evaporation process is rate limited by the diffusive transport of vapour. The one-sided model assumes that the evaporation process is rate limited by the phase change. Both are widely used and they can produce greatly different predictions (Larsson & Kumar Reference Larsson and Kumar2023). Since we are considering a nano liquid thread in its own vapour, we believe it is more appropriate to use the one-sided model (Cazabat & Guéna Reference Cazabat and Guéna2010; Larson Reference Larson2014). The one-sided model provides a constitutive equation that relates the mass flux at the interface to the local surface temperature (Burelbach et al. Reference Burelbach, Bankoff and Davis1988). It is derived by linearising the Hertz–Knudsen relation (Schrage Reference Schrage1953) using the Clausius–Clapeyron law. But since the system we study here could have a high supersaturation ratio, the linearisation becomes inaccurate. It is necessary to directly adopt the Hertz–Knudsen relation to predict the evaporation–condensation flux.

The vapour pressure at the surface of the liquid thread is given by the Kelvin or Thomson–Gibbs equation (Hinds Reference Hinds1999)

(2.3) \begin{equation} P_s = P_0 \exp {\left (\frac {\kappa \gamma V_0}{k_BT}\right)}\approx P_0\left (1+\frac {\kappa \gamma V_0}{k_B T}\right)\!, \end{equation}

where $P_0$ is the saturation vapour pressure at a planar surface and $V_0$ is the molecular volume in liquid phase. Using the Hertz–Knudsen relation, the local volume flux $J_v$ due to evaporation–condensation can be determined

(2.4) \begin{equation} J_v = \frac {\alpha _{c} P_0V_0 \left (1+\dfrac {\kappa \gamma V_0}{k_B T}-S\right)}{\sqrt {2\pi mk_B T}}, \end{equation}

where $m$ is the molecular weight, $S=P_v/P_0$ is the supersaturation ratio of the surrounding vapour, $P_v$ is the realistic surrounding vapour pressure and $\alpha _c$ is the anomalous evaporation/condensation coefficient; $0\lt \alpha _c\lt 1$ , describing the probability of phase change. The exact value of $\alpha _c$ is a matter of debate and is usually used to match experimental results to theoretical results. As the evaporation and condensation processes are separately Poisson processes, $J_v$ should also be a fluctuating quantity. However, in the following analysis we will neglect the fluctuations in $J_v$ . In Appendix A we provide a formal analysis of the influence of these flux fluctuations and show that they are indeed negligible.

To derive the mass conservation equation that accounts for the evaporation–condensation flux we start from the generalised kinematic boundary condition

(2.5) \begin{equation} \boldsymbol{n}\boldsymbol{\cdot }(\boldsymbol{v}|_{r=h}-\dot {\boldsymbol{x}})=J_v ,\end{equation}

where $\boldsymbol{n}$ is the unit normal vector to the interface, $\boldsymbol{v}$ is the fluid velocity and $\dot {\boldsymbol{x}}$ is the interface velocity. The unit normal vector is given by $\boldsymbol{n}=[1,-h^{\prime}]/\sqrt {1+(h^{\prime})^2}$ . The fluid velocity has two components: a radial component $v_r$ and an axial component $v_z$ , thus $\boldsymbol{v}=[v_r,v_z]$ . The interface velocity is given by $\dot {\boldsymbol{x}}=[\partial h/ \partial t,0]$ . After substitution of these parameters into (2.5) we have

(2.6) \begin{equation} -\frac {\partial h}{\partial t} + v_r|_{r=h} - v_z|_{r=h} h^{\prime} = J_v\sqrt {1+\left (h^{\prime}\right)^2}. \end{equation}

Using the radial expansion method (Eggers & Dupont Reference Eggers and Dupont1994), the velocity field can be written as a power series

(2.7) \begin{equation} v_z(r,z,t) = v_0(z,t) + v_2(z,t)r^2 + {\cdots} , \end{equation}

(2.8) \begin{equation} v_r(r,z,t) = -v_0^{\prime}(z,t)\frac {r}{2} - v_2^{\prime}(z,t)\frac {r^3}{4} + {\cdots} . \end{equation}

Inserting (2.8), (2.7) into (2.6), to the lowest order in $r$ we have

(2.9) \begin{equation} -\frac {\partial h}{\partial t}-v_0^{\prime}\frac {h}{2}-v_0h^{\prime} = J_v\sqrt {1+(h^{\prime})^2}. \end{equation}

Since the typical radial scale $l_r$ is much smaller than the axial scale $l_z$ , to the leading order $\sqrt {1+(h^{\prime})^2}=1+O((l_r/l_z)^2)\approx 1$ . Using this approximation, and replacing $v_0$ with $u$ , (2.9) can be transformed into

(2.10) \begin{equation} \frac {\partial h}{\partial t} = -uh^{\prime}-\frac {u^{\prime}h}{2}-J_v .\end{equation}

Substituting (2.4) into (2.10) we obtain the mass conservation equation with evaporation–condensation effects taken into account

(2.11) \begin{equation} \frac {\partial h}{\partial t} = - uh^{\prime} - \frac {u^{\prime}h}{2}-\frac {\alpha _c P_0V_0 \left (1+\dfrac {\kappa \gamma V_0}{k_B T}-S\right)}{\sqrt {2\pi mk_B T}}. \end{equation}

Let $B:={\alpha _c P_0V_0^2}/{\sqrt {{2\pi m (k_BT)}^{3}}}$ , $C:= {\alpha _c P_0V_0(S-1)}/{\sqrt {2\pi mk_B T}}$ , then (2.11) becomes

(2.12) \begin{equation} \frac {\partial h}{\partial t} = - uh^{\prime} - \frac {u^{\prime}h}{2}-B\gamma \kappa + C. \end{equation}

Equations (2.12), (2.2) constitute a new MSLE. We note that our MSLE is different from that derived by Kang & Landman (Reference Kang and Landman2007). In the work of Kang & Landman (Reference Kang and Landman2007) the saturation vapour pressure at a planar surface $P_0$ was treated as the same as the surrounding gas pressure. Kang & Landman (Reference Kang and Landman2007) also assumed that the pressure $P_0$ includes the contribution of inert gas. However, we know that in reality the saturation vapour pressure at a planar surface $P_0$ is an intrinsic property parameter dependent only on temperature.

Using the initial uniform radius $h_0$ , liquid density $\rho$ and surface tension $\gamma$ as basic units, the variables are non-dimensionalised as follows:

(2.13) \begin{equation} \begin{gathered} h=\frac {h^*}{h_0},\quad z=\frac {z^*}{h_0},\quad t = \frac {t^*}{\sqrt {\rho h_0^3/\gamma }},\quad u = \frac {u^*}{\sqrt {\gamma /\rho h_0}}, \quad p=\frac {p^*}{\gamma /h_0}, \quad \mathcal{N}=\frac {\mathcal{N}^*}{\sqrt [4]{\gamma /\rho h_0^5}},\quad \\ B = \frac {B^*}{\sqrt {h_0/\rho \gamma }},\quad C=\frac {C^*}{\sqrt {\gamma /\rho h_0}}, \end{gathered} \end{equation}

where the asterisk denotes dimensional quantities. In the following analysis, we will use dimensionless variables. After the non-dimensionalisation, the MSLE becomes

(2.14) \begin{equation} \frac {\partial h}{\partial t} = - uh^{\prime} - \frac {u^{\prime}h}{2}-B \kappa + C, \end{equation}

(2.15) \begin{equation} \frac {\partial u}{\partial t} = -u u^{\prime} - p^{\prime} + 3 {\textit{Oh}} \frac {\left (h^2u^{\prime}\right)^{\prime}}{h^2} + \sqrt {\frac {6 {\textit{Oh}}}{\pi }}{\textit{Th}}\frac {\left (h\mathcal{N}\right)^{\prime}}{h^2}, \end{equation}

where $\textit{Oh}=\nu \sqrt {\rho /\gamma h_0}$ is the Ohnesorg number and $\textit{Th}=\sqrt {k_BT/\gamma }/h_0$ is a dimensionless number characterising the strength of the thermal fluctuations.

3. Instability in the saturated vapour environment

Let us first consider the case that the two dimensionless parameters $B$ and $C$ are equal ( $B=C$ ), which implies that the surrounding vapour pressure is the same as the saturation vapour pressure at a surface with curvature $\kappa =1/h_0$ , so that $P_v = P_0(1+ {V_0\gamma }/{k_BTh_0})$ . When the liquid thread is placed in a small closed space, the surrounding vapour will quickly reach equilibrium with the liquid phase, leading to $B=C$ . In that case, there is a mass redistribution on the thread as a result of the curvature-dependent evaporation rate, but the total liquid mass is not considerably changed by the evaporation–condensation. The basic state of the system is $\bar {h}=1,\, \bar {u}=0$ so we can assume $h=1+\tilde {h}$ and $u=\tilde {u}$ , with $\tilde {h}$ and $\tilde {u}$ both small perturbative quantities $\ll 1$ . Thus, (2.14), (2.15) are linearised to

(3.1) \begin{equation} \frac {\partial \tilde {h}}{\partial t} = -\frac {\tilde {u}^{\prime}}{2} + B(\tilde {h}+\tilde {h}^{\prime\prime}), \end{equation}

(3.2) \begin{equation} \frac {\partial \tilde {u}}{\partial t} = \tilde {h}^{\prime} + \tilde {h}^{\prime\prime\prime} + 3{\textit{Oh}}\tilde {u}^{\prime\prime} + \sqrt {\frac {6{\textit{Oh}}}{\pi }} {\textit{Th}} \mathcal{N}^{\prime}. \end{equation}

Differentiating (3.1) with respect to $t$ , differentiating (3.2) with respect to $z$ and taking the second derivative of (3.1) with respect to $z$ , yields three equations. Combining these three equations, we can eliminate $\tilde {u}$ and obtain an equation with respect to $\tilde {h}$

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

Following Zhao et al. (Reference Zhao, Sprittles and Lockerby2019), this equation can be solved with finite Fourier transform

(3.4) \begin{equation} H(k,t)=\int _0^L\tilde {h}(z,t)\exp {(-ikz)}{\rm d}z, \end{equation}

(3.5) \begin{equation} N(k,t)=\int _0^L\mathcal{N}(z,t)\exp {(-ikz)}{\rm d}z, \end{equation}

where $k$ is the wavenumber. Finite Fourier transform implicitly forces the radius outside the domain $[0,L]$ to be zero, thus causes a spectral leakage. But the error due to the spectral leakage is expected to be very small when we use a very large $L$ . Applying the finite Fourier transform to (3.3) yields a second-order ordinary differential equation (ODE):

(3.6) \begin{equation} \frac {d^2 H}{{\rm d} t^2} + \left (k^2(B+3{\textit{Oh}})-B\right)\frac {{\rm d}H}{{\rm d}t} + \left (\frac {1}{2}+3{\textit{Oh}}B\right)(k^4-k^2)H=\sqrt {\frac {3{\textit{Oh}}}{2\pi }}{\textit{Th}} k^2 N. \end{equation}

The solution of this equation is the sum of two parts: $H = H_{LE}+H_{fluc}$ , where $H_{LE}$ is the solution of the corresponding homogeneous equation with initial conditions $H(t=0)=H_i$ and $\partial H/\partial t (t=0)= 0$ , and $H_{fluc}$ is the solution of the inhomogeneous equation with zero initial conditions $H(t=0)=0$ and $\partial H/\partial t (t=0)= 0$ . They are given by

(3.7) \begin{equation} H_{LE} = H_i\exp \left (-\frac {\alpha t}{2}\right)\left (\cosh \frac {\beta t}{2} + \frac {\alpha }{\beta }\sinh \frac {\beta t}{2}\right), \end{equation}

(3.8) \begin{equation} H_{fluc} = \sqrt {\frac {3{\textit{Oh}}}{2\pi }}{\textit{Th}} k^2 \int _0^t N(k,t-t^{\prime})G(k,t^{\prime}) {\rm d}t^{\prime}, \end{equation}

where $\alpha = k^2(B+3{\textit{Oh}})-B$ , $\beta = \sqrt {(k^2(B+3{\textit{Oh}})-B)^2 - (2+12{\textit{Oh}}B)(k^4-k^2)}$ . Here, $G(k,t^{\prime})$ is the impulse response (Green’s function) of the linear operator

(3.9) \begin{equation} \mathcal{L}:=\frac {d^2}{{\rm d}t^2} + \alpha \frac {d}{{\rm d}t} + \left (\frac {1}{2}+3{\textit{Oh}} B\right)\big (k^4-k^2\big), \end{equation}

that is, $G(k,t)$ is the solution of

(3.10) \begin{equation} \mathcal{L} G(k,t) = \delta (t) \end{equation}

with initial conditions $G(k,0)=0, \dot {G}(k,0)=1$ . It it found that

(3.11) \begin{equation} G(k,t) = \exp \left (-\frac {\alpha t}{2}\right) \frac {2\sinh \dfrac {\beta t}{2}}{\beta }. \end{equation}

We use $| \boldsymbol{\cdot }|$ to denote the modulus. As in the work of Zhao et al. (Reference Zhao, Sprittles and Lockerby2019) and Zhang et al. (Reference Zhang, Zhao and Si2024), the root mean square (r.m.s.) of $|H|$ is calculated to characterise the strength of the perturbation

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

where we have

(3.13) \begin{equation} \langle |H_{LE}|^2\rangle = |H_i|^2\exp \left (-\alpha t\right)\left (\cosh \frac {\beta t}{2} + \frac {\alpha }{\beta }\sinh \frac {\beta t}{2}\right)^2 \!, \end{equation}

(3.14) \begin{equation} \big \langle |H_{fluc}|^2\big \rangle = \frac {3L{\textit{Oh}}}{\pi }{\textit{Th}}^2k^4 \frac {-\alpha ^2\cosh (\beta t) - \alpha \beta \sinh (\beta t) +\alpha ^2-\beta ^2 + \beta ^2\exp (\alpha t)}{\alpha \beta ^2(\alpha ^2-\beta ^2)\exp (\alpha t)}. \end{equation}

We note that the above expressions for $\langle |H_{LE}|^2\rangle$ and $\langle |H_{fluc}|^2\rangle$ are the same as those in the work of Zhao et al. (Reference Zhao, Sprittles and Lockerby2019) and Zhang et al. (Reference Zhang, Zhao and Si2024), except that in our result the coefficients $\alpha$ and $\beta$ are modified by the parameter $B$ .

As an illustration, we show a comparison of results for two cases: $B=0.0$ and $B=0.1$ . The initial spectrum is set to $H_i=0$ , and $\textit{Oh}=1.0$ , $\textit{Th}=0.1$ , $L=100$ . The evolution of $|H|_{rms}$ in 16 time units is shown in figure 2(a), the corresponding evolution of the surface roughness ( $\varTheta = \sqrt{({ {1}/{\pi L})\int _0^\infty |H|_{rms}^2 {\rm d}k}}$ ) is shown in figure 2(b). It can be seen that the perturbations of small wavenumber ( $k\lt 1$ ) grow faster due to the evaporation–condensation effect. Whereas the perturbations of large wavenumber ( $k\gt 1$ ) grow slower due to the evaporation–condensation effect. The contrasting effect of evaporation–condensation on perturbations of different wavenumber is easy to understand. For a small wavenumber perturbation, the curvature at the trough is larger, hence the evaporation at the trough is enhanced by the Kelvin effect, causing the trough to decrease faster. Whereas for large wavenumber perturbations, the curvature at the crest is larger, hence the evaporation is enhanced at the crest by the Kelvin effect, causing the crest to increase slower. As a result of the dominance of perturbations at $k\lt 1$ the evaporation–condensation leads to an enhancement of perturbation growth. Thus in figure 2(b), the surface roughness grows faster at $B=0.1$ .

Figure 2.

(a) The evolution of $|H|_{rms}$ at $B=0.0$ (dashed lines) and $B=0.1$ (solid lines). The dash-dot line marks the position $k=1$ . (b) The evolution of the surface roughness $\varTheta$ .

4. Instability in the supersaturated vapour environment

Normally, $B$ is very small. For example, for liquid argon at $T=84.09\,\mathrm{K}$ , which will be used in MD simulation below, its density is $\rho =1398\,\mathrm{kg\,m^{-3}}$ , surface tension $\gamma =0.0141\,\mathrm{N\,m^{-1}}$ and the saturation vapour pressure $P_0=56\,\mathrm{kPa}$ (Zhao et al. Reference Zhao, Sprittles and Lockerby2019). If the initial radius of the liquid thread is $h_0=5\,\mathrm{nm}$ , and $\alpha _c$ is assumed to be 1, we can determine the dimensionless parameter $B=( {P_0V_0^2}/{\sqrt {{2\pi m (k_BT)}^{3}}})/\sqrt {h_0/\rho \gamma }\approx 3\times 10^{-4}$ . This $B$ is so small that the effect of evaporation–condensation is completely negligible. We have also examined $B$ for pure water at temperature $0^\circ \mathrm{C} \lt T\lt 100^\circ \mathrm{C}$ . We use the International Association for the Properties of Water and Steam interpolation formula (Wagner & Kretzschmar Reference Wagner and Kretzschmar2008) to determine the surface tension, and the Tetens equation (Monteith & Unsworth Reference Monteith and Unsworth2014) to determine the vapour pressure. If the radius of the liquid thread is $h_0=5\,\mathrm{nm}$ , the liquid density is fixed at $\rho =1000\,\mathrm{kg\,m^{-3}}$ and $\alpha _c$ is assumed to be 1, the dimensionless $B$ can be determined to be monotonically increasing from $6.7\times 10^{-7}$ at $0^\circ \mathrm{C}$ to $6.2\times 10^{-5}$ at $100^\circ \mathrm{C}$ . For liquid propane at $185\,\mathrm{K}$ , if we use the density and surface tension data from Kang (Reference Kang2008) ( $\rho =600\, \mathrm{kg\,m^{-3}}$ , $\gamma =0.01\,\mathrm{N\,m^{-1}}$ ), and vapour pressure determined from our own additional MD simulation ( $P_0\approx 23\,\mathrm{kPa}$ ), it is found that $B\approx 1.3 \times 10^{-4}$ ( $h_0$ is again set to $5\, \mathrm{nm}$ , $\alpha _c$ is assumed to be 1). We can see that $B$ is usually too small to have a considerable effect on the system. Moreover, the dimensional form of $B={\alpha _c P_0V_0^2}/{\sqrt {{2\pi m (k_BT)}^{3}}}$ depends only on the intrinsic properties of the liquid and cannot be easily controlled. Therefore, it is of little practical significance to consider the mass distribution effect caused by the curvature-dependent evaporation in a saturated environment.

On the other hand, the parameter $C$ can be controlled by the environmental vapour pressure $P_v$ ; it can be tuned to be much larger than $B$ and may lead to more significant effects on the instability. Thus, in subsequent discussion we focus on the case $B \ll C \ll 1$ . The effect of the curvature-dependent evaporation is neglected by setting $B=0$ . The basic state ( $\bar {h},\bar {u}$ ) is then determined by setting $B$ and all derivative terms with respect to $z$ in (2.14), (2.15) to zero and solving the two equations with initial conditions $h(z,0)=1,\, u(z,0)=0$ . That is $\bar {h} = 1+Ct$ , $\bar {u}=0$ ( $\bar {h}$ is depicted as the red dashed line in figure 1). In contrast to the previous work of Burelbach et al. (Reference Burelbach, Bankoff and Davis1988) on an evaporating/condensing liquid film, it is straightforward to determine the exact basic state for our problem due to the isothermal assumption. Hence, we can let $h = 1+Ct+\tilde {h}$ , $u=\tilde {u}$ and linearise (2.14), (2.15) to

(4.1) \begin{equation} \frac {\partial \tilde {h}}{\partial t} = -\frac {1+Ct}{2}\tilde {u}^{\prime}, \end{equation}

(4.2) \begin{equation} \frac {\partial \tilde {u}}{\partial t} = \frac {\tilde {h}^{\prime}}{(1+Ct)^2} + \tilde {h}^{\prime\prime\prime} + 3{\textit{Oh}}\tilde {u}^{\prime\prime} + \sqrt {\frac {6{\textit{Oh}}}{\pi }} {\textit{Th}} \frac {\mathcal{N}^{\prime}}{(1+Ct)}. \end{equation}

Eliminating $\tilde {u}$ and applying the finite Fourier transform ((3.4), (3.5)), yields again a second-order ODE

(4.3) \begin{equation} \frac {\partial ^2 H}{\partial t^2} + \left (3{\textit{Oh}} k^2-\frac {C}{1+Ct}\right) \frac {\partial H}{\partial t} + \frac {(1+Ct)^2 k^4-k^2}{2(1+Ct)} H = \sqrt {\frac {3{\textit{Oh}}}{2\pi }}{\textit{Th}} k^2 N. \end{equation}

Using the condition $C\ll 1$ , we let $\varepsilon :=C$ , $\tau := \varepsilon t$ and transform the above equation into

(4.4) \begin{equation} \varepsilon ^2 \frac {\partial ^2 H}{\partial \tau ^2} + \left (\alpha \varepsilon - \frac {\varepsilon ^2}{1+\tau }\right)\frac {\partial H}{\partial \tau } + \frac {(1+\tau)^2k^4-k^2}{2(1+\tau)} H = \sqrt {\frac {3{\textit{Oh}}}{2\pi }}{\textit{Th}} k^2 N, \end{equation}

where $\alpha =3{\textit{Oh}}k^2$ . We first solve the corresponding homogeneous equation

(4.5) \begin{equation} \varepsilon ^2 \frac {\partial ^2 H}{\partial \tau ^2} + \left (\alpha \varepsilon - \frac {\varepsilon ^2}{1+\tau }\right)\frac {\partial H}{\partial \tau } + \frac {(1+\tau)^2 k^4-k^2}{2(1+\tau)} H = 0, \end{equation}

with the Wentzel–Kramers–Brillouin (WKB) theory (Bender & Orszag Reference Bender and Orszag1999). In WKB theory, the solution ansatz is $H(\tau)\sim \exp [S(\tau)/\delta ]$ , $\delta \rightarrow 0$ , where $S(\tau)$ is a phase function that could be imaginary (so both dissipative and dispersive behaviours can be characterised). If the solution is dissipative ( $S(\tau)$ is real), $\delta$ is the thickness of the boundary layer; if the solution is dispersive ( $S(\tau)$ is imaginary), $\delta$ is the characteristic scale of the oscillation wavelength. Here, $S(\tau)$ depends implicitly on $\delta$ , so it is usually expanded as a power series of $\delta$ : $S(\tau) = \sum _n \delta ^n S_n(\tau)$ . Hence, we assume the form of the approximate solution to be

(4.6) \begin{equation} H(\tau)\sim \exp \left [\frac {1}{\delta }\sum _{n=0}^{\infty }\delta ^n S_n(\tau)\right ]\!,\quad \delta \rightarrow 0. \end{equation}

Differentiating $H(\tau)$ twice gives

(4.7) \begin{equation} \frac {\partial H(\tau)}{\partial \tau }\sim \frac {1}{\delta }\sum _{n=0}^{\infty }\delta ^n S^{\prime}_n(\tau)\exp \left [\frac {1}{\delta }\sum _{n=0}^{\infty }\delta ^n S_n(\tau)\right ]\!,\quad \delta \rightarrow 0, \end{equation}

(4.8) \begin{equation} \frac {\partial ^2 H(\tau)}{\partial \tau ^2}\sim \left [\frac {1}{\delta ^2}\left (\sum _{n=0}^{\infty }\delta ^n S^{\prime}_n(\tau)\right)^2 + \frac {1}{\delta }\sum _{n=0}^{\infty }\delta ^nS^{\prime\prime}_n(\tau)\right ]\exp \left [\frac {1}{\delta }\sum _{n=0}^{\infty }\delta ^n S_n(\tau)\right ]\!,\quad \delta \rightarrow 0. \end{equation}

Substituting (4.6), (4.7), (4.8) into (4.5) and dividing off the exponential factors

(4.9) \begin{align} \frac {\varepsilon ^2}{\delta ^2}S_0^{\prime 2}+2\frac {\varepsilon ^2}{\delta }S_0^{\prime}S_1^{\prime} + \frac {\varepsilon ^2}{\delta }S_0^{\prime\prime} + \frac {\varepsilon }{\delta }\alpha S_0^{\prime}+\varepsilon \alpha S_1^{\prime} - \frac {\varepsilon ^2}{\delta }\frac {1}{1+\tau }S_0^{\prime} + {\cdots} = -\frac {(1+\tau)^2 k^4-k^2}{2(1+\tau)}. \end{align}

The largest term on the left-hand side of the above equation is either $( {\varepsilon ^2}/{\delta ^2})S_0^{\prime 2}$ or $( {\varepsilon }/{\delta })\alpha S_0^{\prime}$ , and they must balance the term on the right, so in either case it is required that $\delta =\varepsilon$ . Setting $\delta =\varepsilon$ , and comparing powers of $\varepsilon$ gives a sequence of equations that determines $S_0$ , $S_1$ , ${\cdots}$ . We consider only the first two terms of the sequence: $\varepsilon ^0$

(4.10) \begin{equation} S_0^{\prime 2} + \alpha S_0^{\prime} + \frac {(1+\tau)k^4}{2}-\frac {k^2}{2(1+\tau)} = 0, \end{equation}

and $\varepsilon ^1$

(4.11) \begin{align} 2S_0^{\prime}S_1^{\prime}+S_0^{\prime\prime}+\alpha S_1^{\prime} - \frac {1}{1+\tau }S_0^{\prime} = 0. \end{align}

They are solved order by order. However, $S_0$ cannot be integrated out, so we seek a local approximation around $\tau =0$ and use $1/(1+\tau)= 1-\tau + O(\tau ^2)$ to simplify (4.10). This gives

(4.12) \begin{align} S_0 =-\frac {1}{2}\alpha \tau \pm \frac {\beta ^3}{6k^4} + O(\tau ^2), \end{align}

where $\beta =\sqrt {9{\textit{Oh}}^2k^4 + 2(k^2-k^4)-2(k^2+k^4)\tau }$ . Then $S_1$ can be solved from (4.11) to be

(4.13) \begin{equation} S_1 = \frac {\pm \alpha \log \frac {\zeta -\beta }{\zeta +\beta }+\zeta \log (1+\tau)-\zeta \log \beta }{2\zeta }, \end{equation}

where $\zeta =\sqrt {\alpha ^2+4k^2}$ . Thus, the physical optics approximation of $H$ is given by

(4.14) \begin{equation} \begin{aligned} H_{LE}(t)&\approx A_1\left (\frac {\zeta -\beta }{\zeta +\beta }\right)^{\frac {\alpha }{2\zeta }}\left (1+\varepsilon t\right)^{\frac {1}{2}} \beta ^{-\frac {1}{2}}\exp \left (-\frac {\alpha t}{2}+\frac {\beta ^3}{6\varepsilon (k^4+k^2)}\right)\\ &\quad+ A_2\left (\frac {\zeta -\beta }{\zeta +\beta }\right)^{-\frac {\alpha }{2\zeta }}\left (1+\varepsilon t\right)^{\frac {1}{2}} \beta ^{-\frac {1}{2}}\exp \left (-\frac {\alpha t}{2}-\frac {\beta ^3}{6\varepsilon (k^4+k^2)}\right)\!, \end{aligned} \end{equation}

where $A_1$ and $A_2$ are constants that can be determined with initial conditions.

To evaluate the accuracy of the local approximation we conduct a numerical integration on (4.5) using the Scipy integration function odeint, and compare the result with the prediction of (4.14). The initial spectrum $H_i$ and Oh have been arbitrarily set to $H_i=0.02$ and $\textit{Oh}=1.0$ . The integration is computed up to $Ct=0.4$ , i.e. the average radius of the liquid thread has increased 40 % due to condensation at the end of the integration. The comparison is presented in figure 3(a). The prediction of (4.14) (solid lines) is in good agreement with the numerical integration result (dashed lines). The deviation becomes large when $Ct$ is increased. Figure 3(b) shows that the maximum relative error of (4.14) compared with the numerical integration result is increased up to approximately 25 % at $Ct=0.4$ . But when $Ct\lt 0.2$ the relative error is less than 5 %. Therefore, we expect that our local approximation is accurate when $Ct\lt 0.2$ . The prediction of the classic RP theory is also shown in figure 3 (dotted lines). It can be seen that the inclusion of the condensation effect has significantly impeded the growth of the perturbations. We will discuss this in more detail later.

Figure 3.

(a) The evolution of $|H_{LE}|_{rms}$ . Solid lines: prediction of (4.14); dashed lines: results determined by numerical integration of (4.5); dotted lines: theoretical prediction of the classic RP theory. (b) The maximum relative prediction error of (4.14) compared with the numerical integration result.

To find the solution of the inhomogeneous (4.3), we solve the Green’s function equation

(4.15) \begin{equation} \frac {\partial ^2 G(t,t^{\prime})}{\partial t^2} + \left (\alpha ^2-\frac {\varepsilon }{1+\varepsilon t}\right) \frac {\partial G(t,t^{\prime})}{\partial t} + \frac {(1+\varepsilon t)^2k^4-k^2}{2(1+\varepsilon t)} G(t,t^{\prime}) = \delta (t-t^{\prime}) \end{equation}

with initial conditions:

(4.16) \begin{equation} G(t=0,t^{\prime})=0,\quad \partial G/\partial t (t=0,t^{\prime})=0. \end{equation}

In the regions $t\lt t^{\prime}$ and $t\gt t^{\prime}$ , the solutions have the same form as the solution of the homogeneous equation. In the region $t\lt t^{\prime}$ , the initial conditions (4.16) require that the constants $A_1$ , $A_2$ are both zero. Thus,

(4.17) \begin{equation} G(t,t^{\prime}) = 0, \quad t\lt t^{\prime}. \end{equation}

In the region $t\gt t^{\prime}$ , we write the solution in the form

(4.18) \begin{equation} \begin{aligned} G(t,t^{\prime}) &= A_1(t^{\prime})\left (\frac {\zeta -\beta }{\zeta +\beta }\right)^{\frac {\alpha }{2\zeta }}\left (1+\varepsilon t\right)^{\frac {1}{2}} \beta ^{-\frac {1}{2}} \exp \left (-\frac {\alpha t}{2}+\frac {\beta ^3}{6\varepsilon (k^4+k^2)}\right)\\ &\quad + A_2(t^{\prime}) \left (\frac {\zeta -\beta }{\zeta +\beta }\right)^{-\frac {\alpha }{2\zeta }}\left (1+\varepsilon t\right)^{\frac {1}{2}} \beta ^{-\frac {1}{2}} \exp \left (-\frac {\alpha t}{2}-\frac {\beta ^3}{6\varepsilon (k^4+k^2)}\right)\!,\quad t\gt t^{\prime}. \end{aligned} \end{equation}

The two constants $A_1$ , $A_2$ are determined by patching. At $t=t^{\prime}$ the Green’s function is required to be continuous, and the slope of $G(t,t^{\prime})$ is required to have a finite discontinuity

(4.19) \begin{equation} G(t^{\prime}_{-},t^{\prime})=G(t^{\prime}_{+},t^{\prime}),\quad \frac {\partial G}{\partial t}(t^{\prime}_{+},t^{\prime}) - \frac {\partial G}{\partial t}(t^{\prime}_{-},t^{\prime}) = 1. \end{equation}

Substituting (4.17), (4.18) into (4.19) the constants $A_1(t^{\prime})$ and $A_2(t^{\prime})$ can be determined and we find the solution for $G(t,t^{\prime})$ to be

(4.20) \begin{equation} G(t,t^{\prime}) = \frac {2W(t,t^{\prime})}{D(t^{\prime})\sqrt {\beta }} \exp \left [\frac {-\alpha (t-t^{\prime})}{2}\right ] \sinh \left [\frac {\beta ^3-\beta ^{\prime 3}}{6\varepsilon (k^4+k^2)}-\ln \left (\frac {R(\beta ^{\prime})}{R(\beta)}\right)\right ]\quad t\gt t^{\prime} ,\end{equation}

where

(4.21) \begin{equation} W(t,t^{\prime}) = \beta ^{\prime{\frac {3}{2}}}\sqrt {(1+\varepsilon t)(1+\varepsilon t^{\prime})}, \end{equation}

(4.22) \begin{equation} D(t^{\prime}) = 2k^4(1+\varepsilon t^{\prime})^2-2k^2\left [1-(\varepsilon t^{\prime})^2\right ] - \alpha ^2(1+\varepsilon t^{\prime}) + \alpha \varepsilon , \end{equation}

(4.23) \begin{equation} R(\beta) = \left (\frac {\zeta -\beta }{\zeta +\beta }\right)^{\frac{\alpha }{2\zeta }}, \end{equation}

and $\beta =\sqrt {9{\textit{Oh}}^2k^4 + 2(k^2-k^4)-2(k^2+k^4)\varepsilon t^{\prime}}$ . Then the solution of the inhomogeneous equation is

(4.24) \begin{equation} H_{fluc} = \sqrt {\frac {3{\textit{Oh}}}{2\pi }}{\textit{Th}} k^2 \int _0^t N(k,t-t^{\prime})G(t,t^{\prime}) {\rm d}t^{\prime}. \end{equation}

The r.m.s. of $H_{fluc}$ is given by

(4.25) \begin{equation} |H_{fluc}|_{rms} = \sqrt {\langle |H_{fluc}|^2\rangle }= \sqrt {\frac {3{\textit{Oh}}}{2\pi }{\textit{Th}}^2k^4L\int _0^t G^2(t,t^{\prime}){\rm d}t^{\prime}}. \end{equation}

Substituting (4.14), (4.25) into (3.12) we can obtain the evolution of the r.m.s. of the perturbation spectrum.

5. Instability of a short liquid thread

For a long liquid thread the perturbations develop very fast, the thread may break up before the effects of condensation become very significant. To illustrate the effects of condensation more clearly, we consider a short liquid thread whose dimensionless length is close to the critical length in Rayleigh’s theory: $2\pi$ , so that the breakup time becomes large and the effects of condensation become significant. For a short liquid thread of fixed length $L$ we use a Fourier series to expand the functions $\tilde {h}$ and $\mathcal{N}$

(5.1) \begin{equation} \tilde {h}(z,t) = \sum _{n=-\infty }^{\infty } H^{(n)}(t) \exp \left (nk_0 i z\right), \end{equation}

(5.2) \begin{equation} \mathcal{N}(z,t) = \sum _{n=-\infty }^{\infty } N^{(n)}(t)\exp \left (nk_0 i z\right), \end{equation}

where a periodic boundary is assumed and $k_0=2\pi /L$ ; the expansion coefficients are given by

(5.3) \begin{equation} H^{(n)}(t) = \frac {1}{L}\int _0^L \tilde {h}(z,t)\exp (-nk_0iz) {\rm d}z, \end{equation}

(5.4) \begin{equation} N^{(n)}(t) = \frac {1}{L}\int _0^L \mathcal{N}(z,t)\exp (-nk_0iz) {\rm d}z. \end{equation}

Following the same procedure described in § 4, $H^{(n)}$ is split into two parts: $H^{(n)}=H^{(n)}_{LE}(t)+H^{(n)}_{fluc}(t)$ . And the homogeneous solution $H^{(n)}_{LE}$ is

(5.5) \begin{align} \begin{aligned} H_{LE}^{(n)}(t)&\approx A_1\left (\frac {\zeta -\beta }{\zeta +\beta }\right)^{\frac {\alpha }{2\zeta }}\left (1+\varepsilon t\right)^{\frac {1}{2}} \beta ^{-\frac {1}{2}} \exp \left (-\frac {\alpha t}{2}+\frac {\beta ^3}{6\varepsilon \left [(nk_0)^4+(nk_0)^2\right ]}\right)\\ &\quad + A_2\left (\frac {\zeta -\beta }{\zeta +\beta }\right)^{-\frac {\alpha }{2\zeta }}\left (1+\varepsilon t\right)^{\frac {1}{2}} \beta ^{-\frac {1}{2}} \exp \left (-\frac {\alpha t}{2}-\frac {\beta ^3}{6\varepsilon \left [(nk_0)^4+(nk_0)^2\right ]}\right)\! ,\end{aligned} \end{align}

with $\beta = \sqrt {[9{\textit{Oh}}^2-2(1+\varepsilon ^{\prime} t)] (nk_0)^4 + 2(1-\varepsilon ^{\prime} t)(nk_0)^2}$ , $\zeta =\sqrt {9{\textit{Oh}}^2(nk_0)^4+4(nk_0)^2}$ . The constants $A_1$ and $A_2$ can be determined by the initial conditions $H_{LE}^{(n)}(t=0) = H_i^{(n)}$ , $\partial H_{LE}^{(n)}/\partial t (t=0) = 0$ . The fluctuation-driven part of the Fourier series coefficient is

(5.6) \begin{equation} H^{(n)}_{fluc} = \sqrt {\frac {3{\textit{Oh}}}{2\pi }}{\textit{Th}} (nk_0)^2 \int _0^t N^{(n)}(t-t^{\prime})G(t,t^{\prime}) {\rm d}t^{\prime}, \end{equation}

with the Green’s function $G(t,t^{\prime})$ the same as (4.20) except that all of the $k$ values are replaced by $nk_0$ . For a liquid thread whose dimensionless length is around $2\pi$ , the wavelength with $n=-1$ and $1$ are dominant. So we have

(5.7) \begin{equation} \tilde {h}(z,t) \approx H^{(-1)}(t)\exp (-k_0iz) + H^{(1)}(t)\exp (k_0iz). \end{equation}

Since $H^{(-1)}(t)$ and $H^{(1)}(t)$ are conjugate, the perturbation wave on the liquid thread is of sinusoidal form whose amplitude is

(5.8) \begin{equation} \tilde {h}_{max}(t) = 2|H^{(1)}_{LE}+H^{(1)}_{fluc}|. \end{equation}

Thus, the r.m.s. of $\tilde {h}_{max}(t)$ is

(5.9) \begin{equation} \tilde {h}_{max}(t)|_{rms} = 2\sqrt {\left \langle |H^{(1)}_{LE}|^2\right \rangle + \left \langle |H^{(1)}_{fluc}|^2\right \rangle }. \end{equation}

The first term under the square root in (5.9) is easy to calculate, whereas the second term is given by

(5.10) \begin{equation} \begin{aligned} \left \langle |H^{(1)}_{fluc}|^2 \right \rangle &= \frac {3{\textit{Oh}}}{2\pi L^2}{\textit{Th}}^2 k_0^4 \int _0^t\int _0^t \int _0^L\int _0^L \begin{aligned}[t]\left \langle \mathcal{N}(z,t-t^{\prime})\mathcal{N}(z^{\prime},t-t^{\prime\prime})\right \rangle \exp [-k_0i(z-z^{\prime})]\\\times G(t,t^{\prime})G(t,t^{\prime\prime}){\rm d}z{\rm d}z^{\prime}{\rm d}t^{\prime}{\rm d}t^{\prime\prime} \end{aligned}\\ &= \frac {3{\textit{Oh}}}{2\pi L}{\textit{Th}}^2 k_0^4 \int _0^t G^2(t,t^{\prime}) {\rm d}t^{\prime}, \end{aligned} \end{equation}

where the correlation relation $\langle \mathcal{N}(z,t-t^{\prime}) \mathcal{N}(z^{\prime},t-t^{\prime\prime})\rangle = \delta (z-z^{\prime})\delta (t^{\prime}-t^{\prime\prime})$ has been used for the derivation.

We further try to determine the distribution functions of $|H^{(1)}_{fluc}|^2$ and $\tilde {h}_{max}(t)$ . The expectation of $|H^{(1)}_{fluc}|^2$ is already given by (5.10). To determine the variance of $|H^{(1)}_{fluc}|^2$ , we first calculate the expectation of $|H^{(1)}_{fluc}|^4$

(5.11) \begin{align} \left \langle |H^{(1)}_{fluc}|^4 \right \rangle &=\frac {9{\textit{Oh}}^2}{2\pi ^2L^4}{\textit{Th}}^4 k_0^8 \int _0^t\int _0^t \int _0^t\int _0^t \int _0^L\int _0^L \int _0^L\int _0^L \nonumber\\ &\quad \times \left \langle \mathcal{N}(z,t-t^{\prime})\mathcal{N}(z^{\prime},t-t^{\prime\prime}) \mathcal{N}(z^{\prime\prime},t-t^{\prime\prime\prime})\mathcal{N}(z^{\prime\prime\prime},t-t^{\prime\prime\prime\prime})\right \rangle \nonumber\\ &\quad \times \exp [-k_0i(z-z^{\prime} + z^{\prime\prime}-z^{\prime\prime\prime})] \nonumber\\ &\quad \times G(t,t^{\prime})G(t,t^{\prime\prime})G(t,t^{\prime\prime\prime})G(t,t^{\prime\prime\prime\prime}) \nonumber\\ &\quad \times {\rm d}z{\rm d}z^{\prime}{\rm d}z^{\prime\prime}{\rm d}z^{\prime\prime\prime}{\rm d}t^{\prime}{\rm d}t^{\prime\prime}{\rm d}t^{\prime\prime\prime}{\rm d}t^{\prime\prime\prime\prime} \nonumber\\ &= \frac {9{\textit{Oh}}^2}{2\pi ^2L^2}{\textit{Th}}^4 k_0^8 \left [\int _0^t G^2(t,t^{\prime}) {\rm d}t^{\prime}\right ]^2\! , \end{align}

where the correlation relation (Risken Reference Risken1989)

(5.12) \begin{equation} \begin{aligned} &\left \langle \mathcal{N}(z,t-t^{\prime})\mathcal{N}(z^{\prime},t-t^{\prime\prime})\mathcal{N}(z^{\prime\prime},t-t^{\prime\prime\prime})\mathcal{N}(z^{\prime\prime\prime},t-t^{\prime\prime\prime\prime})\right \rangle \\ &\quad = \delta (z-z^{\prime})\delta (t^{\prime}-t^{\prime\prime})\delta (z^{\prime\prime}-z^{\prime\prime\prime})\delta (t^{\prime\prime\prime}-t^{\prime\prime\prime\prime}) \\ &\quad\quad+ \delta (z-z^{\prime\prime})\delta (t^{\prime}-t^{\prime\prime\prime})\delta (z^{\prime}-z^{\prime\prime\prime})\delta (t^{\prime\prime}-t^{\prime\prime\prime\prime})\\ &\quad\quad + \delta (z-z^{\prime\prime\prime})\delta (t^{\prime}-t^{\prime\prime\prime\prime})\delta (z^{\prime}-z^{\prime\prime})\delta (t^{\prime\prime}-t^{\prime\prime\prime}) \end{aligned} \end{equation}

has been used for the derivation.

Thus, the variance of $|H^{(1)}_{fluc}|^2$ is

(5.13) \begin{equation} \left \langle |H^{(1)}_{fluc}|^4 \right \rangle - \left \langle |H^{(1)}_{fluc}|^2 \right \rangle ^2 = \frac {9{\textit{Oh}}^2}{4\pi ^2L^2}{\textit{Th}}^4 k_0^8 \left [\int _0^t G^2(t,t^{\prime}) {\rm d}t^{\prime}\right ]^2 \!. \end{equation}

It is not a coincidence that the variance of $|H^{(1)}_{fluc}|^2$ is the square of its mean. Because the real and imaginary parts of $H^{(1)}_{fluc}$ are normal random variables with 0 mean and the same variance (let us denote this variance as $\sigma _H^2$ ). We can also prove that any linear combination of them is normally distributed, and the cross-correlation coefficient between them is 0. Hence, they are independent random variables. Then, $|H^{(1)}_{fluc}|^2$ is actually a random variable of scaled Chi-squared distribution with two degrees of freedom: $|H^{(1)}_{fluc}|^2 \sim \sigma _H^2 \chi _2^2$ . Therefore, the mean of $|H^{(1)}_{fluc}|^2$ is $2\sigma _{H}^2$ and its variance is $4\sigma _{H}^4$ . The probability density distribution function of $\psi =|H^{(1)}_{fluc}|^2$ is

(5.14) \begin{equation} P_\psi (\psi ,t) = \frac {1}{2\sigma ^2_{H}}\exp \left (-\frac {\psi }{2\sigma ^2_H}\right)\! ,\end{equation}

with $\sigma ^2_{H}$ determined by comparing $\langle |H^{(1)}_{fluc}|^2\rangle = 2\sigma _H^2$ and (5.10)

(5.15) \begin{equation} \sigma ^2_{H} = \frac {3{\textit{Oh}}}{4\pi L}{\textit{Th}}^2 k_0^4 \int _0^t G^2(t,t^{\prime}) {\rm d}t^{\prime}. \end{equation}

When the initial perturbations are neglected, $\tilde {h}_{max}\approx 2\sqrt {\psi }$ , the probability density function of $\tilde {h}_{max}$ is determined from variable transformation to be

(5.16) \begin{equation} P_h(\tilde {h}_{max},t) \approx \frac {\tilde {h}_{max}}{4\sigma ^2_H}\exp \left (-\frac {\tilde {h}_{max}^2}{8\sigma ^2_H }\right)\!. \end{equation}

6. Molecular dynamics simulation

A molecular dynamics simulation of a liquid argon thread (figure 4) with initial radius $h_0=5\,\mathrm{nm}$ at $T=84.09\, \mathrm{K}$ is conducted to validate our results. The Lennard–Jones 12-6 potential is used to model the interaction among atoms using the energy parameter $\varepsilon _{LJ}=0.9938\,\mathrm{kJ\,mol^{-1}}$ and distance parameter $\sigma _{LJ}=3.405\,Ã$ . The cutoff radius is set to $r_c=15\, Ã$ . The corresponding liquid and saturated vapour density are $\rho =1398\,\mathrm{kg\,m^{-3}}$ and $\rho _g=3.22\,\mathrm{kg\,m^{-3}}$ (Zhao et al. Reference Zhao, Sprittles and Lockerby2019). The kinematic viscosity of the liquid phase is calculated through the Green–Kubo relation (Frenkel & Smit Reference Frenkel and Smit2023) to be: $\nu =1.77\times 10^{-7}\,\mathrm{m^2\,s^{-1}}$ . The surface tension between the liquid and vapour phase is calculated through planar liquid film simulation and the Irving–Kirkwood relation (Kirkwood & Buff Reference Kirkwood and Buff1949; Irving & Kirkwood Reference Irving and Kirkwood1950) to be: $\gamma =1.41\times 10^{-2}\,\mathrm{N\,m^{-1}}$ . The initial configurations of the liquid phase and gas phase are extracted from separate bulk liquid and bulk gas simulations under a Nose–Hoover thermostat. We examine several vapour supersaturation conditions: $S=2.0,\,3.0,\,4.0$ . The corresponding initial vapour densities are $\rho _g = 6.44,\, 9.66,\, 12.88\,\mathrm {kg\,m^{-3}}$ . To sustain the supersaturation condition, a grand canonical Monte Carlo (GCMC) exchange (Frenkel & Smit Reference Frenkel and Smit2023) of atoms with an imaginary ideal gas reservoir with $T=84.09\rm \,K$ , $P_v=1.11,\, 1.67,\, 2.22\,\mathrm{atm}$ is conducted in the marginal regions (thickness $d_{GCMC}=1\,\mathrm{nm}$ ) of the simulation box (see figure 4). The whole system is simulated in canonical ensembles (constant particle number N, volume V, and temperature T; NVT), hence the thermal energy released by condensation is removed from the system and the overall temperature of the system is fixed. We have verified that the supersaturation environment only slightly changes the surface tension, thus, $\gamma =1.41\times 10^{-2}\,\mathrm{N\,m^{-1}}$ will be used throughout in this article. Periodic boundary condition is applied to all 3 directions of the simulation box. Before the NVT simulation, a very short period ( $0.03\,\mathrm{ns}$ ) of Langevin dynamics simulation is also performed. We change the seed of the random number for the Langevin simulation to repeat the whole simulation many times and get an ensemble average for each simulation condition. We use the open source LAMMPS framework (Thompson et al. Reference Thompson2022) to perform all the MD simulations.

Figure 4.

Schematic of the molecular simulation model.

For clarity, we have summarised the physical parameters used in and derived from the MD simulations in table 1. The upper part of the table shows the dimensional parameters ( $T$ , $h_0$ , $\nu$ , $\gamma$ , $\rho$ , $P_v$ ), the lower part shows the corresponding dimensionless parameters ( $\textit{Oh}$ , $\textit{Th}$ , $S$ ). The supersaturation ratio $S$ is changed due to the variation of the vapour pressure $P_v$ . The dimensionless thread length $L$ is also changed to examine both very long and short liquid threads.

Table 1.

Physical parameters used in and derived from the MD simulations.

For analysis, the simulation domain is divided into a series of thin slabs along the $z$ direction. Particle density $\rho _a(z_i,t)$ (including particles in both liquid and vapour phases) in each slab is recorded during the simulation, where $z_i$ denotes the $z$ coordinate of the centre of the $i$ th slab. The vapour density $\rho _g(z_i,t)$ is also recorded. Then the radius of the liquid thread at $z_i$ is calculated by

(6.1) \begin{equation} h(z_i,t)=\sqrt {\frac {L_xL_y(\rho _a-\rho _g)}{\pi (\rho -\rho _g)}}, \end{equation}

where $\rho$ is the density of the liquid phase, which we assume to be constant: $1398\,\mathrm{kg\,m^{-3}}$ . Discrete Fourier transform is then applied to the sequence of $h(z_i,t)$ to get the spectrum $H(k,t)$ .

The LAMMPS input scripts and the pre- and post-processing codes used in this study are publicly available on Zenodo (Mo Reference Mo2025).