1. Introduction
When a bubble is resting close to a liquid–gas interface, its rupture gives rise to the formation of a central jet. This jet breaks up into small droplets, which could transport biological material, toxins, salts, surfactants or dissolved gases (Woodcock et al. Reference Woodcock, Kientzler, Arons and Blanchard1953; MacIntyre Reference MacIntyre1972; Veron Reference Veron2015; Séon & Liger-Belair Reference Séon and Liger-Belair2017; Poulain & Bourouiba Reference Poulain and Bourouiba2018; Zenit & Rodríguez-Rodríguez Reference Zenit and Rodríguez-Rodríguez2018). It is unsurprising therefore that bursting bubble phenomena have received significant interest due to their occurrence in a multitude of natural and industrial applications. In the absence of contaminants, the physical mechanisms of surfactant-free bursting bubbles on the ejection of droplets have been widely studied experimentally (Woodcock et al. Reference Woodcock, Kientzler, Arons and Blanchard1953; Ghabache et al. Reference Ghabache, Antkowiak, Josserand and Séon2014; Ghabache & Séon Reference Ghabache and Séon2016; Séon & Liger-Belair Reference Séon and Liger-Belair2017), numerically (Deike et al. Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Séon2018; Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2019; Singh & Das Reference Singh and Das2019) and through theoretical approaches (Gañán-Calvo Reference Gañán-Calvo2017; Lai, Eggers & Deike Reference Lai, Eggers and Deike2018; Blanco-Rodrìguez & Gordillo Reference Blanco-Rodrìguez and Gordillo2020; Zeff et al. Reference Zeff, Kleber, Fineberg and Lathrop2000).
This previous work has shown that when a nucleated bubble rises through the liquid and then rests close to a free surface, its static resting shape results from a balance between gravitational and surface tension forces. The bubble shape may be characterised by a submerged interface, a liquid layer/cap above the bubble and an interface, which extends away from the bubble cap (Toba Reference Toba1959; Ghabache Reference Ghabache2015). The cap is characterised by a length scale 
$\delta / R_o \sim O (10^{-6})$, where 
$\delta$ and 
$R_o$ refer to the liquid layer thickness and the initial bubble radius, respectively. The layer curvature creates over-pressure relative to the liquid bulk below it. The fluid within the layer drains over time, enabling van der Waals forces to rupture the interface when 
$\delta \rightarrow 0$, forming a hole (Lhuissier & Villermaux Reference Lhuissier and Villermaux2012). The hole leaves an open, unstable cavity that collapses to form a vertical jet whose dynamics is governed by inertial, viscous and capillary forces.
Surfactants can affect the dynamics of surface-tension-driven flows by the reduction of the local surface tension and capillary pressure, and by the formation of Marangoni stresses brought about by gradients in interfacial surfactant concentration. In a recent study, Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020) reviewed the state-of-the-art of the effect of surfactant on the dynamics of capillary singularities due to topological change. Craster, Matar & Papageorgiou (Reference Craster, Matar and Papageorgiou2002) have shown that the presence of insoluble surfactant during the thinning of threads does not alter the well known breakup scalings predicted by Eggers (Reference Eggers1993) as the capillary singularity is such a violent event which convects surfactant away from the pinch-off point. Both McGough & Basaran (Reference McGough and Basaran2006) and Kamat et al. (Reference Kamat, Wagoner, Thete and Basaran2018) have reported that the presence of surfactant changes the fate of the thinning of threads as surfactant enhances the formation of micro-threads, whose thinning dynamics also follows Eggers (Reference Eggers1993). Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020) showed that the presence of surfactant can suppress the ‘end-pinching’ mechanism in the retracting dynamics of ligament threads due to the suppression of stagnation points, and subsequent flow reversal. Most of chemical surfactants are soluble in the bulk liquid, adding a layer of complexity to the phenomena as Marangoni stress can also be regulated by the surfactant sorption kinetics between the bulk and interface. Liao et al. (Reference Liao, Subramani, Franses and Basaran2004) observed experimentally that high concentrations of soluble surfactant favour an asymmetric breakup of a liquid bridge. Jin, Gupta & Stebe (Reference Jin, Gupta and Stebe2006) and Jin & Stebe (Reference Jin and Stebe2007) showed that solubility affects the dynamics of drop formation in terms of thinning. Craster, Matar & Papageorgiou (Reference Craster, Matar and Papageorgiou2009) have shown that, similar to their previous work (Craster et al. Reference Craster, Matar and Papageorgiou2002), the addition of solubility does not affect the scaling predicted by Eggers (Reference Eggers1993), but it does lead to the formation of larger satellite droplets, which was confirmed experimentally by Kovalchuk, Nowak & Simmons (Reference Kovalchuk, Nowak and Simmons2016) and Kovalchuk et al. (Reference Kovalchuk, Jenkinson, Miller and Simmons2018).
As revealed by the foregoing review, studies involving surfactant effects on interfacial flows have received considerable attention, however, to the best of our knowledge, the influence of surfactants on the dynamics of bursting bubbles has not been explored, and this is the subject of the present article. Here, we use three-dimensional numerical simulations to account for surfactant solubility, diffusion, sorption kinetics and Marangoni stresses in bursting phenomena. The rest of this paper is structured as follows: in § 2, the numerical method, governing dimensionless parameters, problem configuration and validation, are introduced. Section 3 presents the results, and concluding remarks are given in § 4.
2. Problem formulation and numerical method
The numerical simulations were performed by solving the two-phase Navier–Stokes equations in a three-dimensional Cartesian domain 
$\boldsymbol {x} = (x, y, z )$ (see figure 1). A hybrid front-tracking/level-set method was used to treat the interface, where surfactant transport was resolved both in the bulk and on the interface (Shin et al. Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018). Here, and in what follows, all variables will be made dimensionless (represented by tildes) using
\begin{equation} \left.\begin{gathered} \tilde{\boldsymbol{x}}=\frac{\boldsymbol{x}}{R_o},\quad \tilde{t}=\frac{t}{T}, \quad \tilde{{\boldsymbol{u}}}=\frac{{\boldsymbol{u}}} {U},\quad \tilde{p}=\frac{p}{\rho U^2}, \quad \tilde{\sigma}=\frac{\sigma}{\sigma_s},\\ \tilde{\varGamma}=\frac{\varGamma}{\varGamma_\infty}, \quad \tilde{C}=\frac{C}{C_{\infty} }, \quad \widetilde{C_s}=\frac{C_s}{C_{\infty} }, \end{gathered}\right\} \end{equation}
where, 
$t$, 
$\boldsymbol {u}$ and 
$p$ stand for time, velocity and pressure, respectively. The physical parameters correspond to the liquid density 
$\rho$, viscosity, 
$\mu$, surface tension, 
$\sigma$, surfactant-free surface tension, 
$\sigma _s$, and gravitational acceleration, 
$g$; 
$T=\sqrt {\rho R_o^3/\sigma _s}$ is the capillary time scale, hence the velocity scale is 
$U=R_o/T= \sqrt {\sigma _s/(\rho R_o)}$. The interfacial surfactant concentration, 
$\varGamma$, is scaled on the saturation interfacial concentration, 
$\varGamma _{\infty }$, whereas the bulk and bulk sub-phase (immediately adjacent to the interface) surfactant concentrations given by 
$C$ and 
$C_s$, respectively, are scaled on the initial bulk surfactant concentration, 
$C_\infty$. As a result of this scaling, the dimensionless equations read
\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} = -Bo \boldsymbol{i}_z + Oh \, \boldsymbol{\nabla} \boldsymbol{\cdot} \left [ \tilde{\mu} (\boldsymbol{\nabla} \tilde{\boldsymbol{u}} +\boldsymbol{\nabla} \tilde{\boldsymbol{u}}^\textrm{T}) \right ] \nonumber\\ \hspace{10pc} +\, \int\limits_{\tilde{A}(\tilde{t})} (\tilde{\sigma} \tilde{\kappa} \boldsymbol{n} + \boldsymbol{\nabla}_s \tilde{\sigma}) \delta (\tilde{\boldsymbol{x}} - \tilde{\boldsymbol{x}}_{f})\mbox{d}\tilde{A}, \end{gather}
\begin{gather} \frac{\partial \tilde{C}} {\partial \tilde{t}}+\tilde{\boldsymbol{u}}\boldsymbol{\cdot} \boldsymbol{\nabla} \tilde{C}= \frac{1}{Pe_b} \nabla^2 \tilde{C}, \end{gather}
\begin{gather} \frac{\partial \tilde{\varGamma}}{\partial \tilde{t}}+\boldsymbol{\nabla}_s \boldsymbol{\cdot} (\tilde{\varGamma}\tilde{\boldsymbol{u}}_{{t}})=\frac{1}{Pe_s} \nabla^2_s \tilde{\varGamma}+ Bi ( k \widetilde{C_s} (1-\tilde{\varGamma})- \tilde{\varGamma} ), \end{gather}
\begin{gather} \tilde{\sigma}=1 + \beta_s \ln{(1 -\tilde{\varGamma})},\end{gather}
where the density and viscosity are given by 
$\tilde {\rho }=\rho _g/\rho + (1 -\rho _g/\rho ) \mathcal {H}(\tilde {\boldsymbol {x}},\tilde {t})$ and 
$\tilde {\mu }=\mu _g/\mu + (1 -\mu _g/\mu ) \mathcal {H}( \tilde {\boldsymbol {x}},\tilde {t})$ wherein 
$\mathcal {H}( \tilde {\boldsymbol {x}},\tilde {t})$ represents a smoothed Heaviside function, which is zero in the gas phase and unity in the liquid phase, while the subscript 
$g$ designates the gas phase; 
$\tilde {\boldsymbol {u}}_{{t}}= ( \tilde {\boldsymbol {u}}_{{s}} \boldsymbol {\cdot } \boldsymbol {t} ) \boldsymbol {t}$ is the tangential velocity at the interface in which 
$\tilde {\boldsymbol {u}}_{{s}}$ represents the interfacial velocity; 
$\kappa$ is the curvature; 
$\boldsymbol {\nabla }_s=({\boldsymbol {I}}-\boldsymbol {n}\boldsymbol {n})\boldsymbol {\cdot } \boldsymbol {\nabla }$ is the interfacial gradient wherein 
$\boldsymbol {I}$ is the identity tensor and 
$\boldsymbol {n}$ is the outward-pointing unit normal; 
$\delta$ is the Dirac delta function, which is equal to unity at the interface, located at 
$\tilde {\boldsymbol {x}}=\tilde {\boldsymbol {x}}_f$, and zero otherwise; 
$\tilde {A} (\tilde {t})$ is the time-dependent interface area.
Figure 1. (a) Initial shape of the bubble resting close to the interface, highlighting the computational domain of size 
$(15 R_o)^3$ (not to scale) in a three-dimensional Cartesian domain 
$\boldsymbol {x} = (x, y, z)$, where the Cartesian resolution is 
$768^3$; (b) comparison of surfactant-free simulations from the current study (squares) with scaling argument of the length of the jet (solid line) and numerical simulations (dots) proposed by Lai et al. (Reference Lai, Eggers and Deike2018).
The dimensionless groups that appear in the above equations are defined as
\begin{gather} Bo=\frac{\rho g R_o^2}{\sigma_s }, \quad La=\frac{1}{Oh^2}= \frac{\rho \sigma_s R_o}{\mu^2 }, \end{gather}
\begin{gather} Bi=\frac{k_d R_o}{U}, \quad k=\frac{k_a C_\infty }{k_d}, \quad Pe_s=\frac{ U R_o}{D_s}, \quad Pe_b=\frac{ U R_o}{D_b}, \quad \beta_s= \frac{\Re T \varGamma_\infty }{\sigma_s }, \end{gather}
where 
$Bo$ stands for the Bond number and represents the ratio of gravitational to capillary forces; 
$Oh=\mu /\sqrt {\rho \sigma _s R_o}$ is the Ohnesorge number that measures the relative importance of viscous to surface tension forces, and 
$La$ is the Laplace number; 
$Bi$ denotes the Biot number representing the ratio of characteristic desorptive to convective time scales; 
$k$ is the ratio of adsorption to desorption time scales; 
$Pe_s$ and 
$Pe_b$ are the interfacial and bulk Péclet numbers, respectively, and represent the ratio of convective to diffusive time scales in the plane of the interface and the bulk, respectively; and 
$\beta _s$ is the surfactant elasticity number, which measures the sensitivity of the surface tension to the surfactant concentration. The parameter 
$\mathcal {R}$ refers to the thermodynamic ideal gas constant value 
$8.314\ \textrm {J}\ \textrm {K}^{-1}\ \textrm {mol}^{-1}$. In addition to the above parameters, we assume an initially uniform surfactant concentration at the interface so that 
$\varGamma (\boldsymbol {x},t=0)=\varGamma _o=\varGamma _{\infty }/4$. The chosen density and viscosity ratios, 
$\rho _g/ \rho =1.2 \times 10^{-3}$ and 
$\mu _g/\mu = 0.018$, respectively, are representative of an air–water system.
The Marangoni stress, 
$\tilde {\tau }$, is expressed as a function of 
$\tilde {\varGamma }$
\begin{equation} \tilde{\tau} \equiv \boldsymbol{\nabla}_s \tilde{\sigma} \boldsymbol{\cdot} \boldsymbol{t} =-\frac{\beta_s}{1 -\tilde{\varGamma}} \boldsymbol{\nabla}_s\tilde{\varGamma} \boldsymbol{\cdot} \boldsymbol{t}, \end{equation}
where 
$\boldsymbol {t}$ is the unit tangent to the interface; tildes are dropped henceforth. From (2.9), the variation of 
$\tau$ can be achieved by altering the elasticity parameter, the magnitude of the interfacial concentration or the interfacial concentration gradients. In the current study, the weakening or strengthening of the Marangoni stresses has been analysed by studying the variation of 
$\beta _s$. The Marangoni time scale, 
$\mu R_o/ \Delta \sigma = {O(10^{-5})}$ s, as compared to the capillary and sorptive time scales, which are of order 
$O(10^{-4})$ s and 
$O(10^{-3}\text {--}10^{-5})$ s, respectively; thus Marangoni stresses will play a key role in the dynamics.
The numerical procedure to solve (2.2)–(2.6) has been presented in detail by Shin, Chergui & Juric (Reference Shin, Chergui and Juric2017) and Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018); only a brief summary is provided here. The Navier–Stokes equations are solved by using a finite-volume method over a staggered grid (Harlow & Welch Reference Harlow and Welch1965). The computational domain comprises of a fixed regular grid (i.e. Eulerian grid) in which the spatial derivatives are approximated by standard centred-difference discretisation, except for the nonlinear term, which is discretised using a second-order essentially non-oscillatory scheme (Shu & Osher Reference Shu and Osher1989). The projection method presented by Chorin (Reference Chorin1968) is used to enforce the incompressibility behaviour.
The interface is tracked by an additional Lagrangian grid by using the front-tracking method and regularly reconstructed by a level-set method. The immersed boundary method of Peskin (Reference Peskin1977) is used for communication between both grids. The temporal integration scheme is based on a second-order Gear method. Thus, even though our method as described in Shin & Juric (Reference Shin and Juric2002, Reference Shin and Juric2009) follows a hybrid level-set/front-tracking approach, incorporating some features of level sets, it is important to note that it fully retains the well-established immersed boundary formulation for surface forces used in the front-tracking method, which requires the surface integral representation in (2.3). The mathematical apparatus is detailed in the pioneering work of Peskin (Reference Peskin1977). Finally, surfactant transport is solved on the interface, where the interfacial surfactant concentration is located on the centre of the triangular front elements. More information can be found in Shin et al. (Reference Shin, Chergui and Juric2017, Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018).
2.1. Numerical set-up, validation and parameters
The dimensionless computational domain size is chosen as 
$(15R_o)^3$, which is found to be sufficiently large to render boundary effects negligible. Hence, a radial component is defined as 
$r=\sqrt {(x-x_o)^2 + (y-y_o)^2}$ where 
$x_o$ and 
$y_o$ are the abscissa and ordinate bubble position, respectively. Solutions are sought subject to Neumann boundary conditions on all variables at the lateral boundaries, 
$p=0$ at the top boundary 
$z=15R_o$, and no-slip at 
$z=0$. At the free surface, we impose 
$\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {\nabla } \tilde {C}=-BiPe_b ( k \widetilde {C_s} (1-\tilde {\varGamma })- \tilde {\varGamma } )$ as a condition on 
$\tilde {C}$ (more information can be found in Shin et al. Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018). Simulations are initialised as a bubble resting immediately beneath the free surface before its rupture. Its initial shape is determined by solving the Young–Laplace equation for 
$Bo=10^{-3}$, which is sufficiently small to avoid the effect of gravity on the bubble bursting process. To initialise bursting, the top spherical cap of the bubble is removed, leaving only a bubble cavity (see figure 1a). The relevant time scale associated with the cap retraction, 
$t_{CR}=R_o/\sqrt {2 \sigma /\rho \delta }= O(10^{-6})$ s, is too short to affect the dynamics. A similar approach has been used by Boulton-Stone & Blake (Reference Boulton-Stone and Blake1993), Garcia-Briones, Brodkey & Chalmers (Reference Garcia-Briones, Brodkey and Chalmers1994), Duchemin et al. (Reference Duchemin, Popinet, Josserand and Zaleski2002) and Deike et al. (Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Séon2018). Our numerical simulations have been validated against the work of Eggers (Reference Eggers1993) and Lai et al. (Reference Lai, Eggers and Deike2018) in terms of liquid thread breakup, and the scaling of the ejected jet length, 
$L_d$, with 
$La$ (see figure 1b), respectively. Solving the small scales of the bursting dynamics is a challenging process. We have assessed the grid dependence of our results ensuring their convergence for grids larger than 
$768^3$.
In the current study, the Laplace number has been fixed to 
$La=2\times 10^{4}$ following the works of Lai et al. (Reference Lai, Eggers and Deike2018) and Gordillo & Rodríguez-Rodríguez (Reference Gordillo and Rodríguez-Rodríguez2019). Additionally, we have set the Bond number 
$Bo$ to be of 
$O(10^{-3})$ to ensure that the shape of the bubble prior to its rupture is spherical even in the presence of surfactants. Operating in this parameter space allows for the elimination of gravitational effects and the isolation of Marangoni-induced dynamics during jet-drop formation. As mentioned earlier, to study the effect of Marangoni stresses on the flow, we have varied the elasticity parameter in the range of 
$0.7<\beta _s<0.9$. In the context of this work, both the superficial and bulk Péclet numbers are set to unity, making neither the diffusion nor the convective terms dominant. With respect to the solubility, we have studied the range 
$10^{-2}<Bi<10^0$, where for the lower end of the parameter scale the sorptive time scales are much larger than those associated with interfacial effects, and the dynamics is expected to be similar to that observed for an insoluble surfactant case, whereas for the upper end, surfactants tend to desorb from the interface and the dynamics begins to resemble that of a surfactant-free case. In summary, we have chosen the values of the surfactant-related parameters to ensure that all of the relevant physical processes associated with surfactant transport such as Marangoni stresses, surface/bulk diffusion and sorption kinetics are represented in the present study.
3. Results
3.1. Bursting bubble dynamics of the surfactant-laden base case
The spatio-temporal evolution of the interfacial dynamics, shown in figure 2, is considered for the case characterised by 
$La=2\times 10^4$, 
$Bo=10^{-3}$, 
$Pe_s=1$, 
$\beta _s =0.9$, 
$Bi=0.01$, 
$k=1$ and 
$\varGamma _o =k/4$. At early times, a large capillary pressure is generated near the nucleated hole due to the high curvature of the interface joining the spherical bubble and the horizontal free surface. This capillary pressure leads to rapid expansion of the hole and the creation of a toroidal capillary wave, which travels toward the bottom of the bubble. The convergence of this wave on the bubble rear leads to the formation of a cusp-like region, which is relieved via the detachment of downward-moving conical bubbles, and the formation of a vertical, upward-directed, high-speed, Worthington jet. Small droplets are subsequently ejected from the tip of the jet, triggered by a Rayleigh–Plateau instability. A pinch-off event is also seen to occur at the base of the Worthington jet that leads to retraction of the emitted ligament into a spheroidally shaped drop.
Figure 2. Spatio-temporal evolution of the dynamics of the interface and of the interfacial surfactant concentration, 
$\varGamma$, with 
$La=2\times 10^4$, 
$Bo=10^{-3}$, 
$Pe_s=1$, 
$\beta _s =0.9$, 
$Bi=0.01$, 
$k=1$ and 
$\varGamma _o =k/4$. The colour bars indicate the magnitude of 
$\varGamma$. (a) Top view of the interface; (b) side view of cavity collapse; (c) Worthington jet (entrapped bubble is not shown).
In figure 3, we show snapshots of the interface, the interfacial concentration, 
$\varGamma$, the Marangoni stress, 
$\tau$, the interfacial tangent velocity component, 
$u_t$, that provides a measure of mobility, and the azimuthal component of the vorticity, 
$\omega _\theta$. These snapshots, which are taken at 
$t=0.360$, 
$0.539$, 
$1.239$ and 
$1.760$, corresponding to the framed panels in figure 2, reflect the strong coupling between these flow variables. Figure 3(a,b) illustrates that the collapsing cavity and the accompanying capillary wave transport surfactant towards the bottom of the bubble, giving rise to a local decrease in surface tension. It is also instructive to separate figure 3(b) into four regions due to the existence of stagnation points. In region ‘
$A$’, 
$\tau < 0$ and 
$u_t < 0$, indicating that the direction of the Marangoni flow is towards the bottom of the bubble, which aids cavity collapse and surfactant transport in this direction. In region ‘
$B$’, 
$\tau >0$ and 
$u_t > 0$, thus the Marangoni flow is directed away from the origin, which retards the collapse process. A similar behaviour is seen in regions ‘
$C$’ and ‘
$D$’, in which the trailing edge of the capillary wave has a similar 
$\varGamma$ distribution, albeit with a smaller magnitude, to its leading edge. As a result, 
$\tau < 0$ and 
$\tau > 0$, and 
$u_t < 0$ and 
$u_t > 0$ in regions ‘
$C$’ and ‘
$D$’, which drives flow towards the bottom of the bubble and its tail, respectively. Figure 3(b,c) also shows the existence of an interfacial stagnation point (
$s \sim 1.4$), where surfactant accumulates, 
$\varGamma$ is highest, and the magnitude of 
$\tau$ is largest. This occurs in the region where 
$\omega _{\theta }$ changes sign as the stagnation point is created (see figure 3d) corresponding to the capillary wave moving towards the bottom of the bubble.
Figure 3. Interface location, 
$\varGamma$ and 
$\tau$, 
$u_t$ and 
$\omega _{\theta }$ are shown in columns one to four, respectively. In columns 1,4 and 2,3, the variation is with respect to the dimensionless radial coordinate, 
$r$, and arc length, 
$s$, respectively. (a–d) and (e–h) show the cavity collapse dynamics 
$t=0.360$ and 
$t=0.539$, respectively, (i–l) the Worthington jet at 
$t=1.239$ and (m–h) the retracting ligament at 
$t=1.760$. The insets in (e,i) respectively focus on the bottom of the collapsing cavity, and the oscillations at the tip of the Worthington jet that will eventually lead to its breakup. A description of regions ‘
$A$’–‘
$D$’ in (b) is provided in the text. The parameter values are the same as in figure 2.
Figure 3(e–h) shows the dynamics at 
$t=0.539$ prior to the cavity collapse. The stagnation point where the surfactant accumulates has moved downward to the bottom of the cavity and the magnitude of 
$\omega _{\theta }$ has increased in comparison to that in figure 3(d). The interfacial surfactant concentration reaches its maximum value as the surfactant-laden capillary wave converges on the flow origin. The Marangoni stresses drive motion from high to low 
$\varGamma$ regions and therefore act to oppose the flow, as indicated by the fact that both 
$\tau >0$ and 
$u_t>0$ over the majority of the spatial domain except in the close vicinity of the bottom of the cavity.
Figure 3(i) shows a snapshot of the jet for 
$t=1.239$ in a situation of pinch-off ‘escape’ with a bulbous region formed at its tip. From figure 3(j), we observe that 
$\varGamma$ has two peaks: one at the jet tip, and a lower one far upstream. Figure 3(l) also shows that the change in the sign of 
$\omega _{\theta }$ is linked to the 
$\varGamma$ distribution in the bulbous region. The associated 
$\tau$ distribution is such that 
$\tau <0$ between the bottom of the jet and close to the bulbous region at the jet tip (that is for 
$0.4 <s< 1.9$), and 
$\tau >0$ elsewhere including in the bulbous region. Flow is driven by capillarity from this region towards the quasi-cylindrical jet body and, from figure 3(k), it is seen that 
$u_t > 0$ over the whole domain. This suggests that Marangoni stresses oppose this capillary-induced motion driving flow from the bottom of the jet towards the bulbous region. When these stresses are not sufficiently strong to overcome the Rayleigh–Plateau instability, an ejection of a droplet occurs after the pinch-off of the jet tip. The first drops detached from the tip are characterised by high interfacial concentration while successive droplet detachments have lower 
$\varGamma$, as demonstrated in figure 2 for 
$t\geqslant 1.114$.
At 
$t=1.680$ (see figure 2), the Worthington jet pinches off from its base, forming an elongated ligament thread. Capillary waves develop on the ligament surface, leading to interfacial oscillations as the detached ligament transitions to a spherical drop. In figure 3(m–p), we have isolated the retracting ligament shown at 
$t=1.760$ from the rest of the flow. In figure 3(o), we observe the formation of four stagnation points at 
$s=0$ 0.55, 0.7 and 0.84, which are connected to the change in the distributions of 
$\varGamma$ and 
$\tau$. From figure 3(n,o), it is also seen that the Marangoni stresses oppose ligament pinch-off since 
$\tau > 0$ and 
$u_t < 0$ (
$\tau <0$ and 
$u_t > 0$) for 
$0 \leqslant s \leqslant 0.55$ (
$0.55 \leqslant s \leqslant 0.7$). For 
$0.7 \leqslant s \leqslant 0.84$ (
$s > 0.84$), the Marangoni stresses oppose (aid) the stretching of the ligament as 
$\tau <0$ and 
$u_t >0$ (
$\tau >0$ and 
$u_t >0$). Close inspection of 
$\omega _{\theta }$ (see figure 3p) reveals that high vorticity production is observed close to the first stagnation point, which corresponds to the ligament ‘neck’. As shown by Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020), the presence of such high vorticity regions near the neck is a requirement for the Marangoni-driven inhibition of capillary-induced ‘end pinching’ of retracting surfactant-laden ligaments.
Figure 4(a) shows the temporal evolution of the jet tip from the moment after the ejection of a jet droplet to its next capillary singularity. The formation of a bulbous edge on the jet tip is driven by capillarity (at 
$t=1.120$). The estimated 
$Oh$ number for the tip, using the jet radius as a reference, is 
$Oh \approx 0.02$. Therefore, we are still in the region of low 
$Oh$ number, and consequently, the inhibition of the ‘end-pinching’ mechanism, identified by Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020) can be invoked in order to rationalise the flow behaviour in this case. By close inspection of 
$\tau$ and 
$u_{tr}$ (the radial component of the tangential velocity 
$u_t$), it is possible to identify that Marangoni-induced flow prevents the capillary breakup and reopens the neck of the jet tip (
$t=1.140$). The reason behind the reopening of the neck is the suppression of the stagnation point close to the neck due to the Marangoni-induced flow brought about by the surfactants (see figure 4d). As was shown by Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020), two stagnation points close to the neck are necessary for capillary breakup to occur (see figure 4e).
Figure 4. (a) Spatio-temporal evolution of the jet tip and its escape from pinch-off for the surfactant-laden base case. (b,c) Represent 
$\varGamma$ and 
$\tau$ on the arc length 
$s$ for time 
$t=1.140$ and 
$1.339$, respectively; (d,e) represent the radial component of the tangential velocity 
$u_{tr}$ as a function of the arc length 
$s$, respectively. The red arrows in (b,c) represent the direction of Marangoni-induced flow, and the blue arrows in (d,e) highlight the direction of the flow driven by capillarity. The value 
$s=0$ corresponds to the tip of the jet at 
$r=0$. The parameter values are the same as in figure 2.
In order to isolate the surface tension-reducing effects of surfactants from those associated with Marangoni stress formation, we show in figure 5(a) snapshots of the interface at 
$t=0.44$, 0.496 and 0.500 for the surfactant-free, surfactant-laden but Marangoni-suppressed and full-Marangoni cases, respectively, prior to cavity collapse for the same parameters as in figure 3. These times were chosen at identical spatial locations of the bubble rear in the axis of symmetry 
$(z \sim 3.28)$. For the Marangoni-suppressed case, the reduced surface tension value is calculated via replacing 
$\varGamma$ by 
$\varGamma _o$ in (2.6). The presence of surfactant has been shown to enhance capillary wave damping by Asaki, Thiessen & Marston (Reference Asaki, Thiessen and Marston1995) due to the interfacial rigidification brought about by 
$\tau$, and this is seen clearly in figure 5(a): the size of the cavity is largest for the full-Marangoni case, at the same stage of the dynamics for the three cases considered.
Figure 5. Interface location, pressure field together with a representation of the streamlines for the surfactant-free and full-Marangoni (
$|\tau |>0$) cases, and pressure on the axis of symmetry, are shown in columns one to four, respectively; (a–d) and (e–h) show the cavity collapse dynamics (at 
$t=0.440$, 0.496 and 0.500, for the surfactant-free, no-Marangoni (
$\tau =0$) and full-Marangoni cases) and the Worthington jet (at 
$t=0.507$, 0.594, and 0.619, for the surfactant-free, 
$\tau =0$, 
$|\tau | >0$ cases), respectively. The parameter values are the same as those used to generate figure 2.
The capillary pressure field shown in figure 5(b,c) for the surfactant-free and full-Marangoni cases is influenced heavily by the capillary pressure and, therefore, the local interfacial curvature and surface tension. From figure 5(b,c), it is seen that the pressure is highest immediately upstream of the capillary wave peak and this pressure gradient drives flow towards the lower-pressure region located at bottom of the cavity that coincides with the axis of symmetry. Furthermore, Marangoni stresses induce a recirculation zone close to the free surface, as shown in figure 5(c). In figure 5(d), we also show the axial distribution of the pressure at the axis of symmetry and this displays a peak at the interface due to capillarity. Owing to the presence of surfactant, the surface tension is reduced, which leads to a concomitant fall in the pressure, as illustrated via comparison of figure 5(b,c). The pronounced reduction in capillary pressure, due primarily to the accumulation of surfactant at the bottom of the cavity, is shown in figure 5(d).
As mentioned above, the convergence of the capillary wave on the bottom of the cavity leads to the formation of a Worthington jet (Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2019), which is shown in figure 5(e) for the surfactant-free, Marangoni-suppressed and full-Marangoni cases for 
$t=0.507$, 
$t=5.94$ and 
$t=0.619$, respectively; again, these times are chosen at nearly identical spatial locations of the jet tip in the axis of symmetry. The larger pressure gradient associated with the surfactant-free case, discussed above, leads to longer jets, with more pronounced bulbous regions at their tips, in comparison to the full-Marangoni case where the retarding Marangoni stresses induce a recirculation zone close to the jet base. Close inspection of figure 4(e) also reveals that the longest jets are associated with the 
$\tau =0$ (rather than the surfactant-free) case wherein there is no Marangoni-induced retardation and for which 
$La$ is lowest. As shown previously (Lai et al. Reference Lai, Eggers and Deike2018), decreasing 
$La$ leads to faster and thinner jets, with a greater propensity for breakup. As a result, a reduction in the number of ejected droplets is observed for the shorter, and slower, full-Marangoni jets: four (see figure 2), seven and nine droplets (the latter two not shown) for the full-Marangoni, surfactant-free and Marangoni-suppressed cases, respectively. The droplet numbers for surfactant-free and Marangoni-suppressed cases agree with Berny et al. (Reference Berny, Deike, Séon and Popinet2020). Finally, as shown in figure 5(h), there is an adverse pressure gradient in all cases, since capillarity drives flow from the bulbous region towards the jet base, which is largest for the surfactant-free case.
The immobilising effect of the Marangoni stresses can be seen in figure 6(a), in which we plot the kinetic energy, defined as 
$E_k=\int _V (\rho \boldsymbol {u}^2/2) \,\textrm {d}V$, where Marangoni stresses reduce the maximal and asymptotic values of 
$E_k$ in comparison to the surfactant-free and no-Marangoni cases. In figure 6(b), we observe further that the motion of the interface is retarded maximally when Marangoni stresses are enabled fully. Inspection of figure 6(c) also reveals that the interfacial area 
$A$, normalised by its initial value, 
$A_o$, reduces over time leading to large 
$\varGamma$ at the moment of interfacial vertical collapse at 
$t \sim 0.550$.
Figure 6. Temporal evolution of kinetic energy, the maximal vertical interface location and interfacial area (normalised by its initial value), (a–c), respectively, for the surfactant-free, full-Marangoni and no-Marangoni cases, for the same parameters as in figure 2.
We have measured the droplet characteristics, in terms of their volume and their average interfacial surfactant concentration, at time 
$t=1.76$, as shown in figure 2, and summarised the results in the table 1. The droplets are numbered depending on the order of ejection from the tip of the jet (i.e. the highest droplet in the computational domain corresponds to the first ejected droplet). The volume of each droplet is normalised with the initial volume of the spherical bubble. This shows the multi-scale nature of the phenomenon. Moreover, the size of the first droplet, 
$r_d/R_o \sim 6\,\%$, which agrees with the experimental observations of Blanchard & Woodcock (Reference Blanchard and Woodcock1957) and Tedesco & Blanchard (Reference Tedesco and Blanchard1979), the empirical power law dependence of Lewis & Schwartz (Reference Lewis and Schwartz2004) and the theoretical scaling relationship proposed by Gañán-Calvo (Reference Gañán-Calvo2017). The predicted surfactant concentration of the ejected droplets is higher in value than the concentration of the liquid reservoir. This agrees with the experimental studies of Blanchard & Syzdek (Reference Blanchard and Syzdek1972) which were carried out within the context of examining bacterial concentrations in bursting bubbles.
Table 1. Droplet characteristics of surfactant-laden case.
3.2. Parametric study
Here, we investigate the effect of system parameters on the dynamics of the bursting phenomenon beginning with the influence of 
$\beta _s$, which controls the relative strength of Marangoni stresses. As shown previously, the presence of surfactant tends to rigidify the interface and delays the cavity collapse (see figure 7a). By increasing 
$\beta _s$, Marangoni stresses are strengthened, leading to greater rigidification of the interface, and consequently, higher retardation of the cavity collapse (see figure 7a). This collapse retardation is also accompanied by a lower surfactant convection towards the point of singularity (see figure 7c). Moreover, the strengthening of Marangoni stresses via increase in 
$\beta _s$ leads to dampening of the tangential velocity along the interface (see figure 7e), to thinner jets and fewer ejected droplets (e.g. seven and four droplets for the lower and upper ends of the elasticity parameter 
$\beta _s$).
Figure 7. Effect of the elasticity number, 
$\beta _s$, on the bursting dynamics for 
$La=2\times 10^4$, 
$Bo=10^{-3}$, 
$Pe_s=1$, 
$Bi =0.01$, 
$k=1$ and 
$\varGamma _o =1/4$. (a,b) Represent the interfacial shape, (c,d) the interfacial surfactant concentration 
$\varGamma$ along the arc length 
$s$ and (e,f) the tangential velocity 
$u_t$ along 
$s$, during the cavity collapse at 
$t=0.500$ and the Worthington jet at 
$t=0.800$, respectively.
We also analyse the effect of surfactant solubility on the bursting dynamics by varying the value of the Biot number, 
$Bi$. For the lower 
$Bi$ values, the sorptive time scales are larger than those associated with the interfacial dynamics and the Marangoni stresses, and consequently, the system behaves as if the surfactant were effectively insoluble. From the representation of the interfacial location of the cavity, we see that the desorptive time scales dominate the system and, consequently, a reduction of 
$\varGamma$ is observed as 
$Bi$ increases (see figure 8a) with a concomitant weakening of the Marangoni stresses which oppose cavity collapse. The increase in interfacial mobility, characterised by a rise in the magnitude of 
$u_t$, is also seen in figure 8(e).
Figure 8. Effect of the Biot number, 
$Bi$, on the bursting dynamics for 
$La=2\times 10^4$, 
$Bo=10^{-3}$, 
$Pe_s=1$, 
$\beta _s =0.9$, 
$k=1$ and 
$\varGamma _o =1/4$. (a,b) Represent the interfacial shape, (c,d) the interfacial surfactant concentration 
$\varGamma$ along the arc length 
$s$ and (e,f) the tangential velocity 
$u_t$ along 
$s$, during the cavity collapse at 
$t=0.500$ and the Worthington jet at 
$t=0.800$, respectively.
We also examine the effect of solubility on the jet dynamics whereupon we observe that, at the lower end of the Biot number, the jet is slowest due to interfacial rigidification brought about by 
$\tau$, as can be seen via the damping of 
$u_t$, shown in figure 8(f). Higher gradients in 
$\varGamma$ are observed for the weakly soluble case (see figure 8d), which leads to the highest strength of 
$\tau$. With increasing 
$Bi$, more interfacial surfactants migrate to the bulk and as a result gradients in 
$\varGamma$ are reduced, leading to an overall weakening of the Marangoni stresses. Finally, we study the effect of solubility on the number of ejected droplets. As before, the desorption of surfactant from the interface as 
$Bi$ increases yields weakening of 
$\tau$ and, subsequently, the surfactant does not prevent ‘end pinching’ of the jet tip from which a larger number of droplets are ejected than in the low 
$Bi$ case; for 
$Bi=1$, seven droplets are ejected in comparison with the four droplets produced by the weakly soluble case characterised by 
$Bi=0.01$.
4. Concluding remarks
The effect of Marangoni-induced flow, brought about by the presence of surfactant, on the dynamics of a bubble bursting through an interface was studied using a hybrid front-tracking/level-set method. Our results indicate that a surfactant-covered toroidal capillary wave forms, following the collapse of the cavity, whose motion is retarded by the surfactant-induced Marangoni stresses; these stresses drive flow from regions of high surfactant concentration (low surface tension) to low concentration (high tension) regions. The immobilising effect of the surfactants due to the Marangoni stresses is also observed via the marked reduction in the system kinetic energy and the generation of shorter, and slower, Worthington jets. The breakup of these jets is accompanied by the formation of fewer droplets in comparison to the surfactant-free case. This behaviour is associated with the ‘end-pinching’ mechanism of the jet tip and the presence of surfactant promotes the neck re-opening through Marangoni flow, induced by the formation of surfactant concentration gradients, and not only via lowering of the mean surface tension value. Finally, we have examined the role of the elasticity number 
$\beta _s$ which is an important parameter that controls the strengthening or weakening of Marangoni stresses along the interface. An increase in 
$\beta _s$ leads to a reduction in tangential velocity, and consequently a retardation in interfacial motion. In a similar way, the surfactant solubility, via the variation of Biot number 
$Bi$, was also examined in this study. An increase in this parameter leads to a rise in the rate of mass transfer from the bulk to the interface, resulting in the overall weakening of Marangoni stresses, and consequently, a larger number of droplets as Marangoni stresses cannot prevent ‘end-pinching’ of the jet tip.
Future research avenues for study are to perform numerical simulations featuring three-dimensional behaviours occurring for large Bond numbers (
$Bo\geqslant 0.5$). It has been shown by Ghabache (Reference Ghabache2015) that, at large values of 
$Bo$, the bursting outcome can be either axisymmetric or purely three-dimensional, producing an oblique Worthington jet. Both surfactant-free and surfactant-laden bubbles should be examined with a similar formulation to that described in § 2. Thus, the bubble initialisation will require taking into account the entirety of the same: the immersed cavity, the cap and the meniscus. If the hole formation is located anywhere else on the cap, a three-dimensional cap retraction will certainly affect the behaviour of the collapsing cavity to finally culminate in an inclined Worthington jet, and complex droplet detachment.
Acknowledgements
The numerical simulations were performed with code BLUE (Shin et al. Reference Shin, Chergui and Juric2017) and the visualisations have been generated using ParaView. We are grateful to the referees for many helpful and constructive comments.
Funding
This work is supported by the Engineering and Physical Sciences Research Council (EPSRC), United Kingdom, through a studentship for C.R.C-A. in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London funded by the EPSRC (EP/L015579/1), and through the EPSRC MEMPHIS (EP/K003976/1) and PREMIERE (EP/T000414/1) Programme Grants. C.R.C-A. also acknowledges the funding and technical support from BP through the BP International Centre for Advanced Materials (BP-ICAM), which made this research possible. O.K.M. also acknowledges funding from PETRONAS and the Royal Academy of Engineering for a Research Chair in Multiphase Fluid Dynamics, and the PETRONAS Centre for Engineering of Multiphase Systems. We also acknowledge HPC facilities provided by the Research Computing Service (RCS) of Imperial College London for the computing time. D.J. and J.C. acknowledge support through computing time at 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. 2020A0082B06721.
Declaration of interests
The authors report no conflict of interest.
