Impact of droplets onto surfactant-laden thin liquid films

Droplets impacting on liquid surfaces are observed in a broad range of natural phenomena and industrial applications. Consequently, its study is driving a significant scientific interest; for comprehensive reviews please see e.g., Yarin (2006); Josserand & Thoroddsen (2016); Cheng et al. (2022). The complex topological features of droplet impact have attracted the fluid mechanics community since the first ground-breaking experiments from Worthington (1908) and Edgerton (1977) 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, 2002; Thoraval et al., 2012; Zhang et al., 2012; Agbaglah & Deegan, 2014; Josserand et al., 2016; Thoroddsen et al., 2012; Fudge et al., 2021).

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 (Josserand et al., 2016; Deegan et al., 2007; Che & Matar, 2017). At high impact speeds, the latter situation results in the formation of an upward Worthington jet (Gekle & Gordillo, 2010), whereas the former case gives rise to the well-known Milkdrop Coronet, which is the subject of this paper (Deegan et al., 2007; Josserand et al., 2016; Thoroddsen, 2002). 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 growth on the rim (see Agbaglah et al. (2013); Zhang et al. (2010)), to then develop non-linearities leading to breakup (through end-pinching, see Wang & Bourouiba (2018)), 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 (2020)). Previous studies have demonstrated the role of surfactants in capillary singularities (see for instance, Ananthakrishnan & Yeung (1994), Craster et al. (2002), Liao et al. (2004), Kamat et al. (2018) and Constante-Amores et al. (2020)); here, we study the role of surfactant-induced flows on the late stages of splashing. Che & Matar (2017) 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_{l}UD/\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 of 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 (2003); Josserand et al. (2016); Agbaglah & Deegan (2014). In contrast, the study of the crown-splash regime requires the use of full three-dimensional (3D) numerical simulations in order to accommodate the natural occurrence of possible symmetry-breaking events, such as, the growth of transversal instabilities in the rim, or rim breakup into droplets. However, 3D simulations still remain a numerical challenge (Gueyffier & Zaleski, 1998) and have not often been reported in the literature, (Liang et al., 2019; Reijers et al., 2019; Chen et al., 2020; Xavier et al., 2020).

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 3D 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 Section 4.

High resolution simulations were performed by solving the two-phase incompressible Navier-Stokes equations with surface tension in a three-dimensional Cartesian domain $\mathbf{x}=\left(x,y,z\right)$ (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., 2018). Here, and in what follows, all variables are made dimensionless (represented by tildes) using

The non-linear Langmuir equation is used to describe $\sigma$ in terms of $\Gamma$, this is

First, we discuss the early-time dynamics of the surfactant-free impact at $Re=1000$ and $We=800$. Figure 2a shows the three-dimensional representation of the interface at early times, i.e. prior transversal instabilities grow 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. (2007) and Josserand et al. (2016). 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 in agreement with the experimental observations by Zhang et al. (2010) and Deegan et al. (2007). Figure 2b below 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 between $t=3-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 Rayleigh-Plateau (RP) instability wavelength, given by $\lambda^{*}=2\pi/k=9.016a$, where $k$ and $a$ are the wavenumber of the fastest-growing instability mode, and the rim radius, respectively (Agbaglah et al., 2013; Zhang et al., 2010). In Figure 2c 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. (2010) and the theoretical and computational work from Agbaglah et al. (2013) and Constante-Amores et al. (2022). A closer inspection of the top panel of figure 2c shows that the interfacial corrugations are not always evenly distributed along the rim, as observed in experiments Thoraval et al. (2012). In particular, we observe large wavelengths, which are the result of interactions between adjacent instability modes, which is consistent with experimental observations.

Figures 3(a-d) show a three-dimensional 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 on figure 3b). The neck thins and stretches, due to a capillary-driven flow, to eventually resulting in the capillary pinch-off of a drop (the well-known ‘end-pinching’ mechanism proposed by Stone & Leal (1989)), see figure 3c. We note that the interfacial dynamics closely resembles experimental observations by Zhang et al. (2010), Josserand & Zaleski (2003), and Deegan et al. (2007), as well as the dynamics of retracting sheets reported by Wang & Bourouiba (2018) and Constante-Amores et al. (2022).

Next, we proceed to study the role of surfactants on the flow dynamics; as illustrated by figures 3(e-p) where we observe the droplet impact at various elasticity parameters $\beta_{s}$. As seen in these figures, 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., $\Gamma$ vanishes on the inner surface, see figure 3e. This observation, that the fluid source 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. (2016). A close inspection of the distribution of $\Gamma$ in the outer sheet shows that $\Gamma$-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 3e, in which larger values of $\Gamma$ are found in the rim than in its adjacent sheet.

