1. Introduction
The investigation of the instability of liquid jets due to surface tension has a long history, tracing back to the pioneering work of Plateau (Reference Plateau1873) and Rayleigh (Reference Rayleigh1878, Reference Rayleigh1892). This instability phenomenon also holds significant importance in understanding the dynamics of liquid films coated on solid fibres, with additional complexities at the liquid–solid interface. This field has attracted significant scientific attention (Quéré Reference Quéré1999), owing to its critical relevance across various technological domains, including additive manufacturing (Deng et al. Reference Deng, Nave, Liang, Johnson and Fink2011; Oliveira, Santos & Miranda Reference Oliveira, Santos and Miranda2020), droplet transport (Lee et al. Reference Lee, Chan, Carlson and Dalnoki-Veress2022), chemical element extraction (Chen et al. Reference Chen, Yang, Zheng, Temprano-Coleto, Dong, Cheng, Yao, Stone, Hu and Ren2023) and water collection through fog harvesting (Chen et al. Reference Chen, Ran, Gan, Zhou, Zhang, Zhang, Zhang and Jiang2018; Zhang et al. Reference Zhang, Zheng, Zhu, Zhu, Si and Xu2022).
The instability of a film coated on a fibre has been extensively studied across various stages. The early time behaviours can be described by the Rayleigh–Plateau instability (Quéré Reference Quéré1999), showing that a film with an outer radius 
$h_0$ becomes unstable to sufficiently long-wavelength disturbances, specifically 
$\lambda > \lambda _{crit} = 2 {\rm \pi}h_0$. Here, 
$\lambda _{crit}$ represents the critical wavelength beyond which the instability ceases to grow. Additionally, the dominant (most unstable/fast growing) modes are influenced by the ratio of 
$h_0$ to the fibre radius 
$a$, validated by experiments (Goren Reference Goren1964). Subsequent nonlinear evolution is modelled using the lubrication approximation (Hammond Reference Hammond1983), resulting in a leading-order lubrication equation that is applicable to films on both the inside and outside of a cylinder. Due to its simplicity, this lubrication equation and its higher-order versions (Craster & Matar Reference Craster and Matar2006; Ruyer-Quil et al. Reference Ruyer-Quil, Treveleyan, Giorgiutti-Dauphiné, Duprat and Kalliadasis2008) have been used in studying various interface dynamics of annular films on fibres. Examples include the transition from absolute unstable regimes to convective ones (Kliakhandler, Davis & Bankoff Reference Kliakhandler, Davis and Bankoff2001; Duprat et al. Reference Duprat, Ruyer-Quil, Kalliadasis and Giorgiutti-Dauphiné2007; Craster & Matar Reference Craster and Matar2009) and capillary drainage involving complex interactions of ‘lobes’ and ‘collars’ formed on the interfaces (Lister et al. Reference Lister, Rallison, King, Cummings and Jensen2006). Recently, these lubrication models have been extended to encompass more complex scenarios by incorporating other physics, such as electric fields (Ding et al. Reference Ding, Xie, Wong and Liu2014), heat transfer (Zeng et al. Reference Zeng, Sadeghpour, Warrier and Ju2017), thermal fluctuations (Zhang, Sprittles & Lockerby Reference Zhang, Sprittles and Lockerby2021) and Van der Waals forces (Tomo, Nag & Takamatsu Reference Tomo, Nag and Takamatsu2022).
One of the important physical factors is liquid–solid slip, which has recently attracted substantial research interest (Secchi et al. Reference Secchi, Marbach, Niguès, Stein, Siria and Bocquet2016; Zhang, Sprittles & Lockerby Reference Zhang, Sprittles and Lockerby2020; Kavokine, Netz & Bocquet Reference Kavokine, Netz and Bocquet2021; Kavokine, Bocquet & Bocquet Reference Kavokine, Bocquet and Bocquet2022) and been found to influence the dynamics of various interfacial flows (Liao et al. Reference Liao, Li, Chang, Huang and Wei2014; Halpern, Li & Wei Reference Halpern, Li and Wei2015; Martínez-Calvo, Moreno-Boza & Sevilla Reference Martínez-Calvo, Moreno-Boza and Sevilla2020a; Zhao et al. Reference Zhao, Liu, Lockerby and Sprittles2022). For the case of cylindrical films, Ding & Liu (Reference Ding and Liu2011) introduced a lubrication equation incorporating slip conditions to investigate the instability of films descending along porous vertical fibres. Their findings revealed that the instability is amplified by the presence of a fluid-porous interface, which is modelled using a slip boundary condition. Regarding annular films within slippery tubes, Liao, Li & Wei (Reference Liao, Li and Wei2013) numerically solved a lubrication equation with leading-order terms, demonstrating that even a fractional amount of wall slip significantly exaggerates the instability, leading to considerably faster drainage compared with the no-slip scenario (Hammond Reference Hammond1983). Haefner et al. (Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015) conducted experimental investigations into the influence of slip on the instability for films coated on horizontal fibres. Similar to the observations of both Ding & Liu (Reference Ding and Liu2011) and Liao et al. (Reference Liao, Li and Wei2013), the wall slip was found to enhance the instability, resulting in increased growth rates of perturbations. The experimental results were also shown to match predictions of a slip-modified lubrication equation. Halpern & Wei (Reference Halpern and Wei2017) subsequently illustrated how wall slip can amplify drop formation in a film descending a vertical fibre. This observation provides a plausible explanation for the discrepancy between experimentally predicted and theoretically derived critical Bond numbers for drop formation. More recent investigations delved into the dynamics of films on slippery fibres within non-isothermal conditions (Chao, Ding & Liu Reference Chao, Ding and Liu2018) and under the influence of intermolecular forces (Ji et al. Reference Ji, Falcon, Sadeghpour, Zeng, Ju and Bertozzi2019), employing more intricate lubrication models. Despite the extensive use of slip-modified lubrication models, their constraints have been exposed by Zhao, Zhang & Si (Reference Zhao, Zhang and Si2023a) through linear instability analysis applied to the axisymmetric Stokes equations. The theoretical framework not only highlights an overestimation of the slip-enhanced perturbation growth rate as compared with classical lubrication models, but also reveals a slip-dependent dominant wavelength, deviating from the constant value posited by prior works using the lubrication method (Liao et al. Reference Liao, Li and Wei2013; Haefner et al. Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015; Halpern & Wei Reference Halpern and Wei2017; Chao et al. Reference Chao, Ding and Liu2018).
Noticeably, in modern applications such as additive manufacturing in space (Reitz et al. Reference Reitz2021; Van Ombergen et al. Reference Van Ombergen2023) and three-dimensional printing with liquid metals (Assael et al. Reference Assael, Kalyva, Antoniadis, Michael Banish, Egry, Wu, Kaschnitz and Wakeham2010; Kondic et al. Reference Kondic, González, Diez, Fowlkes and Rack2020), the role of inertia in governing the dynamics of liquid film interfaces has become notably apparent, contrasting with the predominant neglect of inertia in most preceding investigations. One exception is the work done by Goren (Reference Goren1962), who introduced inertial effects by conducting instability analysis for the full Navier–Stokes (NS) equations. The theoretical findings elucidated distinctions between two limiting cases, i.e. the inviscid case and viscous case without inertia. However, the impact of inertia within the intermediary regime between these two limits remained uncharted. Ding et al. (Reference Ding, Wong, Liu and Liu2013) proposed two coupled equations governing the film thickness and flow rate to study instability and dynamics of a film on a fibre, considering both inertia and slip. Nonetheless, their focus was primarily on scenarios characterized by small to moderate Reynolds numbers, where the influence of inertia on interface dynamics seemed less pronounced. While recent studies have extensively investigated the influence of inertia on the instability of planar films (González, Diez & Sellier Reference González, Diez and Sellier2016; Moreno-Boza, Martínez-Calvo & Sevilla Reference Moreno-Boza, Martínez-Calvo and Sevilla2020a), its effect on the cylindrical films, especially on a fibre with slip, remains unclear, thus motivating the present investigations.
In this work, linear instability analysis of the axisymmetric NS equations is performed to investigate inertia and slip effects on the dynamics of a liquid film on a fibre. Direct numerical simulations of the NS equations are also employed to confirm the theoretical findings and provide more physical insights. The paper is laid out as follows. Non-dimensionalised governing equations for a film on a fibre are introduced in § 2. Linear instability analysis for the governing equations is performed in § 3, where the dispersion relation is derived in § 3.1, followed by two limiting cases: jet flows in § 3.2 and film flows without inertia in § 3.3, respectively. Predictions arising from the theoretical model are presented in § 3.4. Subsequently, direct numerical simulations are performed in § 4. These simulations are compared with the predictions of the theoretical model, specifically concerning the influence on the dominant mode (§ 4.1) and the growth rate (§ 4.2) of perturbations. Nonlinear dynamics extracted from the simulations is also analysed in § 4.2.
2. Model formulation
We consider a Newtonian liquid film on a fibre of the radius 
$a$ with the 
$z$ axis along the centreline (figure 1). The initial radius of the film measured from the 
$z$ axis is 
$r = h_0$. Additionally, gravity is neglected, and we assume uniform external pressure and surface tension.
Schematic of a liquid film on a slippery fibre.
The incompressible NS equations are employed to predict the dynamics of the flow inside the liquid film. To identify the governing dimensionless parameters, we non-dimensionalise the NS equations with the rescaling variables shown below:
\begin{equation} (r,z,h)= \frac{(\tilde{r},\tilde{z},\tilde{h})}{h_0}, \quad (u,v,w) = \frac{(\tilde{u},\tilde{v},\tilde{w})}{U} , \quad t =\frac{U}{h_0} \tilde{t}, \quad p = \frac{h_0}{\gamma}\tilde{p}. \end{equation}
Here 
$\tilde {h}$, 
$\tilde {t}$ and 
$\tilde {p}$ represents the dimensional interface height, time and pressure, respectively (note that the dimensional material parameters are not given tildes); 
$(r,\phi,z)$ are the cylindrical coordinates with the corresponding velocities (
$u,v,w$); 
$U$ represents the characteristic velocity inside the film and 
$\gamma$ is the surface tension of the liquid–gas interface. After eliminating all the derivatives with respect to 
$\phi$ and setting 
$v=0$ in the cylindrical coordinates, the axisymmetric incompressible equations can be written as
$$\begin{gather} \frac{\partial w}{\partial z} + \frac{1}{r} \frac{\partial(ur)}{\partial r} = 0 , \end{gather}$$
$$\begin{gather}\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial r} + w \frac{\partial u}{\partial z} =- \frac{1}{We} \frac{\partial p}{\partial r} + \frac{1}{Re} \left[ \frac{\partial^2 u}{\partial z^2}+ \frac{\partial}{\partial r} \left(\frac{1}{r} \frac{\partial (ur)}{\partial r}\right) \right], \end{gather}$$
$$\begin{gather}\frac{\partial w}{\partial t} + u \frac{\partial w}{\partial r} + w \frac{\partial w}{\partial z} =-\frac{1}{We} \frac{\partial p}{\partial z} + \frac{1}{Re}\left[ \frac{\partial^2 w}{\partial z^2}+ \frac{1}{r}\frac{\partial}{\partial r} \left( r \frac{\partial w}{\partial r} \right) \right] , \end{gather}$$
where the non-dimensional quantity 
$We = \rho U^2 h_0/\gamma$ is the Weber number, which relates the inertial force to the capillary force. Here 
$Re = \rho U h_0/\mu$ is the Reynolds number, showing the ratio between inertial force and the viscous force; 
$\mu$ is the liquid dynamic viscosity and 
$\rho$ is the liquid density.
Since the density of gas around the film is much smaller than that of the liquid, the gas flow outside can be assumed to be dynamically passive to simplify the problem. The liquid–gas interface height 
$h(z, t)$ satisfies the kinematic boundary condition
\begin{equation} \frac{\partial h}{\partial t} + w \frac{\partial h}{\partial z} -u = 0. \end{equation}
The normal stress balance at the interface 
$r=h$ gives
\begin{equation} p- \frac{We}{Re} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{n}, \end{equation}
where 
$\tau$ is the shear stress, which is proportional to the strain rate in Newtonian fluids. Here 
$\boldsymbol {n}$ is the outward normal and 
$\boldsymbol {\nabla } \boldsymbol{\cdot} \boldsymbol {n}$ represents the dimensionless Laplace pressure. The tangential force balance is
\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{t}=0 , \end{equation}
where 
$\boldsymbol {t}$ is the tangential vector. With 
$\boldsymbol {n}$ and 
$\boldsymbol {t}$ expressed in terms of the unit vectors in the 
$z$ direction (
$\hat {\boldsymbol {e}}_z$) and 
$r$ direction (
$\hat {\boldsymbol {e}}_r$),
\begin{align} \boldsymbol{n} =-\frac{\partial_z h}{\sqrt{1+(\partial_z h)^2}} \hat{\boldsymbol{e}}_z + \frac{1}{\sqrt{1+(\partial_z h)^2}} \hat{\boldsymbol{e}}_r \quad\text{and}\quad \boldsymbol{t} = \frac{1}{\sqrt{1+(\partial_z h)^2}} \hat{\boldsymbol{e}}_z + \frac{\partial_z h}{\sqrt{1+(\partial_z h)^2}} \hat{\boldsymbol{e}}_r, \end{align}(2.6) and (2.7) explicitly give
$$\begin{gather} p - \frac{We}{Re}\frac{2}{1+ \left(\partial_z h \right)^2} \left[ \frac{\partial u}{\partial r} - \frac{\partial h}{\partial z} \left( \frac{\partial w}{\partial r} + \frac{\partial u}{\partial z}\right) + \left(\frac{\partial h}{\partial z} \right)^2 \frac{\partial w}{\partial z} \right] \nonumber\\ = \frac{1}{h \sqrt{1+\left(\partial_z h \right)^2}} - \frac{\partial_z^2 h}{\left(1+\left(\partial_z h \right)^2\right)^{{3}/{2}}} \end{gather}$$for the normal forces, and
\begin{equation} 2 \frac{\partial h}{\partial z} \left( \frac{\partial u}{\partial r} - \frac{\partial w}{\partial z} \right) + \left[ 1- \left( \frac{\partial h}{\partial z} \right)^2 \right] \left( \frac{\partial w}{\partial r} + \frac{\partial u}{\partial z} \right) = 0 \end{equation}
for the tangential forces. Here 
$\partial _z$ and 
$\partial _z^2$ refers to the first and second axial derivatives. In terms of the boundary conditions at the fibre surface 
$r = \alpha$ (
$\alpha$ is the dimensionless fibre radius, i.e. 
$\alpha = a/h_0$), we introduce the Navier slip boundary condition (Navier Reference Navier1823) in the tangential direction and the no-penetration boundary condition in the normal direction such that
$$\begin{gather} w = l_s \frac{\partial w}{\partial r}, \end{gather}$$
$$\begin{gather}u = 0. \end{gather}$$
Here, 
$l_s$ represents the dimensionless slip length rescaled by 
$h_0$. The governing equations above can be further simplified under the long-wave approximation to give lubrication equations. Similar to the previous work on planar films (Münch, Wagner & Witelski Reference Münch, Wagner and Witelski2005), different rescaling for the leading orders give different forms of the lubrication model. When inertia is neglected (
$Re \ll 1$), the one-equation model, which has been widely used in the previous works (Craster & Matar Reference Craster and Matar2006; Haefner et al. Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015; Zhang et al. Reference Zhang, Sprittles and Lockerby2020; Zhao et al. Reference Zhao, Zhang and Si2023a), can be obtained from the Stokes equations. The dimensionless format of the no-slip lubrication model is
\begin{equation} \frac{\partial h}{\partial t} = \frac{1}{ h} \frac{\partial}{\partial z} \left[ M(h) \frac{\partial }{\partial z} \left (\frac{1}{h} - \frac{\partial^2 h}{\partial z^2} \right)\right], \end{equation}
where 
$M(h)$ is the mobility term
\begin{equation} M(h) =\frac{1}{16} \left[-3h^4-\alpha^4 + 4 \alpha^2 h^2 + 4 h^4 \ln\left(\frac{h}{\alpha}\right) \right]. \end{equation}
When inertia is not negligible and 
$l_s \gg 1$, we propose another lubrication model consisting of two equations (see Appendix A for the derivation). The dimensionless format of this giant-slip lubrication model is
\begin{cases}
\frac{\partial h^2}{\partial t} + \frac{\partial \left(h^2 w \right)}{\partial z} &= 0, \\
\frac{\partial w}{\partial t}+ w \frac{\partial w}{\partial z} &= -\frac{1}{ We}
\frac{\partial}{\partial z}\left(\frac{1}{h}-\frac{\partial^2h}{\partial z^2} \right)+\frac{1}{ Re}
\frac{3}{h^2-\alpha^2} \frac{\partial \left( h^2 \partial_z w\right)}{\partial z} \nonumber\\
&\quad - \frac{1}{ Re} \frac{2 \alpha^2}{h^2-\alpha^2} \left(\frac{\partial^2 w}{\partial
z^2}-\frac{w}{ \alpha\, l_{s}} \right).
\end{cases}
When 
$\alpha =0$, the lubrication model for the jet flows (Eggers & Dupont Reference Eggers and Dupont1994) is recovered.
3. Instability analysis
In this section, linear instability analysis based on ((2.2)–(2.11)) is performed using the normal mode method, which has been widely used for the instability in different fluid configurations (Rayleigh Reference Rayleigh1878; Tomotika Reference Tomotika1935; Craster & Matar Reference Craster and Matar2006; Li, Yin & Yin Reference Li, Yin and Yin2008; Si et al. Reference Si, Li, Yin and Yin2009; Liang et al. Reference Liang, Deng, Nave and Johnson2011; González et al. Reference González, Diez and Sellier2016).
3.1. Derivation for the dispersion relation
To perform instability analysis, the dimensionless perturbed quantities are set as
\begin{align} u(r,z,t) = \hat{u}(r) {\rm e}^{\omega t+{\rm i}kz},\quad w(r,z,t) = \hat{w}(r) {\rm e}^{\omega t+{\rm i}kz} \quad \text{and} \quad p(r,z,t) = 1 + \hat{p}(r) {\rm e}^{\omega t+{\rm i}kz}, \end{align}
where 
$\omega$ is the growth rate of perturbations and 
$k$ is the wavenumber. Here, we assume that there is no base flow inside the film. The perturbed quantities are linearly decomposed into the pressure term and the viscosity term, 
$\hat {u} = \hat {u}_p + \hat {u}_{\nu }$ and 
$\hat {w} = \hat {w}_p + \hat {w}_{\nu }$.
For the pressure term, the velocity potential 
$\phi$ (i.e. 
$\partial _r \phi = \hat {u}_p\text { and } \partial _z \phi = \hat {w}_p$) is introduce to simplify the problem. The mass equation (2.2) becomes a zero-order Bessel equation
\begin{equation} \frac{\mathrm{d}}{\mathrm{d} r} \left(r \frac{\mathrm{d} \phi}{\mathrm{d} r} \right) - k^2 r \phi= 0, \end{equation}whose solution can be expressed in terms of Bessel functions
\begin{equation} \phi = A_1 \mathrm{I}_0(kr) + B_1 \mathrm{K}_0(kr). \end{equation}
Here 
$\mathrm {I}_0$ and 
$\mathrm {K}_0$ are zero-order modified Bessel functions of the first and second kinds; 
$A_1$ and 
$B_1$ are arbitrary constants awaiting determination. Calculating the derivatives of 
$\phi$ gives the solution of 
$\hat {u}_p$ and 
$\hat {w}_p$,
$$\begin{gather} \hat{u}_p = k\left[ A_1 \mathrm{I}_1(kr) - B_1 \mathrm{K}_1(kr)\right], \end{gather}$$
$$\begin{gather}\hat{w}_p = {\rm i}k\left[ A_1 \mathrm{I}_0(kr) + B_1 \mathrm{K}_0(kr)\right]. \end{gather}$$
Here 
$\mathrm {I}_1$ and 
$\mathrm {K}_1$ are first-order modified Bessel functions of the first and second kinds. Substituting (3.3) into momentum equation (2.4) yields the solution of 
$\hat {p}$, expressed as
\begin{equation} \hat{p} =-We \omega \left[ A_1 \mathrm{I}_0(kr) + B_1 \mathrm{K}_0(kr)\right] . \end{equation}Considering the viscosity term of perturbed quantities, we simplify the momentum equation (2.3) as a first-order Bessel equation,
\begin{equation} \frac{\mathrm{d}^2 \hat{u}_{\nu}}{\mathrm{d} r^2} + \frac{1}{r}\frac{\mathrm{d} \hat{u}_{\nu}}{\mathrm{d} r} - \left(Re \omega + k^2+\frac{1}{r^2} \right) \hat{u}_{\nu} = 0. \end{equation}
So the solution of 
$\hat {u}_{\nu }$ is
\begin{equation} \hat{u}_\nu = A_2 \mathrm{I}_1(lr) + B_2 \mathrm{K}_1(lr), \end{equation}
where 
$l^2 = k^2+Re\,\omega$; 
$A_2$ and 
$B_2$ are another two arbitrary constants. According to (2.2),
\begin{equation} \hat{w}_\nu = \frac{il}{k} \left[ A_2 \mathrm{I}_0(lr) - B_2 \mathrm{K}_0(lr) \right]. \end{equation}Combining the pressure parts and viscosity parts gives us the general solution of perturbed variables, namely
\begin{cases}
\hat{u} &= A_1 k \mathrm{I}_1(kr) + A_2 \mathrm{I}_1(lr) - B_1 k \mathrm{K}_1(kr) + B_2
\mathrm{K}_1(lr) ,\\
\hat{w} &= {\rm i} \left[ A_1 k \mathrm{I}_0(kr) + A_2 l \mathrm{I}_0(lr)/k + B_1 k
\mathrm{K}_0(kr) - B_2 l \mathrm{K}_0(lr)/k \right] ,\\
\hat{p} &= -We\,\omega \left[ A_1 \mathrm{I}_0(kr) + B_1 \mathrm{K}_0(kr)\right].
\end{cases} The dimensionless perturbed quantities for 
$\hat {u}$, 
$\hat {v}$, 
$\hat {p}$, combined with 
$h(x, t) = 1 + \hat {h} {\rm e}^{\omega t + {\rm i}kx}$, are then substituted into the boundary equations (2.5)–(2.11). For the boundary conditions at the interface (
$r = 1$), their linearisation gives
$$\begin{gather} \frac{\mathrm{d} \hat{w} }{\mathrm{d} r }+ {\rm i}k \hat{u} = 0, \end{gather}$$
$$\begin{gather}\hat{p} - 2 \frac{We}{Re}\frac{\mathrm{d} \hat{u}}{ \mathrm{d} r } = \hat{h} \left( k^2-1 \right) , \end{gather}$$
$$\begin{gather}\omega \hat{h} = \hat{u}. \end{gather}$$
Furthermore, for the boundary conditions on the fibre surface (
$r = \alpha$), their linearised forms are
$$\begin{gather} \hat{w} = l_{s}\frac{\mathrm{d} \hat{w}}{\mathrm{d} \hat{r}} , \end{gather}$$
$$\begin{gather}\hat{u} = 0 . \end{gather}$$
According to (3.15), 
$\hat {h}$ in (3.14) can be eliminated to give the final four equations of the boundary conditions, i.e. (3.13), (3.14), (3.16) and (3.17). Substituting the Bessel functions ((3.10)–(3.12)) into these perturbed equations leads to a homogeneous system of linear equations for 
$A_1$, 
$A_2$, 
$B_1$ and 
$B_2$, which has a non-trivial solution only if the determinant of the coefficients vanishes. In this way, we have the final equation
\begin{equation} \left| \renewcommand\arraystretch{1.3} \begin{matrix} k \mathrm{I}_1(k\alpha) & \mathrm{I}_1(l \alpha) & - k \mathrm{K}_1(k\alpha) & \mathrm{K}_1(l \alpha) \\ F_{21} & F_{22} & F_{23} & F_{24} \\ 2 k^3 \mathrm{I}_1(k ) & (k^2+l^2) \mathrm{I}_1(l) & -2 k^3 \mathrm{K}_1(k) & (k^2+l^2) \mathrm{K}_1(l) \\ F_{41} & F_{42} & F_{43} & F_{44} \\ \end{matrix} \right|=0, \end{equation}where
\begin{equation} \left. \begin{aligned} & F_{21} = k^2 \mathrm{I}_0(k \alpha)-l_{s} k^3 \mathrm{I}_1(k\alpha) , \\ & F_{22} = l \mathrm{I}_0(l \alpha)-l_{s} l^2 \mathrm{I}_1(l \alpha) , \\ & F_{23} = k^2 \mathrm{K}_0(k \alpha)+l_{s} k^3 \mathrm{K}_1(k \alpha) ,\\ & F_{24} =-l \mathrm{K}_0(l \alpha) - l_{s} l^2 \mathrm{K}_1(l \alpha) , \\ & F_{41} = \mathrm{I}_0(k) \omega^2 + 2 Re^{-1} k^2 \mathrm{I}'_1(k) \omega + We^{-1} (k^2-1) k \mathrm{I}_1(k) ,\\ & F_{42} = 2 Re^{-1} l \mathrm{I}'_1(l) \omega + We^{-1} (k^2-1) \mathrm{I}_1(l) ,\\ & F_{43} = \mathrm{K}_0(k) \omega^2 - 2 Re^{-1} k^2 \mathrm{K}'_1(k) \omega - We^{-1} (k^2-1) k \mathrm{K}_1(k) ,\\ & F_{44} = 2 Re^{-1} l \mathrm{K}'_1(l) \omega + We^{-1} (k^2-1) \mathrm{K}_1(l). \end{aligned} \right\} \end{equation}
As 
$\omega$ occurs in the argument of some Bessel functions, such as 
$\mathrm {I}_1 (l \alpha )$, (3.18) cannot be solved explicitly for 
$\omega$, except in two limiting cases, which are presented in the following subsections.
3.2. Limiting case of jet flows
When the viscosity of the liquid is neglected (
$Re\rightarrow \infty$), the viscosity terms governing the instability become zero, i.e. 
$u_{\nu }=w_{\nu }=0$. As a result, the dispersion relation is simplified to
\begin{equation} \left| \renewcommand\arraystretch{1.3} \begin{matrix} k \mathrm{I}_1(k \alpha) & -k \mathrm{K}_1(k \alpha) \\ \mathrm{I}_0(k) \omega^2 + We^{-1} (k^2-1) k \mathrm{I}_1(k) & \mathrm{K}_0(k) \omega^2 - We^{-1} (k^2-1) k \mathrm{K}_1(k) \\ \end{matrix} \right|=0. \end{equation}
Here 
$\omega$ can be expressed explicitly as
\begin{equation} \omega = \sqrt{\frac{(1-k^2) k}{We}\frac{ \left[ \mathrm{K}_1(k\alpha) \mathrm{I}_1(k) - \mathrm{I}_1(k\alpha) \mathrm{K}_1(k) \right]}{\mathrm{I}_1(k\alpha) \mathrm{K}_0(k)+ \mathrm{K}_1(k\alpha) \mathrm{I}_0(k)}}. \end{equation}
As the fibre radius approaches infinitesimally small values, the flows inside the film are expected to resemble jet flows. Consequently, substituting 
$\alpha = 0$ into (3.21) yields the dispersion relation for the instability of inviscid jets, originally proposed by Rayleigh (Reference Rayleigh1878), which is expressed as
\begin{equation} \omega = \sqrt{ \frac{(1-k^2) k}{We} \frac{ \mathrm{I}_1(k)}{\mathrm{I}_0(k)}}. \end{equation} When viscosity is taken into account, relying solely on the condition 
$\alpha =0$ is no longer adequate to simplify (3.18) into the dispersion relation of jet flows. Hence, it becomes imperative to introduce ultra-slip boundary conditions (
$l_s \rightarrow \infty$), resulting in
\begin{equation} \left| \renewcommand\arraystretch{1.3} \begin{matrix} 2 k^3 \mathrm{I}_1(k ) & (k^2+l^2) \mathrm{I}_1(l) \\ F_{41} & F_{42} \\ \end{matrix} \right|=0. \end{equation}This relation can be rearranged as
\begin{equation} \omega^2 + \frac{2}{Re}\frac{k^2}{\mathrm{I}_0(k)} \left[\mathrm{I}_1'(k) -\frac{2 kl \mathrm{I}_1(k) \mathrm{I}_1'(l)}{(l^2+k^2) \mathrm{I}_1(l)}\right] \omega = \frac{1}{We} (1- k^2 )k \frac{\mathrm{I}_1(k)}{\mathrm{I}_0(k)} \frac{l^2-k^2}{l^2+k^2}, \end{equation}
which is presented by Goldin et al. (Reference Goldin, Yerushalmi, Pfeffer and Shinnar1969). Using 
$(l^2+k^2)/(l^2-k^2) =1+2 k^2/ ( Re \omega )$ and 
$\mathrm {I}'_1(k) = \mathrm {I}_0(k)-\mathrm {I}_1(k)/k$, we obtain the equivalent representation of (3.24),
\begin{equation} \omega^2 + \frac{2 k^2}{Re} \left[2-\frac{\mathrm{I}_1(k)}{k \mathrm{I}_0(k)}+\frac{2 k^2}{l^2-k^2}\left( 1 -\frac{l\, \mathrm{I}_1(k) \mathrm{I}_0(l)}{k\, \mathrm{I}_0(k) \mathrm{I}_1(l)} \right)\right] \omega = \frac{1}{We} (1- k^2 )k \frac{\mathrm{I}_1(k)}{\mathrm{I}_0(k)} , \end{equation}
which is a widely used form of the dispersion relation for the temporal instability of a viscous Newtonian jet, first proposed by Weber (Reference Weber1931). When 
$Re \rightarrow \infty$, (3.22) is also recovered.
These findings are further confirmed through numerical solutions of (3.18) using the FindRoot function of MATHEMATICA. In the analysis, the capillary velocity is adopted as the characteristic velocity for non-dimensionalisation, namely 
$U=\gamma /\mu$. Consequently, we arrive at 
$Re=We=Oh^{-2}$, where the non-dimensional quantity 
$Oh = \mu /\sqrt {\rho \gamma h_0}$ represents the Ohnesorge number, serving as a linkage between viscous forces, inertial forces and surface-tension forces. For inviscid flows, we set 
$Oh=10^{-3}$. As depicted in figure 2(a), the results from the NS dispersion relation (3.18) gradually converge towards the predictions of Rayleigh's model as the fibre radius diminishes. This outcome is consistent with the theoretical analysis, thus further validating the numerical solutions of (3.18). Although the asymptotic behaviours of inviscid cases are realised on the no-slip boundary condition, for the viscous cases, they only manifest when an ultra-slippery fibre is considered, as elucidated in figure 2(b). Here, the solutions of (3.18) converge towards Goldin's model (3.24) as 
$l_s$ increases. This divergence can be attributed to the differential impacts of shear stresses from the fibre surfaces, influenced by the slip length.
The dispersion relation between the growth rate 
$\omega$ and the wavenumber 
$k$ for the limiting cases of inviscid and viscous fluids. (a) The inviscid liquid film (
$Oh=10^{-3}$) on fibres with different radii: 
$\alpha =0.8$ (green dash-dotted line), 
$0.5$ (blue dashed line), 
$0.2$ (red dotted line), 
$0.01$ (black solid line). (b) The viscous liquid film (
$Oh=0.5$) on extremely thin fibres (
$\alpha =0.01$) with different slip lengths: 
$l_s=0.0$ (green dash-dotted line), 
$0.1$ (blue dashed line), 
$1.0$ (red dotted line), 
$10$ (black solid line). The circles are the predictions of Rayleigh (Reference Rayleigh1878) and Goldin et al. (Reference Goldin, Yerushalmi, Pfeffer and Shinnar1969). The lines are predictions from the NS dispersion relation (3.18).
3.3. Limiting case of film flows without inertia
In the regime where inertia of the liquid film is disregarded, i.e. 
$Re \ll 1$ (or 
$Oh \gg 1$), 
$l$ approximates 
$k$. This leads to the first column in (3.18) coinciding with the second column, and the third with the fourth, resulting in an indeterminate form. To address this issue, we employ the method proposed by Tomotika (Reference Tomotika1935), which involves expanding the Bessel functions in a Taylor series with respect to 
$l$. For instance, 
$\textrm {I}_1(l) = \textrm {I}_1(k) + \textrm {I}_1'(k) (l-k) + \mathrm {O}[(l-k)^2]$. By eliminating the zero-order terms and neglecting higher-order terms (greater than the second order), we arrive at a determinant form, expressed as
\begin{equation} \left| \renewcommand\arraystretch{1.3} \begin{matrix} k \mathrm{I}_1(k\alpha) & k\alpha \mathrm{I}_1'(k\alpha) & - k \mathrm{K}_1(k\alpha) & k\alpha \mathrm{K}_1'(k\alpha) \\ G_{21} & G_{22} & G_{23} & G_{24} \\ 2 k^3 \mathrm{I}_1(k ) & 2 k^2 \left[ \mathrm{I}_1(k) + k \mathrm{I}_1'(k) \right] & -2 k^3 \mathrm{K}_1(k) & 2 k^2 \left[ \mathrm{K}_1(k) + k \mathrm{K}_1'(k) \right] \\ G_{41} & G_{42} & G_{43} & G_{44} \\ \end{matrix} \right|=0. \end{equation}
The definitions of the functions 
$G_{ij}$ can be found in Appendix C. Note that 
$\omega$ appears only linearly in the fourth line of the determinant, the dispersion relation between 
$\omega$ and 
$k$ can be expressed explicitly. After replacing the differentiation of the Bessel functions by 
$\mathrm {I}'_1(k) = \mathrm {I}_0(k)-\mathrm {I}_1(k)/k$ and 
$\mathrm {K}'_1(k) = -\mathrm {K}_0(k)-\mathrm {K}_1(k)/k$, we have
\begin{align} \omega = \frac{k^2-1}{2} \frac{-\mathrm{I}_1(k) \varDelta_1+\mathrm{I}_0(k) \varDelta_2 - \mathrm{K}_1(k) \varDelta_3 + \mathrm{K}_0(k) \varDelta_4 } { \left[ k \mathrm{I}_0(k) -\mathrm{I}_1(k) \right] \varDelta_1 -k \mathrm{I}_1(k) \varDelta_2-\left[ k\mathrm{K}_0(k)+\mathrm{K}_1(k) \right] \varDelta_3 + k \mathrm{K}_1(k) \varDelta_4 } , \end{align}
where details of 
$\varDelta _{ij}$ are shown in Appendix C. This dispersion relation (3.27) is identical to the slip-modified Stokes model proposed by Zhao et al. (Reference Zhao, Zhang and Si2023a), which was derived directly from the Stokes equations (neglecting inertia in the NS equations). Numerical investigations are also performed to further support the theoretical analysis, illustrated in figure 3, where the predictions generated by (3.18) tend to converge towards the Stokes model as 
$Oh$ increases (inertia declines). Remarkably, the rate of convergence for the no-slip cases surpasses that of the slip cases significantly. Specifically, when 
$Oh=3.2 \times 10^{-2}$ (as indicated by the red dotted lines in figure 3a), predictions from the NS dispersion relation closely align with the results of the no-slip Stokes model. However, for the slip cases, 
$Oh \geq 1$ is required for a similar convergence. One plausible explanation for this observation, considering the omission of the base flow, is that the no-slip boundary conditions constrain the fluid motion within the liquid film more effectively than the slip boundary conditions, thereby mitigating the influence of inertia.
The dispersion relation between the growth rate 
$\omega$ and the wavenumber 
$k$ for the limiting cases of thin-film flows (
$\alpha = 0.8$). (a) No-slip cases with different inertial effects, 
$Oh=2 \times 10^{-3}$ (green dash-dotted line), 
$8 \times 10^{-3}$ (blue dashed line), 
$3.2 \times 10^{-2}$ (red dotted line), 
$0.128$ (black solid line). (b) Slip cases (
$l_s=10$) with different inertial effects, 
$Oh=0.01$ (green dash-dotted line), 
$0.1$ (blue dashed line), 
$1$ (red dotted line), 
$10$ (black solid line). The circles are the predictions from the slip-modified Stokes model (Zhao et al. Reference Zhao, Zhang and Si2023a) and the lines are predictions from the NS dispersion relation (3.18).
3.4. Predictions of the dispersion relation
Based on the insights gained from our analysis of limiting cases, we turn to the examination of inertia and slip effects in more general scenarios in this subsection.
In figure 4 we present dispersion relations from (3.18) for various slip lengths, while holding specific values of 
$Oh$ (columns) and 
$\alpha$ (rows). The inertial effects are shown along a given row (fixed 
$\alpha$), with the first column being inertia-dominated flows and the third column corresponding to viscosity-dominated cases. The deviations observed across different 
$l_s$ values indicate that slip predominantly governs the dynamics in viscous flows within the film but has a comparatively minor impact on inertia-dominated cases. This observation is consistent with the outcomes obtained from the limiting case analysis of jet flows (figure 2a). We also investigate the influence of fibre radii (film thickness) on the instability within each column. As film thickness (
$1-\alpha$) increases, the deviations diminish, suggesting that slip effects become less pronounced in thicker films.
The dispersion relation between the growth rate 
$\omega$ and the wavenumber 
$k$ on different boundary conditions of various fibre radii: (a–c) 
$\alpha =0.8$, (d–f) 
$\alpha =0.5$, (g–i) 
$\alpha =0.2$. For the inertial effects: (a,d,g) 
$Oh=10^{-3}$, (b,e,h) 
$Oh=0.1$, (c,f,i) 
$Oh=10$. Line types represent different values of the slip length: 
$l_s=0$ (green dash-dotted line), 
$0.1$ (blue dashed line), 
$1$ (red dotted line), 
$10$ (black solid line).
These two observations can be explained qualitatively by considering variations in velocity profiles within the liquid films, influenced by both slip and inertia. As the slip length increases, the flow field near the solid wall undergoes a transition from parabolic flow with a non-uniform velocity profile to plug flow with a uniform velocity profile (Münch et al. Reference Münch, Wagner and Witelski2005). The parabolic flow field decreases and constitutes only a small fraction of the film thickness with an increase in inertia (Schlichting & Kestin Reference Schlichting and Kestin1961), leading to more uniform velocity profiles. This explains why slip does not significantly impact the instability with 
$Oh=10^{-3}$ (figure 4a,d,g). However, when 
$Oh=10$, most of the flow fields within the films are expected to resemble parabolic profiles, making them more susceptible to the effects of slip. Moreover, as the film thickness increases, the proportion of the parabolic flow field in the films diminishes, and the velocity profiles become more uniform, resulting in weaker influences of slip on the instability.
The critical wavenumbers shown in figure 4 align with the findings of Plateau (Reference Plateau1873), i.e. 
$k_{crit}=2 {\rm \pi}h_0 / \lambda _{crit}=1$, indicating that slip conditions do not impact these values. These values are determined by the interplay between two curvature terms governing the Laplace pressure on the right-hand side of (2.9). The circumferential curvature term, 
$1 / (h \sqrt {1+(\partial _x h )^2} )$, acts as the driving force, while the tangential curvature term, 
$\partial _x^2 h / (1+(\partial _x h )^2 ) ^{3/2}$, acts as the resisting force. The balance of forces outlined in (2.9) yields the term 
$k^2-1$ in 
$F_{4j}$ of the final dispersion relation (3.18), ultimately determining the critical wavenumber as 
$k_{crit}=1$.
Figure 5 further elucidates the relationship between the dominant wavenumber 
$k_{max}$ and 
$Oh$. Remarkably, inertia appears to exert minimal influence on 
$k_{max}$ in no-slip cases. Specifically, for a thin film with 
$\alpha =0.8$ (figure 5a), 
$k_{max}$ for small 
$l_s$ remains unchanged as 
$Oh$ increases. This value is close to the analytical expression 
$k_{max} = \sqrt {1/2}$ derived from (B4) for no-slip cases. This finding offers a plausible explanation for why the no-slip lubrication model (2.13), which neglects inertia, has demonstrated remarkable capabilities in predicting the wavelengths observed in numerous experimental studies (Quéré Reference Quéré1990; Craster & Matar Reference Craster and Matar2006; Duprat et al. Reference Duprat, Ruyer-Quil, Kalliadasis and Giorgiutti-Dauphiné2007; Craster & Matar Reference Craster and Matar2009; Ji et al. Reference Ji, Falcon, Sadeghpour, Zeng, Ju and Bertozzi2019). Conversely, in slip cases, 
$k_{max}$ exhibits a decline with increasing 
$Oh$. This trend holds across different values of slip length (
$l_s$), with more pronounced decreases in 
$k_{max}$ for larger slip lengths. Encouragingly, the predictions of the giant-slip lubrication model (B3) closely align with the results of cases with substantial slip (
$l_s>10$). As 
$Oh$ becomes sufficiently large for all 
$l_s$, 
$k_{max}$ predicted by the NS dispersion relation converges to a constant, whereas in the giant-slip lubrication model, 
$k_{max}$ consistently decreases rapidly. For a thick film with 
$\alpha =0.2$ (figure 5b), though 
$k_{max}$ for the no-slip case decline from 
$0.7$ to 
$0.61$, which cannot be predicted by the no-slip lubrication model (B4), the variation trend of 
$k_{max}$ with 
$Oh$ is similar to that observed in thin films (
$\alpha =0.8$). Furthermore, slip is found not to significantly impact 
$k_{max}$ in film flows dominated by inertia (
$Oh<10^{-2}$), consistent with the findings in figure 4. However, in viscous cases (
$Oh>1$), 
$k_{max}$ decreases significantly as 
$l_s$ increases, corroborating the conclusions drawn in the work of Zhao et al. (Reference Zhao, Zhang and Si2023a).
Influence of inertia (different values of 
$Oh$) on the dominant wavenumber 
$k_{max}$ on fibres of two radii: (a) 
$\alpha =0.8$, (b) 
$\alpha =0.2$. The solid lines are the predictions of the NS dispersion relation (3.18) for different slip lengths: 
$l_s=0$ (green), 
$0.1$ (purple), 
$1.0$ (blue), 
$10$ (black), 
$100$ (red). The dotted lines and dashed lines represent the predictions from the no-slip lubrication model (2.13) and the giant-slip one (2.15).
4. Direct numerical simulations
In this section, direct numerical simulations are performed to corroborate the theoretical findings in § 3 and gain deeper physical understanding of how inertia and slip impact the instability of films on fibres.
The numerical solution of the NS equations is achieved through the finite element method within a computational framework facilitated by COMSOL Multiphysics 6.1. The simulations for the film flows are conducted using the arbitrary Lagrangian–Eulerian (ALE) approach. In this method, free-surface nodes are moved in a Lagrangian manner, deforming the computational domain, while nodes inside the film follow a predefined evolution. This approach offers a distinct advantage over other techniques, such as level set or phase field methods. Unlike these alternatives, the ALE approach ensures an exact capture of the free surface, akin to the Lagrangian approach, while retaining the primary advantage of Eulerian methods that mesh elements are less prone to distortion. Consequently, it has gained widespread use in predicting the dynamics of free-surface flows across various phenomena such as droplet dynamics (Chubynsky et al. Reference Chubynsky, Belousov, Lockerby and Sprittles2020; Chakraborty, Chubynsky & Sprittles Reference Chakraborty, Chubynsky and Sprittles2022), dynamics of a ligament (Wei et al. Reference Wei, Rivero-Rodríguez, Zou and Scheid2021), jet breakup (Martínez-Calvo et al. Reference Martínez-Calvo, Rivero-Rodríguez, Scheid and Sevilla2020b) and the instability of planar films (González et al. Reference González, Diez and Sellier2016; Moreno-Boza, Martínez-Calvo & Sevilla Reference Moreno-Boza, Martínez-Calvo and Sevilla2020b). However, it is worth noting that this method is not suitable for scenarios where the topology of the domain might change, such as in cases involving fluids inside the film after rupture. This limitation arises from the necessity to maintain consistent mesh connectivity throughout the simulation, and the only workaround is to manually change the mesh topology.
The computational domain is a quadrilateral (the section of a hollow fibre in cylindrical coordinates) with a size 
$[\alpha, h_0+\hat {h}] \times [0, L]$, illustrated in figure 6(a). Here 
$L$ is the length of the film/fibre. we assign a value of 
$\alpha =0.5$ to the radius of the fibre, and the initial radius of the film is 
$h_0 = 1$. Small perturbations (
$\hat {h}$) are introduced at the liquid–gas interface. The left and right boundary conditions of the computational domain are considered periodic. The top represents the free surface, following (2.5)–(2.7). The bottom is treated as the slip-wall boundary modelled by (2.11) and (2.12), where the slip length 
$l_s$ serves as an input parameter for this boundary condition. The axial velocity on the boundary 
$w_b = l_s \partial _r w_b$. When 
$l_s=0$, the no-slip boundary condition (
$w_b=0$) is recovered. The computational mesh for the liquid domain utilises non-uniform triangular Lagrange elements, shown in figure 6(b). Special attention is given to placing finer grid elements near the solid boundary to accurately capture the fluid behaviours in the non-uniform velocity profiles. The minimum grid size employed is 
$10^{-3}$. Additionally, the variable-order backward differentiation formula is utilised for the temporal integration of the NS equations. All simulations are conducted using dimensionless units, which are established through rescaling variables as described in (2.1a–d).
Numerical settings: (a) quadrilateral computational domains with different boundary conditions (BCs), (b) non-uniform triangular mesh.
We explore two different configurations in this study, each varying in film length and initial perturbations, allowing us to investigate the combined effects of inertia and slip on the dominant wavelength (§ 4.1) and the evolution of perturbation growth (§ 4.2)
4.1. Dominant wavelengths of perturbations
To investigate the influence of inertia and slip effects on the dominant wavelengths of perturbations, we perform simulations involving long films with a length of 
$L=200$ on various slippery fibres. We explore twenty different cases by considering five values of the Ohnesorge number, specifically 
$Oh=10^{-3}$, 
$10^{-2}$, 
$0.1$, 
$1$ and 
$10$, on fibres characterized by four slip lengths, namely 
$l_{s}=0$, 
$1$, 
$10$ and 
$100$. In these simulations, the system is initiated with random initial perturbations described as 
$h(z,0) = 1+ \varepsilon N(z)$, where 
$\varepsilon = 10^{-3}$ and 
$N(z)$ is a random variable following a normal distribution with a mean of zero and a unit variance. These initial perturbations are designed to replicate the arbitrary disturbances commonly encountered in reality.
Driven by surface tension, small random perturbations gradually evolve over time, giving rise to significant capillary waves, as illustrated in figure 7, which shows the evolution of interface profiles 
$h(z,t)$ for six different cases. To assess the impact of inertial effects, identical initial conditions are assigned for all six cases. For the inertia-dominated cases (
$Oh=10^{-3}$), the evolution of capillary waves in slip cases is nearly indistinguishable from that in no-slip cases, except for the slightly faster growth of perturbations in the slip case compared with the no-slip case (figure 7a). Conversely, when viscosity becomes significant (figure 7b,c), slip noticeably affects 
$h(z,t)$, resulting in faster perturbation growth and longer capillary waves. Furthermore, the wavelengths of the no-slip cases do not appear to be significantly influenced by inertia (see the six waves in figure 7a i,b i,c i), despite the presence of deviations in local interface profiles. All these findings align qualitatively with the theoretical predictions outlined in § 3.4. Additionally, we conduct a comparison between the interface profiles obtained through simulations for the NS equations and those calculated numerically from the lubrication models. The details are presented in Appendix D.
Interface profiles at four time instants, illustrated in different colours, on the fibres of the radius 
$\alpha =0.5$. The inertial effects are presented by different values of 
$Oh$: (a) 
$Oh=10^{-3}$, (b) 
$Oh=0.1$, (c) 
$Oh=10$. Panels (a i,b i,c i) are the predictions of no-slip cases (
$l_s=0$) and (a ii,b ii,c ii) are the results of slip cases (
$l_s=10$).
To quantitatively compare the numerical observations with the theoretical predictions derived from (3.18), we conduct multiple independent simulations (10 for each case) with different initial conditions to collect statistical data of the dominant modes. This statistical approach was proposed by Zhao, Sprittles & Lockerby (Reference Zhao, Sprittles and Lockerby2019) and has been employed in evaluating dominant modes of instability in various films (Zhao et al. Reference Zhao, Zhao, Si and Chen2021; Zhao, Zhang & Si Reference Zhao, Zhang and Si2023b; Zhao et al. Reference Zhao, Zhang and Si2023a). For each simulation, a discrete Fourier transform is applied to the interface position 
$h(z,t)$ to obtain the power spectral density (PSD) of the perturbations. The square root of the ensemble-averaged PSD (
$H_{rms}$) at each time step is depicted in figure 8, with a Gaussian function used to fit the modal distribution (spectrum). The peak of this fitted spectrum corresponds to the dominant wavenumber 
$k_{max}$, as indicated by the black dash-dotted lines. Extracting 
$k_{max}$ from the fitted spectrum at each time instant yields the insets in figure 8. Promisingly, 
$k_{max}$ converges to a constant rapidly in all the cases. Consistent with the findings in figure 7, the two spectra in figure 8(a,c) are nearly identical, suggesting that slip does not significantly impact the instability in the inertia-dominated regime. However, in figure 8(d) the spectrum for the slip case exhibits a smaller dominant wavenumber compared with the no-slip case in figure 8(c). This discrepancy indicates that, in the viscous regime, slip leads to an increase in the wavelength.
The root mean square (r.m.s.) of non-dimensional perturbation amplitude versus non-dimensional wavenumber on fibres with two slip lengths at four time instants: (a) 
$10\,583$ (green), 
$15\,875$ (blue), 
$19\,050$ (red) and 
$21\,166$ (black); (b) 
$264.6$ (green), 
$476.2$ (blue), 
$582.1$ (red) and 
$635.0$ (black); (c) 
$10\,583$ (green), 
$15\,875$ (blue), 
$18\,521$ (red) and 
$20\,632$ (black); (d) 
$26.5$ (green), 
$52.9$ (blue), 
$66.1$ (red) and 
$74.1$ (black). The circles are extracted from numerical simulations fitted by the Gaussian function. The inset shows the time history of the dominant wavenumbers extracted from the spectra.
This statistical analysis is applied for all the cases to generate the symbols (
$\lambda _{max} = 2 {\rm \pi}/k_{max}$) in figure 9, which exhibit a good agreement with theoretical predictions. Consequently, we can draw the conclusion that the dominant modes of the thin-film instability remain largely unaffected by inertia in no-slip cases. However, they become significantly influenced by inertia on slippery surfaces, leading to the formation of longer perturbation waves.
Inertial effects on the dominant wavelengths: a comparison between the theoretical predictions of (3.18) and numerical solutions with different slip lengths.
4.2. Evolutions of perturbation growth
To investigate the impact of inertia and slip on perturbation growth, we conduct simulations of a relatively short film with 
$L=10$ on slippery fibres. The film is initially perturbed with 
$h(z,t) =1+ \varepsilon \cos [2 {\rm \pi}(z/L-1/2)]$, where 
$\varepsilon =0.01$.
Figure 10 presents the time evolution of the film interface on a no-slip fibre with a radius of 
$\alpha =0.5$, with the corresponding fluid structure. The contour in the lower panel of figure 10 depicts the radial velocity 
$u$. The upper half displays the axial velocity 
$w$, showing that opposite fluxes occur, directed towards the left and right boundaries, as the perturbation at the free surface grows due to instability. Notably, while the magnitudes of both 
$u$ and 
$w$ increase during the process, the fluid structure remains similar at the linear stage (the amplitude of disturbances is typically less than 20 % of initial radius, i.e. 
$\hat {h}<0.2$), as evident from the contour distribution in figure 10.
Film thinning of one perturbation wave in an axially symmetric domain. For this case, 
$L=10$, 
$\alpha =0.5$, 
$Oh=0.1$ and 
$l_s=0$. Contours of the axial velocity 
$w(r,z,t)$ (a i,b i,c i) and radial velocity 
$u(r,z,t)$ (a ii,b ii,c ii) inside the film are shown at three time instants: (a) 
$t_1=158$, (b) 
$t_2=264$, (c) 
$t_3=370$. Here 
$h_{min}$ represents the minimum radius of the film.
In figure 11 we present more velocity fields for four distinct cases, considering two different values of 
$Oh$: 
$10^{-3}$ for the inertia-dominated cases and 
$0.1$ for the viscous cases, on both no-slip (left panel) and slip (right panel) fibres. According to the variations of the contours in the upper half of figure 11(b), slip not only alters the velocity distribution near the surface of the fibres but also accelerates instability in the viscous cases, as evidenced by the larger axial velocity component 
$w$. However, in the case of inertia-dominated flow, slip has a negligible impact on instability, corroborating our previous findings. The lower half of figure 11 illustrates the velocity vectors near 
$z=\pm 0.9$ for different cases. It is apparent that slip reduces the velocity gradient 
$\partial _r w$ near the surface. Furthermore, the parabolic flow field in the no-slip case with 
$Oh=10^{-3}$ is observed to be considerably smaller than that in the case with 
$Oh=0.1$. This finding lends additional support to the explanation provided following figure 4. In essence, the more uniform velocity profile in the inertia-dominated case serves to limit the impact of slip, resulting in nearly identical dynamics of the instability.
Inertia and slip effects on axial velocity fields: (a) 
$Oh=10^{-3}$ at 
$t=2027$, (b) 
$Oh=0.1$ at 
$t=106$. Here, the fibre radius 
$\alpha = 0.5$. Panels (a i,a iii) and (b i,b iii) illustrate the contours of the entire configuration (a i,b i) and velocity vectors of the local field near 
$|z |=0.9$ (a iii,b iii) of the no-slip cases. Panels (a ii,a iv) and (b ii,b iv) show the results of the slip case (
$l_s=10$).
Figure 12 illustrates the growth rates of perturbations in the twenty simulated cases. Based on the instability analysis (§ 3.1) that employs the expression 
$h(z,t) = 1+\hat {h}\textrm {e}^{\textrm {i}kz+\omega t}$, the initial perturbation, modelled by a cosine function with a fixed wavenumber 
$k=2 {\rm \pi}/L$, experiences exponential growth. So 
$\ln (1-h_{min})$ is utilised as the 
$y$ coordinate in figure 12 to present the linear growth of perturbations. The numerical results closely align with the theoretical predictions, demonstrating the effectiveness of the NS dispersion relation (3.18) in describing the inertial effects on the instability of films on slippery fibres. Furthermore, it becomes evident that deviations due to slip become more pronounced as inertial effects diminish, thus providing quantitative confirmation of previous findings.
Linear time evolution of the minimum radii of films 
$h_{min}(t)$ with four 
$Oh$ values: (a) 
$Oh=10^{-3}$, (b) 
$Oh=10^{-2}$, (c) 
$Oh=0.1$, (d) 
$Oh=10$. The numerical solutions (dotted lines with symbols) are compared with the predictions of the NS dispersion relation (solid lines) for four lengths: 
$l_{s} = 0$ (green), 
$0.1$ (blue), 
$1$ (red) and 
$10$ (black). The theoretical values of 
$\omega$ are (a) 
$2.81 \times 10^{-4}$ (green), 
$2.91 \times 10^{-4}$ (blue), 
$2.94 \times 10^{-4}$ (red), 
$2.95 \times 10^{-4}$ (black); (b) 
$2.41 \times 10^{-3}$ (green), 
$2.61 \times 10^{-3}$ (blue), 
$2.84 \times 10^{-3}$ (red), 
$2.89 \times 10^{-3}$ (black); (c) 
$8.12 \times 10^{-3}$ (green), 
$1.10 \times 10^{-2}$ (blue), 
$1.99 \times 10^{-2}$ (red), 
$2.34 \times 10^{-2}$ (black); (d) 
$8.90 \times 10^{-3}$ (green), 
$1.35 \times 10^{-2}$ (blue), 
$3.62 \times 10^{-2}$ (red), 
$6.32 \times 10^{-2}$ (black).
In addition to examining the evolution of perturbation growth at the linear stage, we also delve into the nonlinear dynamics, illustrated in figure 13. Due to the dramatic changes in the interface profile at the nonlinear stage, the initial dense mesh experiences significant deformation, resulting in poor grid quality and potential numerical errors. To address this challenge, we implement an approach of remeshing, regenerating the mesh (reducing nodes) in the vicinity near the point of 
$h_{min}$ when 
$h_{min}<0.55$. The simulations encompass both inertia-dominated and viscosity-dominated cases on the no-slip (left panel) and slip (right panel) boundary conditions. Notably, the nonlinear evolution reveals substantial distinctions from the dynamics at the linear stage. The interface shapes are found to deviate from their initial cosine, forming plateau structures at their lowest points. Ultimately, satellite droplets emerge between the two main drops. Additionally, fluid structure within the film no longer exhibits ‘similarity’ at different time instants owing to the drastic changes in interface profiles. In scenarios dominated by inertia, although slip has a minimal effect on the interface shape, it substantially alters the flow structure. Figure 13(c,d) illustrates the evolution of vortical structures within the liquid film on a no-slip wall, resulting in significant oscillations of the interface before rupture. Conversely, no vortices appear within the film on the slippery fibre due to the weak shear forces acting on the fluid near the surface. In viscosity-dominated scenarios, slip not only accelerates perturbation growth significantly, as observed at the linear stage, but also affects the interface profiles near rupture. This alteration results in the formation of filaments, rather than satellite droplets, between the two main drops. One plausible explanation is that the dominant perturbation wavelength of the instability in the slip case (
$\lambda _{max}=15.32$) is notably larger than that in the no-slip case (
$\lambda _{max}=9.38$), resulting in the flatter profile in figure 13(h). The observation also suggests that the volume of satellite droplets is influenced by both inertia and slip.
Evolutions of perturbation growth at the nonlinear stage with two 
$Oh$ values: (a–d) 
$Oh=10^{-3}$, (e–h) 
$Oh=10$. The contours represent the velocity magnitude 
$| U| = \sqrt {u^2+w^2}$ . Panels (a i,b i,c i,d i,e i,f i,g i,h i) illustrate the contours and streamlines of no-slip cases. Panels (a ii,b ii,c ii,d ii,e ii,f ii,g ii,h ii) show the results of the slip cases (
$l_s=10$).
Figure 14 presents variations of satellite droplets, extracted from more than seventy simulations with varying values of 
$Oh$ and 
$l_s$. In figure 14(a) we depict the interface profiles of satellite droplets. These profiles clearly demonstrate that in viscosity-dominated scenarios, slip significantly reduces the volumes of satellite droplets. However, in the case of inertia-dominated scenarios where 
$l_s>1$, slip has no discernible effect on the droplet volumes. Note that these profiles were obtained when 
$h_{min} \leq 0.01\,h_0$. Furthermore, the volume of the satellite droplets is quantified by calculating the area between the two lowest points of the profiles, as shown in figure 14(b). It is evident that in viscous cases the volume of satellite droplets decreases as 
$l_s$ increases, and higher viscosity (larger 
$Oh$) leads to a more rapid rate of decrease. When viscosity dominates the fluids (
$Oh>=10$), the relationship between volume (
$V_{sat}$) and slip length (
$l_s$) approximately follows 
$V_{sat} \sim l_s^{-5}$.
(a) The interface profiles of satellite droplets on different boundary conditions with: 
$l_s=0$ (black dash-dotted lines), 
$1$ (green solid lines), 
$10$ (blue dashed lines) and 
$100$ (red dotted lines). Panels (a i) and (a ii) show the results with 
$Oh=10^{-3}$ and 
$Oh=10$, respectively. (b) Variations of the satellite droplets with slip lengths. The inset provides a schematic of the volume of satellite droplets.
5. Conclusions
In this study a theoretical model is developed based on the linear instability analysis of the axisymmetric NS equations to investigate the influence of inertia and slip on dynamics of liquid films on fibres. The model is verified theoretically by the limiting cases for jet flows and film flows without inertia. The resulting dispersion relation (3.18) reveals some intriguing insights. Firstly, it shows that slip has a relatively minor impact on the instability of flows dominated by inertia, whereas it significantly accelerates the growth of perturbations in viscous film flows. Moreover, the influence of slip appears to be contingent on the thickness of the liquid film, with thinner films exhibiting more pronounced slip effects. We also extract the dominant perturbation modes, denoted as 
$k_{max}$, from the dispersion relation. In cases with no-slip boundary conditions, 
$k_{max}$ remains largely unaffected by inertia. Remarkably, for thin films characterized by 
$\alpha =0.8$, 
$k_{max}$ maintains a nearly constant value, closely aligning with the predictions of the no-slip lubrication model (2.13), even as 
$Oh$ varies. Conversely, when slip is introduced, 
$k_{max}$ exhibits a noticeable decrease as inertia decreases. In the limiting case where 
$l_s \rightarrow \infty$, the giant-slip lubrication model (2.15) offers an approximate prediction for 
$k_{max}$.
To substantiate our theoretical findings, direct numerical simulations of the NS equations are conducted via two distinct fluid configurations: (i) a long film with random initial perturbations to investigate the dominant modes of perturbations, and (ii) a short film with a fixed wavelength to examine the evolution of perturbation growth. The velocity fields extracted from these simulations yield valuable insights. Notably, the parabolic field of the velocity profiles in inertia-dominated cases are observed to be significantly smaller than those in cases governed by viscosity. Given that slip primarily influences non-uniform velocity profiles, this observation provides an explanation for why the impact of slip is dampened by inertia. Furthermore, we delve into the realm of nonlinear dynamics in perturbation growth. We found that nonlinear dynamics, in contrast to the linear stage, leads to the generation of intricate vortical structures near no-slip surfaces in inertia-dominated cases. Interestingly, slip was found to mitigate the impact of these shear stresses, resulting in smoother flows without the presence of vortices. In the context of viscous cases, slip was identified to reduce the volume of satellite droplets (
$V_{sat}$) at the final stages before rupture. When viscosity dominates the fluid dynamics, this reduction in volume followed an approximate power-law relationship, specifically 
$V_{sat} \sim l_s^{-5}$.
Given that experimental evidence has already confirmed the existence of wall slip in films on fibres, as demonstrated in previous studies (Haefner et al. Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015; Ji et al. Reference Ji, Falcon, Sadeghpour, Zeng, Ju and Bertozzi2019), and that slip length can be directly quantified through methods such as those described in Huang, Guasto & Breuer (Reference Huang, Guasto and Breuer2006), Maali & Bhushan (Reference Maali and Bhushan2012) and Maali, Colin & Bhushan (Reference Maali, Colin and Bhushan2016), it is our hope that the wall slip can be controlled experimentally to validate our predictions, especially regarding the impact of inertia on wavelengths (drop size) for various liquids. For instance, experiments involving water can be conducted to examine inertia-dominated scenarios, while experiments with silicone oil can be carried out to investigate viscosity-dominated cases. There are numerous potential extensions to this framework. One avenue is to incorporate the influence of other physical factors, such as electric fields (Ding et al. Reference Ding, Xie, Wong and Liu2014) and intermolecular forces (Ji et al. Reference Ji, Falcon, Sadeghpour, Zeng, Ju and Bertozzi2019; Tomo et al. Reference Tomo, Nag and Takamatsu2022), both of which are known to affect critical wavenumbers (
$k_{crit}$) and dominant modes (
$k_{max}$) in the instability. Additionally, exploring related flow configurations, such as a liquid film flowing down a fibre driven by a body force, as discussed in Liu & Ding (Reference Liu and Ding2021), presents intriguing opportunities. An open problem in this context is to comprehend the impact of inertial effects on the dynamics of travelling waves in various flow regimes.
Acknowledgements
Useful discussions with Dr Y. Zhang, Dr Z. Ding and Mr J. Xiong are gratefully acknowledged.
Funding
This work was supported by the National Natural Science Foundation of China (grant nos. 12202437, 12027801, 12388101, 12272372), the Youth Innovation Promotion Association CAS (grant nos. 2018491, 2023477), the Fundamental Research Funds for the Central Universities, the China Postdoctoral Science Foundation (grant no. 2022M723044) and the China International Postdoctoral Exchange Program Fellowship (grant no. YJ20210177).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation for the giant-slip lubrication model
In this appendix we follow the approach employed by Münch et al. (Reference Münch, Wagner and Witelski2005) to derive the giant-slip lubrication model. To get this lubrication equation from the axisymmetric NS equations, we need to establish the leading-order terms by their asymptotic expansion in 
$\varepsilon$, for which we use the rescaling shown below:
\begin{equation} \tilde{x}= \tilde{\lambda} x,\quad \tilde{r}=\varepsilon \tilde{\lambda} r,\quad \tilde{w} = \tilde{U}w,\quad \tilde{u}=\varepsilon \tilde{U} u,\quad \tilde{t} = \frac{\tilde{\lambda}}{\tilde{U}} t,\quad \tilde{p} = \frac{\mu U }{ \lambda } p. \end{equation}
Here, 
$\tilde {\lambda } = h_0 / \varepsilon$. Substituting all these scalings into the dimensional NS equations yields
\begin{gather} \frac{\partial w}{\partial z} + \frac{1}{r} \frac{\partial(u r) }{\partial r} = 0, \end{gather}
\begin{gather} \varepsilon^2 Re\left( \frac{\partial u}{\partial t}+ w \frac{\partial u}{\partial z} + u \frac{\partial u}{\partial r} \right) =- \frac{\partial p }{\partial r}+ \varepsilon^2 \frac{\partial^2 u}{\partial z^2} +\frac{\partial}{\partial r} \left[\frac{1}{r} \frac{\partial (ur)}{\partial r}\right], \end{gather}
\begin{gather} \varepsilon^2 Re\left( \frac{\partial w}{\partial t}+ w \frac{\partial w}{\partial z} + u \frac{\partial w}{\partial r} \right)=-\varepsilon^2 \frac{\partial p }{\partial z}+ \varepsilon^2 \frac{\partial^2 w}{\partial z^2} + \frac{1}{ r} \frac{\partial}{\partial r}\left(r \frac{\partial w}{\partial r} \right) . \end{gather}
The equations at the liquid–gas interface (
$r=h$) are scaled as
$$\begin{gather} \frac{\partial h}{\partial t} + w \frac{\partial h}{\partial z} -u= 0 ,\nonumber\\ p - \frac{2}{1+ \varepsilon^2 \left(\partial_z h \right)^2} \left[ \frac{\partial u}{\partial r} - \frac{\partial h}{\partial z} \left( \frac{\partial w}{\partial r} + \varepsilon^2 \frac{\partial u}{\partial z}\right) + \varepsilon^2 \left(\frac{\partial h}{\partial z} \right)^2 \frac{\partial w}{\partial z} \right] \end{gather}$$
$$\begin{gather} =\frac{Re}{We} \left[ \frac{1}{\varepsilon^2 h \sqrt{1+ \varepsilon^2 \left(\partial_z h \right)^2}} - \frac{\partial_z^2 h}{\left(1+\varepsilon^2 \left(\partial_z h \right)^2\right)^{{3}/{2}}} \right] , \end{gather}$$
$$\begin{gather} 2 \varepsilon^2 \frac{\partial h}{\partial z} \left( \frac{\partial u}{\partial r} - \frac{\partial w}{\partial z} \right) + \left[ 1- \varepsilon^2 \left(\frac{\partial h}{\partial z} \right)^2 \right] \left( \frac{\partial w}{\partial r} + \varepsilon^2 \frac{\partial u}{\partial z} \right) = 0 . \end{gather}$$
For the boundary conditions at the liquid–solid interface (
$r=\alpha$), the scaled forms are
$$\begin{gather} w = \frac{ l_{s}}{ \varepsilon^2 } \frac{\partial w}{\partial r}, \end{gather}$$
$$\begin{gather}u=0 . \end{gather}$$The rescaled equations can be approximately solved by the perturbation expansion, expressed as
\begin{equation} \left( w,u,p,h\right) = \left( w_0,u_0,p_0,h_0\right) + \varepsilon^2 \left( w_1, u_1, p_1, h_1\right)+ \cdots. \end{equation}
After eliminating all the high-order terms of 
$\varepsilon$ and only keeping the terms of the leading order, we obtain
\begin{gather} \partial_z w_0+ \partial_r(u_0 r)/r = 0 , \end{gather}
\begin{gather} \partial_r p_0 = \partial_r \left[ \partial_r \left(u_0 r \right) / r\right] , \end{gather}
\begin{gather} \partial_r \left(r \partial_r w_0 \right) =0, \end{gather}
\begin{gather} \partial_{t} h_0 + w_0 \partial_z h_0 - u_0 = 0 , \end{gather}
\begin{gather} p_0 -2\left( \partial_r u_0 - \partial_z h_0 \partial_r w_0 \right) = \frac{Re}{ We } \left( 1/h_0 - \partial_z^2 h_0 \right), \quad (r=h) , \end{gather}
\begin{gather} \partial_r w_0 = 0, \quad (r=h) , \end{gather}
\begin{gather} \partial_r w_0 = 0, \quad (r=\alpha) , \end{gather}
\begin{gather} u_0 = 0, \quad (r=\alpha) . \end{gather}
According to (A13), (A17) and (A18), 
$w$ is independent of the 
$r$, i.e. 
$w_0 = w_0(z,t)$. So (A11) is rearranged as
\begin{equation} u_0 ={-} r \partial_z w_0/2 . \end{equation}Substituting (A19) into (A14) and (A15) yields
$$\begin{gather} \partial_{t} h_0 + w_0 \partial_z h_0 + h_0 \partial_z w_0/2= 0, \end{gather}$$
$$\begin{gather}p_0+\partial_z w_0 = \frac{Re}{ We } \left( 1/h_0 - \partial_z^2 h_0 \right) . \end{gather}$$For the terms of the next order, (A4) becomes
\begin{equation} Re\left( \frac{\partial w_0}{\partial t}+ w_0 \frac{\partial w_0}{\partial z} \right)={-} \frac{\partial p_0 }{\partial z}+ \frac{\partial^2 w_0}{\partial z^2} + \frac{1}{ r} \frac{\partial}{\partial r}\left(r \frac{\partial w_1}{\partial r} \right) , \end{equation}and (A7) and (A8) are applied for the boundary conditions
$$\begin{gather} -3 \frac{\partial h_0}{\partial z} \frac{\partial w_0}{\partial z} + \frac{\partial w_1}{\partial r}-\frac{r}{2}\frac{\partial^2 w_0}{\partial z^2} = 0, \quad (r=h) , \end{gather}$$
$$\begin{gather}w_0 - l_{s} \frac{\partial w_1}{\partial r} = 0, \quad (r=\alpha) . \end{gather}$$
Integrating (A22) from 
$\alpha$ to 
$h$ and using boundary conditions (A23) and (A24) leads to
\begin{equation} Re\left( \frac{\partial w_0}{\partial t}+ w_0 \frac{\partial w_0}{\partial z} \right) \int_\alpha^h r\,{\rm d}r= \left.\left( -\frac{\partial p_0 }{\partial z}+ \frac{\partial^2 w_0}{\partial z^2} \right) \int_\alpha^h r\,{\rm d}r +r \frac{\partial w_1}{\partial r} \right|_{r=\alpha}^{r=h} . \end{equation}Rearranging (A25) yields
\begin{align} Re\left( \frac{\partial w_0}{\partial t}+ w_0 \frac{\partial w_0}{\partial z} \right)={-}\frac{\partial p_0 }{\partial z}+ \frac{\partial^2 w_0}{\partial z^2} + \frac{h_0}{h_0^2-\alpha^2}\left( 6 \frac{\partial h_0}{\partial z} \frac{\partial w_0}{\partial z} + h_0 \frac{\partial^2 w_0}{\partial z^2}\right) - \frac{2 \alpha}{h_0^2-\alpha^2} \frac{w_0}{l_{s}}. \end{align}Combining (A20), (A21) and (A26) gives the final lubrication model for the film on a fibre with a giant-slip length, written as
\begin{cases}
\frac{\partial h^2}{\partial t} + \frac{\partial \left(h^2 w \right)}{\partial z} &= 0, \\
\frac{\partial w}{\partial t}+ w \frac{\partial w}{\partial z} &= -\frac{1}{ We}
\frac{\partial}{\partial z}\left(\frac{1}{h}-\frac{\partial^2 h}{\partial z^2} \right)+\frac{1}{
Re} \frac{3}{h^2-\alpha^2} \frac{\partial \left( h^2 \partial_z w\right)}{\partial z} \nonumber\\
&\quad - \frac{1}{ Re} \frac{2 \alpha^2}{h^2-\alpha^2} \left(\frac{\partial^2 w}{\partial
z^2}-\frac{w}{ \alpha\, l_{s}} \right).
\end{cases}
If 
$U = \gamma / \mu$, 
$Re = We = Oh^{-2}$, the momentum equation (A27b) becomes
\begin{equation} \frac{\partial w}{\partial t}+ w \frac{\partial w}{\partial z} ={-} Oh^2 \frac{\partial}{\partial z}\left(\frac{1}{h}-\frac{\partial^2 h}{\partial z^2} \right) + Oh^2 \frac{3}{h^2-\alpha^2} \frac{\partial \left( h^2 \partial_z w\right)}{\partial z} - Oh^2 \frac{2 \alpha^2}{h^2-\alpha^2} \left(\frac{\partial^2 w}{\partial z^2}-\frac{w}{ \alpha l_{s}} \right). \end{equation}It is worth noting that 
$h^2-\alpha ^2$ introduces a singularity as 
$h \rightarrow \alpha$. A similar singularity is also observed in other lubrication equations describing free-surface flows with significant inertial effects, such as liquid jets (Eggers & Dupont Reference Eggers and Dupont1994), liquid sheets (Erneux & Davis Reference Erneux and Davis1993) and planar films on ultra-slip walls (Münch et al. Reference Münch, Wagner and Witelski2005). However, this singularity disappears in lubrication models that describe inertialess flows of bounded thin films (Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997; Kang, Nadim & Chugunova Reference Kang, Nadim and Chugunova2017). The absence is primarily attributed to the distinct scaling used in the lubrication approximation.
Appendix B. Linear instability analysis for lubrication equations
In this appendix the instability analysis is performed for the lubrication equations (2.15) and (2.13) using the normal mode method.
For the giant-slip model (2.15), substituting 
$h(z,t) = 1 + \hat {h} \textrm {e}^{\omega t + \textrm {i} k z}$ and 
$w(z,t) = \hat {w} \textrm {e}^{\omega t + \textrm {i} k z}$ into the linearised lubrication equation gives
$$\begin{gather} \omega \hat{h} ={-} {\rm i}k \hat{w}/2, \end{gather}$$
$$\begin{gather}Oh^{{-}2} \omega \hat{w}= {\rm i}k(1-k^2) \hat{h} - \frac{3 k^2 \hat{w}}{1-\alpha^2}+ \frac{2 \alpha^2}{1-\alpha^2} \left(k^2+ \frac{1}{\alpha\, l_{s}} \right) \hat{w}. \end{gather}$$
Eliminating the 
$\hat {h}$ in (B2) yields the dispersion relation
\begin{equation} \omega^2 + Oh^2 \left[ \frac{3-2\alpha^2}{1-\alpha^2}k^2 - \frac{2 \alpha}{(1-\alpha^2) l_{s}}\right] \omega + Oh^2 \frac{k^2 (k^2-1)}{2}=0. \end{equation}With a similar approach, we can have the dispersion relation from the no-slip lubrication model (2.13), namely
\begin{equation} \omega = (k^2-1) k^2 (3+ \alpha^4 - 4 \alpha^2 + 4\ln\alpha )/16. \end{equation}Appendix C. Auxiliary functions for the dispersion relations
In this appendix we present the definitions of the auxiliary functions for the dispersion relations (3.26) and (3.27).
The 
$G_{ij}$ in (3.26) are
\begin{equation} \left. \begin{aligned} G_{21} & = k^2 \mathrm{I}_0(k \alpha)-l_{s} k^3 \mathrm{I}_1(k\alpha),\\ G_{22} & = k \mathrm{I}_0(k\alpha) + k^2 \alpha \mathrm{I}_0'(k\alpha)-2 l_s k^2 \mathrm{I}_1(k\alpha) - l_s k^3 \alpha \mathrm{I}_1'(k\alpha),\\ G_{23} & = k^2 \mathrm{K}_0(k \alpha)+l_{s} k^3 \mathrm{K}_1(k \alpha) ,\\ G_{24} & ={-} k \mathrm{K}_0(k\alpha) - k^2 \alpha \mathrm{K}_0'(k\alpha)-2 l_s k^2 \mathrm{K}_1(k\alpha) - l_s k^3 \alpha \mathrm{K}_1'(k\alpha), \\ G_{41} & = 2k^2 \mathrm{I}'_1(k) \omega + (k^2-1) k \mathrm{I}_1(k), \\ G_{42} & = 2\left[ k^2 \mathrm{I}_0'(k)- k \mathrm{I}_1'(k) \right] \omega+ (k^2-1)k \mathrm{I}_1'(k),\\ G_{43} & ={-} 2k^2 \mathrm{K}'_1(k) \omega - (k^2-1) k \mathrm{K}_1(k) ,\\ G_{44} & ={-} 2\left[ k^2 \mathrm{K}_0'(k)+ k \mathrm{K}_1'(k) \right] \omega+ (k^2-1)k \mathrm{K}_1'(k). \end{aligned} \right\} \end{equation}For (3.27), we have
$$\begin{gather} \varDelta_1 = \left|\begin{array}{ccc} \alpha \mathrm{I}_0(k \alpha) & \mathrm{K}_1(k\alpha) & \alpha \mathrm{K}_0(k \alpha) \\ H_{22} & H_{23} & H_{24}\\ k \mathrm{I}_0(k) +\mathrm{I}_1(k) & k \mathrm{K}_1(k) & k \mathrm{K}_0(k)-\mathrm{K}_1(k) \\ \end{array}\right|, \end{gather}$$
$$\begin{gather}\varDelta_2 = \left|\begin{array}{ccc} \mathrm{I}_1(k\alpha) & \mathrm{K}_1(k\alpha) & \alpha \mathrm{K}_0(k \alpha) \\ H_{21} & H_{23} & H_{24}\\ k \mathrm{I}_1(k ) & k \mathrm{K}_1(k) & k \mathrm{K}_0(k)-\mathrm{K}_1(k) \\ \end{array}\right|, \end{gather}$$
$$\begin{gather}\varDelta_3 = \left|\begin{array}{ccc} \mathrm{I}_1(k\alpha) & \alpha \mathrm{I}_0(k \alpha) & \alpha \mathrm{K}_0(k \alpha) \\ H_{21} & H_{22} & H_{24}\\ k \mathrm{I}_1(k ) & k \mathrm{I}_0(k) +\mathrm{I}_1(k) & k \mathrm{K}_0(k)-\mathrm{K}_1(k) \\ \end{array}\right|, \end{gather}$$
$$\begin{gather}\varDelta_4 = \left|\begin{array}{ccc} \mathrm{I}_1(k\alpha) & \alpha \mathrm{I}_0(k \alpha) & \mathrm{K}_1(k\alpha) \\ H_{21} & H_{22} & H_{23}\\ k \mathrm{I}_1(k ) & k \mathrm{I}_0(k) +\mathrm{I}_1(k) & k \mathrm{K}_1(k) \\ \end{array}\right|, \end{gather}$$where
\begin{equation} \left. \begin{aligned} H_{21} & = k \mathrm{I}_0(k \alpha)-l_{s} k^2 \mathrm{I}_1(k\alpha) ,\\ H_{22} & = (2- l_s k^2 \alpha ) \mathrm{I}_0(k\alpha) +(\alpha-2 l_s) k \mathrm{I}_1(k\alpha),\\ H_{23} & ={-}k \mathrm{K}_0(k \alpha)-l_{s} k^2 \mathrm{K}_1(k \alpha) ,\\ H_{24} & =(2-l_s k^2 \alpha) \mathrm{K}_0(k\alpha) + (2 l_s-\alpha)k \mathrm{K}_1(k\alpha). \end{aligned} \right\} \end{equation}Appendix D. Comparisons of the interface profiles
In this appendix we compare the interface profiles obtained through simulations for the NS equations with those calculated numerically from the lubrication equations. Here, we focus on inertialess flows (
$Oh \gg 1$) on no-slip fibres. So we perform numerical investigations for the no-slip lubrication equation (2.13) introduced in § 2. Note that Lister et al. (Reference Lister, Rallison, King, Cummings and Jensen2006) used a different lubrication equation (D1) to explore nonlinear dynamics of film interfaces, particularly those exhibiting characteristic collar-and-lobe structures, written as
\begin{equation} \frac{\partial h}{\partial t} ={-}\frac{1}{ 3} \frac{\partial}{\partial z} \left[ (h-\alpha)^3 \left ( \frac{\partial h}{\partial z} + \frac{\partial^3 h}{\partial z^3} \right)\right]. \end{equation}
Here, the assumption 
$h-\alpha \ll \alpha$ simplifies the form of (D1) compared with (2.13). Both of these lubrication equations, with periodic boundary conditions, are solved using a simple second-order finite-difference scheme in both time and space. In the direct numerical simulations of the NS equations, we set 
$Oh=10$ and 
$l_s=0$ to enable a straightforward comparison with the numerical solutions of the lubrication equations. Figure 15(a) illustrates excellent agreement between the numerical predictions of the NS and lubrication equations, starting from the same initial interface profile, for a thin film (
$\alpha =0.8$). However, deviations become more pronounced for a thicker film (
$\alpha =0.5$), consistent with previous theoretical results from instability analysis (Zhao et al. Reference Zhao, Zhang and Si2023a).
Interface profiles of inertialess flows on two fibres of different radii at three time instants: (a) 
$\alpha =0.8$, 
$t_1=0$, 
$t_2=5700$, 
$t_3=8400$; (b) 
$\alpha =0.5$, 
$t_1=0$, 
$t_2=400$, 
$t_3=600$. The blue dotted lines and red dashed lines represent the results predicted by (D1) and (2.13), respectively. The solid lines are obtained from direct numerical simulations for the NS equations.
