1. Introduction
Droplets impacting on liquid surfaces are observed in a broad range of natural phenomena and industrial applications. Consequently, this area of study drives significant scientific interest; for comprehensive reviews please see e.g. Yarin (Reference Yarin2006), Josserand & Thoroddsen (Reference Josserand and Thoroddsen2016) and Cheng, Sun & Gordillo (Reference Cheng, Sun and Gordillo2022). The complex topological features of droplet impact have attracted the fluid mechanics community since the first ground-breaking experiments from Worthington (Reference Worthington1908) and Edgerton (Reference Edgerton1977) over a century ago. However, it was not until the advent of high-speed imaging and high-resolution computational tools that it became possible to capture the early dynamics of the broad range of scales involved in such phenomena (Thoroddsen Reference Thoroddsen2002; Thoraval et al. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012; Thoroddsen et al. Reference Thoroddsen, Thoraval, Takehara and Etoh2012; Zhang et al. Reference Zhang, Toole, Fezzaa and Deegan2012; Agbaglah & Deegan Reference Agbaglah and Deegan2014; Josserand, Ray & Zaleski Reference Josserand, Ray and Zaleski2016; Fudge, Cimpeanu & Castrejón-Pita Reference Fudge, Cimpeanu and Castrejón-Pita2021).
The dynamics of droplet impact on liquids is extremely diverse and is determined by a wide number of factors including the impact speed and the fluid properties, as well as whether the impact occurs on either a deep pool or a thin layer of fluid (Deegan, Brunet & Eggers Reference Deegan, Brunet and Eggers2007; Josserand et al. Reference Josserand, Ray and Zaleski2016; Che & Matar Reference Che and Matar2017). At high impact speeds, the latter situation results in the formation of an upward Worthington jet (Gekle & Gordillo Reference Gekle and Gordillo2010), whereas the former case gives rise to the well-known milkdrop coronet, which is the subject of this paper (Thoroddsen Reference Thoroddsen2002; Deegan et al. Reference Deegan, Brunet and Eggers2007; Josserand et al. Reference Josserand, Ray and Zaleski2016). Instants after impact on to a thin liquid layer, inertia prompts the formation of an ejecta sheet. After some time, its edge forms a capillary-driven cylindrical-shaped rim that evolves through various stages: first, a combination of Rayleigh–Taylor and Rayleigh–Plateau instabilities grow on the rim (see Zhang et al. Reference Zhang, Brunet, Eggers and Deegan2010; Agbaglah, Josserand & Zaleski Reference Agbaglah, Josserand and Zaleski2013), then nonlinearities develop, leading to breakup (through end pinching, see Wang & Bourouiba Reference Wang and Bourouiba2018), and the formation of a cascade of droplet sizes.
Under most natural or industrial conditions, streams are contaminated with surfactants (i.e. surface-active agents), whose concentration variations lead to surface tension gradients, which in turn result in the formation of Marangoni stresses and surface viscosities (see Manikantan & Squires Reference Manikantan and Squires2020). Previous studies have demonstrated the role of surfactants in capillary singularities (see for instance, Ananthakrishnan & Yeung Reference Ananthakrishnan and Yeung1994; Craster, Matar & Papageorgiou Reference Craster, Matar and Papageorgiou2002; Liao et al. Reference Liao, Subramani, Franses and Basaran2004; Kamat et al. Reference Kamat, Wagoner, Thete and Basaran2018; Constante-Amores et al. Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020); here, we study the role of surfactant-induced flows in the late stages of splashing. Che & Matar (Reference Che and Matar2017) demonstrated the role of surfactants on the impact dynamics onto thin liquid films for low and moderate Weber and high Reynolds numbers, e.g. 
$2.6< We<532$ and 
$3000< Re<8654$, respectively (
$Re=\rho _lUD/\mu _l$ and 
${We=\rho _l U^2 D/\sigma }$, where 
$D$, 
$\rho _l$, 
$\mu _l$ 
$\sigma$ and 
$U$ stand for the droplet diameter, and its characteristic density, viscosity, surface tension and velocity, respectively). In this past work, the authors showed that the presence of surfactants affects the post-impact dynamics such as the propagation of capillary waves, the growth of the crown and the generation of secondary droplets. However, Marangoni stresses, surfactant concentration dynamics and other important insights into the impact phenomena were not investigated.
The vast majority of past works have focused on the physical understanding of the early dynamics of the ejected sheet using axisymmetric simulations, which is a valid assumption in the limiting case of 
$t\ll D/U$ (where 
$t$ stands for time) – see for example Josserand & Zaleski (Reference Josserand and Zaleski2003), Josserand et al. (Reference Josserand, Ray and Zaleski2016) and Agbaglah & Deegan (Reference Agbaglah and Deegan2014). In contrast, the study of the crown-splash regime requires the use of full three-dimensional (3-D) numerical simulations in order to accommodate the natural occurrence of possible symmetry-breaking events, such as the growth of transverse instabilities in the rim, or rim breakup into droplets. However, 3-D simulations remain a numerical challenge (Gueyffier & Zaleski Reference Gueyffier and Zaleski1998) and have not often been reported in the literature (Liang et al. Reference Liang, Zhang, Chen, Chen and Shen2019; Reijers et al. Reference Reijers, Liu, Lohse and Gelderblom2019; Chen et al. Reference Chen, Shu, Wang and Yang2020; Xavier et al. Reference Xavier, Zuzio, Averseng and Estivalezes2020).
The present study aims to unravel the role of surfactant-induced Marangoni stresses on the interfacial dynamics during the impact of a surfactant-free droplet on a surfactant-laden thin layer of the same liquid in the crown-splash regime. Here, we perform 3-D calculations of splashing in the presence of surfactants by using a hybrid front-tracking method to track the interface location. Section 2 presents the governing dimensionless parameters, the problem geometry, initial and boundary conditions and the numerical validation. Section 3 provides a discussion of the results and concluding remarks are given in § 4.
2. Problem formulation and numerical method
High-resolution simulations were performed by solving the two-phase incompressible Navier–Stokes equations with surface tension in a 3-D Cartesian domain 
$\boldsymbol {x} = (x, y, z)$ (see figure 1a). The interface between the gas and liquid is described by a hybrid front-tracking/level-set method, where (insoluble) surfactant transport is resolved at the interface (Shin et al. Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018). Here, and in what follows, all variables are made dimensionless (represented by tildes) using
\begin{equation} \tilde{\boldsymbol{x}}=\frac{\boldsymbol{x}}{D_o}, \quad \tilde{t}=\frac{t}{t_{r}}, \quad \tilde{\boldsymbol{u}}=\frac{\boldsymbol{u}} {U}, \quad \tilde{p}=\frac{p}{\rho_{{l}} U^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$, liquid viscosity, 
$\mu _l$, surface tension, 
$\sigma$, surfactant-free surface tension, 
$\sigma _s$, 
$U$ is the drop impact velocity and 
$D_o$ its initial diameter; hence, the characteristic time scale is 
$t_{r}= D_o/U$. The interfacial surfactant concentration, 
$\varGamma$, is scaled by the saturation interfacial concentration, 
$\varGamma _{\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} ={-}\frac{\mathbf{i}_z }{Fr^2}~ + \frac{1}{Re} \boldsymbol{\nabla}\boldsymbol{{\cdot}} [ \tilde{\mu} (\boldsymbol{\nabla} \tilde{\boldsymbol{u}} +\boldsymbol{\nabla} \tilde{\boldsymbol{u}}^T) ] \nonumber\\ \hspace{163pt}+\,\frac{1}{\rm We} \int_{\tilde{A}\tilde{(t)}} (\tilde{\sigma} \tilde{\kappa} \boldsymbol{n} + \nabla_s \tilde{\sigma}) \delta (\tilde{\boldsymbol{x}} - \tilde{\boldsymbol{x}}_{f} )\,{\rm d}\tilde{A}, \end{gather}
\begin{gather}\frac{\partial \tilde{\varGamma}}{\partial \tilde{t}}+\nabla_s \boldsymbol{{\cdot}} (\tilde{\varGamma}\tilde{\boldsymbol{u}}_{t})=\frac{1}{Pe_s} \nabla^2_s \tilde{\varGamma}, \end{gather}
where the density and viscosity are given 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 Heaviside function, which is zero in the gas phase and unity in the liquid phase, while the subscripts ‘
$g$’ and ‘
$l$’ designate the gas and liquid phases, respectively, and 
$\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 and 
$\kappa$ the curvature. The interfacial gradient is given by 
$\nabla _s=({\boldsymbol {I}}-\boldsymbol {n}\boldsymbol {n})\boldsymbol {{\cdot }} \boldsymbol {\nabla }$ wherein 
$\boldsymbol {I}$ is the identity tensor and 
$\boldsymbol {n}$ is the outward-pointing unit normal. In addition, 
$\delta$ is a Dirac delta function, equal to unity at the interface and zero otherwise, and 
$\tilde {A} (\tilde {t})$ is the time-dependent interface area. As justified by the final paragraph of Stone (Reference Stone1990), in our frame of reference 
$\boldsymbol {u} \boldsymbol {{\cdot }} \boldsymbol {n}=0$, which gives rise to 2.4. The dimensionless groups that appear in the governing equations are defined as
\begin{equation} \left.\begin{array}{c} Re =\dfrac{\rho_l U D_o}{\mu_l}, \quad We =\dfrac{\rho_l U^2 D_o}{\sigma_s}, \quad Fr =\dfrac{ U }{\sqrt{g D_o}}, \\ Pe_s=\dfrac{ U D_o}{\mathcal{D}_s},\quad \beta_s= \dfrac{\Re \mathcal{T} \varGamma_\infty}{\sigma_s}, \end{array}\right\} \end{equation}
where 
$Re$, 
$We$, 
$Fr$ and 
$Pe_s$ stand for the Reynolds, Weber, Froude and (interfacial) Péclet numbers. Gravity is negligible during the impact as indicated by the Froude number value, i.e. 
$Fr \sim {O}(10^2)$ (similar assumptions were made by Deegan et al. Reference Deegan, Brunet and Eggers2007). The parameter 
$\beta _s$ is the surfactant elasticity number that is a measure of the sensitivity of the surface tension, 
$\sigma$, to the interface surfactant concentration, 
$\varGamma$. Here, 
$\Re$ is the ideal gas constant value (
$\Re = 8.314$ J K
$^{-1}$ mol
$^{-1}$), 
$\mathcal {T}$ denotes temperature and 
$\mathcal {D}_s$ stands for the interfacial diffusion coefficient.
Figure 1. Schematic representation of the flow configuration, and validation of the numerical procedure: (a) initial configuration, highlighting the computational domain of size (not to scale) 
$(8D_o \times 8D_o \times 4D_o)$ in a 3-D Cartesian domain, 
$\boldsymbol {x} = (x, y, z)$, with a resolution of 
$768^2\times 384$; (b) validation of the numerical results (solid lines) against the experimental data (dashed lines) for the crown rim diameter temporal evolution of a surfactant-free case reported in Che & Matar (Reference Che and Matar2017): effect of varying the dimensionless film depth, 
$h=H/D_o$, while keeping the other parameters constant (
$Re=7514$, 
$We = 249$ and 
$Fr = 13.65$).
The nonlinear Langmuir equation is used to describe 
$\sigma$ in terms of 
$\varGamma$, this is
\begin{equation} \tilde{\sigma}=1 + \beta_s \ln{ (1 -\tilde{\varGamma})}. \end{equation}
The Marangoni stress, 
$\tilde {\tau }$, is expressed as a function of 
$\tilde {\varGamma }$ as
\begin{equation} \frac{1}{\rm We} \nabla_s \tilde{\sigma}\boldsymbol{{\cdot}} {\boldsymbol{t}} \equiv \frac{\tilde{\tau}}{\rm We} ={-} \frac{\rm Ma}{(1-\tilde{\varGamma})} \nabla_s\tilde{\varGamma}\boldsymbol{{\cdot}} {\boldsymbol{t}}, \end{equation}
where 
$Ma=\beta _s/We=Re \mathcal {T} \varGamma _\infty /\rho _l U^2 D_o$ is the Marangoni parameter and 
$\boldsymbol {t}$ is the unit tangent to the interface. Finally, a dimensionless thickness of the liquid layer, 
$h=H/D_o$ is defined as the ratio between the liquid film thickness and the droplet diameter. Tildes are dropped henceforth.
2.1. Numerical set-up, validation and parameters
Figure 1(a) highlights the geometry and the initial and boundary conditions of the problem, which closely follow previous work by Josserand et al. (Reference Josserand, Ray and Zaleski2016), Agbaglah & Deegan (Reference Agbaglah and Deegan2014) and Fudge et al. (Reference Fudge, Cimpeanu and Castrejón-Pita2021). A drop of initial size 
$D_o$ with velocity 
$U$ impacts a uniform layer of thickness 
$H$ of the same liquid. The size of the computational domain is 
$8D_0 \times 8D_0 \times 4D_0$, which is sufficiently large to avoid any effects of artificial reflections from the boundaries. The centre of the drop is located at a small distance above the pool surface (e.g. an initial separation of 
$0.05D_0$). A no-slip boundary and no-penetration conditions are assumed for the bottom wall of the domain, whereas a no-penetration boundary condition is prescribed for the top and lateral domain boundaries (similar to previous work by Batchvarov et al. Reference Batchvarov, Kahouadji, Constante-Amores, Norões Gonçalves, Shin, Chergui, Juric, Craster and Matar2021).
Figure 1(b) illustrates the ability of the numerical technique to reproduce the temporal evolution of the crown experimentally obtained and reported by Che & Matar (Reference Che and Matar2017). After impact, inertia drives the rapid ejection of a vertical sheet from the pool, then capillarity prompts the formation of a rim and then a crown. We observe excellent agreement for 
$h=0.22$, while a small offset is observed for the other two film thicknesses. The discrepancy could be attributed to the experimental conditions not considered here, such as the surface waves produced during drop formation or air-induced flows during droplet fall. Nonetheless, the agreement between the numerical method and experimental results is visible. Additional validation cases are found in Appendix A. The numerical technique has also been validated for drop–interface coalescence, which can be considered as a limiting case of drop impact onto a pool, i.e. 
$U=0$ (see Constante-Amores et al. Reference Constante-Amores, Batchvarov, Kahouadji, Shin, Chergui, Juric and Matar2021a for more details). The nonlinear interfacial dynamics has been validated for the capillary breakup of liquid threads (see Constante-Amores et al. Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020, Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021b for more details). The numerical validation of the surfactant equations has been previously presented in Shin et al. (Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018).
In terms of mesh resolution, we have ensured that our numerical simulations are mesh independent, and as a consequence, the numerical findings do not change by decreasing the cell size for a resolution of 
$768^2 \times 384$ (i.e. 
$D_0/\Delta \boldsymbol {x}=96$ cells, see Appendix B). Additionally, liquid volume and surfactant mass conservation are met under errors of less than 
$10^{-1}\%$ and 
$10^{-2}\%$, respectively. Extensive mesh studies for surface-tension-driven phenomena using the same computational method can be found elsewhere, e.g. Batchvarov et al. (Reference Batchvarov, Kahouadji, Constante-Amores, Norões Gonçalves, Shin, Chergui, Juric, Craster and Matar2021, Reference Batchvarov, Kahouadji, Magnini, Constante-Amores, Craster, Shin, Chergui, Juric and Matar2020) and Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020, Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021b, Reference Constante-Amores, Chergui, Shin, Juric, Castrejón-Pita and Castrejón-Pita2022, Reference Constante-Amores, Abadie, Kahouadji, Shin, Chergui, Juric, Castrejón-Pita and Matar2023).
The dimensionless values for the investigated phenomenon are consistent with experimentally realisable systems. We assume a water–air system, where the density and viscosity ratios are 
$\rho _g / \rho _l = 1.2 \times 10^{-3}$ and 
$\mu _g / \mu _l = 0.018$, respectively. The targeted flow conditions for the surfactant-free base case are based on Deegan et al. (Reference Deegan, Brunet and Eggers2007), who built a phase diagram in the 
$We-Re$ space for 
$h=0.2$. Our study focuses on the ‘crown-splash’ regime (i.e. 
$We>500$ and 
$Re <1100$), as we aim to stay away from the ‘microdroplet-splash’ regime (i.e. 
$We>500$ and 
$Re <1000$). The latter regime is characterised by the formation of a rapid ejecta sheet which undergoes capillary singularities to form microdroplets, requiring an extremely large resolution to capture droplet formation from the ejecta sheet and from the crown. Thus, we have selected 
$Re=1000$, 
$We=800$ and 
$h=0.2$ as the targeted surfactant-free case.
We have examined the case of adding surfactant to the thin liquid film but not the impacting droplet as it is less challenging to control the characteristics of surfactant-free drops experimentally. We have considered insoluble surfactants whose critical micelle concentration (CMC), 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, we have explored the range of 
$0.1 <\beta _s<0.9$ which corresponds to 
${O}(10^{-7})<\textrm {CMC}<{O}(10^{-6})$ mol m
$^{-2}$, for typical values of 
$\sigma _s$. We have set 
$Pe_s=10^2$ following Batchvarov et al. (Reference Batchvarov, Kahouadji, Magnini, Constante-Amores, Craster, Shin, Chergui, Juric and Matar2020) and Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020) who demonstrated that the interfacial dynamics is weakly dependent on 
$Pe_s$ beyond this value.
For a water droplet of a typical size of 
$D_o \sim {O}(10^{-3})$ m, the impact time scale 
$T_{imp}=t_r=D_o/U \sim {O}(10^{-4}\unicode{x2013}10^{-3})$ s; whereas the Marangoni time scale 
$T_{\tau }$ can be estimated by a balance between the Marangoni and the viscous stresses, resulting in 
$T_{\tau } \sim \mu _l D^2 / (h \Delta \sigma ) \sim {O}(10^{-4}\unicode{x2013}10^{-3})$ s (see Che & Matar Reference Che and Matar2017 for more details). Therefore, surfactant-induced Marangoni stresses will play a crucial role in the flow physics. A discussion of the results is presented next.
3. Results
First, we discuss the early-time dynamics of the surfactant-free impact at 
$Re=1000$ and 
$We=800$. Figure 2(a) shows the 3-D representation of the interface at early times, i.e. to show the formation of transverse instabilities growing at the rim. As seen in this figure, instants after impact, we observe the formation of an axisymmetric thin layer of ejecta fluid which draws fluid from the droplet eventually giving rise to a Peregrine sheet, in agreement with Deegan et al. (Reference Deegan, Brunet and Eggers2007) and Josserand et al. (Reference Josserand, Ray and Zaleski2016). As time evolves, the Peregrine sheet grows upwards and expands radially while its film thickness reduces; here, the radial coordinate 
$r$ is defined at the initial point of impact. A closer inspection of the interfacial shape of the ejecta sheet at 
$t=0.5$ demonstrates that its orientation is nearly horizontal, i.e. parallel to the pool. We also see surface tension forces inducing the formation of a hemispherical tip at the edge of the ejecta, which leads to a pressure gradient between the tip and the adjacent sheet, which results in fluid motion from the tip towards the sheet. This behaviour leads to the formation of a sheet perpendicular to the pool at longer times; this is in agreement with the experimental observations by Zhang et al. (Reference Zhang, Brunet, Eggers and Deegan2010) and Deegan et al. (Reference Deegan, Brunet and Eggers2007). Figure 2(b) shows the temporal evolution of the rim size for the surfactant-free case at the crossing of the 
$x$–
$z$ plane (
$y=4.90$). The peak observed in the range 
$t=3\unicode{x2013}5$ is due to the growth of corrugations in the rim. As the rim elongates to form the crown, we observe the growth of an azimuthal instability on the surface. Our results indicate that these naturally occurring rim instabilities select wavelengths that are consistent with the most unstable Rayleigh–Plateau (RP) instability. We have measured the average distance (
$\lambda$) between the ‘undulations’ or crests seen along the rim. We then normalised the result by the most unstable RP instability wavelength, given by 
$\lambda ^*= 2 {\rm \pi}/k= 9.016 a$, where 
$k$ and 
$a$ are the wavenumber of the fastest-growing instability mode and the rim radius, respectively (Zhang et al. Reference Zhang, Brunet, Eggers and Deegan2010; Agbaglah et al. Reference Agbaglah, Josserand and Zaleski2013). In figure 2(c) we plot the selected wavelength as a function of the rim radius; for the cases studied here, 
$\lambda /\lambda ^* = 1.0 \pm 0.1$. These results are in agreement with the experimental results from Zhang et al. (Reference Zhang, Brunet, Eggers and Deegan2010) and the theoretical and computational work from Agbaglah et al. (Reference Agbaglah, Josserand and Zaleski2013) and Constante-Amores et al. (Reference Constante-Amores, Chergui, Shin, Juric, Castrejón-Pita and Castrejón-Pita2022). A closer inspection of the top image of figure 2(c) shows that the interfacial corrugations are not always evenly distributed along the rim, as observed in experiments (Thoraval et al. Reference Thoraval, Takehara, Etoh, Popinet, Ray, Josserand, Zaleski and Thoroddsen2012). In particular, we observe large wavelengths, which are the result of interactions between adjacent instability modes, which is consistent with experimental observations.
Figure 2. Wavelength selection in the crown-splash regime. Panel (a) presents the early interfacial dynamics through a 3-D representation of the predicted interface location for 
$Re=1000$ and 
$We=800$, at 
$t=(0.5, 1.0, 1.5)$ corresponding to columns one to three, respectively. Panel (b) plots the temporal evolution of the rim radius. Panel (c) shows the selection of the wavelength 
$\lambda$ normalised with the theoretical RP instability, 
$\lambda ^*$ as a function of the local rim radius, 
$a$, at early times of the simulation for 
$t<5$. A 3-D representation of the crown for the dimensionless time 
$t=7.5$ is also presented, highlighting the predicted wavelength between crests.
Figure 3(a–d) shows a 3-D representation of the temporal evolution of the interface, from the time instabilities become visible to the time droplets breakup by end pinching. At these longer times (
$t>5$), the rim is weakly affected by the dynamics occurring at the impacting point, and interfacial RP instabilities trigger the development of ligaments (see figure 3a). With increasing time, liquid ligaments grow perpendicular to the rim, resulting in a local pressure-gradient-induced flow from the ligament to the tip, triggering the formation of a capillary-induced blob and a neck (displayed in figure 3b). The neck thins and stretches, due to a capillary-driven flow, eventually resulting in the capillary pinch-off of a drop (the well-known ‘end-pinching’ mechanism proposed by Stone & Leal Reference Stone and Leal1989), see figure 3(c). We note that the interfacial dynamics closely resembles experimental observations by Zhang et al. (Reference Zhang, Brunet, Eggers and Deegan2010), Josserand & Zaleski (Reference Josserand and Zaleski2003) and Deegan et al. (Reference Deegan, Brunet and Eggers2007), as well as the dynamics of retracting sheets reported by Wang & Bourouiba (Reference Wang and Bourouiba2018) and Constante-Amores et al. (Reference Constante-Amores, Chergui, Shin, Juric, Castrejón-Pita and Castrejón-Pita2022).
Figure 3. Effect of 
$\beta _s$ on the drop impact dynamics for insoluble surfactants. Spatio-temporal evolution of the 3-D interface shape for surfactant-free cases, (a–d), and surfactant-laden cases for 
$\beta _s = (0.1, 0.5, 0.7)$ corresponding to panels (e–h), (i–l) and (m–p), respectively. Here, the dimensionless parameters are 
$Re=1000$, 
$We=800$ and 
$h=0.2$. For the surfactant-laden cases, 
$Pe_s=100$ and 
$\varGamma =\varGamma _\infty /2$, the colour indicates the value of 
$\varGamma$ and the legend is shown in panel (e).
Next, we proceed to study the role of surfactants in the flow dynamics; as illustrated by figure 3(e–p) where we observe the droplet impact at various elasticity parameters 
$\beta _s$. As seen in this figure, regardless of the value of 
$\beta _s$, large surfactant concentration gradients are found on the outer interface of the ejecta sheet. This result supports the hypothesis that the ejecta sheet is generated at the interface of the surfactant-laden pool. On the other hand, the inner interface of the ejecta sheet arises from the surfactant-free drop, i.e. 
$\varGamma$ vanishes on the inner surface, see figure 3(e). This observation, that the fluid sources of the inner and outer surfaces of the ejecta sheet are the pool and the droplet, respectively, is in good agreement with the work by Josserand et al. (Reference Josserand, Ray and Zaleski2016). A close inspection of the distribution of 
$\varGamma$ in the outer sheet shows that 
$\varGamma$-values are maximum at the base of the ejecta, and their values reduce upwards (see for example panel (e) of figure 3). This is the result of an increase of the interfacial area, that leads to a dilution of surfactant at the interface, thus the surfactant concentration (surface tension) decreases (increases). In addition, the interfacial expansion results in large convective effects that transport surfactant towards the rim; see figure 3(e), in which larger values of 
$\varGamma$ are found in the rim than in its adjacent sheet.
The increase (decrease) of 
$\varGamma$ (
$\sigma$) implies that a large surface deformation is required to satisfy the normal stress balance at the interface. Similar to the surfactant-free case, capillary-induced flow results in the formation of a rim on the sheet edge, whose thickening leads to the onset of destabilisation. This is in agreement with previous studies from Asaki, Thiessen & Marston (Reference Asaki, Thiessen and Marston1995), who reported a surfactant-driven dampening effect on the capillary waves.
In addition, we expect higher (lower) values of 
$\varGamma$ at the rim wave crests (instability waves) as they belong to a radially converging (diverging) region, which has the effect of locally increasing (reducing) 
$\varGamma$. The gradients in 
$\varGamma$ result in surfactant-induced Marangoni stress flows from the lower-
$\sigma$ radially converging to the higher-
$\sigma$ radially diverging regions. This slows down the interfacial dynamics, and the formation of instabilities at the rim, in agreement with the idealised case presented by Constante-Amores et al. (Reference Constante-Amores, Chergui, Shin, Juric, Castrejón-Pita and Castrejón-Pita2022). By increasing 
$\beta _s$ the surfactant distribution along the interface is enhanced, as displayed in figures 3(i–l) and 3(j–m) for 
$\beta _s=0.5$ and 
$\beta _s=0.7$, respectively.
Surfactant-induced Marangoni stresses only delay the development of the ligament from the rim; these eventually break up via the end-pinching mechanism to form droplets, as shown in figure 3(g,k,o). These results offer an explanation for the experimental observations of Che & Matar (Reference Che and Matar2017). Figure 4 displays the selected instability wavelength in terms of 
$\beta _s$ and 
$a_0$. Indeed, a close inspection indicates that the surfactant does not seem to greatly affect the selection of the wavelength; the predicted 
$\lambda$ is consistent with the most unstable RP instability at a ratio of 
$\lambda /\lambda ^* = 1.0 \pm 0.1$. Consequently, according to our results, the crown dynamics is mostly driven by the RP instability, although longer ligaments developed prior to break up in the presence of surfactants, as observed in figure 3(p). This indicates that Marangoni stresses promote the retardation from end pinching, which, under certain circumstances may even lead to the neck re-opening, further delaying the splash, a phenomenon previously demonstrated by Kamat et al. (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2020) and Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020, Reference Constante-Amores, Chergui, Shin, Juric, Castrejón-Pita and Castrejón-Pita2022).
Figure 4. Effect of the elasticity parameter, 
$\beta _s$, on the selection of wavelength of the undulations normalised with the theoretical RP instability, 
$\lambda ^*$, as a function of the local rim radius, 
$a$.
Two-dimensional projections, in the 
$x$–
$z$ plane (
$y=4)$, of the interfacial shape, 
$\varGamma$, 
$\tau$ and the radial component of the interfacial velocity 
$u_{tz}$ in terms of 
$\beta _s$ and a fixed value of 
$Pe_s$, are displayed in figures 5 and 6 corresponding to 
$t=5$ and 
$t=20$, respectively. As described above, at early times, i.e. 
$t=5$, the drop impact results in a non-uniform interfacial surfactant distribution, with high surface concentrations at the base of the ejected sheet. As seen, the interfaces of the surfactant-free and the 
$\beta _s=0.1$ surfactant-laden cases are practically indistinguishable. However, as the surfactant concentration increases, 
$\varGamma$ distributions trigger a surfactant-driven dynamics from the sheet base to its edge (see figure 5b). As 
$\beta _s$ increases, long ejecta sheets are promoted, resulting in the reduction of the rim and sheet thickness due to mass conservation (displayed in the magnified panel of figure 5a). The increase in 
$\beta _s$ also enhances the 
$\varGamma$ distribution and increases the value of 
$\tau$ (see the increase of 
$\tau$ with 
$\beta _s$ in figure 5c). Furthermore, figure 5(c) highlights the presence of large 
$\tau$ peaks on both sides of the rim due to the accumulation of 
$\varGamma$, and sudden drops at the base of the rim, suggesting a fluid transport from the rim to the sheet. In addition, a larger 
$\tau$ maximum is observed from the side of the inner sheet as the inner sheet is devoid of 
$\varGamma$, leading to an uneven flow motion inside the sheet. The variation of the streamwise interfacial tangential velocity 
$u_{t_{z}}$, along the arc length, is shown in figure 5(d). This figure, at the early stages of the dynamics, exhibits 
$u_t>0$ over the majority of the domain due to the dominance of the vertical radial expansion of the ejecta sheet. A strong variation of 
$u_{t_{z}}$ is observed around the rim with 
$u_{t_{z}}$ tending to zero at a sufficiently large distance away from the impactregion.
Figure 5. Effect of the elasticity parameter, 
$\beta _s$, on the flow dynamics at 
$t=5$. Two-dimensional projections of the interface, 
$\varGamma$, 
$\tau$ and 
$u_{tz}$ in the 
$x$–
$z$ plane (
$y=4$) are shown in (a–d), respectively. In panel (a), a magnified view of the ejecta sheet is also presented. Note that the abscissa in (a) corresponds to the 
$x$ coordinate, and in (b–d) to the arc length, 
$s$. The arc length 
$s$ corresponds to the 
$x$–
$z$ plane (
$y = 4$) intersecting the interface; 
$s$ has been normalised on the full extent of 
$s$ associated with the length of the impact region in each case. The diamond shapes in panels (a,b) indicate the location of the crown. All parameters remain unchanged from figure 3. Surf., surfactant.
Figure 6. Effect of the elasticity parameter, 
$\beta _s$, on the flow dynamics at 
$t=20$. Two-dimensional projections of the interface, 
$\varGamma$, 
$\tau$ and 
$u_{tz}$ in the 
$x$–
$z$ plane (
$y=4$) are shown in (a–d), respectively. Note that the abscissa in (a) corresponds to the 
$x$ coordinate, and in (b–d) to the arc length, 
$s$. The arc length 
$s$ corresponds to the 
$x$–
$z$ plane (
$y = 4$) intersecting the interface; 
$s$ has been normalised on the full extent of 
$s$ associated with the length of the impact region in each case. The diamond shapes in panels (a,b) indicate the location of the crown. All parameters remain unchanged from figure 3.
The influence of 
$\beta _s$ at long times, 
$t=20$, is examined in figure 6. Qualitatively, the surfactant effects are particularly evident as the dynamics has been slowed down by the surface rigidification brought about by the Marangoni stresses, as can be seen in figure 6(d) (the overall magnitude of 
$u_{t_z}$ is smaller with increasing 
$\beta _s$). As 
$\beta _s$ decreases (
$\sigma$ increases), the surfactant is highly convected as 
$\tau$ tends to zero, lowering the interfacial tension and consequently increasing surface deformation. As 
$\beta _s$ increases, the gradient in 
$\varGamma$ (
$\tau$) decreases (increases) near the rim, and the Marangoni stresses predominantly retard the flow and ‘rigidify’ the interface.
As can be observed in the results presented in figure 6(c), an increase in 
$\beta _s$ results in the strengthening of 
$\tau$. Figure 6(d) provides conclusive evidence of the surface rigidification as most of 
$u_{tz}$ is characterised by a negative value due to the decay of the dynamics. The variation of 
$u_{tz}$ decreases as 
$\beta _s$ increases, strengthening the Marangoni stresses arising from the surfactant-induced surface tension gradients. Finally, according to the foregoing results, we conclude that the interfacial dynamics for surfactant-laden cases in the ‘crown regime’ resembles its surfactant-free counterpart except for the retardation due to the surface rigidification brought about by the presence of surfactant-induced Marangoni stresses.
In figure 7 we plot the component of the velocity variation in the sheet-normal direction, 
$u_z$, within the sheet, in a frame of reference moving with the inner sheet at 
$t=5$. For the surfactant-free case, the liquid within the sheet follows a parabolic profile into the rim; this effectively corresponds to a no-slip condition with 
$u_z$ pinned to the interface, due to the absence of surfactants. In contrast, the retardation brought about by the presence of surfactants via Marangoni stresses results in shear flow for the surfactant-laden case. The higher surfactant concentration in the outer rim wall results in 
$\tau$ triggering the deceleration of the fluid entering into the bulbous end near the outer wall.
Figure 7. Profiles of the velocity component, 
$u_z$, in the sheet-normal direction inside the sheet for the surfactant-free and surfactant-laden 
$(\beta _s=0.5)$ cases at 
$t=5$. Here, 
$u_z$ has been calculated in an inner-sheet frame of reference along the cross-stream direction, 
$x$, for the axial location along the rim neck. Note that the axial distance has been normalised with the average sheet thickness, i.e. 
$x/h_s$ (thus 
$x/h_s=0$ and 
$x/h_s=1$ correspond to the inner and outer sheets, respectively). All parameters remain unchanged from figure 3.
Figure 8(a) presents the late stages of a surfactant-laden ligament before the end-pinching mechanism. In the absence of surfactants, the competition between viscosity and capillarity drives the formation of and dynamics of the bulbous edge on the ligament tip. We proceed to calculate the relative importance of viscous to surface tension forces via the local Ohnesorge number for the ligament (i.e. 
$Oh=\mu _l / \sqrt {\rho _l \sigma _s R}$, where 
$R$ stands for the average radius of the ligament) to explain the retardation phenomenon observed in figure 8(a). The resulting local 
$Oh \sim 0.07$ (e.g. low 
$Oh$-regime), and subsequently the local flow dynamics affecting the ligament (i.e. surfactant-driven retardation from end pinching) can be explained by the mechanism proposed by Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020), Kamat et al. (Reference Kamat, Wagoner, Castrejón-Pita, Castrejón-Pita, Anthony and Basaran2020) and Hoepffner & Paré (Reference Hoepffner and Paré2013). Figures 8(b,c) and 8(d,e) show 
$\varGamma$ and 
$\tau$, and the tangential velocity components 
$u_t$ for the framed panels in figure 8(a), respectively. Figure 8(b,d) demonstrates that, near the neck (marked with a blue diamond symbol) 
$\tau$ results in 
$u_{tx}>0$, preventing the formation of a second stagnation point required to trigger the capillary singularity. As shown by Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020), two stagnation points sandwiched between the neck boundaries are needed to induce end pinching, as can be observed in 8(d,e) (‘SP’ stands for stagnation points). Surface tangential stresses prevent capillary break-up for a longer time than in the surfactant-free case, leading to longer ligaments and an increase in surface area. This results in the dilution (increase) of surfactant (surface tension) and the eventual reduction in the strength of surfactant-induced Marangoni stresses thus initiating end pinching.
Figure 8. (a) Spatio-temporal evolution of a surfactant-laden ligament and its retardation from pinch-off driven by surfactant-induced Marangoni flow when 
$\beta _s=0.5$. (b,c) Represent 
$\varGamma$ and 
$\tau$ as a function of the arc length 
$s$ for the framed panels of (a), i.e. 
$t=13.49$ and 
$t=15.49$, respectively; (d,e) represent the tangential velocity 
$u_t$ as a function of the arc length 
$s$. The axial location in (a) has been normalised over the distance of the ligament on the axial direction; the abscissa direction has been normalised with respect to its value at 
$x=0$ for each single panel. The arc lengths have been normalised over the full extent of 
$s$. The diamond markers represent the location of the neck; ‘SP’ in (e) indicates the location of the stagnation points. All parameters remain unchanged from figure 3.
Figure 9 presents the temporal evolution of the neck radius for the surfactant-free and surfactant-laden cases (we cannot report the late stages of the neck due to the lack of snapshots up to its capillary singularity). As seen, at short times, the presence of surfactants at the interface leads to the reopening of the neck (
$\textrm {d}r/\textrm {d}t \sim 0$), and subsequently, to a short period in which surfactants induce retardation from end pinching. At a later time, as surfactants are evacuated from the neck, the dynamics eventually results in the interfacial singularity.
Figure 9. Temporal evolution of the neck radius for the surfactant-free and surfactant-laden cases (
$\beta _s=0.5$), highlighting the surfactant-driven retardation of the interfacial capillarity.
Next, we focus on drop formation and the splash phenomena: we have identified three different modes of droplet shedding. The first mode corresponds to droplet detaching from the ligament via end pinching. The mean lengths of the ligament 
$\langle l\rangle$ which undergoes mode-1 of ejection are: 
$\langle l\rangle =0.775$, 
$\langle l\rangle =0.864$, 
$\langle l\rangle =0.992$ and 
$\langle l\rangle =1.141$, which correspond to the surfactant-free and surfactant-laden cases with 
$\beta _s=(0.1,0.5,0.7)$, respectively. This agrees with a surfactant-driven retardation of the capillary singularity. In the second mode we observe the initial ejection of filaments that are rapidly recombined into a single filament, see figure 11(a–e). The third mode sees the formation of satellite (smaller) droplets during the stretching of the filament neck. These satellite droplets are not observed for large values of the elasticity parameter, in agreement with Kovalchuk et al. (Reference Kovalchuk, Jenkinson, Miller and Simmons2018). For the surfactant-free case, figure 10, ligaments merge due to the presence of local cusps arising from a non-uniform distribution of mass per unit length (in agreement with Gordillo, Lhuissier & Villermaux (Reference Gordillo, Lhuissier and Villermaux2014) and Wang & Bourouiba (Reference Wang and Bourouiba2018)), or by the rim parameter not being a multiple of the RP instability. In the vicinity of a cusp, the rim is not perpendicular to the flow entering the rim; instead, it takes an angle 
$\theta$, as shown in figure 10(a). This, oblique flow induces a local drift velocity along the rim. In a reference frame along the rim, the drift velocity can be calculated from the projection of the incoming velocity 
$u_{in}$ in the longitudinal direction of the rim as 
$u_{drift}=u_{in}(t) \sin (\theta )$ (see Wang & Bourouiba Reference Wang and Bourouiba2018 for more details). Under the conditions of figure 10, we obtain 
$u_{drift} \sim 0.106$ m s
$^{-1}$ and 
$\theta \sim 20^{\circ } \pm 3^{\circ }$. In contrast, from simulations, the velocity field yields 
$u_{in} \sim 0.27$ m s
$^{-1}$ which leads to 
$u_{drift} \sim 0.083$ m s
$^{-1}$, and a good agreement with the direct measurement. The collision of the drifting protrusions results in the formation of a corrugated ligament (see figure 10e), which succumbs to end pinching, resulting in drop formation (not shown). We confirm that we have made the same calculation for the other two flow-induced cuspids which give in similar results.
Figure 10. Spatio-temporal evolution of two adjacent rim protrusions which shift in the spanwise direction induced by local cusps, resulting in their collision; 
$\theta$ denotes the angle made by the left crest with the rim. The time difference between snapshots is 
$\Delta t=1.0$. The parameters remain unchanged from figure 2.
The presence of surfactant enhances the ligament-merging mode due to surfactant-driven Marangoni stresses that induce a local motion along the rim, and a radial shift of the ligament position along the rim. This is highlighted in figure 11, which presents the temporal evolution of two protrusions that eventually become ligaments, and merge prior to pinch-off.
Figure 11. Spatio-temporal evolution of two adjacent rim protrusions shown panels (a–e) which shift in the spanwise direction due to surfactant-induced Marangoni stresses, resulting in their collision and merging. Panel (f) shows a 3-D reconstruction of the 
$|\tau |$ profile with respect to the arc length, 
$s$, across a plane cutting the rim in (a) in half; the arrows represent the direction of 
$|\tau |$. Here, 
$\beta _s=0.5$, and all other parameters remain unchanged from figure 3.
Figure 11(a) demonstrates that the highest 
$\varGamma$ values are found at the protrusions as the fluid motion along the rim to the protrusions acts to accumulate surfactant locally. They are characterised by regions of a radially converging dynamics, and subsequently, increased (reduced) 
$\varGamma$ (
$\sigma$) locally. These results demonstrate that the effect of Marangoni stresses from adjacent protrusions is to bridge their gap, promoting the collision and merging of the ligaments. Figure 11(f) shows a 3-D reconstruction of the 
$|\tau |$ profile with respect to the arc length, 
$s$, across a plane cutting the rim seen in 11(a) in half; the arrows represent the direction of 
$|\tau |$, evidencing the enhancement of the filament merging. To the best of our knowledge, this surfactant-induced mechanism of droplet shedding has not been reported previously. Figure 11 also shows that surfactants are actively convected from the rim to the film. Consequently, the crown is thinner and survives longer in the presence of surfactants, in agreement with the experimental observations reported by Che & Matar (Reference Che and Matar2017).
The total number of ejected droplets for our surfactant-free case is found to be approximately 
$30$ up to times 
$t < 22.5$. In contrast, the number of secondary droplets (at 
$We,Re=800,1000$) as reported by Agbaglah & Deegan (Reference Agbaglah and Deegan2014), corresponds to 30–32 so our numerical predictions are in good agreement. This also verifies that a RP instability is the primary mechanism for wavelength selection.
Next, we turn our attention to the size distribution of the droplets ejected from the impact. Figure 12(a) shows the probability density function (PDF), at 
$t=20$, scaled by the initial droplet diameter. The sampling is conducted after the impact reaches 
$t=20$ and over three snapshots separated by time intervals of 
$\Delta t=0.2$. To get a smoother PDF, more snapshots should be included, but the dynamics is too rapid and due to the extreme computational cost of the current numerical simulation, this task was not undertaken. It is important to note that, as 
$\beta _s$ increases, a higher droplet size is favoured, as seen in figures 12(b–d) and 12(e–g), where the droplet volume, average surfactant concentration and average velocity for each droplet at 
$t=15$ and 
$t=20$, respectively, are shown. The average velocity and surfactant concentration are calculated as 
$\langle |\boldsymbol {u}|\rangle =\int _{\varOmega }|\boldsymbol {u}| \,\textrm {d} \varOmega / \varOmega$ and 
$\langle \varGamma \rangle =\int _{\varOmega }\varGamma \,\textrm {d}\varOmega / \varOmega$, respectively, where 
$\varOmega$ represents the surface area of the droplet. We observe that fewer droplets are produced with increasing 
$\beta _s$; fewer droplets are produced due to the rigidification brought about by surfactant-driven Marangoni stresses (see figure 12b,e). For the high values of 
$\beta _s$, we observe some droplets with large volume corresponding to the merging of ligaments. We also observe that small droplets are produced for the surfactant-free and the 
$\beta _s = 0.1$ surfactant-laden cases, (see figure 12e). These small droplets correspond to satellite droplets, which are characterised by 
$V, \langle |\boldsymbol {u}|\rangle$ and 
$\langle \varGamma \rangle \xrightarrow {} 0$, in good agreement with Kovalchuk & Simmons (Reference Kovalchuk and Simmons2018) and Che & Matar (Reference Che and Matar2017). We conclude that the addition of surfactants leads to a larger droplet size distribution for two reasons: the retardation of the dynamics, and the suppression of end pinching through the reopening of the neck driven by a surfactant-induced flow.
Figure 12. Metrics of splashing as a function of the elasticity parameter 
$\beta _s$. Panel (a) shows the probability density function (PDF) of the droplet size; droplet sizes are made dimensionless by the initial droplet diameter 
$D_0$. Panels (b–d) and (e–g) show the volume, average surfactant concentration and average velocity for every droplet at 
$t=15$ and 
$t=20$, respectively (
$d_i$ stands for the enumeration of droplets). All other parameters remain unchanged from figure 3.
Figure 12(c,f) indicates that the average 
$\langle \varGamma \rangle$ of the droplets is lower than the initial surfactant concentration of the pool (viz. 
$\varGamma =0.5$). This agrees well with the experimental studies of Blanchard & Syzdek (Reference Blanchard and Syzdek1972) and Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021c) on bursting bubbles. Figure 12(c) highlights that 
$\langle \varGamma \rangle$ increases with decreasing 
$\beta _s$ as the surfactant has been highly convected towards the rim. Lastly, figure 12(d,g), which shows 
$\langle |\boldsymbol {u}|\rangle$ at 
$t=15$ and 
$t=20$, respectively, demonstrates the effect of rigidification brought about by surfactant-induced Marangoni flow which corresponds to an overall reduction in 
$\langle |\boldsymbol {u}|\rangle$ for the surfactant-laden droplets.
Lastly, we plot in figure 13 the effect of surfactants on the interfacial area, and kinetic energy defined as 
$E_k=\rho \int _V \boldsymbol {u}^2/2 \,\textrm {d}V$, normalised by their initial values. Inspection of the surface area plot in figure 13(a) reveals a linear growth rate in interfacial area for all cases at short times; at longer times, a monotonic increase of surface area is observed with increasing 
$\beta _s$ due to surfactant-induced effects discussed above. In 13(b), we see that the presence of surfactants by increasing 
$\beta _s$ acts to decrease the overall value of 
$E_k$ (e.g. rigidification brought about surfactants).
Figure 13. Temporal evolution of the total interfacial area (a), and the kinetic energy (b), normalised by their initial values, for the surfactant-free and surfactant-laden cases. All parameters remain unchanged from figure 3.
4. Concluding remarks
We have presented, for the first time, a detailed analysis of the effect of insoluble surfactants on the dynamics of droplet impact on thin-film liquid layers using high-resolution 3-D simulations. This article focuses on the crown-splash regime, in which the transverse destabilisation of the rim leads to the formation of ligaments and droplets. Thus, we have focused on conditions of high Reynolds and Weber numbers. The nature of the hybrid front-tracking level-set numerical simulations enables the coupling, and analysis, of inertia, capillarity, interfacial diffusion and Marangoni stresses owing to the presence of surfactant-induced surface tension gradients. According to our results, the addition of surfactants does not significantly affect the selection of the wavelength of the transverse rim instability; as the predicted wavelength is consistent with the most unstable RP instability. The instability give rise to the formation of ligaments which eventually end pinch into droplets. In the presence of surfactants, we observe a delay in end pinching, driven by surfactant-induced Marangoni stresses, resulting in longer ligaments. We identify three modes of droplet ejection, and demonstrate that surfactant-induced Marangoni stresses result in the bridging of the spanwise surfactant concentration gradients between adjacent ligaments in the rim, eventually leading to their merging. The presence of surfactant-induced Marangoni stresses also leads to interfacial rigidification, which is observed through a reduction of the surface velocity and the kinetic energy, and a larger interface area.
This research is of importance to many applications as droplet impact onto finite-depth liquid layers is found in a wide range of industrial and daily-life applications, from raindrops impinging onto puddles to screen and inkjet printing. This research focuses on the crown-splashing regime but we acknowledge that the presence of surfactants will play a crucial role in other 
$We-Re$ regions, e.g. at the microdroplet-splash regime. Although we have only focused on insoluble surfactants, the presence of surfactant solubility will lead to additional richness and complexity, and we anticipate that the addition of soluble surfactants may result in detrimental effects, and they will be the subject of future work. Additionally, surfactant-induced rheological effects could be important and of great relevance to industry as most industrial fluids are formulated with dispersants (surfactants).
Finally, it is well known that the thickness of the film layer affects the interfacial dynamics. Additionally, most of the literature has focused on the limiting case of 
$t\ll D/U$ in which axisymmetric simulations suffice. The main novelty of this work is the study of the role of the surfactant at longer times after rim formation. At these times, surfactants play a major role in the late stages of droplet detachment, and, to the best of our knowledge, this manuscript is the first computational study on this phenomenon. Consequently, future studies could focus on the surfactant-induced effects of droplet impacting onto various liquid film thicknesses, and the ejecta sheet.
Acknowledgements
C.R.C.-A. and L.K. thank with gratitude Dr A. Batchvarov for the fruitful discussions. We acknowledge HPC facilities provided by the Research Computing Service (RCS) of Imperial College London for the computing time. For the purpose of Open Access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission.
Funding
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). O.K.M. acknowledges funding from PETRONAS and the Royal Academy of Engineering for a Research Chair in Multiphase Fluid Dynamics, and the EPSRC MEMPHIS (EP/K003976/1) and PREMIERE (EP/T000414/1) Programme Grants. 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 2022 A0122B06721. This research was funded in whole or in part by the UK EPSRC grant numbers specified above.
Declaration of interests
The authors report no conflict of interest.
Appendix A
Figure 14(a,b) shows additional numerical results against the experimental data of Che & Matar (Reference Che and Matar2017) for the crown rim evolution. Here, the effect of varying the Weber number, e.g. 
$We=(82,157,249)$, at a constant film depth 
$h=0.22$, 
$Re=(4312,6153,7514)$ and 
$Fr=(7.96,10.80,13.65)$, is seen for various droplet sizes. Due to the difference in droplet size, the dimensionless groups vary between each case: a small droplet with a diameter 
$D=0.00229$ m is described by 
$We=183$, 
$Re=5411$, 
$Fr=16.01$ and 
$h=0.13$; a medium droplet with diameter 
$D=0.00315$ m is characterised by 
$We=249$, 
$Re=6311$, 
$Fr=13.65$ and 
$h=0.095$; a large droplet with 
$D=0.00394$ m is represented by 
$We=315$, 
$Re=9341$, 
$Fr=12.20$ and 
$h=0.076$. We observe good agreement between our numerical predictions and the reported literature; we conclude our numerical approach faithfully captures the rich dynamics of droplet impact onto a thin-film layer of the same liquid.
Figure 14. Additional validation of the numerical framework (solid lines) against the experimental data (dashed lines) of Che & Matar (Reference Che and Matar2017). Panel (a) shows the effect of varying the Weber number and (b) shows the effect of varying droplet size on the crown diameter.
Appendix B
The code has been benchmarked against: (i) simulations of Notz & Basaran (Reference Notz and Basaran2004) on the recoiling of surfactant-free liquid threads, (ii) the RP instability model for liquid jets of Plateau (Reference Plateau1873), (iii) the scaling laws of pinching-off threads by Eggers (Reference Eggers1993) and Lister & Stone (Reference Lister and Stone1998) and (iv) the work of Blanchette & Bigioni (Reference Blanchette and Bigioni2006) on droplet impact at vanishing velocity (droplet coalescence on a flat surface). These validations are found elsewhere, i.e. Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021b, Reference Constante-Amores, Kahouadji, Batchvarov, Seungwon, Chergui, Juric and Matar2020, Reference Constante-Amores, Abadie, Kahouadji, Shin, Chergui, Juric, Castrejón-Pita and Matar2023).
Figure 15 displays the temporal evolution of the kinetic energy, 
$E_k$, and the conservation of liquid volume for the surfactant-free case (
$Re=1000$ and 
$We=800$) for different mesh resolutions. ‘Mesh-1’ and ‘Mesh-2’ refer to different levels of refinement characterised by 
$384^2\times 192$, and 
$768^2\times 384$, respectively. As demonstrated in figure 15, both levels of refinement can reproduce the dynamics with conservation of fluid volume around 
$0.1\,\%$. Based on these results, we conclude that ‘Mesh-2’ is sufficiently refined to ensure mesh-independent results while providing a good compromise with the computational cost. All the results presented in this paper correspond to ‘Mesh-2’.
Figure 15. Mesh study for the surfactant-free case when 
$Re=1000$ and 
$We=800$. The panels highlight the temporal evolution of kinetic energy 
$E_k$, and the relative variation of the liquid volume for two different meshes.
Open access
