Impact of droplets onto surfactant-laden thin liquid films

Abstract We study the effect of insoluble surfactants on the impact of surfactant-free droplets onto surfactant-laden thin liquid films via a fully three-dimensional direct numerical simulation approach that employs a hybrid interface-tracking/level-set method, and by taking into account surfactant-induced Marangoni stresses due to gradients in interfacial surfactant concentration. Our numerical predictions for the temporal evolution of the surfactant-free crown are validated against the experimental work by Che & Matar (Langmuir, vol. 33, 2017, pp. 12140–12148). We focus on the ‘crown-splash regime’, and we observe that the crown dynamics evolves through various stages: from the growth of linear modes (through a Rayleigh–Plateau instability) to the development of nonlinearities leading to primary and secondary breakup events (through droplet shedding modulated by an end-pinching mechanism). We show that the addition of surfactants does not affect the wave selection via the Rayleigh–Plateau instability. However, the presence of surfactants plays a key role in the late stages of the dynamics as soon as the ligaments are driven out from the rim. Surfactant-induced Marangoni stresses delay the end-pinching mechanisms to result in longer ligaments prior to their capillary singularity. Our results indicate that Marangoni stresses bridge the gap between adjacent protrusions promoting the adjacent protrusions' collision and the merging of ligaments. Finally, we demonstrate that the addition of surfactants leads to surface rigidification and consequently to the retardation of the flow dynamics.


Introduction
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 nonlinearities 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 < W e < 532 and 3000 < Re < 8654, respectively (Re = ρ l U D/µ l and W e = ρ l U 2 D/σ, where D, ρ l , µ l σ, 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 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 surfactantladen 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.

Problem formulation and numerical method
High resolution simulations were performed by solving the two-phase incompressible Navier-Stokes equations with surface tension in a three-dimensional Cartesian domain 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. 2018).Here, and in what follows, all variables are made dimensionless (represented by tildes) using where, t, u, and p stand for time, velocity, and pressure, respectively.The physical parameters correspond to the liquid density ρ l , liquid viscosity, µ l , surface tension, σ, surfactant-free surface tension, σ s , and 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, Γ , is scaled by the saturation interfacial concentration, Γ ∞ .As a result of this scaling, the dimensionless equations read where the density and viscosity are given by ρ = ρ g /ρ l + (1 − ρ g /ρ l ) H x, t and μ = µ g /µ l + (1 − µ g /µ l ) H x, t wherein H x, t represents a Heaviside function, which is zero in the gas phase and unity in the liquid phase, while the subscript 'g' 'l' designate the gas, and liquid phases, respectively, and ũt = (ũ s • t) t is the tangential velocity at the interface in which ũs represents the interfacial velocity, and κ is the curvature.The interfacial gradient is given by ∇ s = (I − nn) • ∇ wherein I is the identity tensor and n is the outward-pointing unit normal.In addition, δ is a Dirac delta function, equal to unity at the interface and zero otherwise, and Ã( t) is the time-dependent interface area.
As justified by the final paragraph of Stone (1990), in our frame of reference u • n = 0, which gives rise to equation 2.4.The dimensionless groups that appear in the governing equations are defined as where Re, W e, F r, and P e s stand for the Reynolds, Weber, Froude, and (interfacial) Peclet numbers.Gravity is negligible during the impact as indicated by the Froude number value, i.e.F r ∼ O(10 2 ) (similar assumptions were made by Deegan et al. (2007)).The parameter β s is the surfactant elasticity number that is a measure of the sensitivity of the surface tension, σ, to the interface surfactant concentration, Γ .Here, is the ideal gas constant value ( = 8.314 J K −1 mol −1 ), T denotes temperature, and D s stands for the interfacial diffusion coefficient.The non-linear Langmuir equation is used to describe σ in terms of Γ , this is (2.6) The Marangoni stress, τ , is expressed as a function of Γ as where M a = β s /W e = ReT Γ ∞ /ρ l U 2 D o is the Marangoni parameter and 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.

