1 Introduction
The break-up of free-surface flows has been investigated for a long time. The first quantitative studies of the instability responsible for the spontaneous break-up of cylindrical liquid threads date back to the 19th century; the correct physical description of the instability mechanism was due to Plateau (Reference Plateau1873), who deduced that a small perturbation with a wavelength larger than the circumference of the unperturbed column is unstable, finally breaking up into main drops and smaller satellite droplets in between. A few years later, Rayleigh (Reference Rayleigh1878), Lord Rayleigh (Reference Lord Rayleigh Sec.1892) calculated the most unstable wavelength using a linear temporal stability analysis. The subject experienced a renaissance 50 years ago that has lasted to the present due to its central role in industrial and medical applications such as chemical reactors, ink-jet and three-dimensional printing, additive manufacturing, drug and protein encapsulation, and cytometry, to cite a few; the reader is referred to the reviews of Bogy (Reference Bogy1979), Eggers (Reference Eggers1997), Christopher & Anna (Reference Christopher and Anna2007), Eggers & Villermaux (Reference Eggers and Villermaux2008), Derby (Reference Derby2010) and Anna (Reference Anna2016).
The theoretical approach to the study of the dynamics of jet break-up was first based on the linear stability analysis of infinite liquid threads. As already mentioned, the local temporal approach was pioneered by Lord Rayleigh (Reference Lord Rayleigh Sec.1892). About 80 years later, the local spatial and spatiotemporal problems, in which the liquid jet moves with uniform velocity 
$U$ with respect to the injector, were solved (Keller, Rubinow & Tu Reference Keller, Rubinow and Tu1973; Leib & Goldstein Reference Leib and Goldstein1986a,Reference Leib and Goldsteinb). In particular, it was demonstrated by Keller et al. (Reference Keller, Rubinow and Tu1973) that the spatial and temporal stability analyses are equivalent if 
$U$ is sufficiently larger than the speed of small-amplitude capillary instability waves, 
$U_{\unicode[STIX]{x1D70E}}$. In the spatial setting, the latter condition means that the relative growth of the wave amplitude along one wavelength is small. Thus, in a frame of reference moving with the jet, the amplitude growth is spatially uniform to a first approximation, which explains the equivalence of the temporal and spatial approaches if 
$U\gg U_{\unicode[STIX]{x1D70E}}$. Since the wavelength of the unstable capillary waves is much larger than the unperturbed cylinder radius, 
$\bar{R}$, the scaling of 
$U_{\unicode[STIX]{x1D70E}}$ depends on the value of the associated Reynolds number, 
$Re_{\unicode[STIX]{x1D70E}}=\unicode[STIX]{x1D70C}U_{\unicode[STIX]{x1D70E}}\bar{R}/\unicode[STIX]{x1D707}$, where 
$\unicode[STIX]{x1D70C}$ and 
$\unicode[STIX]{x1D707}$ are the liquid density and viscosity, respectively. In the limit of Euler flow, 
$Re_{\unicode[STIX]{x1D70E}}\gg 1$, the value of 
$U_{\unicode[STIX]{x1D70E}}$ is given by the balance 
$\unicode[STIX]{x1D70E}_{0}/\bar{R}\sim \unicode[STIX]{x1D70C}U_{\unicode[STIX]{x1D70E}}^{2}$, where 
$\unicode[STIX]{x1D70E}_{0}$ is the surface tension, yielding 
$U_{\unicode[STIX]{x1D70E}}\sim \sqrt{\unicode[STIX]{x1D70E}_{0}/(\unicode[STIX]{x1D70C}\bar{R})}$, usually referred to as the capillary velocity, and 
$U/U_{\unicode[STIX]{x1D70E}}\sim \sqrt{We}$, where 
$We=\unicode[STIX]{x1D70C}U^{2}\bar{R}/\unicode[STIX]{x1D70E}_{0}$ is the Weber number. Note that, in this case, 
$Re_{\unicode[STIX]{x1D70E}}=\sqrt{La}\gg 1$, where 
$La=\unicode[STIX]{x1D70C}\bar{R}\unicode[STIX]{x1D70E}_{0}/\unicode[STIX]{x1D707}^{2}$ is the Laplace number, which may also be written as 
$La=Oh^{-2}$ in terms of the usual Ohnesorge number, 
$Oh=\unicode[STIX]{x1D707}/\sqrt{\unicode[STIX]{x1D70C}\bar{R}\unicode[STIX]{x1D70E}_{0}}$. In the opposite limit of Stokes flow, 
$Re_{\unicode[STIX]{x1D70E}}\ll 1$, the appropriate balance is 
$\unicode[STIX]{x1D70E}_{0}/\bar{R}\sim \unicode[STIX]{x1D707}U_{\unicode[STIX]{x1D70E}}/\bar{R}$, whence 
$U_{\unicode[STIX]{x1D70E}}\sim \unicode[STIX]{x1D70E}_{0}/\unicode[STIX]{x1D707}$, usually referred to as the visco-capillary velocity. In this limit, 
$U/U_{\unicode[STIX]{x1D70E}}\sim Ca$, where 
$Ca=\unicode[STIX]{x1D707}U/\unicode[STIX]{x1D70E}_{0}$ is the capillary number, and 
$Re_{\unicode[STIX]{x1D70E}}=La\ll 1$. Therefore, the condition that must be satisfied for the temporal and spatial approaches to be equivalent is that 
$\sqrt{We}\gg 1$ when 
$\sqrt{La}\gg 1$, or that 
$Ca\gg 1$ when 
$La\ll 1$. It is also important to point out that the formation of a slender jet from a nozzle requires that 
$We>We_{c}\sim O(1)$ when 
$\sqrt{La}\gg 1$, or that 
$Ca>Ca_{c}\sim O(1)$ when 
$La\ll 1$, where 
$We_{c}$ and 
$Ca_{c}$ are the critical Weber and capillary numbers for the transition from convective to absolute instability (Leib & Goldstein Reference Leib and Goldstein1986a,Reference Leib and Goldsteinb).
Many experimental studies have been carried out, from the first investigations of Savart (Reference Savart1833), Magnus (Reference Magnus1859), Plateau (Reference Plateau1873), Rayleigh (Reference Rayleigh1882) and Donnelly & Glaberson (Reference Donnelly and Glaberson1966), to the highly accurate measurements of González & García (Reference González and García2009), whose aim was to describe the mechanism of instability and to measure the growth rate of the associated waves in the linear regime. These experiments have shown an excellent agreement with the dispersion relation obtained by Rayleigh (Reference Rayleigh1878), Lord Rayleigh (Reference Lord Rayleigh Sec.1892) and by Chandrasekhar (Reference Chandrasekhar1961). It is important to emphasise that, although linear stability theory cannot describe the final stages of the dynamics prior to pinch-off, it can be used to predict the break-up time 
$\bar{t}_{b}$ with small relative errors, provided that the initial amplitude of the disturbance, 
$\unicode[STIX]{x1D700}$, satisfies 
$\unicode[STIX]{x1D716}=\unicode[STIX]{x1D700}/\bar{R}\ll 1$. In the spatial setting, this fact can be used to estimate the break-up length as 
$U\bar{t}_{b}$, in close agreement with experiments (Kalaaji et al. Reference Kalaaji, Lopez, Attane and Soucemarianadin2003; González & García Reference González and García2009).
However, to describe the satellite formation process, which is the main objective of the present study, a nonlinear approach is needed. In particular, Goedde & Yuen (Reference Goedde and Yuen1970) first investigated such nonlinear effects in detail, comparing their experiments with the weakly nonlinear theory of Yuen (Reference Yuen1968). The satellite drop formation process was first quantified by Rutland & Jameson (Reference Rutland and Jameson1970) and Lafrance (Reference Lafrance1975), while Chaudhary & Maxworthy (Reference Chaudhary and Maxworthy1980) studied how satellite drop formation is affected by forcing the liquid jet with different harmonics, revealing the conditions needed to inhibit their formation. These efforts to control drop formation were mainly motivated by the practical need of increasing the performance of the ink-jet printing devices under development at that time. The increase in computational power finally allowed a fully nonlinear approach by means of direct numerical simulations of the axisymmetric Navier–Stokes equations. In particular Mansour & Lundgren (Reference Mansour and Lundgren1990) and Ashgriz & Mashayek (Reference Ashgriz and Mashayek1995) computed the satellite droplet radii just prior to pinch-off, finding an excellent agreement with the experiments of Rutland & Jameson (Reference Rutland and Jameson1970) and Lafrance (Reference Lafrance1975). Due to the high numerical cost of accurately solving the Navier–Stokes equations with a free boundary up to times close to the break-up singularity, several works have been devoted to develop one-dimensional approximations by expanding the flow variables as powers of the radial coordinate (see e.g. Lee Reference Lee1974; Eggers & Dupont Reference Eggers and Dupont1994; García & Castellanos Reference García and Castellanos1994). These models have been shown to work reasonably well in different configurations (see e.g. Ambravaneswaran et al. Reference Ambravaneswaran, Subramani, Phillips and Basaran2004; Rubio-Rubio, Sevilla & Gordillo Reference Rubio-Rubio, Sevilla and Gordillo2013; Martínez-Calvo, Rubio-Rubio & Sevilla Reference Martínez-Calvo, Rubio-Rubio and Sevilla2018). Of particular importance is the fact that the leading-order model allowed the unravelling of the universal self-similar structure of the local flow close to the singularity (Eggers Reference Eggers1993; Papageorgiou Reference Papageorgiou1995).
The presence of surfactant molecules at an interface induces an effective surface rheology by means of Marangoni stresses and surface viscosities (for reviews, see Fuller & Vermant Reference Fuller and Vermant2012; Langevin Reference Langevin2014), leading to substantial changes in the dynamics with respect to the case of clean interfaces. Indeed, different flow configurations of technological interest are affected by surfactants, for instance liquid bridges (Liao, Franses & Basaran Reference Liao, Franses and Basaran2006), dip coating (Scheid et al. Reference Scheid, Delacotte, Dollet, Rio, Restagno, van Nierop, Cantat, Langevin and Stone2010; Delacotte et al. Reference Delacotte, Montel, Restagno, Scheid, Dollet, Stone, Langevin and Rio2012; Champougny et al. Reference Champougny, Scheid, Restagno, Vermant and Rio2015) or drop break-up (Roché et al. Reference Roché, Aytouna, Bonn and Kellay2009; Ponce-Torres et al. Reference Ponce-Torres, Montanero, Herrada, Vega and Vega2017; Kamat et al. Reference Kamat, Wagoner, Thete and Basaran2018), to cite a few. Regarding liquid threads, the effect of surfactants has been explored by means of theory (Timmermans & Lister Reference Timmermans and Lister2002; Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2018) and numerical simulations (Campana & Saita Reference Campana and Saita2006; Dravid et al. Reference Dravid, Songsermpong, Xue, Corvalan and Sojka2006; McGough & Basaran Reference McGough and Basaran2006; Kamat et al. Reference Kamat, Wagoner, Thete and Basaran2018). In particular, the two latter works focused on the micro-thread cascade that appears close to break-up due to the presence of surfactants. These works also analyse the different scalings close to pinch-off and the evolution of the minimum radius of the thread and its axial position during the unfolding of the micro-cascade. Moreover, Kamat et al. (Reference Kamat, Wagoner, Thete and Basaran2018) revealed that the mechanism responsible for the dynamical surface tension effects induced by surfactants in filament break-up is the action of Marangoni stresses rather than the lowering of surface tension. In the case of a surfactant-free liquid thread, Ashgriz & Mashayek (Reference Ashgriz and Mashayek1995) already reported the axial movement of the location of minimum radius for low-viscosity filaments, 
$La\gg 1$. More recently, Castrejón-Pita et al. (Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, Lister and Basaran2015) generalised this result, showing that this translation occurs for any finite value of 
$La$, leading to the asymptotic inertial–viscous regime (Eggers Reference Eggers1993). Furthermore, by means of experiments and high-precision numerical simulations of the full axisymmetric Navier–Stokes equations, Castrejón-Pita et al. (Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, Lister and Basaran2015) demonstrated that, depending on the value of 
$La$, the thinning of the filament experiences different transitions that delay the occurrence of the universal inertial–viscous regime. In contrast with the latter studies, which focused on a detailed description of the transitions between the different scaling laws prior to pinch-off, the present contribution aims at providing a global parametric description of the satellite formation process when the interface is coated with an insoluble surfactant monolayer.
Given the success of the leading-order one-dimensional approximation in capturing the nonlinear dynamics of clean interfaces under certain configurations and values of 
$La$ (Eggers & Dupont Reference Eggers and Dupont1994; Ambravaneswaran, Wilkes & Basaran Reference Ambravaneswaran, Wilkes and Basaran2002; Notz & Basaran Reference Notz and Basaran2004; Yildirim, Xu & Basaran Reference Yildirim, Xu and Basaran2005; Subramani et al. Reference Subramani, Yeoh, Suryo, Xu, Ambravaneswaran and Basaran2006; Rubio-Rubio et al. Reference Rubio-Rubio, Sevilla and Gordillo2013), similar models have been derived that account for the presence of surfactants, and numerically solved for different flow configurations (Ambravaneswaran & Basaran Reference Ambravaneswaran and Basaran1999; Craster, Matar & Papageorgiou Reference Craster, Matar and Papageorgiou2002, Reference Craster, Matar and Papageorgiou2009; Xu, Liao & Basaran Reference Xu, Liao and Basaran2007). However, as pointed out by Timmermans & Lister (Reference Timmermans and Lister2002) and Martínez-Calvo & Sevilla (Reference Martínez-Calvo and Sevilla2018), a higher-order approximation is needed when the surface viscoelastic stresses are large enough. The failure of the leading-order one-dimensional models to describe the flow for large enough elastic and surface viscous stresses is due to the fact that the velocity profile is uniform in the leading-order equations, and cannot accomodate the shear induced by tangential interfacial stresses. Therefore, following the same strategy as Mansour & Lundgren (Reference Mansour and Lundgren1990) and Ashgriz & Mashayek (Reference Ashgriz and Mashayek1995) for a clean interface, and Dravid et al. (Reference Dravid, Songsermpong, Xue, Corvalan and Sojka2006) and McGough & Basaran (Reference McGough and Basaran2006) for a surfactant-laden interface, in the present contribution our approach is to numerically integrate the Navier–Stokes equations in a temporal setting, thereby avoiding the approximations involved in one-dimensional models. Unlike Dravid et al. (Reference Dravid, Songsermpong, Xue, Corvalan and Sojka2006), we use a nonlinear equation of state to relate the surface tension to the surfactant concentration, derived from first principles, that leads to substantial differences calling out for a careful experimental analysis. Moreover, we perform an exhaustive parametric study, accurately quantifying the volume of the satellite droplet prior to pinch-off and the amount of surfactant trapped at its surface, as a function of the two dimensionless governing parameters, namely the elasticity parameter and the Laplace number.
The remainder of the paper is organised as follows. In § 2 we describe the mathematical model and the numerical method employed for the simulations. In § 3 we first validate the simulations by comparing the initial growth rate of small harmonic disturbances with the results provided by a temporal stability analysis. We then unravel the structure of the parameter plane spanned by the Laplace and elasticity numbers in terms of the satellite formation process, followed by a detailed analysis of the volume of the satellite droplets and their shape at break-up, the mass of surfactant trapped at their surface and the nonlinear correction to the linear break-up time. The detailed time evolution is studied in several representative cases to provide physical explanations of the results obtained. Conclusions are drawn in § 4. Finally, a stringent validation of our numerical technique is presented in the Appendix.
2 Mathematical model and numerical method
We consider the axisymmetric motion of an infinitely long liquid thread of density 
$\unicode[STIX]{x1D70C}$, viscosity 
$\unicode[STIX]{x1D707}$, surface tension 
$\bar{\unicode[STIX]{x1D70E}}$ and unperturbed radius 
$\bar{R}$, which occupies a volume 
${\mathcal{V}}(\bar{t})$ and is embedded in a passive ambient at constant pressure 
$p_{a}$ in the absence of gravity. The interface 
$\unicode[STIX]{x2202}{\mathcal{V}}(\bar{t})$, placed at a radial position 
$\bar{r}=\bar{a}(\bar{z},\bar{t})$, is coated with a superficial concentration 
$\bar{\unicode[STIX]{x1D6E4}}$ of insoluble surfactant molecules (see figure 1a). Note that 
$\bar{r}$, 
$\bar{z}$ and 
$\bar{t}$ stand for the radial and axial coordinates and time, respectively. Henceforth an overbar will denote dimensional variables if not specified otherwise. The effect of the surfactant adsorbed at the interface is to reduce the effective surface tension by an amount that depends on 
$\bar{\unicode[STIX]{x1D6E4}}$, and thus any disturbance of the interface shape generates an imbalance in 
$\bar{\unicode[STIX]{x1D6E4}}$ that produces a surface stress due to gradients of 
$\bar{\unicode[STIX]{x1D70E}}(\bar{\unicode[STIX]{x1D6E4}})$. For simplicity, in the present work we assume that surface viscous stresses can be neglected, thus disregarding the role of the surface shear and dilatational viscosities, 
$\unicode[STIX]{x1D707}_{s}$ and 
$\unicode[STIX]{x1D705}_{s}$, respectively. The latter approximation is accurate provided that the corresponding Boussinesq numbers are small, namely 
$\mathscr{B}_{\unicode[STIX]{x1D707}}=\unicode[STIX]{x1D707}_{s}/(\unicode[STIX]{x1D707}\bar{R})\ll 1$ and 
$\mathscr{B}_{\unicode[STIX]{x1D705}}=\unicode[STIX]{x1D705}_{s}/(\unicode[STIX]{x1D707}\bar{R})\ll 1$ (Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2018). The problem is non-dimensionalised with the visco-capillary time, 
$\unicode[STIX]{x1D707}\bar{R}/\unicode[STIX]{x1D70E}_{0}$, as characteristic time and with 
$\bar{R}$ as characteristic length, 
$\unicode[STIX]{x1D70E}_{0}$ being the surface tension associated with the initial concentration of insoluble surfactant at the interface 
$\bar{\unicode[STIX]{x1D6E4}}(\bar{z},0)=\unicode[STIX]{x1D6E4}_{0}$, which are used to scale the surface tension and the surface concentration, respectively.
Figure 1. (a) Dimensional sketch of the flow configuration. (b) Example of a liquid thread approaching pinch-off for 
$La=0.01$, 
$\unicode[STIX]{x1D6FD}=1$, 
$\unicode[STIX]{x1D716}=10^{-3}$ and 
$k=k_{m}=0.516$ at time 
$t=123$. The contour map represents the dimensionless pressure field 
$p$, and the arrows show the dimensionless velocity field 
$\boldsymbol{u}$, both at the top, while the deformed mesh is shown at the bottom.
The flow is governed by the dimensionless Navier–Stokes equations
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0\quad \text{at }{\mathcal{V}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle La\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\right)=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}\quad \text{at }{\mathcal{V}}, & \displaystyle\end{eqnarray}$$ where 
$\boldsymbol{u}(\boldsymbol{x},t)=u\,\boldsymbol{e}_{r}+w\,\boldsymbol{e}_{z}$ is the velocity field and 
$u$, 
$w$ and 
$\boldsymbol{e}_{r}$, 
$\boldsymbol{e}_{z}$ are the radial and axial velocity components and the corresponding unit vectors, respectively. In (2.2), 
$\unicode[STIX]{x1D64F}=-p\unicode[STIX]{x1D644}+[\unicode[STIX]{x1D735}\boldsymbol{u}+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}}]$ is the stress tensor for an incompressible Newtonian liquid, 
$\unicode[STIX]{x1D644}$ is the standard identity tensor and 
$p(\boldsymbol{x},t)$ is the pressure field. The numerical simulations reported herein were performed using an arbitrary Lagrangian–Eulerian (ALE) method, in which the domain 
$\boldsymbol{x}(\boldsymbol{X},t)\in {\mathcal{V}}(t)$ is parametrised by the initial position 
$\boldsymbol{X}=\boldsymbol{x}(\boldsymbol{X},0)\in {\mathcal{V}}(0)$, defining a time-dependent displacement field, 
$\boldsymbol{x}-\boldsymbol{X}$, which is enforced to satisfy the Laplace equation with proper boundary conditions specified below. The local time derivatives are evaluated in the spatial reference frame as
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}\hat{\boldsymbol{u}}}{\unicode[STIX]{x2202}t}-\frac{\unicode[STIX]{x2202}\boldsymbol{x}}{\unicode[STIX]{x2202}t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\quad \text{at }{\mathcal{V}},\end{eqnarray}$$ where 
$\hat{\boldsymbol{u}}(\boldsymbol{X},t)=\boldsymbol{u}(\boldsymbol{x},t)$ is the velocity in the material reference frame.
Since the interface 
$\unicode[STIX]{x2202}{\mathcal{V}}$ is coated with surfactant, a surface transport equation is needed for 
$\unicode[STIX]{x1D6E4}(\boldsymbol{x},t)$:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}_{\!s}\boldsymbol{\cdot }(\unicode[STIX]{x1D6E4}\,\boldsymbol{u}_{s})=0\quad \text{at }\unicode[STIX]{x2202}{\mathcal{V}},\end{eqnarray}$$ where 
$\boldsymbol{u}_{s}=\boldsymbol{u}(\boldsymbol{x}_{s})$ is the liquid velocity at the interface and 
$\boldsymbol{x}_{s}$ represents any position at the surface 
$\boldsymbol{x}_{s}(\boldsymbol{X}_{s},t)\in \unicode[STIX]{x2202}{\mathcal{V}}(t)$, which is parametrised by its initial position 
$\boldsymbol{X}_{s}=\boldsymbol{x}_{s}(\boldsymbol{X}_{s},0)\in \unicode[STIX]{x2202}{\mathcal{V}}(0)$. Here 
$\unicode[STIX]{x1D735}_{\!s}=\unicode[STIX]{x1D644}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ is the surface gradient operator, where 
$\unicode[STIX]{x1D644}_{s}=\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n}$ is the surface projection tensor and 
$\boldsymbol{n}$ is the outer unit normal vector at the surface. The local time derivatives at the interface are evaluated in the spatial reference frame as
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D6E4}}}{\unicode[STIX]{x2202}t}-\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{s}}{\unicode[STIX]{x2202}t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{s}\unicode[STIX]{x1D6E4}\quad \text{at }\unicode[STIX]{x2202}{\mathcal{V}},\end{eqnarray}$$ where 
$\hat{\unicode[STIX]{x1D6E4}}(\boldsymbol{X}_{s},t)=\unicode[STIX]{x1D6E4}(\boldsymbol{x}_{s},t)$ is the concentration of surfactant in the material frame of reference, which is needed in order to be implemented with the ALE method that is used in the present work. The reader is referred to the works of Stone (Reference Stone1990), Wong, Rumschitzki & Maldarelli (Reference Wong, Rumschitzki and Maldarelli1996) and Pereira & Kalliadasis (Reference Pereira and Kalliadasis2008) for further details of the time derivative of a surface quantity.
Note that the surface diffusion of surfactant has been neglected in the transport equation (2.4). Indeed, in the present work we only consider the limit where the surface Péclet number 
$Pe_{s}=U_{sc}\bar{R}/D_{s}\rightarrow \infty$, where 
$D_{s}$ is the surface diffusion coefficient and 
$U_{sc}$ is the characteristic liquid velocity at the free surface. The correct scaling for 
$U_{sc}$ depends on the value of 
$La$. In the limit of dominant inertia, 
$La\gg 1$, the appropriate velocity scale is the capillary velocity, 
$[\unicode[STIX]{x1D70E}_{0}/(\unicode[STIX]{x1D70C}\bar{R}^{3})]^{1/2}$, so that 
$Pe_{s}=[\unicode[STIX]{x1D70E}_{0}\bar{R}/(\unicode[STIX]{x1D70C}D_{s}^{2})]^{1/2}$. For instance, if we consider a water thread of radius within the range 
$1$–
$100~\unicode[STIX]{x03BC}\text{m}$, the corresponding Laplace numbers lie in the range 
$10^{2}\lesssim La\lesssim 10^{4}$. Typical values of 
$D_{s}$ for SDS, SB12 and other monomers in aqueous solution are within the range 
$10^{-9}\lesssim D_{s}\lesssim 10^{-8}~\text{m}^{2}~\text{s}^{-1}$ when 
$\unicode[STIX]{x1D6E4}$ is below the critical micelle concentration (Siderius, Kehl & Leaist Reference Siderius, Kehl and Leaist2002), providing values of the surface Péclet number in the range 
$10^{4}\lesssim Pe_{s}\lesssim 10^{5}$. Therefore, in configurations where 
$La\gg 1$, it is expected that surface diffusion has a very small effect. In the opposite limit of dominant viscous forces, 
$La\lesssim 1$, the appropriate velocity scale is the visco-capillary velocity, 
$\unicode[STIX]{x1D70E}_{0}/\unicode[STIX]{x1D707}$, leading to 
$Pe_{s}=\unicode[STIX]{x1D70E}_{0}\bar{R}/(\unicode[STIX]{x1D707}D_{s})$. Considering, for instance, a polydimethylsiloxane silicon oil of dynamic viscosity in the range 
$0.1$–
$10$ Pa s, density 
$\unicode[STIX]{x1D70C}\approx 970$ kg m
$^{-3}$ and surface tension 
$\unicode[STIX]{x1D70E}_{0}\approx 21.1$ mN m
$^{-1}$, the Laplace number takes values in the range 
$10^{-4}\lesssim La\lesssim 1$. Although we are not aware of experimental studies reporting typical values of 
$D_{s}$ in highly viscous solutions, if we assume that they are of the same order of magnitude as those of aqueous solutions, the Péclet number lies in the range 
$1\lesssim Pe_{s}\lesssim 10^{6}$. It is thereby deduced that, when 
$La\lesssim 1$, there may be cases where surface diffusion cannot be neglected in the analysis. Therefore, although the influence of surface diffusion on the satellite droplet formation process is not addressed in the present work, it clearly deserves further study, particularly in the case of highly viscous threads.
The presence of surfactant at the interface modifies 
$\unicode[STIX]{x1D70E}$ by decreasing its value as 
$\unicode[STIX]{x1D6E4}$ increases, and thus the stress balance at the interface takes the following form in the limit 
$\mathscr{B}_{\unicode[STIX]{x1D707}}\ll 1$, 
$\mathscr{B}_{\unicode[STIX]{x1D705}}\ll 1$ (Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2018):
$$\begin{eqnarray}\unicode[STIX]{x1D64F}\boldsymbol{\cdot }\boldsymbol{n}=\unicode[STIX]{x1D735}_{\!s}\unicode[STIX]{x1D70E}-\boldsymbol{n}(\unicode[STIX]{x1D735}_{\!s}\boldsymbol{\cdot }\boldsymbol{n})\unicode[STIX]{x1D70E}\quad \text{at }\unicode[STIX]{x2202}{\mathcal{V}},\end{eqnarray}$$ where the viscous stress exerted by the ambient fluid on the interface has been neglected and the ambient pressure 
$p_{a}$ has been set to zero without loss of generality. Additionally, the kinematic condition must also hold at the interface:
$$\begin{eqnarray}\boldsymbol{u}_{s}\boldsymbol{\cdot }\boldsymbol{n}=\frac{\unicode[STIX]{x2202}\boldsymbol{x}_{s}}{\unicode[STIX]{x2202}t}\boldsymbol{\cdot }\boldsymbol{n}\quad \text{at }\unicode[STIX]{x2202}{\mathcal{V}}.\end{eqnarray}$$ Finally, an equation of state that relates the surface tension, 
$\unicode[STIX]{x1D70E}$, to the surface concentration of surfactant, 
$\unicode[STIX]{x1D6E4}$, is also needed. Surface-active molecules at the interface induce a surface pressure 
$\bar{\unicode[STIX]{x1D6F1}}$ which depends on the surfactant concentration, 
$\bar{\unicode[STIX]{x1D6F1}}=\bar{\unicode[STIX]{x1D6F1}}(\bar{\unicode[STIX]{x1D6E4}})$. The surface pressure is defined as the difference in the surface tension due to the presence of surfactant, 
$\bar{\unicode[STIX]{x1D6F1}}(\bar{\unicode[STIX]{x1D6E4}})=\unicode[STIX]{x1D70E}_{clean}-\bar{\unicode[STIX]{x1D70E}}(\bar{\unicode[STIX]{x1D6E4}})$, and thus 
$\unicode[STIX]{x1D735}_{s}\bar{\unicode[STIX]{x1D6F1}}=-\unicode[STIX]{x1D735}_{s}\bar{\unicode[STIX]{x1D70E}}$. In addition, the Gibbs elasticity 
$E$ relates the changes of interface area, 
$\bar{A}$, to the surface pressure through the surface compressibility 
$1/E=-(1/\bar{A})(\unicode[STIX]{x2202}\bar{A}/\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D6F1}})_{\bar{T}}$, where 
$\bar{T}$ is the temperature at the interface, which is assumed to remain constant. Hence,
$$\begin{eqnarray}E=-\bar{A}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D6F1}}}{\unicode[STIX]{x2202}\bar{A}}=\bar{A}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70E}}}{\unicode[STIX]{x2202}\bar{A}}=-\bar{\unicode[STIX]{x1D6E4}}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70E}}}{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D6E4}}},\end{eqnarray}$$ where in the last equation it has been taken into account that, in the insoluble case considered in the present work, the number of surfactant molecules is conserved at the interface. Equation (2.8) can be used to relate 
$\bar{\unicode[STIX]{x1D70E}}$ and 
$\bar{\unicode[STIX]{x1D6E4}}$:
$$\begin{eqnarray}\unicode[STIX]{x1D735}_{s}\bar{\unicode[STIX]{x1D70E}}=\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70E}}}{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D6E4}}}\unicode[STIX]{x1D735}_{s}\bar{\unicode[STIX]{x1D6E4}}=-\frac{E}{\bar{\unicode[STIX]{x1D6E4}}}\unicode[STIX]{x1D735}_{s}\bar{\unicode[STIX]{x1D6E4}}.\end{eqnarray}$$ Making 
$\bar{\unicode[STIX]{x1D70E}}$ and 
$\bar{\unicode[STIX]{x1D6E4}}$ dimensionless with 
$\unicode[STIX]{x1D70E}_{0}$ and 
$\unicode[STIX]{x1D6E4}_{0}$, respectively, equation (2.9) finally yields the dimensionless equation of state
$$\begin{eqnarray}\unicode[STIX]{x1D70E}=1-\unicode[STIX]{x1D6FD}\ln \unicode[STIX]{x1D6E4},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FD}=E/\unicode[STIX]{x1D70E}_{0}$ is the so-called elasticity parameter, also referred to as the Marangoni number (Champougny et al. Reference Champougny, Scheid, Restagno, Vermant and Rio2015). Note that, in the limit of small surface concentration variations around the initial value, 
$\bar{\unicode[STIX]{x1D6E4}}=\unicode[STIX]{x1D6E4}_{0}+\unicode[STIX]{x1D6FF}\bar{\unicode[STIX]{x1D6E4}}$, with 
$\unicode[STIX]{x1D6FF}\bar{\unicode[STIX]{x1D6E4}}\ll \unicode[STIX]{x1D6E4}_{0}$, one has 
$\unicode[STIX]{x1D6E4}=1+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6E4}$ with 
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6E4}\ll 1$, and the equation of state (2.10) can be linearised to yield 
$\unicode[STIX]{x1D70E}=1-\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6E4}$, which is equivalent to the equation of state employed by Dravid et al. (Reference Dravid, Songsermpong, Xue, Corvalan and Sojka2006). However, it is important to emphasise that the relative variations of 
$\unicode[STIX]{x1D6E4}$ during the thread break-up process are not small, as demonstrated in § 3. Therefore, the use of a linearised equation of state introduces considerable errors and is not justified. At this point, the limitations of the nonlinear equation of state (2.10) should be clearly stated. Indeed, the main shortcoming of (2.10) is that 
$\unicode[STIX]{x1D70E}\rightarrow \infty$ as 
$\unicode[STIX]{x1D6E4}\rightarrow 0$, which eventually occurs as the surfactant is depleted from the pinch-off region due to the local advection out of the collapsing neck. Hence, a different equation of state is required to properly model the dynamics of the smallest scales close to break-up, which properly captures the saturation of 
$\unicode[STIX]{x1D70E}$ to the clean interface value as 
$\unicode[STIX]{x1D6E4}\rightarrow 0$ (see e.g. McGough & Basaran Reference McGough and Basaran2006; Kamat et al. Reference Kamat, Wagoner, Thete and Basaran2018). Nevertheless, for the purposes of the present contribution, the equation of state (2.10) is perfectly valid. In effect, at the smallest scales that need to be resolved to provide robust measures of 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$, the relative variations of 
$\unicode[STIX]{x1D6E4}$ and 
$\unicode[STIX]{x1D70E}$ are small enough as to guarantee the validity of (2.10) in all the results reported herein.
The surface stress balance (2.6) can be rewritten in terms of 
$\unicode[STIX]{x1D6E4}$ as
$$\begin{eqnarray}\unicode[STIX]{x1D64F}\boldsymbol{\cdot }\boldsymbol{n}=-\frac{\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6E4}}\,\unicode[STIX]{x1D735}_{s}\unicode[STIX]{x1D6E4}-\boldsymbol{n}\left(\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }\boldsymbol{n}\right)\left(1-\unicode[STIX]{x1D6FD}\ln \unicode[STIX]{x1D6E4}\right)\quad \text{at }\unicode[STIX]{x2202}{\mathcal{V}},\end{eqnarray}$$ which, together with (2.1)–(2.5) and (2.7), form a closed system to determine 
$\boldsymbol{u}$, 
$p$, 
$\unicode[STIX]{x1D6E4}$ and 
$\boldsymbol{x}_{s}$.
Concerning the computational domain and the corresponding boundary conditions, in the temporal approach adopted herein we only consider half a perturbation wavelength subjected to the following symmetry conditions:
$$\begin{eqnarray}w=0,\quad \frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}=0\quad \text{and}\quad \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6E4}}{\unicode[STIX]{x2202}z}=0\quad \text{at }z=0,\unicode[STIX]{x03C0}/k,\end{eqnarray}$$ where 
$k$ is the dimensionless axial wavenumber, together with the axisymmetry condition
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}w}{\unicode[STIX]{x2202}r}=0\quad \text{and}\quad u=0\quad \text{at }r=0.\end{eqnarray}$$ Finally, regarding the initial conditions imposed at 
$t=0$, we perturb the position of the liquid cylinder with a harmonic disturbance of amplitude 
$\unicode[STIX]{x1D716}$:
$$\begin{eqnarray}\boldsymbol{x}_{s}=z\boldsymbol{e}_{z}+[R-\unicode[STIX]{x1D716}\cos (kz)]\boldsymbol{e}_{r},\end{eqnarray}$$ where 
$R=(1-\unicode[STIX]{x1D716}^{2}/2)^{1/2}$ is a dimensionless radius defined in terms of 
$\unicode[STIX]{x1D716}$, such that the liquid volume remains constant as 
$\unicode[STIX]{x1D716}$ varies (Ashgriz & Mashayek Reference Ashgriz and Mashayek1995). We also assume that the liquid thread is initially at rest and that the surfactant concentration is uniform:
$$\begin{eqnarray}\boldsymbol{u}(\boldsymbol{x},0)=\mathbf{0},\quad \unicode[STIX]{x1D6E4}(\boldsymbol{x}_{s},0)=1.\end{eqnarray}$$ Note that the assumption of a uniform initial concentration of surfactant is a good approximation, since our main results have been obtained in the limit 
$\unicode[STIX]{x1D716}\ll 1$ in which the deviations from a uniform concentration can be neglected. As explained in § 1, our results can also be applied to describe the spatial instability and subsequent downstream break-up of liquid jets moving with uniform velocity 
$U$ with respect to the injector reference frame, provided that 
$U\gg U_{\unicode[STIX]{x1D70E}}$, where 
$U_{\unicode[STIX]{x1D70E}}$ is the speed of small-amplitude capillary waves. If the latter condition is satisfied, the spatial evolution of the jet is obtained by the downstream advection of the temporal results presented herein with a uniform velocity 
$U$. In particular, the jet break-up length is given by 
$U\bar{t}_{b}$ to a first approximation.
The problem depends on four dimensionless parameters, namely the Laplace number 
$La$, the elasticity parameter 
$\unicode[STIX]{x1D6FD}$, the axial wavenumber 
$k$ and the amplitude of the initial perturbation 
$\unicode[STIX]{x1D716}$. However, in the present work we are concerned with the unforced break-up of cylindrical threads due to small-radius disturbances. Therefore, all the results were obtained by setting 
$k=k_{m}$, where 
$k_{m}(La,\unicode[STIX]{x1D6FD})$ is the most unstable wavenumber (see § 3.1). Moreover, it will be shown that, in the small-disturbance limit, 
$\unicode[STIX]{x1D716}\ll 1$, the only result that depends on 
$\unicode[STIX]{x1D716}$ is the break-up time of the thread, 
$t_{b}(La,\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D716})$. However, our results have revealed that the functional dependence of 
$t_{b}$ can be split into a contribution predicted by linear theory in explicit form, 
$t_{b,L}(La,\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D716})$, plus a nonlinear correction, 
$\unicode[STIX]{x0394}t_{NL}(La,\unicode[STIX]{x1D6FD})$, which does not depend on 
$\unicode[STIX]{x1D716}$. Consequently, only two independent dimensionless parameters appear in our formulation, namely 
$La$ and 
$\unicode[STIX]{x1D6FD}$.
To perform the numerical simulations, the liquid domain 
$0\leqslant r\leqslant a(z,t),0\leqslant z\leqslant \unicode[STIX]{x03C0}/k$ is partitioned into a rectangular or triangular finite-element mesh which is dynamically deformed using the ALE method. In particular, the displacement field, 
$\boldsymbol{x}-\boldsymbol{X}$, is enforced to satisfy the Laplace equation, and the normal mesh velocity, 
$\boldsymbol{n}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{u}$, solves the kinematic condition (2.7). To that end, equations (2.1)–(2.4), together with the boundary and initial conditions (2.12)–(2.15), are written in weak form following the methodology described by Rivero-Rodríguez & Scheid (Reference Rivero-Rodríguez and Scheid2018a,Reference Rivero-Rodríguez and Scheidb), and the spatial discretisation is carried out using the finite-element method provided by COMSOL, where Lagrange linear (P1) elements are used for 
$p$ and quadratic (P2) elements are used for 
$\boldsymbol{x}$, 
$\boldsymbol{u}$ and 
$\unicode[STIX]{x1D6E4}$. The time discretisation was performed using the first-order backward Euler method with adaptive time stepping. Figure 1(b) shows a representative deformed mesh for a simulation with 
$La=0.01$, 
$\unicode[STIX]{x1D6FD}=1$, 
$\unicode[STIX]{x1D716}=10^{-3}$ and 
$k=k_{m}=0.516$ at time 
$t=123$, together with the pressure field as a contour plot and the velocity field represented by arrows. All the results reported were carefully checked as being mesh-independent, with an integration tolerance of the order of 
$10^{-6}$–
$10^{-7}$. In addition, it was checked that the relative variations of liquid volume and surfactant mass where smaller than 
$10^{-5}$ during each simulation. The numerical code has been validated with the linear theory in § 3.1. In the nonlinear regime, the validation was performed by comparing our results with those of Ashgriz & Mashayek (Reference Ashgriz and Mashayek1995) for a clean interface and with those of McGough & Basaran (Reference McGough and Basaran2006) and Kamat et al. (Reference Kamat, Wagoner, Thete and Basaran2018) for a surfactant-laden thread (not shown). In particular, the Appendix is devoted to show the performance of our numerical framework close to pinch-off, comparing our results with the different theoretical scalings of the minimum radius as a function of time to break-up.
3 Results and discussion
Since we are interested in the spontaneous break-up of the surfactant-laden thread, all the results were computed from an initial condition where the liquid cylinder is perturbed with the wavenumber of maximum amplification, 
$k_{m}(La,\unicode[STIX]{x1D6FD})$. Hence, the results of a linear stability analysis are first summarised in § 3.1 to obtain 
$k_{m}$ and 
$\unicode[STIX]{x1D714}_{m}$, the latter being the maximum temporal growth rate. Note that 
$k_{m}$ is needed to define the initial geometry and the initial condition (2.14), while 
$\unicode[STIX]{x1D714}_{m}$ is used to compute the nonlinear correction to the linear break-up time, which is defined in § 3.2. In addition, the linear theory has also been used to validate the numerical code by comparing the associated maximum temporal growth rate, 
$\unicode[STIX]{x1D714}_{m}$, with the results extracted from the numerical simulations during the initial transient of exponential amplitude growth. Sections 3.2 and 3.3 are devoted to the analysis of the nonlinear break-up and the satellite formation dynamics, separating the weak-elasticity limit and the surfactant-laden case. To that end, we have performed direct numerical simulations of (2.1)–(2.15) until times very close to pinch-off. In particular, we report a parametric study for different values of 
$La$ and 
$\unicode[STIX]{x1D6FD}$, computing the volume of the satellite droplet, the mass of surfactant trapped at its interface, the satellite shape at pinch-off and the break-up time.
At this point, it has to be pointed out that a similar phenomenology was previously reported by Dravid et al. (Reference Dravid, Songsermpong, Xue, Corvalan and Sojka2006) for 
$La=0.01$ and 
$100$, although using the linearised equation of state 
$\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D6E4})=1-\unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D6E4}-1)$. In addition, those authors did not consider the natural break-up of the thread, since the disturbance wavenumber 
$k$ was restricted to fixed values different from the most amplified one, 
$k_{m}$.
3.1 Linear stability analysis
To obtain the dispersion relation 
$D(\unicode[STIX]{x1D714},k)=0$ relating the temporal growth rate 
$\unicode[STIX]{x1D714}$ and the axial wavenumber 
$k$, all the flow variables are slightly perturbed around a uniform stationary state and decomposed as temporal normal modes:
$$\begin{eqnarray}(u,w,p,a,\unicode[STIX]{x1D70E},\unicode[STIX]{x1D6E4})=(0,0,1,1,1,1)+\unicode[STIX]{x1D716}(\hat{u} ,\hat{v},\hat{p},\hat{a},\hat{\unicode[STIX]{x1D70E}},\hat{\unicode[STIX]{x1D6E4}})\exp (\text{i}kz+\unicode[STIX]{x1D714}t).\end{eqnarray}$$ Introducing (3.1) into the system (2.1)–(2.4) and keeping terms proportional to 
$\unicode[STIX]{x1D716}$, the following dispersion relation is obtained:
$$\begin{eqnarray}\displaystyle & & \displaystyle La\,\unicode[STIX]{x1D714}^{2}F(k)-k^{2}(1-k^{2})+\unicode[STIX]{x1D6FD}k^{2}[1+F(k)(F(\tilde{k})-2)]\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{k^{4}}{La}\left[4-\frac{\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D714}}\left(2-\frac{1-k^{2}}{\unicode[STIX]{x1D714}}\right)\right][F(k)-F(\tilde{k})]+2\unicode[STIX]{x1D714}k^{2}(2F(k)-1)=0,\end{eqnarray}$$ where 
$\tilde{k}=\sqrt{k^{2}+La\,\unicode[STIX]{x1D714}}$ and 
$F(x)=xI_{0}(x)/I_{1}(x)$. Here, 
$I_{n}(x)$ denotes the 
$n$th-order modified Bessel function of the first kind. Note that dispersion relation (3.2) is exactly the same as the one deduced by Timmermans & Lister (Reference Timmermans and Lister2002), and is also a particular case of the one provided by Martínez-Calvo & Sevilla (Reference Martínez-Calvo and Sevilla2018) in the limit of negligible surface viscosities. The Rayleigh–Chandrasekhar dispersion relation is recovered when 
$\unicode[STIX]{x1D6FD}\rightarrow 0$ (Lord Rayleigh Reference Lord Rayleigh Sec.1892; Chandrasekhar Reference Chandrasekhar1961).
Figure 2. (a) Semi-logarithmic plot of the radius amplitude 
$A(t)$ as a function of time, extracted from two numerical simulations for 
$\unicode[STIX]{x1D716}=10^{-4}$, 
$\unicode[STIX]{x1D6FD}=1$ and two values of the Laplace number, namely 
$La=(0.01,100)$. The corresponding optimal wavenumbers, 
$k_{m}(La,\unicode[STIX]{x1D6FD})$, highlighted in (b) with stars, are used to build the initial conditions, and their values are indicated near each curve together with the associated linear temporal growth rates, 
$\unicode[STIX]{x1D714}_{m}(La,\unicode[STIX]{x1D6FD})$ and 
$La$. (b) Temporal growth rate 
$\unicode[STIX]{x1D714}$ as a function of the axial wavenumber 
$k$, computed with the dispersion relation (3.2) (solid lines) and with the numerical simulations (circles), for 
$\unicode[STIX]{x1D6FD}=1$ and two different values of 
$La=(0.01,100)$, indicated near each curve. The maximum growth rates 
$\unicode[STIX]{x1D714}_{m}$ computed in (a) are marked with stars. (c) Isocontours of the most amplified wavenumber 
$k_{m}(La,\unicode[STIX]{x1D6FD})$ and (d) its corresponding growth rate 
$\unicode[STIX]{x1D714}_{m}(La,\unicode[STIX]{x1D6FD})$.
As shown experimentally by Goedde & Yuen (Reference Goedde and Yuen1970), and numerically by Mansour & Lundgren (Reference Mansour and Lundgren1990) and Ashgriz & Mashayek (Reference Ashgriz and Mashayek1995), a convenient way to compute the temporal growth rate of small disturbances is through the radius amplitude, extracted from the present simulations as 
$A(t)=(\max _{z}[a(z,t)]-\min _{z}[a(z,t)])/2$. Figure 2(a) shows the temporal evolution of 
$A(t)$ in semi-logarithmic scale, extracted from two numerical simulations for an initial perturbation amplitude 
$\unicode[STIX]{x1D716}=10^{-4}$, an elasticity parameter 
$\unicode[STIX]{x1D6FD}=1$ and two values of the Laplace number, 
$La=0.01$ and 
$La=100$, close to the Stokes and Euler regimes, respectively. In each case, the most amplified wavenumber, 
$k_{m}(La,\unicode[STIX]{x1D6FD})$, is used to build the initial condition. As expected due to the smallness of 
$\unicode[STIX]{x1D716}$, figure 2(a) shows that during most of the time the amplitude grows exponentially, i.e. 
$A\propto \exp (\unicode[STIX]{x1D714}_{m}t)$, and thus the maximum temporal growth rate, 
$\unicode[STIX]{x1D714}_{m}(La,\unicode[STIX]{x1D6FD})$, can be easily computed as the slope of the linear region in the semi-logarithmic plot, 
$\unicode[STIX]{x1D714}_{m}=\text{d}\ln (A)/\text{d}t$. It can also be deduced from figure 2(a) that there is an initial transient during which the growth of 
$A(t)$ is not exponential, which can be explained by the fact that the initial conditions in the numerical simulations are imposed on the shape of the interface, but disregard the associated disturbances in the velocity, pressure and surfactant concentration fields. As shown in figure 2(b), this procedure was used to obtain 
$\unicode[STIX]{x1D714}$ for different values of 
$k$ (symbols), and the results are compared with the amplification curves 
$\unicode[STIX]{x1D714}(k)$ computed from the dispersion relation (3.2) (solid lines), affording an excellent agreement that validates the numerical code in the linear regime. Finally, figure 2(c,d) show the isocontours of 
$k_{m}$ and 
$\unicode[STIX]{x1D714}_{m}$, respectively, as a function of 
$La$ and 
$\unicode[STIX]{x1D6FD}$ extracted from (3.2), whose values will be used hereafter.
Figure 3. The structure of the 
$(La,\unicode[STIX]{x1D6FD})$ parameter plane. An abrupt transition takes place along the solid line, 
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{c}(La)$, across which both the satellite volume and the entrapped mass of surfactant experience a discontinuous jump, such that both magnitudes are larger above the solid line. The inset shows the jumps in the satellite volume, 
$\unicode[STIX]{x0394}V_{sat}(La)=V_{sat}(\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D6FD}_{c}\rightarrow 0^{+})-V_{sat}(\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D6FD}_{c}\rightarrow 0^{-})$, and in the associated entrapped mass of surfactant, 
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}_{sat}$. Both jumps, together with 
$\unicode[STIX]{x1D6FD}_{c}$, increase monotonically as 
$La$ decreases, and reach respective Stokes asymptotes as 
$La\rightarrow 0$, namely 
$\unicode[STIX]{x0394}V_{sat}\rightarrow 0.022$, 
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}_{sat}\rightarrow 0.045$ and 
$\unicode[STIX]{x1D6FD}_{c}\rightarrow 0.98$. The filled circle indicates the origin of the discontinuous transition, 
$(La,\unicode[STIX]{x1D6FD}_{c})=(7.5,0.55)$, at which both jumps become zero. For 
$La>7.5$, the satellite volume is a continuous function of 
$\unicode[STIX]{x1D6FD}$. The open circles correspond to the values of 
$La$ and 
$\unicode[STIX]{x1D6FD}$ of the shapes just before pinch-off shown in figure 4.
Figure 4. The satellite shapes just prior to pinch-off in the 
$(La,\unicode[STIX]{x1D6FD})$ parameter plane (see open circles in figure 3). The vertical lines indicate the axial positions, 
$z_{min}$, of the minimum thread radii, 
$a_{min}$.
3.2 Satellite formation regimes and transitions in the 
$(La,\unicode[STIX]{x1D6FD})$ parameter plane
Let us first present the structure of the 
$(La,\unicode[STIX]{x1D6FD})$ parameter plane in terms of the satellite formation process. To that end, we conducted an exhaustive parametric study in which the Laplace and elasticity parameters were varied in small steps within wide ranges, namely 
$0.01\leqslant La\leqslant 100$ and 
$0\leqslant \unicode[STIX]{x1D6FD}\leqslant 1$. Thus, for each pair of values of 
$La$ and 
$\unicode[STIX]{x1D6FD}$, we simulated the instability-driven time evolution of the thread from an initial condition with 
$\unicode[STIX]{x1D716}\ll 1$ until a time 
$t_{b}$ very close to break-up. In total, around 
$10^{4}$ time-dependent simulations were carried out to characterise the 
$(La,\unicode[STIX]{x1D6FD})$ parameter plane shown in figures 3 and 5.
Figure 5. Isocontours in the 
$(La,\unicode[STIX]{x1D6FD})$ parameter plane of (a) the normalised satellite volume 
$V_{sat}$, (b) the normalised mass of surfactant trapped at its interface 
$\unicode[STIX]{x1D6F4}_{sat}$, (c) the nonlinear correction to the break-up time 
$\unicode[STIX]{x0394}t_{NL}$ and (d) the sphericity of the satellite droplet 
${\mathcal{S}}$.
At this point, it is important to emphasise that the fate of the main and satellite drops after pinch-off is outside the scope of the present work, and therefore we do not explore the possible successive break-up events that may take place and lead to the formation of sub-satellites. Keeping this in mind, we have extracted the satellite volume at the last numerical step, 
$t=t_{b}$. Normalising its value with the total volume provides the definition
$$\begin{eqnarray}V_{sat}=\frac{\displaystyle \int _{0}^{z_{min}}a^{2}\,\text{d}z}{\displaystyle \int _{0}^{\unicode[STIX]{x03C0}/k_{m}}a^{2}\,\text{d}z},\end{eqnarray}$$ where 
$z_{min}$ is the axial position where the liquid column reaches its minimum radius, 
$a_{min}$, at 
$t=t_{b}$. A more common measure of the satellite size is its equivalent radius, 
$R_{sat}$, which is the radius of a spherical drop of the same volume as the satellite (Rutland & Jameson Reference Rutland and Jameson1971; Mansour & Lundgren Reference Mansour and Lundgren1990; Ashgriz & Mashayek Reference Ashgriz and Mashayek1995; Mashayek & Ashgriz Reference Mashayek and Ashgriz1995). All the results reported herein in terms of 
$V_{sat}$ can be easily converted to 
$R_{sat}$ through the equation 
$R_{sat}=[3\unicode[STIX]{x03C0}V_{sat}/(2k_{m})]^{1/3}$. Following the same procedure, we have also computed the mass of surfactant trapped at the satellite surface which, normalised with the total mass of surfactant, provides the definition
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{sat}=\frac{\displaystyle \int _{0}^{z_{min}}a\unicode[STIX]{x1D6E4}\sqrt{1+\left(\frac{\unicode[STIX]{x2202}a}{\unicode[STIX]{x2202}z}\right)^{2}}\,\text{d}z}{\displaystyle \int _{0}^{\unicode[STIX]{x03C0}/k_{m}}a\unicode[STIX]{x1D6E4}\sqrt{1+\left(\frac{\unicode[STIX]{x2202}a}{\unicode[STIX]{x2202}z}\right)^{2}}\,\text{d}z}.\end{eqnarray}$$ We point out that, since 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ are always obtained when 
$a_{min}$ is within the range 
$a_{min}\sim 10^{-4}$–
$8\times 10^{-3}$, the sensitivity of these magnitudes to the exact value of 
$a_{min}$ is negligible, such that both represent very robust measures. Similarly, the corresponding break-up time 
$t_{b}$ is barely sensitive to the value of 
$a_{min}$.
In contrast with 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$, which do not depend on the initial amplitude in the limit 
$\unicode[STIX]{x1D716}\ll 1$, the break-up time is a function of the form 
$t_{b}(La,\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D716})$ such that 
$t_{b}\rightarrow \infty$ as 
$\unicode[STIX]{x1D716}\rightarrow 0$. Indeed, the break-up time can be easily estimated from linear theory through the equation 
$a_{min}(t)\sim 1-\unicode[STIX]{x1D716}\exp (\unicode[STIX]{x1D714}_{m}t)$, where 
$\unicode[STIX]{x1D714}_{m}$ is the growth rate associated with the most amplified wavenumber 
$k_{m}$ shown in figure 2(c,d), leading to the estimation 
$t_{b}\sim \ln (\unicode[STIX]{x1D716}^{-1})/\unicode[STIX]{x1D714}_{m}$. Based on the latter result, we define the nonlinear correction to the linear break-up time as
$$\begin{eqnarray}\unicode[STIX]{x0394}t_{NL}=t_{b}-\frac{\ln (\unicode[STIX]{x1D716}^{-1})}{\unicode[STIX]{x1D714}_{m}},\end{eqnarray}$$ where 
$t_{b}$ is obtained by extrapolating 
$a_{min}$ to zero using the last few computed time steps. Unlike 
$t_{b}$, 
$\unicode[STIX]{x0394}t_{NL}$ only depends on 
$La$ and 
$\unicode[STIX]{x1D6FD}$, but not on 
$\unicode[STIX]{x1D716}$, provided only that 
$\unicode[STIX]{x1D716}\ll 1$. The latter fact is demonstrated in § 3.3. Finally, we have also computed the sphericity of the satellite droplet at pinch-off as
$$\begin{eqnarray}{\mathcal{S}}=\frac{2\left(\frac{3}{4}\displaystyle \int _{0}^{z_{min}}a^{2}\,\text{d}z\right)^{2/3}}{\displaystyle \int _{0}^{z_{min}}a\sqrt{1+\left(\frac{\unicode[STIX]{x2202}a}{\unicode[STIX]{x2202}z}\right)^{2}}\,\text{d}z},\end{eqnarray}$$ which is the ratio between the surface of a sphere of the same volume as the satellite and its actual surface. The quantification of the satellite formation process will be based on the four functions 
$V_{sat}$, 
$\unicode[STIX]{x1D6F4}_{sat}$, 
$\unicode[STIX]{x0394}t_{NL}$ and 
${\mathcal{S}}$, extracted from the numerical simulations. These four functions only depend on 
$La$ and 
$\unicode[STIX]{x1D6FD}$ when 
$\unicode[STIX]{x1D716}$ is sufficiently small, as is demonstrated in § 3.3. Thus, the main results reported herein have been computed in the limit 
$\unicode[STIX]{x1D716}\rightarrow 0$.
The structure of the (
$La$, 
$\unicode[STIX]{x1D6FD}$) parameter plane is summarised in figures 3 and 4 in terms of the satellite formation process. In particular, figure 3 depicts the most salient features of the parameter plane, and figure 4 displays several satellite shapes at the last computed numerical step just prior to pinch-off, whose associated values of 
$La$ and 
$\unicode[STIX]{x1D6FD}$ are indicated with circles in figure 3. The most important feature of the parameter plane is the solid line shown in figure 3, which represents a discontinuous transition that takes place for a critical elasticity, 
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{c}(La)$ for 
$La<7.5$. In particular, both the satellite volume and the associated entrapped mass of surfactant experience sudden jumps from certain values 
$V_{sat}(\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D6FD}_{c}\rightarrow 0^{-})$ and 
$\unicode[STIX]{x1D6F4}_{sat}(\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D6FD}_{c}\rightarrow 0^{-})$ to larger values 
$V_{sat}(\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D6FD}_{c}\rightarrow 0^{+})$ and 
$\unicode[STIX]{x1D6F4}_{sat}(\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D6FD}_{c}\rightarrow 0^{+})$. Indeed, the inset of figure 3 shows the jumps experienced by the satellite volume, 
$\unicode[STIX]{x0394}V_{sat}(La)=V_{sat}(\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D6FD}_{c}\rightarrow 0^{+})-V_{sat}(\unicode[STIX]{x1D6FD}-\unicode[STIX]{x1D6FD}_{c}\rightarrow 0^{-})$, and by the associated entrapped mass of surfactant, 
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}_{sat}$. Both jumps and 
$\unicode[STIX]{x1D6FD}_{c}$ increase monotonically as 
$La$ decreases, and reach respective Stokes asymptotes as 
$La\rightarrow 0$, namely 
$\unicode[STIX]{x0394}V_{sat}\rightarrow 0.022$, 
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F4}_{sat}\rightarrow 0.045$ and 
$\unicode[STIX]{x1D6FD}_{c}\rightarrow 0.98$. The filled circle in figure 3 indicates the origin of the discontinuous transition, 
$(La,\unicode[STIX]{x1D6FD}_{c})=(7.5,0.55)$, at which both jumps become zero. For values of 
$La>7.5$, 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ are continuous functions of 
$La$ and 
$\unicode[STIX]{x1D6FD}$.
As shown in figure 4, for values of 
$\unicode[STIX]{x1D6FD}=0<\unicode[STIX]{x1D6FD}_{c}$ and 
$\unicode[STIX]{x1D6FD}=0.5<\unicode[STIX]{x1D6FD}_{c}$ the sequence of interface shapes at pinch-off depends continuously on 
$La$, with the trend that larger satellites are formed as 
$La$ increases, reaching the regular limit of inviscid flow as 
$La\rightarrow \infty$. For 
$\unicode[STIX]{x1D6FD}<\unicode[STIX]{x1D6FD}_{c}$ and small values of 
$La$, figure 4 reveals that the main drops are separated by very thin threads of tiny volume whose break-up behaviour has been characterised in previous studies (see e.g. Kowalewski Reference Kowalewski1996). For 
$\unicode[STIX]{x1D6FD}<\unicode[STIX]{x1D6FD}_{c}$ and intermediate values of 
$La$, the main drops are separated by a satellite centred at 
$z=0$ that is connected to the main drops by very thin threads (see e.g. the case for 
$La=1$ and 
$\unicode[STIX]{x1D6FD}=0$ in figure 4). Finally, for 
$\unicode[STIX]{x1D6FD}<\unicode[STIX]{x1D6FD}_{c}$ and large values of 
$La$, the satellite drop is directly connected to the main drops. In contrast, when 
$\unicode[STIX]{x1D6FD}=1>\unicode[STIX]{x1D6FD}_{c}$, figure 4 shows a different picture, where large satellites are formed for all values of 
$La$. These results have also been analysed quantitatively, and are discussed in detail below.
From figures 3 and 4 it is deduced that, although the physical mechanisms are different, both the liquid inertia and the interfacial elastic stress favour the formation of satellites. In particular, surface elasticity tends to form spherical-shaped satellites at pinch-off, whereas the increase of the liquid inertia generates oval-shaped satellites. In the set of shapes close to pinch-off shown in figure 4, a discontinuous transition is observed for 
$La=0.01$ and 1, as 
$\unicode[STIX]{x1D6FD}$ increases. However, for 
$La=10>7.5$ a continuous transition of the thread shape is observed as 
$\unicode[STIX]{x1D6FD}$ increases. Finally, for 
$La=100$, the upper row evidences that the influence of the elastic stress on the shape of the thread is much weaker when the value of 
$La$ is large enough. The physics underlying these transitions can be explained in terms of the coupling between the liquid inertia, the viscous stress, the surface tension and the interfacial elastic stress. The competition between these forces is discussed in § 3.3, based on the trends exhibited by the functions 
$V_{sat}$, 
$\unicode[STIX]{x1D6F4}_{sat}$, 
$\unicode[STIX]{x0394}t_{NL}$ and 
${\mathcal{S}}$, and also by analysing the temporal evolution of the interface shapes starting from small disturbances, depending on the values of 
$La$ and 
$\unicode[STIX]{x1D6FD}$.
3.3 Nonlinear dynamics of a surfactant-laden interface: satellite drop formation
To unveil the effect of liquid inertia, viscous stresses and surface elasticity on the satellite droplet formation regimes, here we present and discuss the quantitative results of the detailed numerical analysis that has been carried out in the present work.
Figure 5 shows the isocontours of 
$V_{sat}$, 
$\unicode[STIX]{x1D6F4}_{sat}$, 
$\unicode[STIX]{x0394}t_{NL}$ and 
${\mathcal{S}}$ in the 
$(La,\unicode[STIX]{x1D6FD})$ parameter plane. We first observe that, at the discontinuous transition that occurs for 
$La<7.5$, the value of 
$V_{sat}$ increases from 
$10^{-3}$–1.5 % to 2–2.3 %, whereas 
$\unicode[STIX]{x1D6F4}_{sat}$ increases from 10
$^{-3}$–1.5 % to 3.5–4.7 %. The exact value of both jumps as functions of 
$La$ can be seen in the inset of figure 3. In contrast, for 
$La>7.5$ or 
$\unicode[STIX]{x1D6FD}>\unicode[STIX]{x1D6FD}_{c}(La)$, the values of 
$V_{sat}$, 
$\unicode[STIX]{x1D6F4}_{sat}$, 
$\unicode[STIX]{x0394}t_{NL}$ and 
${\mathcal{S}}$ vary continuously.
As a first general observation, it is deduced from figure 5 that the linear theory may either underestimate or overestimate the break-up time, in a way that does not necessarily coincide with the transitions in the satellite formation process. Indeed, 
$t_{b}$ is underestimated for 
$La\gg 1$ independently of the value of 
$\unicode[STIX]{x1D6FD}$. However, for 
$La\ll 1$, 
$t_{b}$ is overestimated for 
$0.28\lesssim \unicode[STIX]{x1D6FD}\lesssim 1$, while it is underestimated outside this range. Regarding the sphericity 
${\mathcal{S}}$, figure 5 confirms the trend deduced from figure 4: the most spherical satellite shapes, with 
${\mathcal{S}}\gtrsim 0.9$, take place for 
$\unicode[STIX]{x1D6FD}\gtrsim \unicode[STIX]{x1D6FD}_{c}$ and 
$La\lesssim 10$. In contrast, the shapes become most elongated, with 
${\mathcal{S}}\lesssim 0.2$, when 
$\unicode[STIX]{x1D6FD}<\unicode[STIX]{x1D6FD}_{c}$ and 
$La\lesssim 0.1$ (grey area in figure 5d).
Figure 6. Temporal evolution of the liquid thread radius a (a–d), of the axial velocities at the free surface 
$w_{s}$ and at the axis 
$w_{a}$ (e–h) and the radial surface velocity 
$u_{s}$ (i), for 
$\unicode[STIX]{x1D716}=10^{-3}$, 
$\unicode[STIX]{x1D6FD}=0$, (a–d) 
$La=0.01$, 
$k=k_{m}=0.150$, and (e–h) 
$La=100$, 
$k=k_{m}=0.635$. The vertical lines in each last snapshot indicate the axial position 
$z_{min}$ of minimum radii 
$a_{min}$, being 
$z_{min}=1.49$ and 
$a_{min}=3.63\times 10^{-5}$ for 
$La=0.01$ and 
$z_{min}=3.12$ and 
$a_{min}=1.29\times 10^{-4}$ for 
$La=100$. (i) Zoomed region close to the neck at the instant shown in (h).
3.3.1 Analysis of the temporal evolution of clean interfaces
To present the dynamics of satellite droplet formation, we take as reference cases the two canonical temporal evolutions of clean interfaces (
$\unicode[STIX]{x1D6FD}=0$) illustrated in figure 6, for 
$La=0.01$ in figure 6(a–d), close to the Stokes limit, and for 
$La=100$ in figure 6(e–h), an almost inviscid case close to the Euler limit (as shown in § 3.3.4). Specifically, we plot snapshots at different times, indicated in the labels, of the jet radius 
$a$ (upper rows), the axial surface velocity 
$w_{s}$ (middle rows, black lines), the axial velocity at the centreline 
$w_{a}$ (middle rows, green lines) and the radial surface velocity 
$u_{s}$ (bottom rows). In both cases the initial disturbance amplitude is very small, 
$\unicode[STIX]{x1D716}=10^{-3}$, and thus the initial evolution is triggered by the Plateau–Rayleigh instability mechanism, and can be described with linearised theory. This initial stage is not shown in figure 6 for conciseness, but it can be appreciated in figure 2(a). The initial disturbance, of most amplified wavelength 
$k_{m}$, creates an axial capillary pressure gradient that induces a flow from the valley to the crest of the wave. The latter mechanism finally leads to the break-up of the liquid thread and the formation of two main drops with either a liquid thread or a satellite droplet in between.
A key feature that determines the nonlinear evolution of the destabilised thread is the fact that the axial curvature makes the capillary pressure gradient to be locally larger in the regions that connect the central part of the thread with the growing crests, as evidenced by the surface and axis velocities in the snapshot of figure 6(f). This enhanced pressure gradient drives liquid towards the crests faster in the nearby regions than in the central part, and explains the appearance of two local minima in the jet radius for large enough values of 
$La$, as can be clearly appreciated in snapshots of figure 6(f,g) for 
$La=100$. In addition, the axial position of the minimum radii 
$z_{min}$ is advected with the flow along with the maximum pressure gradient, i.e. towards higher values of 
$z$ as time advances (Ashgriz & Mashayek Reference Ashgriz and Mashayek1995; Castrejón-Pita et al. Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, Lister and Basaran2015). These two local minima become the two neck regions where pinch-off takes place, leading to the formation of an oval-shaped satellite droplet, as can be observed in the snapshot of figure 6(h). This scenario applies to cases where 
$La\gg 1$ (figure 6e–h), for which the viscous stress is negligible, and the capillary pressure gradient is entirely transferred to liquid inertia leading to a self-accelerated process.
In contrast, when 
$La\ll 1$ (figure 6a–d), the viscous dissipation inhibits the growth of higher harmonics, and larger pressure gradients are needed to overcome the viscous damping, as has already been pointed out by Ashgriz & Mashayek (Reference Ashgriz and Mashayek1995). Hence, the axial movement of the minimum radius is delayed by the viscous stress, since it weakens the capillary pressure and the concomitant liquid advection. Consequently, the central region shrinks almost uniformly until the last instants before break-up, giving rise to long and thin filaments without the formation of appreciable satellite droplets before detachment. Notice also that, for 
$La\ll 1$, the axial velocities at the centreline, 
$w_{a}$, and at the interface, 
$w_{s}$, are almost equal (green and black lines, respectively, in figures 6a–d), indicating that the radial profile of axial velocity inside the thread is nearly uniform at low Laplace numbers.
Figure 7. Values of 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ as a function of 
$\unicode[STIX]{x1D6FD}$ for (a,b) 
$La=0.01$ and (c,d) 
$La=100$, and for different values of the initial perturbation amplitude 
$\unicode[STIX]{x1D716}$ indicated in the legend. The dotted line represents the results of using the linear equation of state 
$\unicode[STIX]{x1D70E}=1-\unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D6E4}-1)$ (Dravid et al. Reference Dravid, Songsermpong, Xue, Corvalan and Sojka2006). The insets show 
$V_{sat}$ as a function of 
$\unicode[STIX]{x1D716}$ in logarithmic scale for (a) 
$\unicode[STIX]{x1D6FD}=1.3$ and (c) 
$\unicode[STIX]{x1D6FD}=0.7$, demonstrating that the final stage of the liquid thread just before pinch-off becomes independent of 
$\unicode[STIX]{x1D716}$ when its value is sufficiently small. The insets in (b,d) show the nonlinear correction to the break-up time 
$\unicode[STIX]{x0394}t_{NL}$ as a function of 
$\unicode[STIX]{x1D6FD}$.
3.3.2 Analysis of the temporal evolution of surfactant-laden interfaces
To explain the different trends and transitions observed in figure 5, let us first focus on the effect of 
$\unicode[STIX]{x1D6FD}$ for the particular cases of 
$La=0.01$ and 
$La=100$. Figure 7 shows 
$V_{sat}$ (figure 7a,c), 
$\unicode[STIX]{x1D6F4}_{sat}$ (figure 7b,d) and 
$\unicode[STIX]{x0394}t_{NL}$ (insets in figure 7b,d) as functions of 
$\unicode[STIX]{x1D6FD}$ for 
$La=0.01$ in figure 7(a,b) and for 
$La=100$ in figure 7(c,d). In addition, we have computed the results for several values of 
$\unicode[STIX]{x1D716}$ indicated in the legend of figure 7(a), with the aim of clearly establishing the limit of infinitesimal disturbances. In particular, figure 7 shows that 
$V_{sat}$, 
$\unicode[STIX]{x1D6F4}_{sat}$ and 
$\unicode[STIX]{x0394}t_{NL}$ become independent of 
$\unicode[STIX]{x1D716}$ provided that 
$\unicode[STIX]{x1D716}$ is small enough, as stated before. Indeed, the insets in figure 7(a,c), which show the dependence of 
$V_{sat}$ on 
$\unicode[STIX]{x1D716}$, clearly demonstrate that the value of 
$V_{sat}$ reaches the infinitesimal-disturbance plateau when 
$\unicode[STIX]{x1D716}\lesssim 0.1$. Figure 7 also displays the results obtained with the linear equation of state 
$\unicode[STIX]{x1D70E}=1-\unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D6E4}-1)$ (dotted line), instead of the nonlinear one (2.10). It is important to note that the use of the linear equation of state leads to substantial quantitative differences with respect to the nonlinear one (2.10). In particular, the linear equation underestimates the values of 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ considerably. We note also that we compared our numerical results using the linear equation of state with those reported by Dravid et al. (Reference Dravid, Songsermpong, Xue, Corvalan and Sojka2006), finding very good agreement. However, their results were calculated for wavenumbers 
$k\neq k_{m}$, and the satellite droplet was measured by those authors by means of the thread radius at 
$z=0$ close to pinch-off, instead of using either 
$V_{sat}$ or 
$R_{sat}$.
Figure 7(a,b) shows the discontinuous transition in 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ that occurs when 
$\unicode[STIX]{x1D6FD}$ is increased above the critical value 
$\unicode[STIX]{x1D6FD}_{c}(La=0.01)=0.978\pm 0.0003$. For 
$\unicode[STIX]{x1D6FD}>\unicode[STIX]{x1D6FD}_{c}$ (
$La=0.01$) a satellite drop centred at 
$z=0$ is formed, trapping approximately 
$2.1\,\%$ of the total volume of liquid and 
$4.5\,\%$ of the total mass of surfactant.
Figure 8. Temporal evolution of the liquid thread radius 
$a$ (first row), surfactant concentration 
$\unicode[STIX]{x1D6E4}$ (second row, black lines), surface tension 
$\unicode[STIX]{x1D70E}$ (second row, blue lines), axial velocity at the interface 
$w_{s}$ (third row, black lines) and at the centreline 
$w_{a}$ (third row, green lines) and radial surface velocity 
$u_{s}$ (fourth row), for 
$La=0.01$, 
$\unicode[STIX]{x1D716}=10^{-3}$ and 
$\unicode[STIX]{x1D6FD}=0.960<\unicode[STIX]{x1D6FD}_{c}(La=0.01)$, with 
$k=k_{m}=0.508$. The vertical line in the last snapshot of 
$a$ indicates the position of 
$z_{min}$. Here 
$z_{min}=0.33$ and 
$a_{min}=7.29\times 10^{-4}$.
3.3.3 Physical explanation of the discontinuous transition
To explain the abrupt transition induced by the presence of surfactants, figures 8 and 9 show the temporal evolution of the liquid thread for 
$La=0.01$ and two different values of 
$\unicode[STIX]{x1D6FD}$, namely 
$\unicode[STIX]{x1D6FD}=0.960<\unicode[STIX]{x1D6FD}_{c}(La=0.01)$, with 
$k=k_{m}=0.508$, and 
$\unicode[STIX]{x1D6FD}=0.979>\unicode[STIX]{x1D6FD}_{c}(La=0.01)$, with 
$k=k_{m}=0.512$, respectively. In both cases, we have computed the thread radius 
$a$ (first row), the surfactant concentration 
$\unicode[STIX]{x1D6E4}$ together with the surface tension 
$\unicode[STIX]{x1D70E}$ (second row), the axial velocity at the interface 
$w_{s}$ and at the centreline 
$w_{a}$ (third row) and the radial surface velocity 
$u_{s}$ (fourth row). Time is indicated in the labels.
The presence of surfactants introduces two main effects. The advection of surfactant molecules outside the central region of the thread increases the local surface tension in this region, as can be observed in figures 8(a) and 9(a). This surfactant depletion generates two opposed effects. First, the axial capillary pressure gradient is enhanced, since the value of 
$\unicode[STIX]{x1D70E}$ becomes larger in the central region, where 
$\unicode[STIX]{x1D6E4}$ is smaller, while 
$\unicode[STIX]{x1D70E}$ becomes smaller away from the centre, where 
$\unicode[STIX]{x1D6E4}$ is larger. Second, there is a stabilising effect induced by the elastic or Marangoni stress, which competes with the destabilising Plateau–Rayleigh mechanism enhanced by the first effect. Actually, the gradient of 
$\unicode[STIX]{x1D70E}$ generates a tangential stress at the interface directed towards increasing values of 
$\unicode[STIX]{x1D70E}$, which opposes the drainage flow and tends to replenish the central zone with surfactant.
In the case of 
$\unicode[STIX]{x1D6FD}<\unicode[STIX]{x1D6FD}_{c}(La=0.01)$, figure 8(b,c) shows that the Marangoni stress reduces the axial surface velocity, 
$w_{s}$, compared with the centreline velocity, 
$w_{a}$, the difference between both velocities being larger in the region where 
$\unicode[STIX]{x1D735}_{s}\unicode[STIX]{x1D70E}$ is higher. As the fluid is drained from the centre for increasing times, 
$\unicode[STIX]{x1D735}_{s}\unicode[STIX]{x1D70E}$ becomes larger. When 
$\unicode[STIX]{x1D6FD}<\unicode[STIX]{x1D6FD}_{c}(La=0.01)$ the capillary pressure gradient is able to remove most of the liquid from the centre. Eventually, close to pinch-off, inertia becomes important and the flow is reverted close to 
$z=0.33$, so that the rate of thinning increases in this region and 
$z_{min}$ moves towards the latter axial position where the liquid thread finally detaches forming a tiny satellite droplet with 
$V_{sat}<10^{-5}$, as evidenced by figure 8(d). Note that, during thread evolution, two bulges connecting the central and outer regions grow due to the reduction of the surface velocity, and are finally connected by a thin liquid thread close to pinch-off.
When 
$\unicode[STIX]{x1D6FD}>\unicode[STIX]{x1D6FD}_{c}(La=0.01)$ the foregoing explanation still holds, but the elastic stress is large enough to revert the flow near the interface at early times far from break-up, as shown in figure 9(a). The associated stagnation point diffuses radially inwards, and leads to a counterflow separating a region where liquid flows towards the centre and induces the formation of a satellite from another region where the incipient main drop is fed with liquid. Consequently, the thread detaches in between these two regions. If 
$\unicode[STIX]{x1D6FD}$ increases further, the break-up time increases and the flow reversal occurs at earlier stages, so that 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ increase monotonically, as shown in figure 7(a,b).
Figure 9. Same as figure 8 but for 
$\unicode[STIX]{x1D6FD}=0.979>\unicode[STIX]{x1D6FD}_{c}(La=0.01)$, with 
$k=k_{m}=0.512$. The insets are zooms showing the normalised velocity vector field and isocontours of the pressure field. Here 
$z_{min}=1.93$ and 
$a_{min}=2.5\times 10^{-4}$.
When 
$La=100$, figure 7(c,d) shows that the effect of surface elasticity is much weaker in the case of dominant inertia, as was anticipated both in figure 5 and also by the shapes shown in the upper row of figure 4. The small influence of insoluble surfactants in the inviscid limit, 
$La\gg 1$, had been already noted in the linear stability analyses of Whitaker (Reference Whitaker1976), Hansen, Peters & Meijer (Reference Hansen, Peters and Meijer1999) and Timmermans & Lister (Reference Timmermans and Lister2002). Indeed, the effect of Marangoni stresses is confined to a thin boundary layer at the free surface, where the viscous stress rapidly restores any imbalance of 
$\unicode[STIX]{x1D70E}$, and which does not have any influence in the bulk liquid motion. Consequently, for 
$La=100$, the satellite volume 
$V_{sat}$ varies only slightly with respect to the value of a clean liquid thread, 
$V_{sat}(\unicode[STIX]{x1D6FD}=0,La=100)\simeq 0.03$, with a minimum at 
$\unicode[STIX]{x1D6FD}\simeq 0.203$, whereas 
$\unicode[STIX]{x1D6F4}_{sat}$ increases monotonically as 
$\unicode[STIX]{x1D6FD}$ increases. To explain this result, figures 10 and 11 show two sets of snapshots of 
$a$, 
$\unicode[STIX]{x1D6E4}$, 
$\unicode[STIX]{x1D70E}$, 
$w_{s}$, 
$w_{a}$ and 
$u_{s}$ for 
$\unicode[STIX]{x1D6FD}=0.203$, at which 
$V_{sat}$ is minimum, and for 
$\unicode[STIX]{x1D6FD}=1$, respectively.
In the weak-elastic limit, 
$\unicode[STIX]{x1D6FD}\rightarrow 0$, a satellite droplet with volume 
$V_{sat}\simeq 3\%$ is formed at pinch-off, as already shown in figures 5 and 7(c). The satellite volume decreases as 
$\unicode[STIX]{x1D6FD}$ increases in the range 
$0<\unicode[STIX]{x1D6FD}\lesssim 0.203$. Indeed, when 
$\unicode[STIX]{x1D6FD}$ increases, the Marangoni rigidification of the interface slows down the pinch-off process by decreasing the interfacial velocities, as evidenced by the time evolution of 
$w_{s}$, 
$w_{a}$ and 
$u_{s}$ in figure 10 with respect to figure 6(e–h). The latter behaviour, together with the fact that the pressure gradient is locally enhanced due to the variations of 
$\unicode[STIX]{x1D70E}$, explain why a larger volume is drained out of the satellite droplet compared to the case of a clean interface. However, the Marangoni stress that opposes the drainage flow away from the centre reduces the advection of surfactant towards the main drops, and thus the value of 
$\unicode[STIX]{x1D6F4}_{sat}$ increases, as shown in figure 7(d). The snapshot in figure 10(c) shows that the flow is reversed near the neck region, as happens for a clean interface (see e.g. figure 6h,i). However, in the elastic regime the flow reversal takes place earlier than in the clean interface limit. This behaviour at high values of 
$La$ and low values of 
$\unicode[STIX]{x1D6FD}$ was previously noticed by Kamat et al. (Reference Kamat, Wagoner, Thete and Basaran2018), who showed that the stagnation point occurs at earlier stages in surfactant-laden interfaces compared with clean interfaces, due to the strong Marangoni stress in the neck region.
Figure 10. Same as figure 8 but for 
$La=100$ and 
$\unicode[STIX]{x1D6FD}=0.203$, with 
$k=k_{m}=0.625$. Here 
$z_{min}=3.09$ and 
$a_{min}=7.89\times 10^{-4}$.
A representative case of 
$La=100$ and 
$\unicode[STIX]{x1D6FD}>0.203$ is shown in the snapshots of figure 11 for 
$\unicode[STIX]{x1D6FD}=1$. The main change with respect to the preceding case is the fact that for 
$\unicode[STIX]{x1D6FD}=1$ the Marangoni stress is strong enough to revert the surface flow at earlier stages, as shown in figure 11(b,d). Therefore, the stagnation point appears earlier than in the case of figure 10, and diffuses almost instantaneously in the radial direction, leading to a satellite droplet with larger values of the normalised volume and of the surfactant mass. It can thus be deduced that the minimum value of 
$V_{sat}$ displayed in figure 7(c) appears due to a competition between the two aforementioned opposite effects induced by the presence of surfactants.
For 
$La<7.5$, the two effects described previously coexist when 
$\unicode[STIX]{x1D6FD}$ is increased, as shown by the isocontours of 
$V_{sat}$ in figure 5. For instance, when 
$La=1$, 
$V_{sat}$ first decreases as 
$\unicode[STIX]{x1D6FD}$ increases, and when the elastic stress is strong enough, the flow is reversed and the discontinuous transition occurs. Note that, in the latter case, inertia is important since 
$La$ is of order unity, and a small but finite satellite droplet exists in the clean limit, 
$\unicode[STIX]{x1D6FD}\rightarrow 0$ (see e.g. the second row of figure 4), where 
$V_{sat}=0.394\,\%$ (a value significantly larger than in the limit 
$La\ll 1$, as shown in the isocontours of figure 5). Hence, the main difference with respect to the limit 
$La\gg 1$ is that in this case, since 
$V_{sat}(\unicode[STIX]{x1D6FD}\rightarrow 0)$ is small, the increase of 
$\unicode[STIX]{x1D6FD}$ reduces the satellite volume and may even make it negligible. For 
$La<7.5$, 
$\unicode[STIX]{x1D6F4}_{sat}$ also decreases monotonically together with 
$V_{sat}$ when 
$\unicode[STIX]{x1D6FD}<\unicode[STIX]{x1D6FD}_{c}$, which can be explained by the fact that 
$V_{sat}$ is already small when 
$\unicode[STIX]{x1D6FD}=0$, so that 
$\unicode[STIX]{x1D6F4}_{sat}$ necessarily decreases when 
$\unicode[STIX]{x1D6FD}$ is increased.
Let us recall at this point that the critical elasticity, 
$\unicode[STIX]{x1D6FD}_{c}(La)$, decreases as 
$La$ increases within the range 
$0<La<7.5$, as shown in figures 3 and 5. The reason for the latter trend is the fact that the advection of surfactant away from the central region is enhanced by the liquid inertia, so that 
$\unicode[STIX]{x1D735}_{s}\unicode[STIX]{x1D70E}$ also increases, and thus the value of 
$\unicode[STIX]{x1D6FD}$ for which the elastic stress reverts the flow is smaller. Furthermore, the value of 
$V_{sat}(La,\unicode[STIX]{x1D6FD}\rightarrow 0)$ increases as 
$La$ becomes larger, and therefore the jumps experienced by 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ at the discontinuous transition, 
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{c}$, decrease, as deduced from the inset of figure 3. Finally, for 
$La>7.5$, the discontinuous transition disappears.
Figure 11. Same as figure 10 but for 
$\unicode[STIX]{x1D6FD}=1$ with 
$k=k_{m}=0.647$. Here 
$z_{min}=2.94$ and 
$a_{min}=5.14\times 10^{-4}$.
Figure 12. (a) Normalised satellite volume 
$V_{sat}$ and (b) normalised mass of surfactant trapped at its interface 
$\unicode[STIX]{x1D6F4}_{sat}$ as a function of the Laplace number 
$La$ on log–log scales for different values of 
$\unicode[STIX]{x1D6FD}$ indicated in the legend. The inset shows the dependence of the maximum amplification wavenumber 
$k_{m}$ with respect to 
$La$ on a log–log scale. The circle with error bars corresponds to the experiments of the natural break-up of a liquid jet of water performed by Rutland & Jameson (Reference Rutland and Jameson1970).
3.3.4 Scaling laws for $V_{sat}$ and $\unicode[STIX]{x1D6F4}_{sat}$ as functions of $La$
Figure 12 shows 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ as functions of 
$La$ for different values of 
$\unicode[STIX]{x1D6FD}$ indicated in the legend. The circle with error bars corresponds to the experiment of Rutland & Jameson (Reference Rutland and Jameson1971) of the natural break-up of a liquid jet of clean water, which is in close agreement with our numerical result for 
$\unicode[STIX]{x1D6FD}=0$. The inset displays the most unstable wavenumber, 
$k_{m}$, as a function of 
$La$, showing the inviscid plateau 
$k_{m}\simeq 0.697$ for 
$La\gg 1$ (Rayleigh Reference Rayleigh1878), as well as the power-law dependence for small values of 
$La$. The latter power law can be deduced from the long-wave approximation of the dispersion relation (3.2) or, equivalently, from the leading-order one-dimensional model deduced by Eggers & Dupont (Reference Eggers and Dupont1994) and García & Castellanos (Reference García and Castellanos1994). In the clean case, 
$\unicode[STIX]{x1D6FD}=0$, the leading-order one-dimensional results are 
$k_{m}\sim (2+3\sqrt{2}La^{-1/2})^{-1/2}$ and 
$\unicode[STIX]{x1D714}_{m}\sim (2\sqrt{2}+6La^{-1/2})^{-1}$ (Eggers & Villermaux Reference Eggers and Villermaux2008). The latter long-wave result provides very accurate results in the whole range of 
$La$, since 
$k\in (0,1)$ accomplishes the slenderness assumption. In the inviscid limit, 
$La\rightarrow \infty$, both 
$\unicode[STIX]{x1D714}_{m}$ and 
$k_{m}$ are slightly overestimated by the one-dimensional model, namely 
$\unicode[STIX]{x1D714}_{m}\rightarrow 2^{-3/2}$ and 
$k_{m}\rightarrow 2^{-1/2}$. However, in the Stokes limit, 
$La\ll 1$, the values of 
$\unicode[STIX]{x1D714}_{m}\rightarrow 1/6$ and 
$k_{m}=3^{-1/2}2^{-1/4}La^{1/4}$ are in excellent agreement with the exact linear theory. When 
$\unicode[STIX]{x1D6FD}>1/2$, the elastic stress regularises 
$k_{m}$ in the limit of 
$La\rightarrow 0$, as analysed in detail by Timmermans & Lister (Reference Timmermans and Lister2002) (see also the isocontours of 
$k_{m}$ in figure 2c).
In the limit of a clean interface, 
$\unicode[STIX]{x1D6FD}=0$, 
$V_{sat}$ increases monotonically with 
$La$, as previously shown in figures 4 and 5, and explained in figure 6. In particular, our numerical results reveals that the satellite volume scales as 
$V_{sat}=0.00421La^{1.64}$ when 
$La\lesssim 2$, and thus 
$V_{sat}\rightarrow 0$ and 
$\unicode[STIX]{x1D6F4}_{sat}\rightarrow 0$ as 
$La\rightarrow 0$. When 
$La$ is finite, a satellite drop is always formed, since the liquid thread always experiences a transition to the inertial–viscous regime (Eggers Reference Eggers1993; Castrejón-Pita et al. Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, Lister and Basaran2015) and thus 
$z_{min}$ moves from 
$z=0$ towards higher values when 
$t$ is close enough to 
$t_{b}$. The elongated satellite droplet formed when 
$La\ll 1$ can break up into more droplets after pinch-off as it relaxes, depending on the value of 
$La$ (Notz & Basaran Reference Notz and Basaran2004; Castrejon-Pita, Castrejon-Pita & Hutchings Reference Castrejon-Pita, Castrejon-Pita and Hutchings2012; Anthony et al. Reference Anthony, Kamat, Harris and Basaran2019; Wang et al. Reference Wang, Contò, Naz, Castrejón-Pita, Castrejón-Pita, Bailey, Wang, Feng and Sui2019), unless 
$La\rightarrow 0$ (Eggers & Fontelos Reference Eggers and Fontelos2005). Alternatively, using the expression for the equivalent radius 
$R_{sat}$ developed in § 3.2, which depends on 
$k_{m}$, and since 
$k_{m}=3^{-1/2}2^{-1/4}La^{1/4}$ within the range of 
$La$ for which 
$V_{sat}$ exhibits power-law scaling, it is deduced that 
$R_{sat}=0.34La^{0.463}$.
When 
$La\gtrsim 10$, the value of 
$V_{sat}$ reaches a plateau of about 3 %, as already discussed in the context of figures 5 and 7(c). Equivalently, since 
$k_{m}\simeq 0.697$ when 
$La\gg 1$, 
$R_{sat}\simeq 0.588$ in the inviscid limit, in excellent agreement with the experiments of Rutland & Jameson (Reference Rutland and Jameson1971) (circle with error bars in figure 12a), and also with the numerical simulations of Ashgriz & Mashayek (Reference Ashgriz and Mashayek1995). In the weak-elasticity limit, 
$\unicode[STIX]{x1D6FD}<0.05$, the behaviour of 
$\unicode[STIX]{x1D6F4}_{sat}$ is identical to that of 
$V_{sat}$, displaying the same scaling law within the same range in 
$La$, and also reaching an inviscid plateau of about 2.9 % when 
$La\gtrsim 10$. This scaling law for 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ prevails when 
$\unicode[STIX]{x1D6FD}<\unicode[STIX]{x1D6FD}_{c}$ and 
$La<La_{c}$, although the prefactor changes with 
$\unicode[STIX]{x1D6FD}$ as shown in figure 12. In particular, when 
$\unicode[STIX]{x1D6FD}$ increases the prefactor is smaller and thus 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ reach smaller values as 
$La\rightarrow 0$. This can be explained by the translation of 
$z_{min}$, which is inhibited as 
$\unicode[STIX]{x1D6FD}$ becomes higher and thus the surface stress exerted at the interface increases.
Figure 12 also shows that, when 
$La\gg 1$, the 3 % plateau reached by 
$V_{sat}$ is barely affected by 
$\unicode[STIX]{x1D6FD}$ since, as explained previously, inertia dominates and the elastic stress cannot induce any substantial change in the bulk motion (Whitaker Reference Whitaker1976; Hansen et al. Reference Hansen, Peters and Meijer1999; Timmermans & Lister Reference Timmermans and Lister2002). As inertia increases, the influence of 
$\unicode[STIX]{x1D6FD}$ on 
$V_{sat}$ becomes even weaker than in the case of 
$La=100$ displayed in figure 7(c). Although the satellite shape at pinch-off is the same because viscosity cannot balance the elastic stress and transmit it to the bulk, 
$\unicode[STIX]{x1D6F4}_{sat}$ reaches different inviscid limits as 
$\unicode[STIX]{x1D6FD}$ increases. In particular, the 2.9 % inviscid plateau reached by 
$\unicode[STIX]{x1D6F4}_{sat}$ in the weak-elasticity limit increases with 
$\unicode[STIX]{x1D6FD}$, the reason being the same as in the case of figure 7(d).
When 
$\unicode[STIX]{x1D6FD}=1>\unicode[STIX]{x1D6FD}_{c}(La)$ for arbitrary values of 
$La$, a satellite is always formed with 
$V_{sat}\gtrsim 2\,\%$ and 
$\unicode[STIX]{x1D6F4}_{sat}\gtrsim 4\,\%$, and both increase smoothly with 
$La$ as inertia becomes more important. In fact, 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ reach respective plateaus in both the Stokes and Euler limits. For 
$La\rightarrow 0$, the values of 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ are about 2.30 % and 4.67 %, respectively, whereas in the limit 
$La\rightarrow \infty$, their values are 3 % and 12 %, approximately. This trend prevails provided that 
$\unicode[STIX]{x1D6FD}>\unicode[STIX]{x1D6FD}_{c}(La\rightarrow 0)=0.978$, as shown in the isocontours of 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ in figure 5. Hence, in the elasticity-dominated regime where 
$\unicode[STIX]{x1D6FD}>\unicode[STIX]{x1D6FD}_{c}$, the effect of inertia on 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ is weaker, although the column for 
$\unicode[STIX]{x1D6FD}=1$ in figure 4 reveals that the shape of the thread at pinch-off changes substantially. As already mentioned, inertia tends to form oval-shaped satellites, which are more likely to break up in the relaxation process after pinch-off, whereas the surface elasticity tends to form spherical satellites which will not experience secondary break-up events.
4 Conclusions
In this paper we have reported an exhaustive numerical study of the unforced break-up of free axisymmetric threads of a Newtonian liquid whose interface is coated with insoluble surfactants. Our main objective was to describe and explain how the presence of these molecules affects the nonlinear dynamics of the liquid thread and the satellite drop formation regimes when the dynamics is triggered by the most dangerous initial disturbance. Under these conditions we have shown that, when the initial perturbation amplitude is sufficiently small, the flow depends on two dimensionless parameters, namely the Laplace number 
$La$ and the elasticity parameter 
$\unicode[STIX]{x1D6FD}$. Our numerical simulations have allowed us to characterise the influence of these two parameters on the satellite volume 
$V_{sat}$, the mass of surfactant trapped at its interface 
$\unicode[STIX]{x1D6F4}_{sat}$, the nonlinear correction to the linear break-up time 
$\unicode[STIX]{x0394}t_{NL}$ and the satellite sphericity 
${\mathcal{S}}$, all of them computed at times very close to break-up. It is important to emphasise that our numerical simulations do not contemplate the post break-up behaviour of the threads and satellites, including their relaxation or eventual secondary break-up events. Indeed, an accurate analysis of the dynamics beyond break-up is an important though technically challenging task, which is out of the scope of the present study. Clearly, a future task to be pursued would be to extend the present results by performing numerical simulations that are able to compute the post pinch-off dynamics to reveal the ultimate state of the unstable liquid thread.
We have found a discontinuous transition at a critical elasticity number 
$\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{c}(La)$ within the range 
$0<La<7.5$, at which 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ change abruptly. We have explained this behaviour in terms of a competition between the Plateau–Rayleigh instability mechanism and the elastic or Marangoni stresses that arise due to interfacial surface tension gradients. When 
$\unicode[STIX]{x1D6FD}$ is high enough, the elastic stress that opposes the flow induced by the capillary pressure gradient is able to revert it at the interface. Afterwards, the surface stagnation point diffuses radially inwards, and finally a net flux of liquid swells the central region forming a satellite droplet prior to pinch-off.
When 
$La<7.5$, 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ increase from a non-zero satellite droplet for 
$\unicode[STIX]{x1D6FD}<\unicode[STIX]{x1D6FD}_{c}$ to a larger value when 
$\unicode[STIX]{x1D6FD}>\unicode[STIX]{x1D6FD}_{c}$. When 
$La\lesssim 0.2$, the critical elasticity number reaches a plateau, 
$\unicode[STIX]{x1D6FD}_{c}=0.978$. Finally when 
$La>7.5$ the abrupt transition disappears. In between, 
$\unicode[STIX]{x1D6FD}_{c}$ decreases monotonically with 
$La$, since inertia enhances the gradients of surface tension.
For a clean liquid thread, 
$\unicode[STIX]{x1D6FD}\rightarrow 0$, we have provided a new scaling law for the normalised satellite volume, namely 
$V_{sat}=0.00421La^{1.64}$, which is valid for 
$La\lesssim 1$. We have shown the existence of a regular weak-elasticity limit, 
$\unicode[STIX]{x1D6FD}<0.05$, for which the latter scaling law holds, and for which the normalised mass of surfactant carried by the satellite, 
$\unicode[STIX]{x1D6F4}_{sat}$, exhibits the same scaling law as 
$V_{sat}$. In this limit, when inertia is sufficiently dominant, namely 
$La\gtrsim 10$, both 
$V_{sat}$ and 
$\unicode[STIX]{x1D6F4}_{sat}$ reach respective limits of about 3 % and 2.9 %, the value of 3 % being in close agreement with previous experiments (Rutland & Jameson Reference Rutland and Jameson1971) and numerical simulations (Ashgriz & Mashayek Reference Ashgriz and Mashayek1995).
When 
$La=100$ the 3 % inviscid plateau in 
$V_{sat}$ varies slightly with 
$\unicode[STIX]{x1D6FD}$, and displays a minimum within the range 
$0.2\lesssim \unicode[STIX]{x1D6FD}\lesssim 0.4$, whereas the 2.9 % inviscid plateau of 
$\unicode[STIX]{x1D6F4}_{sat}$ increases monotonically. The existence of this minimum has been explained by the competition between two opposed effects induced by the presence of surfactants: (I) the reduction of the surface tension 
$\unicode[STIX]{x1D70E}$ when 
$\unicode[STIX]{x1D6E4}$ increases, which enhances the capillary pressure gradient, and (II) the Marangoni stress exerted at the interface due to the gradients of 
$\unicode[STIX]{x1D70E}$. The initial decrease of 
$V_{sat}$ when 
$\unicode[STIX]{x1D6FD}$ grows is due to (I), whereas the increase above the minimum value is due to (II), which is able to revert the flow at earlier stages of the thread evolution when 
$\unicode[STIX]{x1D6FD}$ is sufficiently high. The decrease of 
$V_{sat}$ with 
$\unicode[STIX]{x1D6FD}$ also coexists with the discontinuous transition for 
$La<7.5$. Additionally, the increase of 
$\unicode[STIX]{x1D6F4}_{sat}$ when 
$La>7.5$ is explained by the reduction of the interfacial velocity due to (II), which tends to accumulate surfactant molecules at the satellite.
When 
$La\gg 1$, the effect of surface elasticity is very weak and the 3 % plateau of 
$V_{sat}$ does not vary with 
$\unicode[STIX]{x1D6FD}$, since its effect is confined to a thin Marangoni boundary layer at the interface, where viscous dissipation tends to restore a modified but constant value of 
$\unicode[STIX]{x1D70E}$. The most important effect of 
$\unicode[STIX]{x1D6FD}$ in the inviscid limit is the fact that 
$\unicode[STIX]{x1D6F4}_{sat}$ increases with 
$\unicode[STIX]{x1D6FD}$ for the reason explained in the previous paragraph.
Here, we have considered a nonlinear equation of state for 
$\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D6E4})$ that is deduced from the equilibrium thermodynamics of the interface together with the conservation of molecules in the insoluble limit. We have shown that using this nonlinear equation leads to substantial quantitative differences with respect to the use of its linearised version (Dravid et al. Reference Dravid, Songsermpong, Xue, Corvalan and Sojka2006). These differences call for a careful experimental study of the present jet flow configuration or a similar one, e.g. a cylindrical liquid bridge between two static discs whose length is above the critical one for spontaneous break-up. The latter configuration has been recently studied experimentally by Kovalchuk et al. (Reference Kovalchuk, Jenkinson, Miller and Simmons2018) for concentrations above the critical micelle concentration. An experimental campaign would also be needed to probe the validity of the insoluble approximation. In fact, it would be interesting to extend the present numerical study to the soluble case, contemplating both bulk diffusion and sorption kinetics. To that end, the bulk diffusion equation together with appropriate adsorption and desorption kinetic equations should be coupled to the equations integrated in the present work (Karapetsas & Bontozoglou Reference Karapetsas and Bontozoglou2013). We believe that the numerical techniques employed herein should be able to properly tackle the soluble problem with minor modifications.
A natural and important extension of the present work is tackling the forced jet problem, in which the wavenumber 
$k$ is not restricted to the most unstable one, and the amplitude 
$\unicode[STIX]{x1D716}$ is not necessarily small. Another feature that deserves future work is the effect of surface diffusion on the satellite drop formation regimes described herein, especially in cases where 
$La\lesssim O(1)$, for which the surface diffusion time could be of the order of the thread break-up time. Similarly, for small-scale threads, the surface shear and dilational viscosities could also play an important role (Boussinesq Reference Boussinesq1913; Scriven Reference Scriven1960; Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2018). The present numerical analysis should also be extended together with experiments to reveal the conditions under which diffusive and surface viscous effects become relevant, and how they affect the transitions described in the present work.
Acknowledgements
A.M.-C. and A.S. thank the Spanish MINECO, Subdirección General de Gestión de Ayudas a la Investigación, for its support through project DPI2015-71901-REDT, and the Spanish MCIU-Agencia Estatal de Investigación through project DPI2017-88201-C3-3-R. These research projects have been partly financed through FEDER European funds. A.M.-C. also acknowledges support from the Spanish MECD through grant FPU16/02562, and its associated programme Ayudas a la Movilidad 2017 during his stay at TIPs–ULB, Brussels. J.R.-R. and B.S. thank the FRS-FNRS for financial support, in particular under the umbrella of the Wolflow project. The authors wish to express their deep gratitude to one anonymous reviewer for making insightful comments, which led to a significant improvement of the present work.
Appendix. Validation of the numerical method
Figure 13. (a) Minimum thread radius 
$a_{min}$ as a function of the time to break-up 
$\unicode[STIX]{x1D70F}$ for two different values of the Laplace number, namely 
$La=18.9$ and 
$100$, and 
$\unicode[STIX]{x1D6FD}=0$. The dashed lines indicate the scaling laws in the different regimes, and the symbols correspond to the results extracted from the numerical simulations of Castrejón-Pita et al. (Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, Lister and Basaran2015). (b) Shape of the thread for the case 
$La=18.9$ at 
$t=138.017$, where 
$a_{min}=1.43\times 10^{-4}$ and 
$z_{min}=3.65$, and which corresponds to the star symbol in (a) for 
$\unicode[STIX]{x1D70F}=1.08\times 10^{-3}$. The zoomed region shows the micro-filament formed just prior to pinch-off. (c) Local mesh in the micro-filament region.
To demonstrate the performance of our numerical technique, in this appendix we report numerical simulations aimed at comparing our results with the well-known scaling laws of 
$a_{min}$ as a function of the time to break-up, 
$\unicode[STIX]{x1D70F}=t_{b}-t$, for two different values of the Laplace number, namely 
$La=18.9$ and 
$100$, in the case of a clean interface, 
$\unicode[STIX]{x1D6FD}=0$. Figure 13 shows 
$a_{min}$ as a function of 
$\unicode[STIX]{x1D70F}$, where the dashed lines represent the different scaling laws, and symbols have been extracted from the results of Castrejón-Pita et al. (Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, Lister and Basaran2015) for the particular case of 
$La=18.9$. Since 
$La$ is moderately high in both cases, intertial and capillary forces balance initially, providing 
$a_{min}\sim \unicode[STIX]{x1D70F}^{2/3}$ (Keller & Miksis Reference Keller and Miksis1983; Day, Hinch & Lister Reference Day, Hinch and Lister1998; Eggers & Fontelos Reference Eggers and Fontelos2015), a regime usually referred to as the inertial (I) regime. However, as the thread thins, viscous forces come into play, as shown numerically and experimentally by Castrejón-Pita et al. (Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, Lister and Basaran2015), leading to the linear behaviour 
$a_{min}=0.0709\unicode[STIX]{x1D70F}$ (Papageorgiou Reference Papageorgiou1995), which is known as the viscous (V) regime. Finally, when 
$a_{min}$ is sufficiently small, inertial, capillary and viscous forces balance, leading to what is usually known as the inertial–viscous (IV) regime, in which 
$a_{min}=0.0304\unicode[STIX]{x1D70F}$ (Eggers Reference Eggers1993). As revealed by figure 13(a), our numerical method is in excellent agreement with these scaling laws close to pinch-off, and with the numerical computations of Castrejón-Pita et al. (Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, Lister and Basaran2015), thereby validating our numerical framework. Finally, although not shown here for conciseness, we have checked that the unphysical singularity of the equation of state (2.10) as 
$\unicode[STIX]{x1D6E4}\rightarrow 0$ leads to a spurious deviation from the asymptotic IV regime, which precludes its use in correctly predicting the smallest scales prior to pinch-off for 
$\unicode[STIX]{x1D6FD}\neq 0$. To that end, a different equation of state that provides the clean-interface constant value of 
$\unicode[STIX]{x1D70E}$ as 
$\unicode[STIX]{x1D6E4}\rightarrow 0$ must be used (McGough & Basaran Reference McGough and Basaran2006; Kamat et al. Reference Kamat, Wagoner, Thete and Basaran2018).
