1. Introduction
The capillary breakup of liquid films into droplets plays a major role in various unsteady fluid fragmentation phenomena, such as atomisation or droplet splashing (Villermaux Reference Villermaux2020). It is therefore unsurprising that the sheet retraction dynamics has received significant interest since the initial experimental observations by Dupre (Reference Dupre1867) and Rayleigh (Reference Rayleigh1879). Capillary film retraction follows the puncturing of a static liquid film, where a capillary-induced flow drives the opening of the interface (a hole) at a constant speed, i.e. the Taylor–Culick velocity (Taylor Reference Taylor1959; Culick Reference Culick1960). Early pictures by Rayleigh (Reference Rayleigh1891) and Ranz (Reference Ranz1959) show that the opening of the sheet results in the formation of a rim with roughly cylindrical caps, as the film moves away from the puncture. The desire to understand the fundamental mechanism underlying the development of the hole-driven expansion has led to numerous studies in this field, see for instance Debrégeas, Martin & Brochard-Wyart (Reference Debrégeas, Martin and Brochard-Wyart1995); Keller, King & Ting (Reference Keller, King and Ting1995); Brenner & Gueyffier (Reference Brenner and Gueyffier1999); Fullana & Zaleski (Reference Fullana and Zaleski1999); Sünderhauf, Raszillier & Durst (Reference Sünderhauf, Raszillier and Durst2002); Roisman, Horvat & Tropea (Reference Roisman, Horvat and Tropea2006); Savva & Bush (Reference Savva and Bush2009); Krechetnikov (Reference Krechetnikov2010); Gordillo et al. (Reference Gordillo, Agbaglah, Duchemin and Josserand2011); Lhuissier & Villermaux (Reference Lhuissier and Villermaux2011); Villermaux & Bossa (Reference Villermaux and Bossa2011) and Agbaglah, Josserand & Zaleski (Reference Agbaglah, Josserand and Zaleski2013); Agbaglah (Reference Agbaglah2021). Debrégeas et al. (Reference Debrégeas, Martin and Brochard-Wyart1995) and Savva & Bush (Reference Savva and Bush2009) concluded that the retracting dynamics of a liquid sheet is governed by a balance between inertial, viscous and surface-tension forces. Accordingly, the Ohnesorge number, 
$Oh$, (i.e. ratio of viscous to capillary forces) becomes the most appropriate parameter to parametrise the flow dynamics. Different regimes have been identified that depend on 
$Oh$: (i) in the 
$Oh< 1$ regime, which is characterised by a nearly inviscid flow, the dynamics is driven by surface tension, leading to the formation of capillary waves ahead of the roughly cylindrical rim; (ii) in the 
$Oh>10$ regime, viscous forces dominate the dynamics causing the suppression of the rim; and (iii) an intermediate regime bridges the two previous regimes.
The ground-breaking experiments of McEntee & Mysels (Reference McEntee and Mysels1969) proposed that a surface instability along the rim drives the formation of ligaments, and eventually, the ejection of droplets. Several studies have suggested that different physical mechanisms prompt the detachment of these droplets: a Rayleigh–Plateau (RP) instability (discussed below), a nonlinear amplification mechanism (Yarin & Weiss Reference Yarin and Weiss1995), a Richtmyer–Meshkov instability and a competition/collaboration between the Richtmyer–Meshkov and Rayleigh–Taylor (RT) instabilities (Krechetnikov Reference Krechetnikov2010).
However, Bremond & Villermaux (Reference Bremond and Villermaux2006) (for the fragmentation of the lamella of colliding jets), Rieber & Frohn (Reference Rieber and Frohn1999) and Zhang et al. (Reference Zhang, Brunet, Eggers and Deegan2010) (for the fragmentation of the sheet crown from impacts of drops onto a thin layer of liquid) concluded that the RP instability is responsible for the onset of drop detachment. These findings contrast with the previous study by Fullana & Zaleski (Reference Fullana and Zaleski1999), that suggests that the growth of the rim during the retraction prevents the onset of a RP instability. Finally, Agbaglah et al. (Reference Agbaglah, Josserand and Zaleski2013) (for a numerical sheet retraction in a frame of reference moving at the Taylor–Culick velocity) and Wang et al. (Reference Wang, Dandekar, Bustos, Poulain and Bourouiba2018) (for drop impact onto a solid surface) showed that both the RT and RP instabilities are responsible for droplet generation from the rim. The former is important at early times during the deceleration of the ejected sheet, whereas the latter is predominant at longer times once the rim has already formed. According to recent studies by Wang & Bourouiba (Reference Wang and Bourouiba2018, Reference Wang and Bourouiba2021), the ligament dynamics can result in end pinching, break up into satellite droplets or recombination into a single fluid volume.
McEntee & Mysels (Reference McEntee and Mysels1969) used surfactants in their experimental work to study the bursting of liquid films without directly endorsing the effects of Marangoni-induced flows. In contrast, recent studies have demonstrated the crucial role of surfactants on the dynamics of the capillary singularity (Ambravaneswaran, Phillips & Basaran Reference Ambravaneswaran, Phillips and Basaran2000; Craster, Matar & Papageorgiou Reference Craster, Matar and Papageorgiou2002; Timmermans & Lister Reference Timmermans and Lister2002; Liao, Franses & Basaran Reference Liao, Franses and Basaran2006; Kamat et al. Reference Kamat, Wagoner, Thete and Basaran2018). Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020) and Kamat et al. (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2020) showed that surfactant-induced Marangoni stresses inhibit end pinching for retracting liquid threads via the suppression of stagnation points, which in turn leads to flow reversal near the vicinity of the neck.
Three-dimensional numerical simulations of retracting liquid sheets are scarce in the literature due to the need for large aspect ratios to induce both the growth of instabilities at the rim, and the development of the nonlinear flow dynamics. In this work, the role of surfactant-induced Marangoni stresses on the retraction of thin liquid sheets and the detachment of droplets is studied in a three-dimensional, nonlinear framework. This paper is structured as follows: § 2 introduces the numerical method, governing dimensionless parameters, problem configuration and validation; § 3 provides the results from the simulations; and concluding remarks are given in § 4.
2. Problem formulation and numerical method
Numerical simulations were performed by solving the transient two-phase Navier–Stokes equations in a three-dimensional Cartesian domain 
$\boldsymbol {x} = (x, y, z )$ (see figure 1). The interface is captured via a hybrid front-tracking/level-set method and a convective–diffusion equation is solved for the surfactant transport along the deforming interface (Shin et al. Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018). All variables are made dimensionless following:
\begin{equation} \tilde{\boldsymbol{x}}=\frac{\boldsymbol{x}}{h_0}, \quad \tilde{t}=\frac{t}{t_{inv}}, \quad \tilde{\boldsymbol{u}}=\frac{\boldsymbol{u}} {u_{TC}}, \quad \tilde{p}=\frac{p}{\rho_l u_{TC}^{2}}, \quad \tilde{\sigma}=\frac{\sigma}{\sigma_s}, \quad \tilde{\varGamma}=\frac{\varGamma}{\varGamma_\infty}, \end{equation}
where, 
$t$, 
$\boldsymbol {u}$ and 
$p$ stand for time, velocity and pressure, respectively. The physical parameters correspond to the liquid density 
$\rho _l$, viscosity 
$\mu _l$, surface tension 
$\sigma$, surfactant-free surface tension 
$\sigma _s$ and 
$u_{TC}=\sqrt {2 \sigma _s/ (\rho h_0)}$, is the well-known Taylor–Culick speed; here, 
$h_0$ stands for the sheet thickness. In this study we focus on the low-
$Oh$ regime; hence, the characteristic time scale is given by the capillary time 
$t_{inv}= h_0/u_{TC}=\sqrt {\rho _l h_0^{3}/ (2 \sigma _s)}$. As seen, the interfacial surfactant concentration, 
$\varGamma$, is also made dimensionless through the saturation interfacial concentration, 
$\varGamma _{\infty }$. As a result of the scaling in (2.1a–f), the dimensionless forms of the governing equations for the flow and the surfactant transport are respectively expressed as
\begin{gather} {\boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}}=0,} \end{gather}
\begin{gather} {\tilde{\rho} \left(\frac{\partial \tilde{\boldsymbol{u}}}{\partial \tilde{t}}+\tilde{\boldsymbol{u}} \boldsymbol{\cdot}\boldsymbol{\nabla} \tilde{\boldsymbol{u}}\right) + \boldsymbol{\nabla} \tilde{p} = Oh ~ \boldsymbol{\nabla}\boldsymbol{\cdot} \left [ \tilde{\mu} (\boldsymbol{\nabla} \tilde{\boldsymbol{u}} +\boldsymbol{\nabla} \tilde{\boldsymbol{u}}^{\rm T}) \right ] + \int _{\tilde{A}\widetilde{(t)}} \left(\tilde{\sigma} \tilde{\kappa} \boldsymbol{n} + \boldsymbol{\nabla}_s \tilde{\sigma} \right) \delta \left(\tilde{\boldsymbol{x}} - \tilde{\boldsymbol{x}}_{_f} \right)\,\text{d}\tilde{A},} \end{gather}
\begin{gather} {\frac{\partial \tilde{\varGamma}}{\partial \tilde{t}}+\boldsymbol{\nabla}_s \boldsymbol{\cdot} (\tilde{\varGamma} \tilde{\boldsymbol{u}} _{\rm t}) =\frac{1}{Pe_s} \nabla^{2}_s \tilde{\varGamma},} \end{gather}
which correspond to the equations of mass and momentum conservation and the convective–diffusion equations for the interfacial concentration, respectively. Here, the density 
$\tilde {\rho }$ and viscosity 
$\tilde {\mu }$ are defined by 
$\tilde {\rho }=\rho _g/\rho _l + (1 -\rho _g/\rho _l) \mathcal {H}(\tilde {\boldsymbol {x}},\tilde {t})$ and 
$\tilde {\mu }=\mu _g/\mu _l+ (1 -\mu _g/\mu _l) \mathcal {H}( \tilde {\boldsymbol {x}},\tilde {t})$ wherein 
$\mathcal {H}(\tilde {\boldsymbol {x}},\tilde {t})$ represents a smoothed Heaviside function. In this work, 
$\mathcal {H}(\tilde {\boldsymbol {x}},\tilde {t})$ is zero in the gas phase and unity in the liquid. The subscript 
$g$ designates the gas phase, 
$\tilde {\boldsymbol {u}}_{{t}}= (\tilde {\boldsymbol {u}}_{{s}} \boldsymbol {\cdot} \boldsymbol {t} ) \boldsymbol {t}$ stands for the velocity vector tangential to the interface in which 
$\tilde {\boldsymbol {u}}_{{s}}$ corresponds to the interfacial velocity, 
$\kappa$ stands for the interfacial curvature obtained from the Lagrangian interface structure and 
$\boldsymbol {\nabla }_s=({\boldsymbol{\mathsf{I}}}-\boldsymbol {n}\boldsymbol {n})\boldsymbol {\cdot} \boldsymbol {\nabla }$ stands for the surface gradient operator wherein 
$\boldsymbol{\mathsf{I}}$ is the identity tensor and 
$\boldsymbol {n}$ is the outward-pointing unit normal to the interface. Finally, 
$\tilde {\boldsymbol {x}}_f$ is the parameterisation of the interface 
$\tilde {A} (\tilde {t})$, and 
$\delta$ represents a Dirac delta function that is non-zero when 
$\tilde {\boldsymbol {x}}=\tilde {\boldsymbol {x}}_f$ only.
(a) Schematic view of the sheet and rim geometry: the initial shape of the retracting sheet (not to scale) is characterised by a cylindrical rim with 
$e=0.2$ and an initial instability amplitude of 
$\epsilon =0.25$. The black line represents the arclength, 
$s$, across the 
$x$–
$z$ plane 
$(y=\lambda /2)$. (b) Presents the numerical results for the retraction velocity in terms of the Taylor–Culick value 
$u_{TC}$ for a surfactant-free liquid sheet with 
$Oh=0.1$, 
$\epsilon =0$ and 
$e=2$.
The dimensionless groups that govern the dynamics are defined as
\begin{equation} Oh=\frac{\mu_l}{\sqrt{\sigma_s h_0 \rho_l}}, \quad Pe_s=\frac{ u_{TC} h_0}{D_s}, \quad \beta_s= \frac{\mathrm{Re} T \varGamma_\infty}{\sigma_s}, \end{equation}
where 
$Oh$ and 
$Pe_s$ stand for the Ohnesorge and (interfacial) Péclet numbers, respectively, while 
$\beta _s$ is the surfactant elasticity number which provides a measure of the sensitivity of 
$\sigma$ to 
$\varGamma$; here, 
$\mathrm {Re}$ is the ideal gas constant value (
$\mathrm {Re} = 8.314$ J K
$^{-1}$ mol
$^{-1}$), 
$T$ denotes temperature and 
$D_s$ stands for the diffusion coefficient. The nonlinear Langmuir equation expresses 
$\sigma$ in terms of 
$\varGamma$, i.e. 
$\tilde {\sigma }=1 + \beta _s \ln { (1 -\tilde {\varGamma })}$, and the Marangoni stress, 
$\tilde {\tau }$, is given as a function of 
$\tilde {\varGamma }$ as 
$\tilde {\tau }= \boldsymbol {\nabla }_s \tilde {\sigma }\boldsymbol {\cdot} {\boldsymbol {t}} =-\beta _s/ (1-\tilde {\varGamma })\boldsymbol {\nabla }_s\tilde {\varGamma }\boldsymbol {\cdot} {\boldsymbol {t}}$. Tildes are dropped henceforth.
We present a brief summary of the numerical method used in this study. The Navier–Stokes equations are solved by a finite volume method on a staggered grid (Harlow & Welch Reference Harlow and Welch1965). The computational domain is discretised by a fixed regular grid (i.e. Eulerian grid) and the spatial derivatives are approximated by standard centred difference discretisation, except for the nonlinear term, which makes use of a second-order essentially non-oscillatory scheme (Sussman, Smereka & Osher Reference Sussman, Smereka and Osher1994). The interface is tracked explicitly by an additional Lagrangian grid by using the front-tracking method, and the interface is reconstructed by using the level contour reconstruction method Shin & Juric (Reference Shin and Juric2002, Reference Shin and Juric2009). Communication between the Eulerian and Lagrangian grids (for the transfer of the geometric information of the interface) is done by using the discrete delta function and the immersed boundary method of Peskin (Reference Peskin1977). The advection of the Lagrangian interface is done by integrating 
${\rm d}\boldsymbol {x}_f /{\rm d}t = \boldsymbol {V}$ with a second-order Runge–Kutta method, where 
$\boldsymbol {V}$ stands for the interfacial velocity which has been calculated by interpolation from the Eulerian velocity. The code is parallelised through an algebraic domain-decomposition technique and the communication between subdomains for data exchange is managed by the message passing interface protocol. More details of the numerical method used to solve the above equations is described in detail in Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018)
2.1. Numerical set-up and validation
Figure 1(a) illustrates the geometry used in this study; we have initialised the sheet film dynamics using a linearly unstable configuration, as proposed by Agbaglah et al. (Reference Agbaglah, Josserand and Zaleski2013). Thus, we have considered a thin liquid sheet of thickness 
$h_0$, and initial length 
$L_0$, connected to a cylindrical rim of radius 
$a_0$. The initial shape of the rim is given by 
$a(\omega,\varepsilon )=a_0 [1+\varepsilon \cos ( \omega z ) ]$, where 
$\varepsilon$ and 
$\omega$ stand for the perturbation amplitude and the growth rate, respectively. The dispersion relation for the RP instability suggests that the largest growth rate, 
$\omega$, is obtained for a wavenumber 
${k} \sim 0.6970$, which corresponds to 
$\lambda _{max}=2{\rm \pi} /{k} =9.016 a_0$ (where 
$\lambda _{max}$ is the most unstable wavelength). All simulations were run using an initial perturbation amplitude of 
$\varepsilon =0.25$ with 
${k} \sim 0.6970$, unless stated otherwise in the text. See the supplementary material (SM) available at https://doi.org/10.1017/jfm.2022.768 for further discussion regarding the selection of 
$\varepsilon$.
Additionally to the dimensionless flow parameters, the flow dynamics is also controlled by the ratio of the film thickness to the radius of the rim, i.e. 
$e=h_0/a_0$. In this study, we consider the zero acceleration limit, and 
$e=0.2$, following the work by Agbaglah et al. (Reference Agbaglah, Josserand and Zaleski2013), where their stability analysis and dispersion relationship depended on the initial acceleration (see SM for additional information about the choice of 
$e$).
As the liquid sheet retracts, the thickness of the rim grows over time as the rim engulfs the sheet, 
$e$ decreases and the rim deceleration vanishes, resulting in the rim dynamics being fully driven by the RP instability. Therefore, the flow dynamics depends on two competing time scales: the time scale given by the RP instability in the (spanwise) 
$z$-direction (
$T_{RP}$), and the time scales of the rim growth (
$T_{ret}$). The time scale of the action of surface tension and droplet detachment is given by the capillary time, 
$T_{RP}=\sqrt {\rho _l a_0^{3}/ \sigma _s}$. On the other hand, mass conservation results in the scaling 
$a ({\rm d}a/{\rm d}t) \sim h u_{TC}$, leading to 
$T_{ret}=a/({\rm d}a/{\rm d}t) \sim \sqrt {\rho _l a_0^{4}/(2 \sigma _s h_0)}$; for 
$T_{RP} \ll T_{ret}$, then, the liquid sheet must satisfy 
$\sqrt {h_0/a_0}<<1$, as 
$a_0^{2} \sim h_0L_0$, which leads to the next condition 
$\sqrt [4]{h_0/L_0}\ll 1$ that is accomplished at very high aspect ratios (see Mirjalili, Chan & Mani (Reference Mirjalili, Chan and Mani2018) for more details). Finally, we justify the initial sinusoidal perturbation by looking into the breakup time 
$T_B$ as a function of 
$T_B=f(Oh,\varepsilon )$ (see expression 
$5$ from Driessen et al. Reference Driessen, Jeurissen, Wijshoff, Toschi and Lohse2013); for 
$T_B=f(Oh,\varepsilon =0.25) \sim T_{RP}$.
The computational domain corresponds to 
$6.75 \lambda _{max} \times \lambda _{max} \times \lambda _{max}$ (i.e. 
$304.29 h_0 \times 45.08h_0 \times 45.08h_0$) where the 
$x$-coordinate is aligned with the length of the sheet. The domain is sufficiently large in the streamwise direction to avoid the effect of artificial reflections from the boundary in this direction. The simulations are initialised with fluids at rest in the absence of gravity. The periodic boundary condition is imposed on all the variables in the (spanwise) 
$z$-direction of the domain. A no-penetration boundary condition is prescribed for the bottom and top domains. The domain is discretised with a regular cubic mesh of size 
$\Delta \boldsymbol {x}= h_0/6$, ensuring that our numerical predictions are mesh independent; our results do not vary with decreasing grid size. Note that similar mesh sizes have been used in previous studies (Agbaglah Reference Agbaglah2021; Gordillo et al. Reference Gordillo, Agbaglah, Duchemin and Josserand2011) for three- and two-dimensional simulations of surfactant-free retracting liquid sheets. Additionally, liquid volume and surfactant mass conservation are satisfied with errors of under 
$0.001\,\%$. Extensive mesh studies for surface-tension-driven phenomena with the same numerical methodology have been presented in Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020,Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021c, Reference Constante-Amores, Batchvarov, Kahouadji, Shin, Chergui, Juric and Matara). At the early stages of the simulation, capillary action drives the formation of a smooth cylindrical rim at the free end of the liquid sheet. After a short transient stage, the dynamics reaches the Taylor–Culick regime, characterised by a constant retraction velocity. Figure 1(b) contrasts our numerical results against the expected Taylor–Culick retraction velocity for 
$Oh=0.1$, 
$\epsilon =0$ and 
$e=2$.
As seen, the constant retraction velocity is 
${\sim }8\,\%$ smaller than the predicted Taylor–Culick velocity. This difference is consistent with the observations by Agbaglah (Reference Agbaglah2021) and Song & Tryggvason (Reference Song and Tryggvason1999) and is explained by the large density and viscosity ratios. The retraction dynamics of a liquid thread and the nonlinear dynamics before pinch-off have been validated previously by Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020). We refer to Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018) for the accuracy of the surfactant equations.
2.2. Parameters
Our conditions assume a retracting water liquid sheet in air, i.e. 
$\rho _g / \rho _l = 1.2 \times 10^{-3}$ and 
$\mu _g / \mu _l = 0.018$, as in the previous works of Agbaglah et al. (Reference Agbaglah, Josserand and Zaleski2013) and Agbaglah (Reference Agbaglah2021). We assume a thickness of a soap film of 
$h_0 = 2.0 \ \mu$m, i.e. 
$Oh=0.0833$ (in agreement with Meister & Scheele Reference Meister and Scheele1969). The parameter 
$\beta _s$ depends on 
$\varGamma _\infty$ and therefore on the critical micelle concentration, i.e. 
$\varGamma _\infty \sim {O}(10^{-6})$ mol m
$^{-2}$ for NBD-PC (1-palmitoyl-2-12-[(7-nitro-2-1,3-benzoxadiazol-4- yl)amino]dodecanoyl-sn-glycero-3 -phosphocholine); thus, here, we explore the range of 
$0.1 < \beta _s < 0.5$. We have set 
$Pe_s = 10^{2}$ following Batchvarov et al. (Reference Batchvarov, Kahouadji, Magnini, Constante-Amores, Craster, Shin, Chergui, Juric and Matar2020), who suggested that the interfacial dynamics is weakly dependent on 
$Pe_s$ beyond this value i.e. convective effects, driven by surface-tension gradients (Marangoni stresses), rather than surface diffusion, dominate the interfacial distribution of surfactant when 
$Pe_s >10^{2}$. Hence, for a liquid sheet characterised by 
$h_0 = 2.0\,\mathrm {\mu }$m, the relevant time scales are given by 
$T_{RP} \sim {O}(10^{-8})$ s, 
$T_{ret} \sim {O}(10^{-6})$ s and the Marangoni time scale, 
$T_{\tau } \sim \mu _l h_0/ \Delta \sigma \sim {O}(10^{-8})$ s. Consequently, Marangoni stresses play a key role in the flow dynamics.
3. Results
We start the discussion of the results by presenting a phenomenological picture of the interfacial dynamics in a phase diagram in a 
$\beta _s{-}Oh$ space. This dynamics ranges from a nearly inviscid sheet (
$Oh \sim 10^{-3}$) to typical soap films (
$Oh \sim 10^{-1}$). Figure 2 shows the regime map in terms of the interfacial dynamics predicted from our numerical simulations. Our results show that two distinct regimes can be identified based on the outcome of the retracting dynamics. At low 
$\beta _s$ and 
$Oh$, surfactant-driven Marangoni stresses do not promote the reopening of the adjacent sheet connected to the rim; thus, capillary waves result in the formation of a hole behind the rim. This hole grows in the spanwise direction and separates the rim from the sheet. In fact, this dynamics has previously been observed by Mirjalili et al. (Reference Mirjalili, Chan and Mani2018) and Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021b). The rest of the 
$\beta _s\unicode{x2013}Oh$ space is dominated by surfactant-induced Marangoni stresses which reopen the adjacent sheet to promote droplet detachment via RP instability. Below, we focus on this regime, and provide an extensive explanation of the interfacial dynamics.
Regime map of the phenomenological interfacial dynamics in the 
$\beta _s$–
$Oh$ space. Circles represent numerical simulations in which capillary waves induce the formation of a hole behind the rim; hollow square symbols represent a droplet shedding from the rim via a RP instability. Here, 
$Pe_s=100$, 
$\varGamma =\varGamma _\infty /2$, 
$e=0.2$ and 
$\epsilon =0.25$.
The spatio-temporal interfacial dynamics for the surfactant-free case, with 
$Oh=0.0833$, is shown in figure 3(a–e). As seen, the initial rim instability grows to develop nonlinearities that eventually lead to a capillary singularity (break up). At early stages of the simulation, the fluid is pulled from the rim to form an elongated ligament (see figure 3c). Over time, the ligament undergoes end pinching, where capillarity overcomes viscosity to form a bulbous end (see figure 3d) terminating in the breakup of a droplet (see figure 3e). The interfacial shape prior to the capillary singularity closely resembles experimental observations by Zhang et al. (Reference Zhang, Brunet, Eggers and Deegan2010) and Wang & Bourouiba (Reference Wang and Bourouiba2018, Reference Wang and Bourouiba2021). Additionally, we observe surface waves resulting from the development of capillary instabilities at the upstream of the rim.
Effect of surfactant elasticity (
$\beta _s$) on the sheet retraction dynamics for insoluble surfactants. Here, the dimensionless parameters are 
$Oh = 0.0833$, 
$e=0.2$ and 
$\epsilon =0.25$; and for the surfactant-laden cases, 
$Pe_s=100$ and 
$\varGamma =\varGamma _\infty /2$. For the surfactant-laden cases, the colour indicates the surfactant concentration 
$\varGamma$.
Attention is now turned to the effect of surface-active agents by varying the surfactant elasticity parameter, 
$\beta _s$, with 
$Oh=0.0833$, 
$Pe_s=100$ and 
$\varGamma =\varGamma _\infty /2$. For 
$\beta _s=0.1$, the spatio-temporal interfacial dynamics resembles that of a surfactant-free case, from the formation of the ligament to the eventual end pinching: the dynamics results in droplet shedding (see figure 3f–j). In this case, surfactant is accumulated at the base of the cylindrical rim as it gradually increases its size. Once the ligament is formed, 
$\tau$ induces a flow towards the ligament (see figure 3h) resulting in a longer ligament prior to pinch-off. However, as surfactant-induced Marangoni stresses are not sufficiently strong to evacuate 
$\varGamma$, the dynamics results in modest surfactant concentration gradients, and, in agreement with Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020), the ligaments break up. By increasing 
$\beta _s$, 
$\tau$ acts to reduce 
$\varGamma$ gradients by inducing a flow from the rim towards the sheet, and from the rim towards the ligament (regions of low surfactant concentration), resulting in a severe reduction of the rim radius (displayed in figure 3). The surfactant-driven flow leads to the reopening of the fluid sheet adjacent to the retracting rim, thereby increasing its film thickness (see figure 4i); the Marangoni-induced flow also acts to suppress the capillary waves ahead of the rim (see figures 3(k–o) and 3(j–n) for 
$\beta _s=0.3$ and 
$\beta _s=0.5$, respectively).
Panels (a–h) show the effect of the elasticity parameter, 
$\beta _s$, on the flow and surfactant concentration during sheet retraction. Two-dimensional projection of the interface, 
$\varGamma$ and 
$\tau$ and 
$u_{tx}$ and 
$u_{tz}$ in the 
$x\unicode{x2013}z$ plane (
$y=\lambda /2$) are shown in (a–d) and (e–h) for 
$t=141.42$ and 
$t=189.50$, respectively. Note that the abscissa in (a,e) corresponds to the 
$z$ coordinate, and in (b–d) and (f–h) to the arclength, 
$s$. The diamond shapes show the location of the bulbous necks. Panel (i) shows the effect of 
$\beta _s$ on the film thickness of the sheet in the 
$x\unicode{x2013}y$ plane (
$z=\lambda /2$). Panels (j–k) show the effect of surfactants on the vorticity, 
$\omega _z$, for the surfactant-free case (left images), and the surfactant-laden case (right images), for 
$\beta _s=0.5$ at 
$t=141.42$ and 
$t=189.50$, respectively.
We now turn our attention to the two-dimensional interfacial shape represented by 
$\varGamma$ and spanwise Marangoni stresses, the tangential component of the surface velocity 
$u_{t_{z}}$ and the streamwise surface velocity 
$u_{t_{x}}$ (where the velocity is given in the reference frame of 
$u_x$ at 
$z=0$) presented in figure 4 at 
$t=141.42$ and 
$t=189.50$, respectively. As observed, the surfactant-driven flow, in the spanwise direction, results in the retardation of the dynamics brought about by surfactant-induced interfacial rigidification (see the dampening of 
$u_{t_{x}}$ in the frame of difference in figure 4(b,f); as 
$\beta _s>0$, 
$u_{t_{x}} < 0$). We observe that surfactant concentration gradients at 
$t=141.42$ give rise to Marangoni stresses which drive a flow from the horizontal part of the rim towards the base of the ligament, and from the tip of the ligament to its base (for high 
$\beta _s$, the ligament tip has a nearly uniform distribution of 
$\varGamma$). The combined effect of Marangoni stresses is to bridge the surfactant gradient between the tip of the ligament and the rim (see figure 4h). This effect is sufficiently strong to promote the development of a thinner bulbous end, and cause a surfactant-driven escape from end pinching to a longer ligament (in agreement with Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020) and Kamat et al. (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2020)). As expected, at a later stage, the nearly uniform distribution of surfactant concentration results in the elimination of Marangoni stresses (see figure 4h). Figure 4(c,g) shows that 
$u_{t_{z}}> 0$ upstream and 
$u_{t_{z}} < 0$ downstream of the neck, i.e. a stagnation point is found in between. As the neck narrows, capillary-driven flow causes further neck thinning and the creation of a velocity maximum at 
$u_{t_{z}}$, as illustrated in figure 4(c,g), which is indicative of singularity formation. The 
$u_{t_{z}}$ profile for the surfactant-free case is characterised by the presence of a large velocity maximum and two stagnation points, with the neck located in between. By close inspection of the 
$u_{t_{z}}$ profile of the surfactant-laden cases, it becomes clear that only one stagnation point is found near the neck for high 
$\beta _s$. The surfactant-induced Marangoni stresses result in the suppression of one of the stagnation points – in agreement with Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020, Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021c).
Figure 4(j,k) shows the vorticity field, i.e. 
$\omega =\triangledown \times \boldsymbol {u}$, in terms of the interface location for the surfactant-free and surfactant-laden cases for 
$t=141.42$ and 
$t=189.50$ (the location is given in a reference frame where the rim position is found at 
$z=0$). For all cases, a primary vortex (with a positive sign) is formed at the front region of the bulbous edge and is shed into the ambient fluid. For the surfactant-free case, we observe a negative vortex ring near the bulbous neck resulting from a change of curvature; this vortex ring grows over time, ending in breakup.
For the surfactant-laden case at 
$t=141.42$, 
$\tau$ induces vorticity from the bulbous tip to the neck and from the rim to the ligament. This flow causes the neck to reopen; after escaping the singularity, the flow through the neck triggers the formation of a jet towards the centre of the bulbous region (these findings are consistent with Constante-Amores et al. Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020). Additionally, the escape from pinch-off causes the length of the ligament to grow over time in the 
$(x)$-direction, forming a secondary bulbous and a secondary neck as well as a shift in the axial curvature signs (see figure 4k).
To provide conclusive evidence that Marangoni stresses drive the reopening of the neck, we run an additional simulation in which Marangoni stresses, or surface-tension gradients, are deactivated, i.e. 
$\boldsymbol {\nabla }_s \sigma =0$, to isolate the effects of a reduction in surface tension from those arising from Marangoni stresses. The Marangoni-suppressed case resembles the surfactant-free dynamics culminating in end pinching (see Supplementary animation). Additionally, we performed an additional surfactant-free simulation at the ‘effective’ Ohnesorge number, obtained by using the surface tension achieved by the surfactant, i.e. 
$Oh_{eff}=\mu /\sqrt {\rho h_0 \sigma _0}=0.103$ (for 
$\beta _s=0.5$ and 
$\varGamma _0=0.5$), finding similar results. Consequently, we claim Marangoni stresses are responsible of the change of the interfacial dynamics discussed above.
Finally, figure 5 shows the effect of the surfactant dynamics on the filament's tip location, 
$L$, the kinetic energy, 
$E_k= \int _{\mathcal {V}} (\rho \boldsymbol {u}^{2})/2 \,{\rm d}{\mathcal {V}}$, and the length of the ligament (defined from the tip of the bulbous edge to its rim). Here, the kinetic energy has been normalised by the surface energy 
$E_s = A_0\sigma _s$, where 
$A_0$ is the initial surface of the sheet. The evidence for a surfactant-driven retardation of the dynamics can be seen by inspection of 
$E_k$ and the streamwise location of the tip, which reveals that increasing 
$\beta _s$ monotonically decreases the overall value of both 
$E_k$ and 
$L$. Thus, the interfacial flow dynamics is retarded by Marangoni stresses resulting in a reduction in the retraction velocity. In figure 5(c), we report the length of the ligament until pinch-off, the addition of surfactant increases the ligament length prior to its breakup, becoming particularly pronounced at high 
$\beta _s$.
Temporal evolution of the ligament's tip position, the kinetic energy and the length of the ligament are shown in (a–c), respectively, for the surfactant-free and surfactant-laden cases. Simulations are stopped when the sheet sheds its first drop. All parameters remain unchanged from figure 3.
4. Concluding remarks
Results from three-dimensional numerical simulations of the retraction dynamics of thin liquid sheets in the presence of insoluble surfactants were presented. The liquid properties in the simulations were chosen to be consistent with a realistic air/water system, and the retraction velocity of the rim formed at the end is in good agreement with the Taylor–Culick speed. For the surfactant-laden cases, we have demonstrated that surfactant-induced Marangoni stresses drive a flow from the high surfactant concentration regions to the low concentration ones, reducing the kinetic energy, affecting the location of the bulbous tip and suppressing end pinching. With increasing elasticity number, the Marangoni stresses play a major role in the interfacial dynamics with the progressive elimination of the capillary wave structures upstream of the rim, where a complete interfacial rigidification is observed for large elasticity numbers. Additionally, Marangoni stresses drive flow towards the sheet, resulting in the reopening of the film thickness adjacent to the rim. Future research avenues are related to examining the role of solubility on the flow dynamics.
Supplementary material and movie
Supplementary material and movie are available at https://doi.org/10.1017/jfm.2022.768.
Funding
C.R.C.-A. and A.A.C.-P. acknowledge the support from the Royal Society through a University Research Fellowship (URF/R/180016), an Enhancement Grant (RGF/EA/181002) and two NSF/CBET-EPSRC grants (Grant Nos. EP/S029966/1 and EP/W016036/1). All authors are grateful by the computing time granted by the Institut du Developpement et des Ressources en Informatique Scientifique (IDRIS) of the Centre National de la Recherche Scientifique (CNRS), coordinated by GENCI (Grand Equipement National de Calcul Intensif) Grant No. 2022A0122B06721. Simulations were performed using code BLUE (Shin, Chergui & Juric Reference Shin, Chergui and Juric2017) and the visualisations were generated using Paraview.
Declaration of interests
The authors report no conflict of interest.
Open access
