1. Introduction
1.1. Surfactants and surface viscosity
Complex interfaces are those whose mechanics cannot be described solely in terms of the interfacial tension. Capsules, vesicles, polymer blends and lipid bilayers are examples of fluid particles delimited by those interfaces. These fluidic entities arise both in numerous natural processes and in multiple biological and industrial applications. Therefore, understanding the factors involved in the dynamics of complex interfaces is of paramount importance at the fundamental level and at the practical one.
Surfactants produce many beneficial effects. For instance, they maintain desired wetting conditions and stabilize emulsions and foams by hindering the coalescence of droplets and bubbles. The major mechanical consequence of the presence of surface-active molecules is probably the local reduction of the capillary pressure (the so-called soluto-capillarity effect). However, when a surfactant monolayer is adsorbed onto an interface, both Marangoni and viscous surface stresses may play a significant role too. Marangoni stress arises due to the interfacial tension (surfactant concentration) gradient, while viscous surface stress is associated with the variation of the surface velocity. While Marangoni stress tends to eliminate inhomogeneities of surfactant concentration, viscous surface stress opposes surface velocity gradient.
The viscous surface stress obeys different constitutive relationships depending on the surfactant molecule nature (Fuller & Vermant Reference Fuller and Vermant2012). In fact, adsorbed surfactant monolayers at fluid surfaces usually exhibit rheological properties (Choi et al. Reference Choi, Steltenkamp, Zasadzinski and Squires2011; Langevin Reference Langevin2014). For a Newtonian interface (Scriven Reference Scriven1960; Langevin Reference Langevin2014), the viscous surface stress can be calculated from the Boussinesq–Scriven constitutive equation in terms of the shear and dilatational surface viscosities, which depend on the surfactant surface concentration (Manikantan & Squires Reference Manikantan and Squires2017; Luo, Shang & Bai Reference Luo, Shang and Bai2019).
Surface viscosity is not frequently accounted for in microfluidics, probably due to the considerable uncertainty about the values taken by this property (Zell et al. Reference Zell, Nowbahar, Mansard, Leal, Deshmukh, Mecca, Tucker and Squires2014; Ponce-Torres et al. Reference Ponce-Torres, Rubio, Herrada, Eggers and Montanero2020). However, it can significantly affect the dynamics of interfaces even for values much smaller than the bulk ones. For instance, it is believed that shear surface viscosity can stabilize foams and emulsions by increasing the drainage time during the coalescence of two bubbles/droplets (Fischer & Erni Reference Fischer and Erni2007; Ozan & Jakobsen Reference Ozan and Jakobsen2019). This effect is similar to that produced by Marangoni stresses, which may mask the role played by surface viscosity in many experiments. Therefore, elucidating the influence of surface viscosity on the dynamics of surfactant-covered interfaces is of great relevance from both fundamental and practical points of view. It is worth mentioning that surface viscosity may become relevant because of its contribution to the balance of fluid momentum and its effect on the surfactant distribution over the interface, which determines the capillary pressure profile and Marangoni stress (Ponce-Torres et al. Reference Ponce-Torres, Montanero, Herrada, Vega and Vega2017, Reference Ponce-Torres, Rubio, Herrada, Eggers and Montanero2020).
1.2. Surfactant-free droplets in linear flows
The deformation and stability of a small drop immersed in a viscous flow is a paradigmatic problem that can be used to evaluate the role of different factors involved in interfacial dynamics. If the droplet size is much smaller than the scale of variation of the imposed flow, then the droplet deforms and bursts due to the local velocity gradient at the droplet location, thus the outer flow can be regarded as linear. Hyperbolic, simple shear and extensional (straining) flows are examples of this type of fluid motion. While a simple shear flow can be easily produced experimentally (Bentley & Leal Reference Bentley and Leal1986), a uniaxial extensional flow facilitates the theoretical analysis owing to the existence of a symmetry axis (Rallison Reference Rallison1984; Stone Reference Stone1994). Inertia is almost always neglected in both the droplet and the outer bath (creeping or Stokes flow). In many cases, the inner phase is a bubble much less viscous than the bath. Taylor was one of the first to analyse the dynamics of a clean droplet in a linear flow both theoretically (Taylor Reference Taylor1932, Reference Taylor1964) and experimentally (Taylor Reference Taylor1934), showing how the droplet can develop conical ends before its breakup. Subsequent studies have examined the deformation and instability of bubbles and drops in both shear and extensional flows (see, for example, Rallison Reference Rallison1984; Stone Reference Stone1994).
The deformation and breakup of surfactant-free bubbles and drops in an axisymmetric extensional flow under different approximations have been studied in the seminal works of Taylor (Reference Taylor1964), Buckmaster (Reference Buckmaster1972), Acrivos and collaborators (Barthls-Biesel & Acrivos Reference Barthls-Biesel and Acrivos1973; Youngren & Acrivos Reference Youngren and Acrivos1976; Acrivos & Lo Reference Acrivos and Lo1978; Rallison & Acrivos Reference Rallison and Acrivos1978; Brady & Acrivos Reference Brady and Acrivos1982), Hinch (Reference Hinch1980), Sherwood (Reference Sherwood1984), Li, Barthes-Biesel & Helmy (Reference Li, Barthes-Biesel and Helmy1988) and Stone & Leal (Reference Stone and Leal1989), and, more recently, by Howell & Siegel (Reference Howell and Siegel2004) and Favelukis (Reference Favelukis2016). Both analytical slender-body solutions and numerical simulations show how the droplet deformation increases with the capillary number (the strain rate in terms of the inverse of the visco-capillary time) until this parameter reaches a critical value beyond which the steady flow ceases to be stable, and the droplet breaks up (Stone Reference Stone1994). Courrech du Pont & Eggers (Reference Courrech du Pont and Eggers2020) and Eggers (Reference Eggers2021) have recently shown that the drop tip always remains rounded, but its curvature increases exponentially with the square of the flow strength. The shape of the tip is described by a similarity solution, which is independent of the outer flow and system geometry, in agreement with earlier results of Eggers & Courrech du Pont (Reference Eggers and Courrech du Pont2009).
The stability analysis allows one to determine which solutions are stable (Taylor Reference Taylor1964; Acrivos & Lo Reference Acrivos and Lo1978; Hinch Reference Hinch1980) and the conditions for the onset of bursting or fluid ejection (Barthls-Biesel & Acrivos Reference Barthls-Biesel and Acrivos1973). Rallison & Acrivos (Reference Rallison and Acrivos1978) considered the creeping motion condition in both the droplet and the outer bath, and showed that steady shapes could be obtained only for capillary numbers below a certain threshold. Hinch & Acrivos (Reference Hinch and Acrivos1979) described the mechanism for the instability: the increased inner pressure pushes the drop ends outwards, increasing the drop length. Due to the volume conservation, the drop has to become narrower, which means that the interior lubrication pressure has to increase even more, thus promoting instability. After instability, the drop can either be torn apart and break up or develop into a quasi-stationary state in which fluid is drained gradually through a thin thread (tip streaming).
The evolution of the shape of a slender inviscid drop can be calculated approximately by using polynomials with time-dependent coefficients (Hinch Reference Hinch1980). Brady & Acrivos (Reference Brady and Acrivos1982) considered the droplet inertia and concluded that its stabilizing effect is so weak that it can be ignored for all practical purposes. The governing equations have been extended to analyse non-axisymmetric flows resulting from the combination of two-dimensional and extensional flows (Howell & Siegel Reference Howell and Siegel2004). The end-pinching breakup mechanism has also been numerically studied in the context of droplets deformed by an extensional flow (Stone & Leal Reference Stone and Leal1989).
One of the most interesting phenomena arising when a droplet is submerged in an extensional flow is tip streaming (Anna Reference Anna2016; Montanero & Gañán-Calvo Reference Montanero and Gañán-Calvo2020), i.e. the ejection of a very thin fluid thread from the tip of the stretched droplet. Zhang (Reference Zhang2004) claimed that the steady recirculating stream arising in a droplet attached to a capillary submerged in an extensional flow evolves toward tip streaming when the capillary number exceeds a critical value. This seems to indicate that steady tip streaming can be obtained by purely hydrodynamic means not only in confined geometries (Suryo & Basaran Reference Suryo and Basaran2006; Gañán-Calvo et al. Reference Gañán-Calvo, González-Prieto, Riesco-Chueca, Herrada and Flores- Mosquera2007; Evangelio, Campo-Cortés & Gordillo Reference Evangelio, Campo-Cortés and Gordillo2016) but also in the presence of an unbounded outer flow. However, these results were obtained with uncontrolled approximations, which should be checked compared to full numerical simulations.
1.3. Surfactant-covered droplets in linear flows
The presence of an insoluble surfactant monolayer significantly affects the deformation and stability of droplets subject to both non-axisymmetric shear flows (Li & Pozrikidis Reference Li and Pozrikidis1997; Bazhlekov, Anderson & Meijer Reference Bazhlekov, Anderson and Meijer2006; Lee & Pozrikidis Reference Lee and Pozrikidis2006; Feigl et al. Reference Feigl, Megias-Alguacil, Fischer and Windhab2007; Mandal, Das & Chakraborty Reference Mandal, Das and Chakraborty2017) and axisymmetric extensional flows (Stone & Leal Reference Stone and Leal1990; Pawar & Stebe Reference Pawar and Stebe1996; Eggleton, Pawar & Stebe Reference Eggleton, Pawar and Stebe1999; Eggleton, Tsai & Stebe Reference Eggleton, Tsai and Stebe2001; Booty & Siegel Reference Booty and Siegel2005; Feigl et al. Reference Feigl, Megias-Alguacil, Fischer and Windhab2007; Vlahovska, Lawzdziewicz & Loewenberg Reference Vlahovska, Lawzdziewicz and Loewenberg2009; Liu et al. Reference Liu, Ba, Wu, Li, Xi and Zhang2018). In the two cases, the surfactant is swept towards the droplet end, where the interfacial tension decreases. This convective effect competes with depletion of surfactant due to interfacial dilatation in that region. At high capillary numbers, the first effect becomes more important than the last one, and the local curvature at the apex increases for the capillary pressure to balance normal stresses. As a consequence, the deformation in the presence of surfactant is larger than that of a drop with a constant interfacial tension equal to the initial equilibrium value (Stone & Leal Reference Stone and Leal1990; Li & Pozrikidis Reference Li and Pozrikidis1997; Feigl et al. Reference Feigl, Megias-Alguacil, Fischer and Windhab2007), and the critical capillary number decreases. The interface is immobilized by the surfactant monolayer (Lee & Pozrikidis Reference Lee and Pozrikidis2006; Bazhlekov et al. Reference Bazhlekov, Anderson and Meijer2006) because Marangoni convection counteracts the external flow (Milliken, Stone & Leal Reference Milliken, Stone and Leal1993; Vlahovska et al. Reference Vlahovska, Lawzdziewicz and Loewenberg2009). Liu et al. (Reference Liu, Ba, Wu, Li, Xi and Zhang2018) conducted numerical simulations to analyse the droplet deformation in the presence of a three-dimensional shear flow. For surfactant-laden droplets, the critical capillary number decreases as the Reynolds number increases, as occurs with clean droplets.
In his pioneering work, De Bruijn (Reference De Bruijn1993) described the tip streaming occurring in a surfactant-laden droplet submerged in a simple shear flow. In this case, the reduction of the interfacial tension in the drop pole resulted in the ejection of a fluid thread much smaller than the drop size. This was one of the first mechanisms used to produce tip streaming in both droplets (De Bruijn Reference De Bruijn1993; Eggleton et al. Reference Eggleton, Tsai and Stebe2001) and bubbles (Booty & Siegel Reference Booty and Siegel2005) in a clear and reproducible way. De Bruijn (Reference De Bruijn1993) concluded that surfactants trigger tip streaming, which then disappears as surfactants are convected away from the tip. In fact, there is a widespread belief that surfactant is a necessary element to produce tip streaming in unbounded extensional flows. Tip streaming has been described from simulations of such systems (Eggleton et al. Reference Eggleton, Tsai and Stebe2001; Booty & Siegel Reference Booty and Siegel2005; Bazhlekov et al. Reference Bazhlekov, Anderson and Meijer2006; Wang, Siegel & Booty Reference Wang, Siegel and Booty2014).
The deformation and stability of a droplet covered with an insoluble surfactant in an extensional flow have been studied theoretically on several occasions. Milliken et al. (Reference Milliken, Stone and Leal1993) showed that Marangoni stresses cause the drop to behave as if it were more viscous, and the surfactant monolayer facilitates the formation of pointed ends over the drop stretching. Pawar & Stebe (Reference Pawar and Stebe1996) examined the effect of both surface saturation and non-ideal interactions among the surfactant molecules. When surface diffusion is negligible, the interface can be essentially split into mobile regions near the drop equator and motionless caps near the drop poles (Eggleton et al. Reference Eggleton, Pawar and Stebe1999). The transfer of soluble surfactant molecules between the interface and the bulk reduces the surfactant effects mentioned above (Milliken & Leal Reference Milliken and Leal1994). All the above-mentioned studies were conducted for zero Reynolds number. Liu et al. (Reference Liu, Ba, Wu, Li, Xi and Zhang2018) analysed confinement and inertia effects. The surfactant monolayer destabilizes the droplet, i.e. it reduces the critical value of the capillary number for all the confinements analysed. On the contrary, that value increases with the Reynolds number for both clean and surfactant-laden interfaces.
The simulations of Eggleton et al. (Reference Eggleton, Tsai and Stebe2001) showed how the reduction of the interfacial tension in the droplet apex due to the surfactant accumulation results in tip streaming, as occurs in a simple shear flow (De Bruijn Reference De Bruijn1993). The size of the ejected thread increases with the initial surfactant surface concentration. Wang et al. (Reference Wang, Siegel and Booty2014) studied the tip streaming numerically in a droplet covered with a soluble surfactant, placed in a uniaxial extension flow. They showed that the adsorption of a small amount of surfactant onto the interface produces tip streaming. The size of the emitted droplet decreases as the Biot number (the ratio of the flow time scale to that for surfactant desorption) increases. Wrobel et al. (Reference Wrobel, Booty, Siegel and Wang2018) extended the analysis by combining the extensional flow at infinity with flow focusing provided by two annular baffles. Eggleton & Stebe (Reference Eggleton and Stebe1998) considered the adsorption-desorption-limited case.
1.4. Surfactant-covered droplets in linear flows: the surface viscosity
The influence of interfacial rheology on the dynamics of drops submerged in linear flows has been considered in several works. Flumerfelt (Reference Flumerfelt1980) conducted a first-order perturbation analysis of the deformation and orientation of drops in shear and extensional flows involving the shear and dilatational surface viscosities. He showed that these dynamical interfacial properties play a critical role in the droplet dynamics and may explain the discrepancies between experimental droplet deformations and theoretical predictions. Numerical results have shown that surface viscosity can suppress the interfacial motion and reduce the magnitude of the deformation of the drop in a simple shear flow (Pozrikidis Reference Pozrikidis1994). Gounley et al. (Reference Gounley, Boedec, Jaeger and Leonetti2016) found that shear (dilatational) surface viscosity stabilizes (destabilizes) the drops under shear flow. Both shear and dilatational surface viscosities reduce the degree of non-uniformity of the surfactant concentration over the drop surface by inhibiting local convection and dilution of surfactant, respectively (Luo et al. Reference Luo, Shang and Bai2019). Narsimhan (Reference Narsimhan2019) developed a perturbation theory to describe the behaviour of droplets in both shear and extensional flows. Similar to what occurs in simple shear flow, the shear (dilatational) surface viscosity stabilizes (destabilizes) the drops under uniaxial extensional flow (Singh & Narsimhan Reference Singh and Narsimhan2020). It should be noted that these studies, which account for viscous surface stresses, assume uniform interfacial tension, and therefore they neglect both soluto-capillarity and Marangoni convection.
1.5. The goal of this paper
A coherent picture of both the instability of a drop, and the subsequent nonlinear evolution, is still missing for both clean and surfactant-laden drops. For that reason, we will conduct numerical simulations to calculate the steady droplet deformation under subcritical conditions and the breakup process for capillary numbers larger than the critical one. We will consider the full hydrodynamic model, which includes inertia in both the inner and outer phases, soluto-capillarity and Marangoni convection due to inhomogeneities in the surfactant distribution, as well as pressure-dependent shear and dilatational surface viscosities. Therefore, our description does not invoke approximations considered in previous works, such as the limit of small droplet deformation (Flumerfelt Reference Flumerfelt1980; Vlahovska, Loewenberg & Blawzdziewicz Reference Vlahovska, Loewenberg and Blawzdziewicz2005; Vlahovska et al. Reference Vlahovska, Lawzdziewicz and Loewenberg2009; Narsimhan Reference Narsimhan2019; Singh & Narsimhan Reference Singh and Narsimhan2020).
The inclusion of gradients of surfactant concentration and surface tension will allow us to study the interplay among soluto-capillarity, Marangoni stress and viscous surface stress. This study could not be carried out in previous works (Pozrikidis Reference Pozrikidis1994; Gounley et al. Reference Gounley, Boedec, Jaeger and Leonetti2016; Narsimhan Reference Narsimhan2019; Singh & Narsimhan Reference Singh and Narsimhan2020), in which the inhomogeneity in surfactant concentration and surface tension was neglected. Besides, we will take into account nonlinear effects in both the surface equation of state (surface saturation and non-ideal interactions among the surfactant molecules) and the dependence of the surface viscosities on the surfactant concentration.
Our analysis extends that of Vlahovska et al. (Reference Vlahovska, Loewenberg and Blawzdziewicz2005, Reference Vlahovska, Lawzdziewicz and Loewenberg2009), who calculated the deformation of a surfactant-covered droplet in an extensional flow up to third order in the capillary number, neglecting the surfactant viscosity. The numerical simulations in the present work include several effects not accounted for by Milliken et al. (Reference Milliken, Stone and Leal1993): inertia, nonlinear effects in the surface equation of state, and surfactant viscosity. Our model is similar to that solved by Luo et al. (Reference Luo, Shang and Bai2019) for a droplet in a simple shear flow. They also included inertia, although the Reynolds number was set to a small value to eliminate its effects. They considered a nonlinear dependency of the shear and dilatational surface viscosities on the surfactant density but a linearized Langmuir equation of state.
We will calculate the steady solutions and determine their stability by conducting a global linear stability analysis. Direct (time-dependent) numerical simulations will be conducted for unstable configurations to investigate the droplet breakup mode. The calculation of the eigenmodes, combined with direct numerical simulations, provides a consistent and comprehensive picture of the droplet dynamics. This twofold approach has not been carried out in previous studies. Attention will be paid to the influence of the surface viscosity on both the droplet stability and the breaking mode. Viscous surface stress is expected to change the shape of the tip of the ejected fluid thread by altering the balance of surface stresses and the distribution of surfactant over the interface.
2. Governing equations
Consider a droplet of radius 
$a$, density 
$\rho _i$ and viscosity 
$\mu _i$ placed in the centre of a linear uniaxial extensional flow of a liquid of density 
$\rho _o$ and viscosity 
$\mu _o$ (figure 1). The interface is loaded with an insoluble surfactant. The interfacial tension 
$\sigma$ and surfactant surface concentration (surface coverage) 
$\varGamma$ at equilibrium are 
$\sigma _{{eq}}$ and 
$\varGamma _{{eq}}$, respectively. The variation of the interfacial tension with the surfactant surface concentration is approximated by the Langmuir equation of state (Tricot Reference Tricot1997)
\begin{equation} \sigma=\sigma_0+\varGamma_{\infty} R_g T\, \ln\left(1-\varGamma/\varGamma_{\infty}\right), \end{equation}
where 
$\sigma _0$ is the interfacial tension of the clean surface, 
$\varGamma _{\infty }$ is the maximum packing density, 
$R_g$ is the gas constant and 
$T$ is the temperature. The viscous surface stress is given by the Boussinesq–Scriven law in terms of the shear and dilatational surface viscosities 
$\mu ^S_{s}$ and 
$\mu ^S_{d}$ (Scriven Reference Scriven1960; Langevin Reference Langevin2014). Some studies (Kim et al. Reference Kim, Choi, Zell, Squires and Zasadzinski2013; Samaniuk & Mermant Reference Samaniuk and Mermant2014) have shown that the surface viscosity depends exponentially on the surface pressure 
$\varPi =\sigma _0-\sigma$ for some typical surfactants. In our study (Manikantan & Squires Reference Manikantan and Squires2017),
\begin{equation} \mu^S_s=\mu^S_{s,{eq}} \exp({(\varPi-\varPi_{{eq}})/\varPi_c}), \end{equation}
where 
$\mu ^S_{s,{eq}}$ and 
$\varPi _{{eq}}$ are the values taking by the corresponding quantities at equilibrium, while 
$\varPi _c$ is the surface pressure scale over which the surface viscosities change. For many insoluble surfactants, 
$\varPi _c$ is only a few milliNewtons per metre (Manikantan & Squires Reference Manikantan and Squires2017). Positive/negative values of this last quantity correspond to 
$\varPi$-thickening/thinning surfactants. In this work, we will assume that the ratio 
$\lambda ^S=\mu ^S_d/\mu ^S_s$ takes a constant value, independent of the surface pressure. Therefore, all the above comments on the shear surface viscosity also apply to the dilatational one.
Sketch of the fluid domain. The dashed line indicates the computational domain.
The droplet is placed in a flow given by the equations
\begin{equation} u^{({o})}={-}Gr/2,\quad w^{(o)}=G\, z, \end{equation}
where 
$u^{({o})}$ and 
$w^{(o)}$ are the radial and axial components of the outer velocity field, and 
$G$ is the intensity of the extensional flow.
Hereafter, all the variables are made dimensionless with the characteristic length 
$a$, velocity 
$\sigma _{{eq}}/\mu _0$, time 
$\mu _0 a/\sigma _{{eq}}$ and stress 
$\sigma _{{eq}}/a$ (Stone & Leal Reference Stone and Leal1989; Milliken et al. Reference Milliken, Stone and Leal1993). The dimensionless axisymmetric incompressible Navier–Stokes equations for the velocity 
$\boldsymbol {v}^{(k)}(r,z;t)$ and pressure 
$p^{(k)}(r,z;t)$ fields are
\begin{gather} [ru^{(k)}]_r+rw^{(k)}_z =0, \end{gather}
\begin{gather}{\rm Re}\,\rho^{\delta_{ik}}(u^{(k)}_t + u^{(k)} u^{(k)}_r+ w^{(k)} u^{(k)}_z)={-}p^{(k)}_r+\mu^{\delta_{ik}}[ u^{(k)}_{rr}+(u^{(k)}/r)_r+u^{(k)}_{zz}], \end{gather}
\begin{gather}{\rm Re}\,\rho^{\delta_{ik}}(w^{(k)}_t + u^{(k)} w^{(k)}_r+ w^{(k)}w^{(k)}_z) ={-} p^{(k)}_z+\mu^{\delta_{ik}}[w^{(k)}_{rr}+w^{(k)}_r/r+w^{(k)}_{zz}], \end{gather}
where 
$t$ is the time, 
$r$ (
$z$) is the radial (axial) coordinate, 
$u^{(k)}$ (
$w^{(k)}$) is the radial (axial) velocity component, 
${\rm Re}=\rho _o \sigma _{{eq}} a/\mu _o^2$ is the Reynolds number, 
$\rho =\rho _i/\rho _o$ and 
$\lambda =\mu _i/\mu _o$ are the density and viscosity ratios, respectively, and 
$\delta _{ij}$ is the Kronecker delta. In the above equations and henceforth, the superscripts 
$k=i$ and 
$k=o$ refer to the inner and outer phases, respectively, while the subscripts 
$t$, 
$r$ and 
$z$ denote the partial derivatives with respect to the corresponding variables. The action of the gravitational field has been neglected due to the smallness of the fluid configuration.
The interface equations are solved using the intrinsic surface coordinate 
$s$. The kinematic compatibility condition reads
\begin{equation} f_t+ \boldsymbol{v}^{(k)}\boldsymbol{\cdot} {\boldsymbol\nabla} f=0, \end{equation}
where 
$f(\boldsymbol {r}_s,t)=0$ determines the interface position 
$\boldsymbol {r}_s$. This position can also be defined in terms of the function 
$F(z,t)$, which represents the distance of an interface element to the symmetry axis 
$z$. The equilibrium of normal and tangential stresses at the interface yields
\begin{equation} \|\tau^{(k)}_n\|=\tau^S_n, \quad \|\tau^{(k)}_t\|=\tau^S_t, \end{equation}
where 
$\|A^{(k)}\|$ denotes the difference 
$A^{(i)}-A^{(o)}$ between the values taken by the quantity 
$A$ on the two sides of the interface, and 
$\tau ^{(k)}_n$ and 
$\tau ^{(k)}_t$ represent the bulk stresses normal and tangential to the interface, respectively. For a Newtonian inertialess interface (Scriven Reference Scriven1960), and using intrinsic surface coordinates, the normal and tangential surface stresses 
$\tau ^S_n$ and 
$\tau ^S_t$ read (see Appendix A) (Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2018)
\begin{gather} \tau^S_n=\kappa\hat{\sigma}+ B_s \bar{\kappa}\left[{-}F\left(\frac{v_t}{F}\right)' + \bar{\kappa} v_n \right] +\lambda^S B_s\kappa\left[\frac{(Fv_t)'}{F} + \kappa v_n\right], \end{gather}
\begin{gather}\tau^S_t=\hat{\sigma}' + \left[B_s(1+\lambda^S)\left(\frac{(Fv_t)'}{F} + \kappa v_n\right)\right]' +2\left(B_s K - B_s'\frac{F'}{F}\right)v_t - 2\kappa_1(B_sv_n)', \end{gather}
where 
$\hat {\sigma }=\sigma /\sigma _{{eq}}$ is the ratio of the local interfacial tension to its equilibrium value, 
$B_s=\mu ^S_s/(\mu _oa)$ is the Boussinesq number, 
$\kappa =\kappa _1+\kappa _2$ is (twice) the mean curvature of the interface, 
$\kappa _1$ and 
$\kappa _2$ are the curvatures along the meridians and parallels in the normal direction, respectively, 
$\bar {\kappa } = \kappa _1 - \kappa _2$, 
$K = \kappa _1\kappa _2$, 
$v_n$ and 
$v_t$ are the normal and tangential components of the velocity on the interface, and the prime denotes the derivative with respect to the intrinsic surface coordinate. Equations (2.9) and (2.10) reduce to those derived by Martínez-Calvo (Reference Martínez-Calvo2020) in cylindrical coordinates.
The surfactant surface concentration is calculated by integrating the conservation equation
\begin{equation} \hat{\varGamma}_t+{\boldsymbol\nabla}^S\boldsymbol{\cdot} (\hat{\varGamma}\boldsymbol{v}^S)+\hat{\varGamma} ({\boldsymbol \nabla}^{S}\boldsymbol{\cdot} {\boldsymbol n})\,{\boldsymbol n}\boldsymbol{\cdot}{\boldsymbol v}=\frac{1}{{Pe}^S}\, {\nabla}^{S2}\hat{\varGamma}, \end{equation}
where 
$\hat {\varGamma }=\varGamma /\varGamma _{\infty }$ is the surface coverage defined as the ratio of the surfactant surface concentration 
$\varGamma$ to the maximum packing density 
$\varGamma _{\infty }$, 
${\boldsymbol n}$ is the unit outward normal vector, 
${Pe}^S=a\sigma _{{eq}}/(\mu _o D^S)$ is the surface Péclet number and 
$D^S$ is the surface diffusion coefficient. This coefficient must also depend on the surface concentration, consistently with the surface viscosities. However, we follow previous works and neglect this effect because it has no significant influence on our results, given the large value taken by the surface Péclet number in our simulations.
The variation of the interfacial tension with the surface coverage is approximated by the Langmuir equation of state (2.1), whose dimensionless form is
\begin{equation} \hat{\sigma}=1+{Ma} \ln\left(\frac{1-\hat{\varGamma}}{1-\hat{\varGamma}_{{eq}}}\right), \end{equation}
where 
$Ma=\varGamma _{\infty } R_g T/\sigma _{{eq}}$ is the Marangoni (elasticity) number and 
$\hat {\varGamma }_{{eq}}=\varGamma _{{eq}}/\varGamma _{\infty }$. The values of 
$\hat {\varGamma }_{{eq}}$ and 
$Ma$ must be chosen so that the interfacial tension does not become negative at any point of the interface in the course of the simulation. Taking into account the exponential dependence (2.2) and the equation of state (2.12), the Boussinesq number can be calculated as
\begin{equation} B_s=B_{s,{eq}} \left(\frac{1-\hat{\varGamma}_{{eq}}}{1-\hat{\varGamma}}\right)^{{Ma}/\hat{\varPi}_c}, \end{equation}
where 
$B_{s,{eq}}$ is its value at equilibrium and 
$\hat {\varPi }_c=\varPi _c/\sigma _{{eq}}$. For 
${Ma}\ll \hat {\varPi }_c$, the Boussinesq number (i.e. the surface viscosity) becomes constant. Since 
$\hat {\varPi }_c$ takes values much smaller than unity for many surfactants (Manikantan & Squires Reference Manikantan and Squires2017), this simplification is valid only for very small values of 
$Ma$.
The conservation of the droplet and surfactant mass yields
\begin{equation} \int_{0}^{\hat{a}} F^2\, {\rm d}z=\frac{2}{3},\quad \int_{0}^{\hat{a}_s} \hat{\varGamma}\, F\, {\rm d}s=2\, \hat{\varGamma}_{{eq}}, \end{equation}
respectively, where 
$\hat {a}$ and 
$\hat {a}_s$ are the (dimensionless) half-lengths of the cross-sectional shape and the interface, respectively. Equations (2.14a,b) must be considered to close the steady problem and to set the initial conditions for the transient numerical simulations. The droplet mass must remain constant during the droplet breakup due to the liquid incompressibility and kinematic boundary conditions. Analogously, (2.11) ensures the conservation of the surfactant mass over time in the transient problem. However, and to reduce the numerical errors, we also enforce (2.14a,b) at any instant of the transient simulations.
The 
$r$ and 
$z$ axes delimiting the computational domain are symmetry axes (figure 1). We impose the velocity field 
$u^{(o)}=-Cr/2$ and 
$w^{(o)}=Cz$ at the two other boundaries of the computational domain (figure 1), where 
$C=Ga\mu _o/\sigma _{{eq}}$ is the capillary number. Transient simulations of the droplet breakup start from the steady solution for a capillary number just below the critical one, which allows us to reduce the computing time considerably.
As can be seen, the problem is formulated in terms of the set of dimensionless numbers 
$\{\rho$, 
${\rm Re}$; 
$\lambda$, 
$C$; 
${Pe}^S$, 
$Ma$, 
$\hat {\varGamma }_{{eq}}$; 
$B_{s,{eq}}$, 
$\lambda ^S$, 
$\hat {\varPi }_c\}$ (table 1). In the absence of surfactant, the parameter space reduces to 
$\{\rho, {\rm Re}; \lambda, C\}$. Besides, the problem is formulated only in terms of 
$\lambda$ and 
$C$ when inertia is neglected. It is worth mentioning that the capillary number 
$C$ indicates the dimensionless characteristic velocity, 
$G a/(\sigma _{{eq}}/\mu _o)$, of the outer fluid. The Reynolds number is defined in terms of material properties and the droplet radius, which allows us to fix its value in an experimental run in which the strain rate (the capillary number) is increased.
Dimensionless parameters, their physical meanings and the values considered in the simulations. The numbers in bold type correspond to most of the simulations.
To calculate the linear global modes of the steady solutions, we assume the temporal dependence
\begin{equation} \left. \begin{aligned} U & =U_0+\delta U\, {\rm e}^{-{\rm i}\omega t} \quad (|\delta U|\ll |U|), \\ (r_s,z_s) & =(r_{s0},z_{s0})+(\delta r_s,\delta z_s)\, {\rm e}^{-{\rm i}\omega t} \quad (|\delta r_s|\ll r_s, |\delta z_s|\ll z_s ), \end{aligned} \right\} \end{equation}
where 
$U$ represents any unknowns of the problem, and 
$U_0$ and 
$\delta U$ stand for the corresponding base flow (steady) solution and the spatial dependence of the eigenmode, respectively. In addition, 
$(r_s,z_s)$ denotes the interface position, 
$(r_{s0},z_{s0})$ denotes the interface position in the base flow, 
$(\delta r_s,\delta z_s)$ denotes the perturbation, and 
$\omega =\omega _r+\textrm {i}\omega _i$ is the eigenfrequency characterizing the perturbation evolution. If the growth rate 
$\omega _i$ of the dominant mode (i.e. that with the largest 
$\omega _i$) is positive, then the base flow is asymptotically unstable under small-amplitude perturbations (Theofilis Reference Theofilis2011). In this work, we will determine the critical value of the capillary number 
$C$ for which the base flow becomes asymptotically unstable as a function of the rest of the governing parameters.
3. Numerical method
We used the numerical method proposed by Herrada & Montanero (Reference Herrada and Montanero2016) to solve the theoretical model described in the previous section. Here, we summarize the main characteristics of this method. The inner and outer fluid domains are mapped onto two quadrangular domains through a non-singular mapping. A quasi-elliptic transformation (Dimakopoulos & Tsamopoulos Reference Dimakopoulos and Tsamopoulos2003) is applied in the outer bath. All the derivatives appearing in the governing equations are expressed in terms of 
$t$ and the spatial coordinates resulting from the mapping. These equations are discretized in the (mapped) radial direction with 
$n_{\chi }^{(i)}$ and 
$n_{\chi }^{(o)}$ Chebyshev spectral collocation points (Khorrami, Malik & Ash Reference Khorrami, Malik and Ash1989) in the inner and outer regions, respectively. We use fourth-order finite differences with 
$n_{\xi }^{(i)}$ and 
$n_{\xi }^{(o)}$ equally spaced points to discretize the (mapped) axial direction in the inner and outer regions, respectively. As can be seen in figure 2, the grid points are equally spaced along the interface. In simulations for very small values of the viscosity ratio 
$\lambda$, the curvature of the droplet tip increases considerably. In those particular cases, we accumulated the grid points in the vicinity of the interface tip. We introduced a stretching function so that the distance decreases linearly with the distance to the droplet apex. The distance between two consecutive points near the apex is around seven times smaller than the distance near the equator.
Details of the grid used in the simulations.
In the transient numerical simulations, second-order backward finite differences are used to discretize the time domain. The time step is adapted in the course of the simulation according to the formula 
$\Delta t=\Delta t_0/v_{{tip}}$, where 
$\Delta t_0$ is the time step at the initial instant, and 
$v_{{tip}}$ is the droplet tip velocity. The time-dependent mapping of the physical domain does not allow the algorithm to surpass the interface pinch-off, therefore the evolution of the emitted droplet cannot be analysed.
To calculate the base flow and its eigenmodes, we progressively increased the capillary number, using the solution obtained for the previous case as an initial guess. We do not calculate the base flow from a dynamical simulation, starting from the previous shape. Our method solves the set of nonlinear equations to find the stationary solution corresponding to a subcritical capillary number, which constitutes a major contribution compared to previous works. To simulate the droplet breakup, we considered as initial condition the steady solution for a subcritical case close to the stability limit and then increased the capillary number up to its prescribed value. The end of the transient simulation is the last time step for which the numerical method converged.
The results presented in this work for 
$\lambda =0.1$ (liquid droplet) were calculated with 
$\{n_{\chi }^{(i)}=21$, 
$n_{\xi }^{(i)}=701$, 
$n_{\chi }^{(o)}=41$, 
$n_{\xi }^{(o)}=1403\}$. In the transient simulations, we set 
$\Delta t_0=0.02$. Figure 3 shows the velocity 
$v_{{tip}}$ of the droplet tip as a function of time 
$t$ for four spatiotemporal discretizations. As can be observed, the results are practically insensitive to the grid refinement for the discretization used in our simulations.
Velocity 
$v_{{tip}}$ of the droplet tip as a function of time 
$t$ for four spatial and temporal discretizations: 
$\{n_{\chi }^{(i)}=18$, 
$n_{\xi }^{(i)}=501$, 
$n_{\chi }^{(o)}=31$, 
$n_{\xi }^{(o)}=1003$, 
$\Delta t_0=0.02\}$ (green diamonds), 
$\{n_{\chi }^{(i)}=21$, 
$n_{\xi }^{(i)}=701$, 
$n_{\chi }^{(o)}=41$, 
$n_{\xi }^{(o)}=1403$, 
$\Delta t_0=0.02\}$ (black circles), 
$\{n_{\chi }^{(i)}=21$, 
$n_{\xi }^{(i)}=701$, 
$n_{\chi }^{(o)}=41$, 
$n_{\xi }^{(o)}=1403$, 
$\Delta t_0=0.01\}$ (blue squares) and 
$\{n_{\chi }^{(i)}=25$, 
$n_{\xi }^{(i)}=801$, 
$n_{\chi }^{(o)}=45$, 
$n_{\xi }^{(o)}=1603$, 
$\Delta t_0=0.02\}$ (red triangles). The values of the governing parameters are 
$\{{\rm Re}=0$, 
$\lambda =0.1$, 
$C=0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$B_{s,{eq}}=0\}$.
As mentioned above, bubbles develop sharp tips even for capillary numbers smaller than unity, which demands a higher spatial resolution along the 
$\xi$ axis. For this reason, the steady simulations for 
$\lambda =10^{-3}$ were conducted for 
$\{n_{\chi }^{(i)}=5$, 
$n_{\xi }^{(i)}=1301$, 
$n_{\chi }^{(o)}=21$, 
$n_{\xi }^{(o)}=2603\}$. To show the accuracy of these simulations, we compared our numerical solution for 
$\lambda =10^{-7}$ in the absence of surfactant with that obtained by Eggers & Courrech du Pont (Reference Eggers and Courrech du Pont2009) for 
$\lambda =0$ (figure 4). As can be observed, an excellent agreement is obtained for the mean curvature 
$\kappa _{{tip}}$ at the bubble tip even for curvature radii of the order of 
$10^{-3}$. This is a very stringent test, which indicates that the tip shape (whose curvature is at least two orders of magnitude smaller in the presence of surface viscosity) is well-resolved.
(a) Mean curvature 
$\kappa _{{tip}}$ at the bubble tip as a function of the capillary number 
$C$ for 
$\{{\rm Re}=0$, 
$\lambda =10^{-7}\}$ (symbols). The solid line corresponds to the numerical result of Eggers & Courrech du Pont (Reference Eggers and Courrech du Pont2009) for the inviscid case. (b) Bubble shape for 
$C=0.33$.
The droplet steady shape is characterized by the deformation
\begin{equation} D=\frac{\hat{a}-\hat{b}}{\hat{a}+\hat{b}}, \end{equation}
where 
$\hat {a}$ and 
$\hat {b}$ are the half-length and half-breadth of the cross-sectional shape, respectively. Figure 5 shows a comparison between our numerical results and those of Eggers & Courrech du Pont (Reference Eggers and Courrech du Pont2009). The figure also shows the perturbation theory to first and second order in 
$Ca$ of Barthls-Biesel & Acrivos (Reference Barthls-Biesel and Acrivos1973).
Droplet deformation 
$D$ as a function of the capillary number 
$C$ for 
$\{{\rm Re}=0$, 
$\lambda =10^{-7}\}$ in the absence of surfactant (symbols). The solid line corresponds to the numerical result of Eggers & Courrech du Pont (Reference Eggers and Courrech du Pont2009) for 
$\lambda =0$. The dotted and dashed lines correspond to the perturbation theory to first and second order in 
$Ca$ of Barthls-Biesel & Acrivos (Reference Barthls-Biesel and Acrivos1973).
4. Results
Before presenting the numerical results, we here justify our choice for the values of the governing parameters. Given the large dimension of the parameter space, we will consider two values of the bulk viscosity ratio, 
$\lambda =10^{-1}$ and 
$10^{-3}$, which represent a liquid droplet and a bubble suspended in a viscous liquid bath, respectively. The intermediate value 
$\lambda =10^{-2}$ will be considered to analyse the dependence of the interface perturbation on 
$\lambda$ at the marginal stability. We will analyse the effect of inertia by considering 
$\rho =1$ (density-matched liquids) and 
${\rm Re}=1$ and 10. We will restrict ourselves to 
$Ma=0.2$, which corresponds to a strong surfactant (Eggleton et al. Reference Eggleton, Pawar and Stebe1999, Reference Eggleton, Tsai and Stebe2001; Wang et al. Reference Wang, Siegel and Booty2014). In most cases, we will consider a moderately dense surfactant monolayer characterized by a surface coverage 
$\hat {\varGamma }_{{eq}}=0.5$. We will decrease the surface coverage down to 0.1 to examine the effect of the amount of surfactant on the droplet stability.
Given the small values taken by the surface diffusion coefficient of most surfactants (Tricot Reference Tricot1997), the surface Péclet number is set to 
$Pe^S=10^3$ in all the simulations. The shear surface viscosity of moderately viscous surfactants can take values of the order of 
$10^{-6}$ Pa s m (Ponce-Torres et al. Reference Ponce-Torres, Montanero, Herrada, Vega and Vega2017), while this property decreases down to values as small as 
$10^{-10}$ Pa s m for surfactants commonly used in experiments, such as sodium dodecyl sulfate (SDS) (Zell et al. Reference Zell, Nowbahar, Mansard, Leal, Deshmukh, Mecca, Tucker and Squires2014; Ponce-Torres et al. Reference Ponce-Torres, Rubio, Herrada, Eggers and Montanero2020). To obtain realistic values for the Boussinesq number, consider, for instance, a water (
$\mu _i=1$ mPa s) droplet 
$a=1$ mm in radius submerged in a liquid bath with viscosity 
$\mu _0=10$ mPa s. In this case, the above-mentioned values of the shear surface viscosity lead to Boussinesq numbers of the order of 
$10^{-1}$ and 
$10^{-5}$, respectively. We will consider values of this parameter in the interval 
$10^{-6}$–1. In some cases, we will consider higher values of the surface viscosity to highlight its effect.
The case 
$B_{s,{eq}}=1$ corresponds to a viscous surfactant. In most cases, we will analyse the effect of the interfacial rheology for equal surface viscosities (
$\lambda ^S=1$). To elucidate the role of those viscosities separately, we will also consider the cases 
$\lambda ^S=0$ and 
$10^3$. As mentioned in § 2, most surfactants exhibit a thickening behaviour (the surface viscosities increase with the surfactant concentration), 
$\varPi _c$ being only a few milliNewtons per metre (Manikantan & Squires Reference Manikantan and Squires2017). For this reason, we will take 
$\hat {\varPi }_c=0.1$ in all our simulations. Table 1 displays the values of the dimensionless governing parameters in our simulations.
4.1. Droplet shape and stability
This subsection examines both the steady deformation and stability of droplets submerged in an extensional flow. The stability is determined from the spectrum of eigenvalues obtained for a given base flow. For the sake of illustration, figure 6 shows the spectrum of eigenvalues with 
$\omega _i>-4.66$ for an inertialess (
${\rm Re}=0$) drop close to the stability limit. The dominant eigenvalue is an imaginary number that becomes positive for 
$C\gtrsim 0.0986$. At marginal stability, 
$C\simeq 0.0986$, both the frequency and damping rate (the growth rate with the sign reversed) of the dominant eigenvalue vanish. This means that the flow becomes unstable under stationary linear perturbations, contrary to what happens in, for example, the jetting mode of flow focusing, in which instability is caused by a supercritical Hopf bifurcation (Cruz-Mazo et al. Reference Cruz-Mazo, Herrada, Gañán-Calvo and Montanero2017; Cabezas et al. Reference Cabezas, Rubio, Rebollo-Muñoz, Herrada and Montanero2021). The results shown in figure 6 are qualitatively the same as those obtained for the rest of the parameter configurations analysed in this work. The flow becomes unstable under stationary linear perturbations, preventing the numerical method from converging to the base flow solution when the capillary number is fixed in the supercritical regime. In other words, the stability limit corresponds to the capillary number for which the numerical method ceases to converge to a steady solution, a correspondence implicitly assumed in previous works.
Spectrum of eigenvalues with 
$\omega _i>-4.66$ for 
$\{{\rm Re}=0$, 
$\lambda =0.1$, 
$C=0.09855$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$B_{s,{eq}}=0\}$. The arrow indicates the eigenvalue of the critical mode.
The droplet deformation 
$D$ monotonically increases with the increasing capillary number and increases sharply near the saddle node bifurcation (figure 7), which corresponds to a turning point (Taylor Reference Taylor1964; Acrivos & Lo Reference Acrivos and Lo1978). As can be observed, the presence of a surfactant monolayer considerably reduces the critical capillary number and the maximum deformation reached by the droplet. In fact, the critical capillary number for the configuration considered in figure 7 increases up to 0.171 in the absence of surfactant. The surfactant added to the interface is driven towards the droplet tip by the outer stream. The droplet tip weakens due to the resulting reduction of the interfacial tension, and therefore the droplet breaks up for lower values of the capillary number. As can be observed in figure 7, the damping rate of the dominant mode for a surfactant-covered droplet, 
$-\omega _i^*$, monotonically decreases as the capillary number approaches its critical value. The curve 
$\omega _i^*(C)$ for a clean interface (
$\hat {\varGamma }_{{eq}}=0$) shows the crossover of two modes.
Droplet deformation 
$D$ (a), and imaginary part of the dominant eigenvalue 
$\omega _i^*$ (b), as functions of the capillary number 
$C$ for 
$\{{\rm Re}=0$, 
$\lambda =0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$B_{s,{eq}}=0\}$. The figure also shows the results for a clean interface (
$\hat {\varGamma }_{{eq}}=0$). The dotted lines are guides for the eye.
The drag force exerted by the outer flow produces a recirculation pattern in the droplet (figure 8). The fluid particles next to the interface are driven towards the apex. The hydrostatic pressure increases there, and the induced pressure gradient makes the liquid flow back next to the droplet symmetry axis. This flow pattern resembles those produced in other microfluidic configurations such as flow focusing (Cruz-Mazo et al. Reference Cruz-Mazo, Herrada, Gañán-Calvo and Montanero2017; Cabezas et al. Reference Cabezas, Rubio, Rebollo-Muñoz, Herrada and Montanero2021) and electrospray (Ponce-Torres et al. Reference Ponce-Torres, Rebollo-Muñoz, Herrada, Gañán-Calvo and Montanero2018). However, the intensity of this recirculation is very small owing to the interface immobilization caused by the surfactant monolayer. As pointed out by Milliken et al. (Reference Milliken, Stone and Leal1993), Marangoni force practically counteracts the imposed external flow and the velocity on the surface almost vanishes, thus the interior fluid is practically motionless, regardless of the viscosity ratio. In fact, the velocity field inside the droplet of figure 8 is two orders of magnitude smaller than that of the imposed external flow. To gain insight into the physical mechanisms responsible for the global instability of this flow, we consider the perturbation 
$\delta p^{(j)}(r,z)$ of the pressure field. Figure 8 shows the isocontours of 
$|\delta p^{(j)}(r,z)|$ for the eigenmode causing the base flow instability. As can be observed, the perturbation of the pressure field increases sharply right in front of the droplet apex, which indicates that the steady flow destabilization originates at that point.
Streamlines of the base flow and isocontours of the magnitude of the perturbed pressure field 
$|\delta p^{(j)}(r,z)|$ for 
$\{{\rm Re}=0$, 
$\lambda =0.1$, 
$C=0.0986$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$B_{s,{eq}}=0\}$. The colour scale shows the values of 
$|\delta p^{(j)}(r,z)|$ normalized with its maximum value. The red line is the interface location.
Figure 9 shows the interface displacement due to the growth of the critical eigenmode at the quasi-marginally stable state for different surfactant concentrations. This displacement corresponds to the interface deformation at the initial (linear) phase of the droplet breakup. As will be shown in § 4.3, the instability described by this eigenmode leads to tip streaming beyond the linear regime in the presence of surfactant. Interestingly, the perturbation affects most of the interface, not only the droplet tip, which indicates that the small scale characterizing the tip streaming is fixed in the nonlinear phase of the droplet deformation. In fact, the interface perturbation in the linear regime is qualitatively the same as that of a droplet in the absence of surfactant (
$\hat {\varGamma }_{{eq}}=0$), which breaks up following the central pinching mode, as will be shown in § 4.3. The perturbation also affects a considerable portion of the interface when the viscosity ratio is reduced down to 
$\lambda =10^{-2}$ (figure 10), although it becomes more localized around the droplet tip. This behaviour is different from that observed in selective withdrawal, where the eigenmode is highly localized even for larger viscosity ratios (Eggers & Courrech du Pont Reference Eggers and Courrech du Pont2010). This difference may be the result of the volume constraint in our problem.
Droplet shape in the base flow (shaded area) and interface displacement due to the eigenmode, 
$(r_s,z_s)=(r_{s0},z_{s0})+\phi \, ({\rm Re}[\delta r_s],{\rm Re}[\delta z_s])$ (dashed lines) (see (2.15)), for 
$\hat {\varGamma }_{{eq}}=0$ and 
$C=0.171$, 
$\hat {\varGamma }_{{eq}}=0.1$ and 
$C=0.13$, 
$\hat {\varGamma }_{{eq}}=0.25$ and 
$C=0.106$, and 
$\hat {\varGamma }_{{eq}}=0.5$ and 
$C=0.09855$. The values of the rest of governing parameters are 
$\{{Re}=0$, 
$\lambda =0.1$, 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$B_{s,{eq}}=0\}$. The value of the arbitrary constant 
$\phi$ in the linear analysis has been chosen to appreciate the interface deformation.
Droplet shape in the base flow (shaded area) and interface displacement due to the eigenmode, 
$(r_s,z_s)=(r_{s0},z_{s0})+\phi \, ({\rm Re}[\delta r_s],{\rm Re}[\delta z_s])$ (dashed line) (see (2.15)), for 
$\lambda =10^{-2}$ and 
$C=0.2505$ in the absence of surfactant. The value of the arbitrary constant 
$\phi$ in the linear analysis has been chosen to appreciate the interface deformation.
We selected in figure 11 a viscosity ratio two orders of magnitude smaller than that considered in figure 7, which corresponds approximately to replacing a water droplet with an air bubble, both submerged in the same liquid bath. As can be observed, the deformation, damping rate and critical capillary number are hardly changed. This shows that the instability mechanism for a surfactant-laden drop is completely different from that without surfactant since in the latter case, instability is controlled by the value of 
$\lambda$ alone. The fact that the damping rate is hardly affected by the viscosity ratio indicates that viscous dissipation in the droplet bulk barely contributes to the damping of the droplet oscillations. This occurs due to the interface immobilization caused by the Marangoni stress, as explained above. The comparison with the clean interface case (
$\hat {\varGamma }_{{eq}}=0$) shows how the interval of stable capillary numbers significantly reduces due to the presence of the surfactant monolayer. As occurs for 
$\lambda =0.1$, the curve 
$\omega _i^*(C)$ for a clean interface shows the crossover of two aperiodic modes, as indicated by the solid symbols in figure 11(b). We could not determine the stability limit for 
$\lambda =10^{-3}$ and 
$\hat {\varGamma }_{{eq}}=0$ due to spatial discretization errors caused by the pointed shape of the droplet tip in that case.
Droplet deformation 
$D$ (a), and imaginary part of the dominant eigenvalue, 
$\omega _i^*$ (b), as functions of the capillary number 
$C$ for 
$\{{\rm Re}=0$, 
$\lambda =10^{-3}$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$B_{s,{eq}}=0\}$. The figure also shows the results for a clean interface (
$\hat {\varGamma }_{{eq}}=0$). The dotted lines are guides for the eye.
A dense and strong surfactant monolayer (
$\hat {\varGamma }_{{eq}}=0.5$, 
$Ma=0.2$) greatly destabilizes the droplet, considerably reducing the critical capillary number. The comparison between figures 7 and 11 shows that in the presence of the surfactant monolayer, the shape and stability of the droplet are hardly affected by the viscosity ratio for the reduced interval of stable capillary numbers. In the absence of surfactant, the critical capillary number increases up to 
$C=0.171$ for 
$\lambda =0.1$ and 
$C\simeq 0.32$ for 
$\lambda =10^{-3}$, there is more flow inside the droplet, and the interior viscosity becomes relevant.
The effect of the droplet and fluid bath inertia on both the droplet deformation and stability is analysed in figure 12. For 
${\rm Re}=1$, the deformation and critical capillary number are hardly affected by inertia. However, the deformation significantly increases, and the critical capillary number considerably decreases, for 
${\rm Re}=10$. The destabilizing effect of inertia has also been observed in an extensional outer flow with non-zero Reynolds number (Acrivos & Lo Reference Acrivos and Lo1978) and for a droplet covered with an inviscid surfactant in the presence of a three-dimensional shear flow (Liu et al. Reference Liu, Ba, Wu, Li, Xi and Zhang2018). However, it is the opposite effect to that found by Brady & Acrivos (Reference Brady and Acrivos1982) when the drop inertia is taken into account. Interestingly, the eigenmode responsible for the instability is aperiodic (
$\omega _r= 0$) even for large inertia. This contrasts with what occurs in other microfluidic extensional flows such as the jetting mode of flow focusing (Cruz-Mazo et al. Reference Cruz-Mazo, Herrada, Gañán-Calvo and Montanero2017; Cabezas et al. Reference Cabezas, Rubio, Rebollo-Muñoz, Herrada and Montanero2021) and electrospray (Ponce-Torres et al. Reference Ponce-Torres, Rebollo-Muñoz, Herrada, Gañán-Calvo and Montanero2018), in which the loss of stability is caused by an oscillatory mode (
$\omega _r\neq 0$, supercritical Hopf bifurcation). The curve 
$\omega _i(C)$ for 
${\rm Re}=10$ shows the crossover of different modes. The solid circles correspond to a stable (damped) oscillatory eigenmode, which becomes the dominant one for 
$0.03\lesssim C\lesssim 0.07$.
Droplet deformation 
$D$ and imaginary part of the dominant eigenvalue 
$\omega _i^*$ as functions of the capillary number 
$C$ for different values of the Reynolds number number, 
${\rm Re}=0$, 1 and 10. The values of the rest of the governing parameters are 
$\{\rho =1$, 
$\lambda =0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$B_{s,{eq}}=0\}$. The solid circles correspond to damped oscillatory perturbations. The vertical dotted lines indicate the values of the critical capillary number.
Vlahovska et al. (Reference Vlahovska, Lawzdziewicz and Loewenberg2009) calculated analytically the steady deformation of a droplet covered with an inviscid surfactant up to third order in 
$C$. In this solution, the Marangoni convection counteracts the external flow and completely suppresses the flow inside the droplet, no matter how small the Marangoni number is. As a consequence, the droplet deformation for 
${Ma}>0$ does not depend on the viscosity ratio 
$\lambda$. Figure 13 compares the droplet deformation 
$D$ obtained in our simulation for zero surface viscosity and the first-order and third-order perturbation theory of Vlahovska et al. (Reference Vlahovska, Lawzdziewicz and Loewenberg2009). As can be observed, the third-order approximation significantly improves the linear theory, although it considerably underestimates the droplet deformation next to the critical capillary number.
Droplet deformation 
$D$ as a function of the capillary number 
$C$ for 
$\lambda =0.1$ and 
$10^{-3}$. The values of the rest of the governing parameters are 
$\{{\rm Re}=0$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$B_{s,{eq}}=0\}$. The symbols are the simulation results, while the solid and dashed lines are the first-order and third-order perturbation theory of Vlahovska et al. (Reference Vlahovska, Lawzdziewicz and Loewenberg2009).
4.2. Droplet shape and stability. Influence of surface viscosity
In this subsection, we examine the influence of the surface viscosity on the droplet deformation 
$D$ and the imaginary part of the dominant eigenvalue, 
$\omega _i^*$. As can be observed in figure 14, the viscous surface stress does not significantly affect the droplet deformation even for the viscous surfactant (
$B_{s,{eq}}=1$). The competition between the capillary and viscous surface stresses can be measured in terms of the product 
$C\, B_{s,{eq}}=G\,\mu _{s,{eq}}^S/\sigma _{{eq}}$. This product takes values smaller than 
$10^{-1}$ in all the cases considered. Normal stresses control the droplet shape, which essentially determines the critical capillary number. This partially explains why this critical value is practically the same for all the values of the Boussinesq number. More importantly, Marangoni stress suppresses the interface motion, which renders the surface viscosity stress negligible.
Droplet deformation 
$D$, and imaginary part of the dominant eigenvalue 
$\omega _i^*$, as functions of the capillary number 
$C$ for different values of the Boussinesq number 
$B_{s,{eq}}=0$, 0.01 and 1. The values of the rest of the governing parameters are 
$\{{\rm Re}=0$, 
$\lambda =0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$\lambda ^S=1$, 
$\hat {\varPi }_c=0.1\}$. The dotted lines are guides for the eye.
The damping rate of the dominant mode is significantly influenced by surface viscosity only for 
$B_{s,{eq}}=1$ (viscous surfactant). It is somewhat counter-intuitive that the damping rate decreases as the surface viscosity increases. In fact, one would expect a dissipative factor, such as viscous surface stresses, to stabilize the system under small-amplitude perturbations, increasing the rate at which those perturbations are damped. This result resembles what occurs in other capillary systems, such as flowing rivulets (Herrada et al. Reference Herrada, Mohamed, Montanero and Gañán-Calvo2015), where the same basic flow can become unstable under infinitesimal perturbations when viscosity exceeds a certain critical value. The decrease of the damping rate as the surface viscosity increases may be explained as follows. Surface viscosity contributes to the interface immobilization, reducing the fluid speed over that surface and therefore inside the droplet. The decrease of the flow in the droplet reduces the viscous damping of the perturbation. To show the interface immobilization due to surface viscosity, we run simulations for 
$B_{s,{eq}}=0$ and 10, and uniform surfactant concentration, i.e. in the absence of soluto-capillarity and Marangoni stress. Once this stress has been removed, surface viscosity becomes noticeable and suppresses the interface motion (figure 15). While for 
$B_{s,{eq}}=0$ the interface velocity is of the order of the outer velocity (
$v_t\sim C$), it practically vanishes in the case 
$B_{s,{eq}}=10$. The damping rates for 
$B_{s,{eq}}=0$ and 10 are 
$-0.656$ and 
$-0.138$, respectively.
Surface velocity 
$v_t$ on the interface for 
$B_{s,{eq}}=0$ and 10. The values of the rest of the governing parameters are 
$\{{\rm Re}=0$, 
$\lambda =0.1$, 
$C=0.04$; 
$\lambda ^S=1$, 
$\hat {\varGamma }_{{eq}}=0.5\}$. We imposed the condition 
$\hat {\varGamma }=\hat {\varGamma }_{{eq}}$ to suppress soluto-capillarity and Marangoni stress.
As mentioned above, the viscous surface stress is much smaller than the capillary pressure because the product 
$C\,B_{s,{eq}}$ takes values much lower than unity. However, surface viscosity could still significantly affect the droplet shape by changing the surfactant distribution over the interface and therefore the capillary pressure profile. However, this is not the case, as demonstrated in figure 16, where the distributions of interfacial tension, surfactant density and tangential surface velocity at the critical point are shown both in the absence of surface viscosity and for the viscous case, 
$B_{s,{eq}}=1$. The figure also shows the corresponding droplet shapes. The surface velocity takes very small values between the stagnation points located at 
$z=0$ and the half-length of the cross-sectional shape 
$z=\hat {a}$. In fact, those values are two orders of magnitude smaller than the characteristic outer fluid velocity 
$C$, and much smaller than the values reached at the droplet tip during tip streaming in the supercritical case (see § 4.3). As a consequence, the viscous surface stress takes small values and hardly modifies the distribution of surfactant over the interface. For that reason, the interfacial tension profile and therefore the droplet shape are essentially the same as those obtained in the inviscid case at the corresponding stability limits. The droplet loaded with a viscous surfactant is slightly more stretched by the outer flow for the same capillary number. Previous results have shown that shear surface viscosity immobilizes the interface when a droplet is submerged in a simple shear flow (Pozrikidis Reference Pozrikidis1994; Luo et al. Reference Luo, Shang and Bai2019), while dilatational surface viscosity reduces the surfactant concentration gradient (Luo et al. Reference Luo, Shang and Bai2019). As can be seen in figure 16(b,c), these effects are also observed in our simulations for uniaxial extensional flow, although their magnitude is small for 
$B_{s,{eq}}=1$. The above results allow us to conclude that the effects of surface viscosity are greatly diminished by Marangoni convection. In the absence of this convection, surface viscosity may play a significant role in both the droplet deformation and stability (Narsimhan Reference Narsimhan2019; Singh & Narsimhan Reference Singh and Narsimhan2020). This condition can be reached for 
$B_{s,{eq}}/({Pe}^S\,{Ma})$ of at least order unity.
Surface tension 
$\hat {\sigma }$ (a), surface coverage 
$\hat {\varGamma }$ (b), tangential surface velocity 
$v_t$ (c), and interface position 
$F$ (d), as functions of the axial coordinate 
$z$ for 
$B_{s,{eq}}=0$ and 
$C=0.09855$ (blue lines) and 
$B_{s,{eq}}=1$ and 
$C=0.09835$ (black lines). The values of the rest of the governing parameters are 
$\{{\rm Re}=0$, 
$\lambda =0.1$, 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$\lambda ^S=1$, 
$\hat {\varPi }_c=0.1\}$.
Figure 17 compares the bubble deformation and the damping rate of the dominant mode for an inviscid and viscous surfactant monolayer. As can be observed, the damping rate is significantly influenced by the Boussinesq number, which means that the surface viscosity plays a more important role than the inner one.
Droplet deformation 
$D$ and imaginary part of the dominant eigenvalue 
$\omega _i^*$ as functions of the capillary number 
$C$ for different values of the Boussinesq number 
$B_{s,{eq}}=0$ and 1. The values of the rest of the governing parameters are 
$\{{\rm Re}=0$, 
$\lambda =10^{-3}$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$\lambda ^S=1$, 
$\hat {\varPi }_c=0.1\}$. The dotted lines are guides for the eye.
We now focus on the role played by the shear and dilatational viscosities separately. The influence of the shear and dilatational surface viscosities on the droplet stability has been studied recently by Singh & Narsimhan (Reference Singh and Narsimhan2020). They concluded from a perturbation theory that for capillary and Boussinesq numbers of order unity, dilatational surface viscosity increases the droplet deformation, while shear surface viscosity produces the opposite effect. Consequently, the dilatational viscosity is found to have a destabilizing impact, while the shear viscosity increases the droplet stability. A similar conclusion was obtained previously for a simple shear flow with both small (Narsimhan Reference Narsimhan2019) and arbitrary (Gounley et al. Reference Gounley, Boedec, Jaeger and Leonetti2016) droplet deformations. These studies did not consider the variation of the interfacial tension over the interface, therefore neither soluto-capillarity nor Marangoni convection was taken into account. To determine the combined influence on the droplet stability of the two effects mentioned above and shear/dilatational viscous stresses, we conducted simulations for 
$\lambda ^S=0$, 1 and 10
$^3$, i.e. when only shear viscosity is considered, dilatational viscosity is taken into account as well, and dilatational viscosity becomes dominant, respectively (figure 18). Both viscosities contribute to decreasing the rate at which oscillations are damped in the stable regime, but they hardly affect the critical capillary number. In fact, their effect is much smaller than that observed in the absence of soluto-capillarity and Marangoni convection (Singh & Narsimhan Reference Singh and Narsimhan2020). We observe a little stabilizing effect of the dilatational viscosity only for the case 
$\lambda ^S=10^3$.
Droplet deformation 
$D$ and imaginary part of the dominant eigenvalue 
$\omega _i^*$ as functions of the capillary number 
$C$ for 
$\lambda ^S=0$, 1 and 
$10^3$. The values of the rest of the governing parameters are 
$\{{\rm Re}=0$, 
$\lambda =0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$B_{s,{eq}}=1$, 
$\hat {\varPi }_c=0.1\}$. The dotted lines are guides for the eye.
Figure 19 compares the droplet shapes and surfactant distributions of two marginally stable droplets when the equilibrium surface coverage is reduced while keeping constant both the Marangoni number and the equilibrium Boussinesq number. For a thickening behaviour (the surface viscosities increase with the surfactant concentration), this is equivalent to reducing the amount the surfactant adsorbed onto the interface but increasing the strength and viscosities of the surfactant monolayer. In this case, the surfactant molecules accumulate in the droplet tip to a greater extent, leaving practically empty most of the droplet surface. In this way, the interfacial tension reduction is more localized in the droplet tip, which helps the droplet to adopt a pointed shape. This deformation is enhanced by the increase of the critical capillary number, which can be explained as follows. The droplet breaks up when the amount of surfactant convected to its tip is sufficiently large for the interfacial tension to fall below a certain value. As the surfactant surface concentration decreases, the capillary number needed to produce such convection increases, and so does its critical value. For instance, the critical capillary number for the configuration considered in figure 19 increases from 
$C=0.098$ to 0.113 when the surfactant concentration decreases from 
$\hat {\varGamma }_{{eq}}=0.5$ to 0.1. One can conclude that larger droplet deformations can be reached if the surface coverage is reduced while keeping the Marangoni and Boussinesq numbers constant. In the example mentioned above, the critical capillary number increases by around 15 % when the surfactant density is decreased by a factor of five. This indicates that the surfactant density cannot be considered as a dominant factor in the droplet stability for 
$\lambda =0.1$ when the Marangoni number and the equilibrium Boussinesq number are fixed. In fact, the droplet deformation calculated from the perturbation theory for inviscid surfactants (Vlahovska et al. Reference Vlahovska, Lawzdziewicz and Loewenberg2009) is independent of the surfactant density. The comparison between the results for 
$\hat {\varGamma }_{{eq}}=0.1$, 
$B_{s,{eq}}=0$ (figure 9) and 
$\hat {\varGamma }_{{eq}}=0.1$, 
$B_{s,{eq}}=1$ (figure 19) indicates that surface viscosity makes the linear interface displacement increase near the droplet tip.
(a,b) Surfactant distribution 
$\hat {\varGamma }$ (a) and droplet shape 
$F$ (b) for 
$\{\hat {\varGamma }_{{eq}}=0.1$, 
$C=0.113\}$ (black) and 
$\{\hat {\varGamma }_{{eq}}=0.5$, 
$C=0.0983\}$ (blue). (c,d) Droplet shape in the base flow (shaded area) and interface displacement due to the critical eigenmode, 
$(r_s,z_s)=(r_{s0},z_{s0})+\phi \, ({\rm Re}[\delta r_s],{\rm Re}[\delta z_s])$ (dashed lines). The values of the rest of the governing parameters are 
$\{Re=0$, 
$\lambda =0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$; 
$B_{s,{eq}}=1$, 
$\lambda ^S=1$, 
$\hat {\varPi }_c=0.1\}$.
Overall, surface viscosity is overshadowed by interface elasticity and has little effect on droplet deformation and its stability for the conditions considered in this work. For instance, we have verified that the velocity and pressure perturbation fields for 
$B_{s,{eq}}=1$ are very similar to those of the inviscid case. The influence of the surface viscosity on the interface perturbation is very small as well. Surface viscosity significantly affects only the damping rate of the dominant mode.
4.3. Tip streaming
We devote the rest of this paper to analysing the tip streaming arising for supercritical capillary numbers. To this end, we conducted direct numerical simulations that show the breakup process leading to the ejection of tiny droplets from the droplet tips. As mentioned in § 3, the simulations start from the steady solution for a subcritical case close to the stability limit.
Before considering the effect of a surfactant monolayer on the droplet's breakup mode, we present the evolution of a droplet with a clean interface. The capillary number is just above the critical one. As can be observed in figure 20, the breakup mechanism selected by the droplet is the so-called centre pinching mode, in which the mother droplet shrinks in its central part and ultimately produces large daughter droplets. This simulation indicates that the uniaxial extensional flow demands the presence of surfactant to produce tip streaming, at least for the parameter conditions considered in this figure. This result differs from that obtained by Zhang (Reference Zhang2004), who claimed that tip streaming could be produced with the uniaxial extensional flow in the absence of surfactant.
Shape of an inertialess droplet in the absence of surfactant at different instants for 
$\lambda =0.1$ and 
$C=0.175$.
Figure 21 shows the breakup of a surfactant-loaded droplet in the absence of viscous surface stresses for a capillary number slightly larger than the critical value 
$C\simeq 0.0986$ (figure 7). As observed by Eggleton et al. (Reference Eggleton, Tsai and Stebe2001), tip streaming arises at the poles of the deformed droplet during the last phase of its breakup. The surfactant convection caused by the outer flow overcomes the Marangoni convection produced by the gradient of surfactant concentration. Consequently, the surfactant accumulates in the droplet tip, which favours the fast ejection of an ultra-thin fluid thread (the diameter of the ejected droplet at 
$t=133.9$ is 0.03). As can be seen in the left-hand column of figure 21, the surfactant concentration increases at the droplet equator due to the surface compression in that region. The surfactant convection does not compensate for this effect because that parallel is a stagnation line.
Droplet shapes for 
$t=60$, 120, 130 and 133.9 (from top to bottom). The colours indicate the surfactant surface concentration 
$\hat {\varGamma }$ (a,d,g,j), interfacial tension 
$\hat {\sigma }$ (b,e,h,k), and tangential surface velocity 
$v_t$ (c,f,i,l). The values of the governing parameters are 
$\{{Re}=0$, 
$\lambda =0.1$, 
$C=0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$B_{s,{eq}}=0\}$.
Figure 22 compares the velocity 
$v_{{tip}}$ of the droplet tip calculated with two initial conditions: (i) a spherical droplet with a uniform surfactant distribution, and (ii) the steady solution for a subcritical case close to the stability limit (also considered in the rest of the simulations in this work). The droplet tip moves at practically the same speed regardless of the initial condition, even though the simulation corresponds to the highest Reynolds number 
${\rm Re}=10$ considered in this work. When the initial droplet shape is spherical, the tip is dragged by the outer flow, and its velocity takes values similar to those of that flow at the symmetry axis; i.e. 
$v_{{tip}}\sim C$. Then the tip slows down and the droplet approaches quasi-statically its maximum stable deformation 
$D=0.37$. It eventually exceeds this value, and then the droplet tip accelerates. The marginally stable shape can be regarded as an intermediate state adopted by the droplet before tip streaming arises. Comparison with the results shown in figure 23 indicates that the droplet tip remains practically still for much longer in the inertialess case.
Velocity 
$v_{{tip}}$ of the droplet tip as a function of the droplet deformation 
$D$ for 
$\{\rho =1$, 
${\rm Re}=10$, 
$\lambda =0.1$, 
$C=0.084$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$B_{s,{eq}}=0\}$. The circles and triangles correspond to simulations starting from a spherical droplet with a uniform surfactant distribution and from the steady solution for a subcritical case close to the stability limit, respectively.
Velocity 
$v_{{tip}}$ of the droplet tip as a function of the time 
$t$ (a) and the droplet deformation 
$D$ (b) for 
$B_{s,{eq}}=0$, 
$10^{-5}$, 
$10^{-4}$ and 
$10^{-3}$. The values of the rest of the governing parameters are 
$\{{\rm Re}=0$, 
$\lambda =0.1$, 
$C=0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$\lambda ^S=1$, 
$\hat {\varPi }_c=0.1\}$. The colours indicate the surfactant surface concentration 
$\hat {\varGamma }$.
4.4. Tip streaming: influence of surface viscosity
Surface viscosity opposes the large surface velocity gradient demanded by tip streaming. On the other hand, the accumulation of surfactant in the droplet tip increases the surface viscosity in that region. A natural question is whether the surfactant viscosity can inhibit or suppress the tip streaming phenomenon described in figure 21. To answer this question, we analyse in figures 23–28 four time-dependent simulation runs conducted for different values of the surface viscosities.
Figure 23 shows the velocity 
$v_{{tip}}$ of the droplet tip as a function of time 
$t$ in the absence of viscous surface stresses, for three values of the Boussinesq number 
$B_{s,{eq}}$. The function 
$v_{{tip}}(t)$ and therefore the tip position are essentially the same in all the cases analysed. The viscous surface stresses slightly delay the ejection process. The tip velocity becomes of the order of the characteristic outer fluid velocity 
$C$ once the ejected fluid thread has formed.
Figure 24 shows snapshots obtained from the four transient simulations. The surface viscosity affects the growth of the ejected liquid thread from deformations 
$D\gtrsim 0.6$, increasing the thread diameter during the rest of the liquid emission. Interestingly, the size of the ejected droplet depends significantly on the Boussinesq number: the higher the surface viscosity, the larger the droplet formed at the end of the ejected thread. Specifically, the ejected droplet diameter in figure 24(e) (measured as twice the maximum of the interface radius at the last simulation instant) is approximately 0.030, 0.039, 0.12, 0.22 and 0.32 for 
$B_{s,{eq}}=0$, 
$10^{-6}$, 
$10^{-5}$, 
$10^{-4}$ and 
$10^{-3}$, respectively. These values are plotted in figure 25. The liquid volume accumulated in the droplet tip increases up to three orders of magnitude for 
$B_{s,{eq}}=10^{-3}$. In fact, this fluid ejection may not be qualified as true tip streaming. It must be noted that this value of the Boussinesq number corresponds to relatively small surface viscosity. In fact, Boussinesq numbers of that order of magnitude can be obtained when a water drop 1 mm in radius is submerged in a liquid bath 10 mPa s in viscosity, and the interface is loaded with a surfactant of shear viscosity 
$10^{-8}$ Pa s m.
Droplet shape in the course (a–d) and at the end (e) of the transient simulation for 
$B_{s,{eq}}=0$, 
$10^{-5}$, 
$10^{-4}$ and 
$10^{-3}$. The values of the rest of the governing parameters are 
$\{{\rm Re}=0$, 
$\lambda =0.1$, 
$C=0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$\lambda ^S=1$, 
$\hat {\varPi }_c=0.1\}$. The colours in graph (e) indicate the surfactant surface concentration 
$\hat {\varGamma }$.
Droplet diameter 
$d$ at the end of the transient simulation for 
$B_{s,{eq}}=10^{-6}$, 
$10^{-5}$, 
$10^{-4}$ and 
$10^{-3}$. The values of the rest of the governing parameters are 
$\{{\rm Re}=0$; 
$\lambda =0.1$, 
$C=0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$\lambda ^S=1$, 
$\hat {\varPi }_c=0.1\}$.
As mentioned above, figure 24(a–d) shows that the droplet evolution is affected by the surface viscosity for 
$D\gtrsim 0.6$. A similar conclusion can be obtained from the evolution of the tip curvature (figure 26). This quantity remains practically independent from the surface viscosity until 
$D\simeq 0.55$, and then increases more rapidly with decreasing surface viscosity. In all the cases, the curvature reaches a maximum whilst the tip is still accelerating (see figure 23), and then decreases due to the inflation of the ejected droplet. The tip keeps on accelerating, dragged by the outer flow even after the inflation of the ejected droplet has begun.
Curvature 
$\kappa _{{tip}}$ of the droplet tip as a function of the droplet deformation 
$D$ for 
$B_{s,{eq}}=0$, 
$10^{-5}$, 
$10^{-4}$ and 
$10^{-3}$. The values of the rest of the governing parameters are 
$\{{\rm Re}=0$, 
$\lambda =0.1$, 
$C=0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$\lambda ^S=1$, 
$\hat {\varPi }_c=0.1\}$.
Now we analyse the effect of the dilatational surface viscosity. The droplet shape and surfactant distribution obtained for 
$\lambda ^S=0$ and 1 (figure 27a,b) are practically the same, which indicates that the shear surface viscosity has a much more noticeable influence on the tip streaming phenomenon than that produced by the dilatational one. In fact, the droplet behaviour is very similar to that of the inviscid surfactant if only the dilatational viscosity is accounted for (
$\lambda ^S=10^3$ and 
$B_{s,{eq}}=10^{-6}$). In other words, the tip streaming arising for an inviscid monolayer is suppressed by the shear viscosity. The qualitatively different roles played by the shear and dilatational viscosity are in contrast to what occurs in the dynamics of slender threads, where that difference is merely quantitative (Martínez-Calvo & Sevilla Reference Martínez-Calvo and Sevilla2018; Wee et al. Reference Wee, Wagoner, Kamat and Basaran2020).
Shape and surfactant distribution at 
$t=132.3$ for 
$\lambda ^S=0$ and 
$B_{s,{eq}}=10^{-3}$ (a), 
$\lambda ^S=1$ and 
$B_{s,{eq}}=10^{-3}$ (b), and 
$\lambda ^S=10^3$ and 
$B_{s,{eq}}=10^{-6}$ (c). The values of the rest of the governing parameters are 
$\{{\rm Re}=0$, 
$\lambda =0.1$; 
$C=0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$, 
$\hat {\varPi }_c=0.1\}$. The colours indicate the surfactant surface concentration 
$\hat {\varGamma }$.
Finally, we examine the effect of the surface stresses due to shear viscosity on the ejected droplet diameter. We have verified that the droplet diameter hardly changes when the term
\begin{equation} \left[B_s\left(\frac{(Fv_t)'}{F} + \kappa v_n\right)\right]' \end{equation}is ‘turned off’ in (2.10). On the contrary, switching off the terms
\begin{equation} B_s \bar{\kappa}\left[{-}F\left(\frac{v_t}{F}\right)'+ \bar{\kappa} v_n \right] \quad \text{or}\quad 2\left(B_s K - B_s'\frac{F'}{F}\right)v_t - 2\kappa_1(B_sv_n)' \end{equation}in (2.9) or (2.10), respectively, significantly affects the droplet size. When the terms (4.2) are turned off simultaneously, the resulting tip streaming is very similar to that occurring in the inviscid case. Therefore, these normal and tangential terms are responsible for the increase in the droplet diameter due to shear surface viscosity.
Figure 28 allows us to understand why the size of the droplet tip is affected by tangential stress caused by the shear surface viscosity despite the small value of the Boussinesq number. The tangential surface stress (2.10) takes very small values over the interface except at the droplet apex and the two ends of the liquid thread connecting the parent and ejected droplet. In fact, the magnitude of that stress reaches its maximum value right in front of the ejected droplet. This stress distribution slightly reduces the surfactant accumulated in the droplet tip. However, this small reduction leads to a large increase in the interfacial tension because the concentration in that region is close to the maximum packing density. The increase of the interfacial tension in the droplet tip makes the interface curvature decrease in that region. Therefore, viscous surface stress changes the shape of the droplet tip by altering the balance of surface stresses and the distribution of surfactant over the interface. This phenomenon resembles what occurs during the pinch-off of an interface covered with a surfactant monolayer, where surface viscosity becomes relevant because of its considerable influence on the transport of surfactants over the interface (Ponce-Torres et al. Reference Ponce-Torres, Montanero, Herrada, Vega and Vega2017).
Surface tension 
$\hat {\sigma }$ (a), surface coverage 
$\hat {\varGamma }$ (b), surface shear stress 
$\tau _{ts}^S$ (c), surface velocity 
$v_t$ (d), and interface position 
$F$ (e), as functions of the axial coordinate 
$z$ for the Boussinesq numbers at equilibrium 
$B_{s,{eq}}=0$ and 
$10^{-3}$. The values of the rest of the governing parameters are 
$\{{\rm Re}=0$, 
$\lambda =0.1$, 
$C=0.1$; 
${Pe}^S=10^3$, 
$Ma=0.2$, 
$\hat {\varGamma }_{{eq}}=0.5$; 
$\lambda ^S=1$, 
$\hat {\varPi }_c=0.1\}$.
5. Conclusions
We studied numerically the steady deformation and breakup of a droplet covered with an insoluble surfactant in a uniaxial extensional flow, focusing our attention on the role of the surface viscosities in both the subcritical and supercritical regimes. We considered the full hydrodynamic model, which comprises arbitrary large droplet deformations, a variation of the interfacial tension over the droplet surface (soluto-capillarity and Marangoni convection effects), and both the droplet and outer bath inertia.
In all the cases analysed, the frequency and damping rate of the linear mode responsible for instability vanish at the critical capillary number. This implies that a stationary linear perturbation destabilizes the flow. The surfactant monolayer considerably reduces the interval of the capillary number for which a steady droplet deformation can be produced. Neither the droplet deformation nor the stability is significantly affected by the droplet viscosity within that interval. In the absence of surfactant, the critical capillary number increases, and the effect of the droplet viscosity becomes noticeable close to the marginal stability. In addition, inertia significantly increases the droplet deformation and decreases the critical capillary number for Reynolds numbers of the order of 10.
Our results consistently show that neither the droplet deformation nor the stability is significantly affected by the surface viscosities, even for moderately viscous (
$B_{s,{eq}}=1$) and very viscous (
$B_{s,{eq}}=10$) surfactants. However, viscous surface stresses considerably reduce the damping rate of the dominant mode. This reduction does not essentially alter the critical capillary number. We conclude that the soluto-capillarity and Marangoni convection affect the droplet stability more significantly than the surface viscosities. Marangoni convection opposes the external flow. The velocity considerably decreases on the interface and therefore inside the droplet. This essentially explains why both the droplet and monolayer viscosities play a secondary role in the droplet's steady deformation and stability.
Our simulations confirm that surfactant is a crucial ingredient to produce tip streaming. The interface deformation in the linear regime affects the entire drop and is qualitatively the same as that observed when a clean droplet breaks up following the central pinching mode. However, the nonlinear evolution is drastically altered by the accumulation of surfactant molecules in the droplet tip. The small local scale characterizing tip streaming is fixed in this phase of the droplet bursting.
One of the main conclusions of this work is that shear surface viscosity considerably changes the shape and size of the droplet tip in the final stage of tip streaming. Interestingly, this occurs even for very small surface viscosities, such as that of SDS, a surfactant commonly used in experiments. We have explained this phenomenon in terms of the influence of the viscous surface stress on both the balance of normal interfacial stresses and the surfactant transport over the interface when the ejected droplet is formed. Specifically, the shear viscosity slightly reduces the amount of surfactant in the droplet tip, which entails a sharp increase in the interfacial tension for the parameter conditions considered in our simulations. The effect of the dilatational viscosity is much less noticeable.
The present work can be extended in several directions. Probably the most interesting one is to study the effects on tip streaming of factors not considered in § 4.3, such as inertia, bulk viscosity ratio, Marangoni number and surfactant concentration. A systematic analysis of these effects requires an enormous computational effort, given the computing time consumed by each transient simulation. Our simulations correspond to configurations with 
${Pe}^S\,{Ma}/B_{s,{eq}}=a^2/(\mu _s D^S)\gg 1$, which explains why the surface viscosity plays a subdominant role vs Marangoni convection in both the droplet deformation and stability. One can expect surface viscosity to become relevant as the droplet diameter decreases. The analysis of the behaviour of micrometre droplets constitutes another interesting extension of the present study.
The Langmuir equation of state (2.1) gives an unrealistic sharp reduction of the interfacial tension for surfactant concentrations close to the maximum packing density 
$\varGamma _{\infty }$. In fact, the interfacial tension becomes negative for 
$(\varGamma _{\infty }-\varGamma )/\varGamma _{\infty }<\textrm {e}^{-1/{Ma}^*}$, where 
$Ma^*=\varGamma _{\infty }RT/\sigma _0$ is the Marangoni number defined in terms of the surface tension 
$\sigma _0$ of the clean interface. However, experiments show that the surface tension reaches a plateau as 
$\varGamma \to \varGamma _{\infty }$, which is commensurate with 
$\sigma _0$. The pointed shape of the droplet tip for 
$B_{s,{eq}}=0$ can be attributed partially to the unrealistic behaviour of the Langmuir equation in this limit. It may be of particular interest to determine the influence of the surface tension minimum value on tip streaming for supercritical capillary numbers. This element may help us to understand the discrepancies between theoretical predictions and experiments (Eggleton et al. Reference Eggleton, Tsai and Stebe2001).
Funding
This research has been supported by the Spanish Ministry of Economy, Industry and Competitiveness under grants DPI2016-78887 and PID2019-108278RB, by Junta de Extremadura under grant GR18175, and by Junta de Andalucía under grant P18-FR-3623.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Calculation of viscous surface stresses
In this Appendix, we derive the interfacial stresses produced by the surfactant monolayer using the interface intrinsic (local) coordinate system. We consider an axisymmetric column of fluid with axis 
$z$ and radius 
$h(z)$. We choose the cylindrical coordinate system 
$(x^1,x^2,x^3) = (z,\theta,r)$ with scale factors 
$h_1=h_3=1$, 
$h_2 = r$. The metric tensor is diagonal with 
$g_{jj} = h_j^2$ and 
$g^{jj} = h_j^{-2}$. The determinant is 
$g = g_{ii} = r^2$. The non-vanishing Cristoffel symbols are 
$\varGamma ^3_{22} = -x^3 \equiv -r$ and 
$\varGamma ^2_{32} = \varGamma ^2_{23} = 1/x^3 \equiv 1/r$.
As the local coordinates for the surface we choose 
$(u^1,u^2) = (s,\theta )$, where 
$s$ is the arc length. Let 
$z = g(s)$; it then follows that 
$g'^2 + h'^2 = 1$. We will always denote the derivative with respect to the arc length by a prime. Now the Cartesian coordinates of the surface are 
${\boldsymbol y} = {\boldsymbol y}(z,\theta ) = (h(z)\cos \theta,h(z)\sin \theta,z)$, and in cylindrical components, 
$x^i(s,\theta ) = (g(s),\theta,h(s))$. As a result, the tangent vectors are 
$t_1^i = (g',0,h')$ and 
$t_2^i = (0,1,0)$. Since 
$a_{\alpha \beta } = g_{ij}t_{\alpha }^it_{\beta }^j$, we find that 
$a_{\alpha \beta }$ is diagonal with 
$a_{11} = 1$ and 
$a_{22} = h^2$, and determinant 
$a = h^2$, so that 
$a^{11} = 1$ and 
$a^{22} = 1/h^2$.
The 
$\epsilon$-tensors are 
$\varepsilon _{ijk} = \sqrt {g}\epsilon _{ijk}$ and 
$\varepsilon ^{\alpha \beta } = \epsilon ^{\alpha \beta }/\sqrt {a}$, so that
\begin{equation} n_i = \tfrac{1}{2}\varepsilon^{\alpha\beta}\varepsilon_{ijk} t_{\alpha}^j t_{\beta}^k = ({-}h',0,g') \end{equation}
is the normal. The second fundamental form is diagonal with 
$b_{11} = g'h''- h'g''$ and 
$b_{22} = -h g'$. The eigenvalues of 
$b_{\alpha }^{\beta } = b_{\alpha \gamma }a^{\gamma \beta }$ are the principal curvatures. Putting
\begin{equation} b_{\alpha}^{\beta} = \begin{pmatrix} g' h'' - h' g'' & 0 \\ 0 & -\dfrac{g'}{h} \end{pmatrix} \equiv \begin{pmatrix} -\kappa_2 & 0 \\ 0 & -\kappa_1 \end{pmatrix}, \end{equation}we have
\begin{equation} \kappa_1 = \frac{g'}{h}, \quad \kappa_2 = h'g'' - g' h'' \end{equation}and
\begin{equation} 2H ={-}\kappa_1-\kappa_2 = g'h'' - h'g'' - \frac{g'}{h}, \quad K = \kappa_1\kappa_2 = \frac{g'(h'g'' - g'h'')}{h}, \end{equation}
where 
$H$ and 
$K$ are the mean and Gaussian curvatures, respectively.
Let 
${\boldsymbol n}$ be the normal that points from inside the fluid to the outer phase. Let 
${\boldsymbol F}$ be the force on the interface (Scriven Reference Scriven1960; Aris Reference Aris1962), so 
${\boldsymbol F} = - {\boldsymbol n}\boldsymbol {\cdot }\boldsymbol {\sigma }$. Assuming negligible inertia of the interface, the covariant version for the force 
${\boldsymbol F}$ is of the form
\begin{equation} -F^i = t_{\alpha}^i T_t^{\alpha} + n^i T_n = t_1^i T_t^1 + n^i T_n, \end{equation}
using the fact that 
$F^i$ can have tangential components only in the 
$t_1^i$-direction. Thus, using 
${\boldsymbol t} = (1,0,h')/\sqrt {1+h_z^2} = t_1^i$, we obtain
\begin{equation} {\boldsymbol n}\boldsymbol{\cdot}\boldsymbol{\sigma}\boldsymbol{\cdot}{\boldsymbol n} = T_n, \quad {\boldsymbol n}\boldsymbol{\cdot}\boldsymbol{\sigma}\boldsymbol{\cdot}{\boldsymbol t} = T_t^1, \end{equation}where (Aris Reference Aris1962)
\begin{equation} T_n = N_1 + N_2 + N_3 \equiv 2 H \sigma + 2H(\mu_1^s + \mu_2^s)t_{\lambda}^j a^{\lambda\mu}W_{j,\mu} + 2\mu_1^s t_{\lambda}^j\varepsilon^{\lambda\alpha}\varepsilon^{\beta\mu} b_{\alpha\beta}W_{j,\mu} \end{equation}and
\begin{align} T_t^1&=T_1 + T_2 + T_3 + T_4 + T_5 + T_6\nonumber\\ &\equiv a^{1\beta}\sigma_{,\beta} + a^{1\beta} \left((\mu_1^s + \mu_2^s)a^{\lambda\mu}t_{\lambda}^j W_{j,\mu}\right)_{,\beta} + 2\mu_1^sK a^{1\beta} t_{\beta}^j W_j - \mu_1^s\varepsilon^{1\beta} \left(\varepsilon^{\lambda\mu}(t_{\mu}^jW_j)_{,\lambda}\right)_{,\beta}\nonumber\\ &\quad-2\mu_1^s\varepsilon^{1\lambda}b_{\lambda\beta}\varepsilon^{\beta\mu} (n^jW_j)_{,\mu} - \mu_{1,\beta}^s\epsilon^{\alpha\lambda}\epsilon^{\beta\mu} t_{\alpha}^i\left(t_{\lambda}^j W_{j,\mu}+t_{\mu}^j W_{j,\lambda}\right). \end{align}
Here 
$\mu _1^s$ and 
$\mu _2^s$ are the shear and dilatational surface viscosities, respectively (called 
$\epsilon$ and 
$\kappa$ by Scriven (Reference Scriven1960)), and 
$\sigma$ is the surface tension. We have generalized the expressions to allow surface viscosities to depend on space. As a result, there are extra terms in 
$T^{\alpha \beta }_{,\beta }$, so 
$T_2$ is modified, and there is an additional term 
$T_6$.
The velocity in cylindrical coordinates is 
$W_i = (v_z,0,v_r)$ and 
$W_{j,\mu } = \partial W_j/\partial u^{\mu }-\varGamma ^l_{jk} W_lt^k_{\mu }$, resulting in 
$W_{j,1} = (v_z',0,v_r'),$ 
$W_{j,2} = (0,\varGamma _{22}^3 W_3 t_2^2,0) = (0,h v_r,0)$. Beginning with the normal balance, 
$N_1 = 2H \sigma$ is the usual surface tension term. Now
\begin{equation} t_{\lambda}^j a^{\lambda\mu}W_{j,\mu} = t_1^jW_{j,1}+\frac{t_2^j}{h^2}\,W_{j,2} = g'v_z' + h' v_r' + \frac{v_r}{h}, \end{equation}so
\begin{equation} N_2 = 2H(\mu_1^s+\mu_2^s)\left[g'v_z' + h' v_r' + \frac{v_r}{h}\right]. \end{equation}
Further, to calculate 
$N_3$, note that
\begin{equation} \varepsilon^{\lambda\alpha}\varepsilon^{\beta\mu}b_{\alpha\beta} = \begin{pmatrix} -b_{22}/a & 0 \\ 0 & -b_{11}/a \end{pmatrix} = \begin{pmatrix} \kappa_1 & 0 \\ 0 & \dfrac{\kappa_2}{h^2} \end{pmatrix} \end{equation}and
\begin{equation} t_1^j W_{j,1} = g'v_z' + h'v_r', \quad t_2^j W_{j,2} = hv_r, \end{equation}so that
\begin{equation} N_3 = 2\mu_1^s\left[\kappa_1(g'v_z' + h'v_r') + \frac{\kappa_2}{h}v_r\right] . \end{equation}The terms from the tangential balance are
\begin{equation} T_1 = a^{1\beta}\sigma_{,\beta} = \sigma' \end{equation}and
\begin{align} T_2& = a^{11}\left[(\mu_1^s+\mu_2^s) \left(a^{\lambda\mu}t_{\lambda}^jW_{j,\mu}\right)\right]_{,1} = \left[(\mu_1^s+\mu_2^s)\left(a^{11}t_1^jW_{j,1} + a^{22}t_2^jW_{j,2}\right)\right]'\nonumber\\ &= \left[(\mu_1^s+\mu_2^s)\left(g'v_z' + h'v_r' + \frac{v_r}{h}\right)\right]'. \end{align}Further,
\begin{equation} T_3 = 2\mu_1 K a^{11} t_1^jW_j = 2\mu_1 K a^{11} t_1^jW_j = 2\mu_1 K(g'v_z + h' v_r). \end{equation} The curl vanishes, and we have 
$T_4 = 0$. Next, we have 
$n^j = n_j$, so
\begin{equation} v_n = n^j W_j ={-}h' v_z + g'v_r. \end{equation}We have
\begin{equation} \varepsilon^{1\lambda}\varepsilon^{\beta\mu}b_{\lambda\beta} = \varepsilon^{12}\varepsilon^{21}b_{22} = \frac{g'}{h}, \end{equation}so
\begin{equation} T_5 ={-}2\mu_1^s\,\frac{g'}{h}\,v_n'. \end{equation}Finally, the new term is
\begin{equation} T_6 ={-}(\mu_1^s)'\varepsilon^{1\lambda}\varepsilon^{1\mu} \left(t_{\lambda}^jW_{j,\mu}+t_{\mu}^jW_{j,\lambda}\right) ={-}2(\mu_1^s)'\,\frac{t_2^j}{h^2}\,W_{j,2} ={-}2(\mu_1^s)'\,\frac{v_r}{h}. \end{equation}Taking into account the above results, the boundary condition for the normal component is
\begin{equation} -{\boldsymbol n}\boldsymbol{\cdot}\boldsymbol{\sigma}\boldsymbol{\cdot}{\boldsymbol n} = \kappa\sigma + \mu_1^s \bar{\kappa}\left[{-}g'v_z' - h'v_r' + \frac{v_r}{h}\right] + \mu_2^s\kappa\left[g'v_z'+ h'v_r' + \frac{v_r}{h}\right], \end{equation}
where 
$\kappa = \kappa _1 + \kappa _2$ and 
$\bar {\kappa } = \kappa _1 - \kappa _2$.
We have
\begin{equation} v_n ={-}h'v_z + g' v_r, \quad v_t = g'v_z + h' v_r \end{equation}for the normal and tangential velocities, respectively, so
\begin{equation} v_z = g'v_t - h' v_n, \quad v_r = h'v_t + g' v_n. \end{equation}From this, it follows that
\begin{equation} v_t' = g' v_z' + h' v_r' - \kappa_2 v_n, \end{equation}so that
\begin{equation} {-{\boldsymbol n}}\boldsymbol{\cdot}\boldsymbol{\sigma}\boldsymbol{\cdot}{\boldsymbol n} = \kappa\sigma + \mu_1^s \bar{\kappa}\left[{-}h\left(\frac{v_t}{h}\right)' + \bar{\kappa} v_n \right] + \mu_2^s\kappa\left[\frac{(hv_t)'}{h} + \kappa v_n\right]. \end{equation}For the tangential component, we have
\begin{align} {\boldsymbol n}\boldsymbol{\cdot}\boldsymbol{\sigma}\boldsymbol{\cdot}{\boldsymbol t}&=\sigma' + \left[\mu_1^s\left(\frac{(hv_t)'}{h} + \kappa v_n\right)\right]' + 2\left(\mu_1^s K - (\mu_1^s)'\,\frac{h'}{h}\right)v_t - 2\kappa_1(\mu_1^sv_n)'\nonumber\\ &\quad +\left[\mu_2^s\left(\frac{(hv_t)'}{h} + \kappa v_n\right)\right]'. \end{align}The dimensionless forms of (A25) and (A26) are (2.9) and (2.10), respectively.
Open access