Numerical Setup, validation and parameters
Figure 1a highlights the geometry and the initial and boundary conditions of the problem, which closely follows previous work by Josserand et al. (2016); Agbaglah &Deegan (2014), andFudge et al. (2021).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 × 8D 0 × 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. (2021)).
Figure 1b illustrates the ability of the numerical technique to reproduce the temporal evolution of the crown experimentally obtained and reported by Che & Matar (2017).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 other two film thicknesses.The discrepancy could be attributed to conditions found in experimental conditions but 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. (2021a) for more details).The nonlinear interfacial dynamics have been validated for the capillary breakup of liquid threads (see Constante-Amores et al. (2020, 2021b) for more details).The numerical validation of the surfactant equations has been previously presented in Shin et al. (2018).
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 × 384 (i.e., D 0 /∆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-tensiondriven phenomena using the same computational method can be found elsewhere, e.g.(Batchvarov et al. 2021(Batchvarov et al. , 2020;;Constante-Amores et al. 2020, 2021b, 2022, 2023).
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 ρ g /ρ l = 1.2 × 10 −3 and µ g /µ l = 0.018, respectively.The targeted flow conditions for the surfactant-free base case are based on Deegan et al. (2007) who built a phase diagram in the W e − Re space for h = 0.2.Our study focuses on the 'crown-splash' regime (i.e., W e > 500 and Re < 1100), as we aim to stay away from the 'microdropletsplash' regime (i.e., W e > 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, W e = 800 and h = 0.2 as the targeted surfactant-free case.
For a water-droplet of a typical size of D o ∼ O(10 −3 )m, the impact time scale )s; whereas the Marangoni time-scale T τ can be estimated by a balance between the Marangoni and the viscous stresses, resulting in T τ ∼ µ l D 2 /(h∆σ) ∼ O(10 −4 − 10 −3 )s (see Che & Matar (2017) 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.

Results
First, we discuss the early-time dynamics of the surfactant-free impact at Re = 1000 and W e = 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 (λ) 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 λ * = 2π/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, λ/λ * = 1.0 ± 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 β s .As seen in these figures, regardless of the value of β 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., Γ 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, A three-dimensional representation of the crown for the dimensionless time t = 7.5 is also presented, highlighting the predicted wavelength between crests.respectively, is in good agreement with the work by Josserand et al. (2016).A close inspection of the distribution of Γ in the outer sheet shows that Γ -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 Γ are found in the rim than in its adjacent sheet.
The increase (decrease) of Γ (σ) 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 Γ at the rim wave crests (instability waves) as they belong to a radially converging (diverging) region, which have the effect of locally increasing (reducing) Γ .The gradients in Γ result in surfactant- induced Marangoni stress flows from the lower-σ radially-converging to the higher-σ 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 β s the surfactant distribution along the interface is enhanced, as displayed in figure 3i-l and figure 3j-m for β s = 0.5 and β 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 β 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 λ is consistent with the most unstable RP instability at a ratio of λ/λ * = 1.0 ± 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 endpinching, 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, Γ , τ , and the radial component of the interfacial velocity u tz in terms of β s and a fixed value of P e 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 β s = 0.1 surfactantladen cases are practically indistinguishable.However, as the surfactant concentration increases, Γ distributions trigger surfactant-driven dynamics from the sheet-base to its edge (see figure 5b).As β 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 β s also enhances the Γ distribution and increases the value of τ (see the increase of τ with β s in figure 5c).Furthermore, Figure 5c highlights the presence of large τ peaks at both sides of the rim due to the accumulation of Γ , and sudden drops at the base of the rim, suggesting a fluid transport from the rim to the sheet.In addition, a larger τ maximum is observed from the side of the inner sheet as the inner sheet is devoid of Γ leading to an uneven flow motion inside the sheet.The variation of the streamwise interfacial tangential velocity u tz , 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 tz is observed around the rim with u tz tending to zero at a sufficiently large distance away from the impact region.Figure 7: Profiles of the velocity component, u z , in the sheet-normal direction inside the sheet for the surfactant-free and surfactant-laden (β 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 sheet, respectively).All parameters remain unchanged from figure 3.
The influence of β 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 tz is smaller with increasing β s ).As β s decreases (σ increases), the surfactant is highly convected as τ tends to zero, lowering the interfacial tension and consequently increasing surface deformation.As β s increases, the gradient in Γ (τ ) 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 β s results in the strengthening of τ .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 β 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 τ 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 endpinching 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 = µ l / √ ρ l σ s R, here R stands for the   (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 surfactantinduced 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 ∼ 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 β 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 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, θ, 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 in in the longitudinal direction of the rim as u drift = u in (t) sin(θ) (see Wang & Bourouiba (2018) for more details).Under the conditions of figure 10, we obtain u drift ∼ 0.106 m/s and θ ∼ 20 • ± 3 • .In contrast, from simulations, the velocity field yields u in ∼ 0.27 m/s which leads to u drift ∼ 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 surfactantdriven 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 Γ 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) Γ (σ) 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 |τ | 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 |τ | 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 W e, 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 (P DF ), 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 ∆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 β 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 < |u| >= Ω |u|dΩ/Ω and < Γ >= Ω Γ dΩ/Ω, respectively, where Ω represents the surface area of the droplet.We observe that fewer droplets are produced with increasing β 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 β 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 β s = 0.1 surfactant-laden cases, (see figure 12e).These small droplets correspond to satellite droplets, which are characterised by V, < |u| > and < Γ >− → 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 < Γ > of the droplets is lower than the initial surfactant concentration of the pool (viz.Γ = 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 < Γ > increases with decreasing β s as the surfactant Lastly, we plot in figure 13 the effect of surfactants on the interfacial area, and kinetic energy defined as E k = ρ V 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 β s due to surfactant-induced effects discussed above.In 13b, we see that the presence of surfactants by increasing β s acts to decrease the overall value of E k (e.g., rigidification brought about surfactants).

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 highresolution 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 finitedepth 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 W e − 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 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.

Appendix A
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.W e = (82, 157, 249), at a constant film depth h = 0.22, Re = (4312, 6153, 7514) and F r = (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 W e = 183, Re = 5411, F r = 16.01, and h = 0.13; a medium droplet with diameter D = 0.00315 m is characterised by W e = 249, Re = 6311, F r = 13.65, and h = 0.095; a large droplet with D = 0.00394 m is represented by W e = 315, Re = 9341, F r = 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 1 :
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 × 8D o × 4D o ) in a three-dimensional Cartesian domain, x = (x, y, z), with a resolution of 768 2 × 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 (2017): effect of varying the dimensionless film depth, h = H/D o , while keeping the other parameters constant (Re = 7514, W e = 249, and F r = 13.65).

Figure 2 :
Figure 2: Wavelength selection in the crown splash regime.Panel (a) presents the early interfacial dynamics through a three-dimensional representation of the predicted interface location for Re = 1000 and W e = 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 λ normalized with the theoretical RP-instability, λ * as a function of the local rim radius, a, at early times times of the simulation for t < 5.A three-dimensional representation of the crown for the dimensionless time t = 7.5 is also presented, highlighting the predicted wavelength between crests.

Figure 4 :
Figure 4: Effect of the elasticity parameter, β s in the selection of wavelength of the undulations normalized with the theoretical RP instability, λ * , as a function of the local rim radius, a.

Figure 5 :
Figure 5: Effect of the elasticity parameter, β s , on the flow dynamics at t = 5.Twodimensional projections of the interface, Γ , τ , 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.

Figure 6 :
Figure 6: Effect of the elasticity parameter, β s , on the flow dynamics at t = 20.Twodimensional projections of the interface, Γ , τ , 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.

Figure 8 :
Figure 8: (a) Spatio-temporal evolution of a surfactant-laden ligament and its retardation from pinchoff driven by surfactant-induced Marangoni flow when β s = 0.5.(b) and (c) represent Γ and τ as a function of the arc length s for the framed panels of (a), i.e., t = 13.49and t = 15.49,respectively; (d) and (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 :
Figure 9: Temporal evolution of the neck radius for the surfactant-free and surfactantladen case (β s = 0.5), highlighting the surfactant-driven retardation of the interfacial capillarity.
(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 withKovalchuk et al. (2018).For the surfactant-free case, figure10, ligaments

Figure 10 :
Figure 10: Spatio-temporal evolution of two adjacent rim protrusions which shift in the spanwise direction induced by local cusps, resulting in their collision; θ denotes the angle made by the left crest with the rim.The time difference between snapshots is ∆t = 1.0.The parameters remain unchanged from figure 2.

Figure 11 :
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 three-dimensional reconstruction of the |τ | profile with respect to the arc-length, s, across a plane cutting the rim in (a) in half; the arrows represent the direction of |τ |.Here, β s = 0.5, and all other parameters remain unchanged from figure 3.

Figure 12 :
Figure 12: Metrics of splashing as a function of the elasticity parameter β s .Panel (a) shows the probability density function (P DF ) 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 13 :
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.

Figure 14 :
Figure 14: Additional validation of the numerical framework (solid lines) against the experimental data (dashed lines) of Che & Matar (2017).(a) sees the effect of varying the Weber number and (b) shows the effect of varying droplet size on the crown diameter.

Figure 15 :
Figure 15: Mesh study for the surfactant-free case when Re = 1000 and W e = 800.The panels highlight the temporal evolution of kinetic energy E k , and the relative variation of the liquid volume for two different meshes