The increase (decrease) of $\Gamma$ ($\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 et al. (1995), who reported a surfactant-driven dampening effect on the capillary waves.

In addition, we expect higher (lower) values of $\Gamma$ at the rim wave crests (instability waves) as they belong to a radially converging (diverging) region, which have the effect of locally increasing (reducing) $\Gamma$. The gradients in $\Gamma$ 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. (2022). By increasing $\beta_{s}$ the surfactant distribution along the interface is enhanced, as displayed in figure 3i-l and figure 3j-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 end-pinching mechanism to form droplets, as shown in figure 3g, k,o. These results offer an explanation to the experimental observations of Che & Matar (2017). 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 are mostly driven by the RP instability, although longer ligaments developed prior to break up in the presence of surfactants, as observed in figure 3p. 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. (2020); Constante-Amores et al. (2020, 2022)

Two-dimensional projections, in the $x-z$ plane ($y=4)$, of the interfacial shape, $\Gamma$, $\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, $\Gamma$ distributions trigger 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 $\Gamma$ distribution and increases the value of $\tau$ (see the increase of $\tau$ with $\beta_{s}$ in figure 5c). Furthermore, Figure 5c highlights the presence of large $\tau$ peaks at both sides of the rim due to the accumulation of $\Gamma$, 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 $\Gamma$ 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 5d. 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 impact region.

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 have been slowed down by the surface rigidification brought by the Marangoni stresses, as can be seen in Figure 6d (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 $\Gamma$ ($\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 6c, an increase in $\beta_{s}$ results in the strengthening of $\tau$. Figure 6d provides conclusive evidence of the surface rigidification as most of the $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’ resemble 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 near the outer wall.

Figure 8a 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 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}$, here $R$ stands for the average radius of the ligament) to explain the retardation phenomenon observed in figure 8a. 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. (2020), Kamat et al. (2020) and Hoepffner & Paré (2013). Figures 8b,c and 8d,e show $\Gamma$ and $\tau$, and the tangential velocity components $u_{t}$ for the framed panels in figure 8a, respectively. Figure 8b,d demonstrate 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. (2020), two stagnation points sandwiched between the neck boundaries are needed to induce end-pinching, as can be observed in 8d-e (e.g ‘SP’ stands for stagnation points). Surface tangential stresses prevent capillary breakup 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 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 ($dr/dt\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.

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 $<l>$ which undergoes mode-1 of ejection are: $<l>=0.775$, $<l>=0.864$, $<l>=0.992$, and $<l>=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) to (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. (2018). 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 et al. (2014) and Wang & Bourouiba (2018)), 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 10a. 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_{\rm{in}}$ in the longitudinal direction of the rim as $u_{\rm{drift}}=u_{\rm{in}}(t)\sin(\theta)$ (see Wang & Bourouiba (2018) for more details). Under the conditions of figure 10, we obtain $u_{\rm{drift}}\sim 0.106$ m/s and $\theta\sim 20^{\circ}\pm 3^{\circ}$. In contrast, from simulations, the velocity field yields $u_{\rm{in}}\sim 0.27$ m/s which leads to $u_{\rm{drift}}\sim 0.083$ m/s, 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 resulting in similar results.

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 pinchoff.

Figure 11a demonstrates that the highest $\Gamma$ values are found at the protrusions as the fluid motion along the rim to the protrusions acts to the accumulation of surfactant locally. They are characterised by regions of a radially converging dynamics, and subsequently, increased (reduced) $\Gamma$ ($\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 11f shows a three-dimensional reconstruction of the $|\tau|$ profile with respect to the arc-length, $s$, across a plane cutting the rim seen in 11a 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 surfactant is highly 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 (2017).

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 (2014), correspond 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 impact. Figure 12a 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 are 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 12b-d and 12e-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 $<|\bf{u}|>=\int_{\Omega}|\bf{u}|d\Omega/\Omega$ and $<\Gamma>=\int_{\Omega}\Gamma d\Omega/\Omega$, respectively, where $\Omega$ represents the surface area of the droplet. We observe that fewer droplets are produced with increasing $\beta_{s}$ is, fewer droplets are produced due to the rigidification brought 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,<|{\bf u}|>$ and $<\Gamma>\xrightarrow{}0$, in good agreement with Kovalchuk & Simmons (2018) and Che & Matar (2017). 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.

Figures 12(c,f) indicate that the average $<\Gamma>$ of the droplets is lower than the initial surfactant concentration of the pool (viz. $\Gamma=0.5$). This agrees well with the experimental studies of Blanchard & Syzdek (1972) and Constante-Amores et al. (2021c) on bursting bubbles. Figure 12c highlights that $<\Gamma>$ increases with decreasing $\beta_{s}$ as the surfactant has been highly convected towards the rim. Lastly, figures 12(d,g), which show $<|{\bf u}|>$ at $t=15$ and $t=20$, respectively, demonstrate the effect of rigidification brought about by surfactant-induced Marangoni flow which corresponds to an overall reduction in $<|{\bf u}|>$ 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}{\bf u}^{2}/2dV$, normalised by their initial values. Inspection of the surface area plot in figure 13a 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 13b, 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).

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 three-dimensional 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 transversal rim instability; as the predicted wavelength is consistent with the most unstable Rayleigh-Plateau 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 microdoplet-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 and 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 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.

Declaration of Interests. The authors report no conflict of interest.

CRC-A and LK thank with gratitude Dr A. Batchvarov for the fruitful discussions. AAC-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. We 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 2022 A0122B06721. This research was funded in whole or in part by the UK EPSRC grant numbers specified above. 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.

Figure 14 a,b show additional numerical results against the experimental data of Che & Matar (2017) 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.

The code has been benchmarked against: i) simulations of Notz & Basaran (2004) on the recoiling of surfactant-free liquid threads, ii) the Rayleigh-Plateau instability model for liquid jets of Plateau (1873), iii) the scaling laws of pinching off threads by Eggers (1993) and Lister & Stone (1998), and iv) the work of Blanchette & Bigioni (2006) on droplet impact at vanishing velocity (droplet coalescence on a flat surface). These validations are found elsewhere, i.e. Constante-Amores et al. (2021b, 2020, 2023).

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 a 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 a ‘Mesh-2’ mesh.