domain: spi nlin

text:
C R Constante-Amores

L Kahouadji

S Shin

J Chergui

D Juric

J R Castrejón-Pita

O K Matar

A A Castrejón-Pita

Impact of droplets onto surfactant-laden thin liquid films

. We focus on the 'crown-splash regime', and we observe that the crown dynamics evolves through various stages: from the 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, while also promoting the spanwise merging between ligaments. Finally, we show 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., [START_REF] Yarin | Drop impact dynamics: Splashing, spreading, receding, bouncing[END_REF]; [START_REF] Josserand | Drop impact on a solid surface[END_REF]; [START_REF] Cheng | Drop impact dynamics: Impact force and stress distributions[END_REF]. The complex topological features of droplet impact have attracted the fluid mechanics community since the first ground-breaking experiments from [START_REF] Worthington | A study of splashes[END_REF] and [START_REF] Edgerton | Stopping Time: The Photographs of Harold Edgerton[END_REF] 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 [START_REF] Thoroddsen | The ejecta sheet generated by the impact of a drop[END_REF][START_REF] Thoraval | von kármán vortex street within an impacting drop[END_REF][START_REF] Zhang | Evolution of the ejecta sheet from the impact of a drop with a deep pool[END_REF][START_REF] Agbaglah | Growth and instability of the liquid rim in the crown splash regime[END_REF]Josserand et al. 2016;[START_REF] Thoroddsen | Micro-bubble morphologies following drop impacts onto a pool surface[END_REF][START_REF] Fudge | Dipping into a new pool: The interface dynamics of drops impacting onto a different liquid[END_REF].

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;[START_REF] Deegan | Complexities of splashing[END_REF][START_REF] Che | Impact of droplets on liquid films in the presence of surfactant[END_REF]. At high impact speeds, the latter situation results in the formation of an upward Worthington jet [START_REF] Gekle | Generation and breakup of worthington jets after cavity collapse. part 1. jet formation[END_REF], whereas the former case gives rise to the well-known Milkdrop Coronet, which is the subject of this paper [START_REF] Deegan | Complexities of splashing[END_REF]Josserand et al. 2016;[START_REF] Thoroddsen | The ejecta sheet generated by the impact of a drop[END_REF]. 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 [START_REF] Agbaglah | Longitudinal instability of a liquid rim[END_REF]; [START_REF] Zhang | Wavelength selection in the crown splash[END_REF]), to then develop nonlinearities leading to breakup (through end-pinching, see [START_REF] Wang | Unsteady sheet fragmentation: droplet sizes and speeds[END_REF]), and the formation of a cascade of droplet sizes.

In fact, [START_REF] Deegan | Complexities of splashing[END_REF] presented a phenomenological picture of the interfacial dynamics in terms of a phase diagram in the W e -Re space (here, W e and Re stand for the Weber and Reynolds numbers, which represent the ratio of inertia to surface tension forces and inertia to viscous forces, respectively). In the 'Peregrine-sheet regime', i.e. low Re and moderate W e, and high Re and low W e, the dynamics are characterised by the generation of an axisymmetric ejecta sheet and rim. In the 'crown-splash regime', i.e. for low Re and high W e, the rim undergoes longitudinal instabilities to trigger the formation of droplets. In the 'microdoplet-splash regime', i.e. for high Re and moderate W e, fast and diminutive droplets are formed from the ejecta sheet immediately upon impact. Microdoplet-splash and crown-splash regimes are expected for high Re and high W e.

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 [START_REF] Manikantan | Surfactant dynamics: hidden variables controlling fluid flows[END_REF]). Previous studies have demonstrated the role of surfactants in capillary singularities (see for instance, [START_REF] Ananthakrishnan | Nonlinear interaction of a vortex pair with clean and surfactant-covered free surfaces[END_REF], [START_REF] Craster | Pinchoff and satellite formation in surfactant covered viscous threads[END_REF], [START_REF] Liao | Effects of soluble surfactants on the deformation and breakup of stretching liquid bridges[END_REF], [START_REF] Kamat | Role of marangoni stress during breakup of surfactant-covered liquid threads: Reduced rates of thinning and microthread cascades[END_REF] and Constante-Amores et al. (2020a)); here, we study the role of surfactant-induced flows on the late stages of splashing. [START_REF] Che | Impact of droplets on liquid films in the presence of surfactant[END_REF] 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. In this past work, the authors showed that the presence of surfactants affect 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, D, and U stand for the time, the droplet diameter, and its characteristic velocity, respectively) -see for example [START_REF] Josserand | Droplet splashing on a thin liquid film[END_REF]; Josserand et al. (2016); [START_REF] Agbaglah | Growth and instability of the liquid rim in the crown splash regime[END_REF]. 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 [START_REF] Gueyffier | Formation de digitations lors de l'impact d'une goutte sur un film liquide[END_REF] and have not often been reported in the literature, [START_REF] Liang | Non-simultaneous impact of multiple droplets on liquid film[END_REF][START_REF] Reijers | Oblique droplet impact onto a deep liquid pool[END_REF][START_REF] Chen | Oblique drop impact on thin film: Splashing dynamics at moderate impingement angles[END_REF][START_REF] Xavier | Toward direct numerical simulation of high speed droplet[END_REF].

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 surfactant transport is resolved at the interface [START_REF] Shin | A hybrid interface tracking -level set technique for multiphase flow with soluble surfactant[END_REF]. Here, and in what follows, all variables are made dimensionless (represented by tildes) using

x = x D o , t = t t r , ũ = u U , p = p ρU 2 , σ = σ σ s , Γ = Γ Γ ∞ , (2.1)
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

∇ • ũ = 0, (2.2) ρ( ∂ ũ ∂ t + ũ • ∇ũ) + ∇p = - i z F r 2 + 1 Re ∇ • μ(∇ũ + ∇ũ T ) + + 1 We à ( t) (σκn + ∇ s σ)δ x -xf d Ã, (2.3) ∂ Γ ∂ t + ∇ s • ( Γ ũt ) = 1 P e s ∇ 2 s Γ , (2.4)
where the density and viscosity are given by ρ = ρ g /ρ + (1 -ρ g /ρ) H x, t and μ = µ g /µ + (1 -µ g /µ) 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' designates the gas phase, 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. 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 [START_REF] Deegan | Complexities of splashing[END_REF]). The parameter β s is the surfactant elasticity number that is a measure of the sensitivity of the surface tension, σ, to the interface surfactant concentration, Γ . The non-linear Langmuir equation is used to describe σ in terms of Γ , this is σ = 1 + β s ln (1 -Γ ).

Re = ρ l U D o µ l , W e = ρ l U 2 D o σ s , F r = U √ gD o , P e s = U D o D s , β s = ℜT Γ ∞ σ s , (2.5) (a) (b)
(2.6)

The Marangoni stress, τ , is expressed as a function of Γ as

1 We ∇ s σ • t ≡ τ We = - Ma (1 -Γ ) ∇ s Γ • t, (2.7)
where M a = β s /W e = ReT Γ ∞ /ρU 2 D o is the Marangoni parameter and t is the unit tangent to the interface. 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. 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 shows geometry and the initial and boundary conditions of the problem which closely follows previous work by Josserand et al. (2016); Agbaglah &[START_REF] Agbaglah | Growth and instability of the liquid rim in the crown splash regime[END_REF][START_REF] Fudge | Dipping into a new pool: The interface dynamics of drops impacting onto a different liquid[END_REF]. 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 is 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. (2020a)).

Figure 1b illustrates the ability of the numerical technique to reproduce the temporal evolution of the crown experimentally obtained and reported by [START_REF] Che | Impact of droplets on liquid films in the presence of surfactant[END_REF]. 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. The numerical technique has also been validated for dropinterface coalescence, which can be considered as a limiting case of drop impact onto a pool, i.e., U = 0 (see [START_REF] Constante-Amores | Direct numerical simulations of transient turbulent jets: vortexinterface interactions[END_REF] for more details). The nonlinear interfacial dynamics have been validated for the capillary breakup of liquid threads (see Constante-Amores et al. (2020a, 2021) for more details). The numerical validation of the surfactant equations has been previously presented in [START_REF] Shin | A hybrid interface tracking -level set technique for multiphase flow with soluble surfactant[END_REF].

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 (772) 3 . 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. 2020a,b;Constante-Amores et al. 2020a, 2021, 2022).

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 [START_REF] Deegan | Complexities of splashing[END_REF] 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 and W e = 800 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. Γ ∞ ∼ O(10 -6 ) mol m -2 for NBD-PC (1-palmitoyl-2-12-[(7nitro-2-1,3-benzoxadiazol-4-yl)amino]dodecanoyl-sn-glycero-3 -phosphocholine); thus, we have explored the range of 0.1 < β s < 0.9 which corresponds to O(10 -7 ) <CMC< O(10 -6 ) mol m -2 , for typical values of σ s . We have set P e s = 10 2 following Batchvarov et al. (2020b) and Constante-Amores et al. (2020a) who showed that the interfacial dynamics are weakly-dependent on P e s beyond this value.

For a water-droplet of a typical size of D o ∼ O(10 -3 )m, the impact time scale T imp = t r = D o /U ∼ O(10 -4 -10 -3 )s; whereas the Marangoni time-scale T τ can be estimated by a balance between the Marangoni and the viscous stresses, resulting in T τ ∼ µD 2 /(h∆σ) ∼ O(10 -4 -10 -3 )s (see [START_REF] Che | Impact of droplets on liquid films in the presence of surfactant[END_REF] 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 in the x -z plane (y = 4) 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 [START_REF] Deegan | Complexities of splashing[END_REF] 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 shows 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 [START_REF] Zhang | Wavelength selection in the crown splash[END_REF] and [START_REF] Deegan | Complexities of splashing[END_REF]. 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 growth rate of the fastest-growing instability mode, and the rim radius, respectively. Figure 2b shows 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 [START_REF] Zhang | Wavelength selection in the crown splash[END_REF] and the theoretical and computational work from [START_REF] Agbaglah | Longitudinal instability of a liquid rim[END_REF] and Constante-Amores et al. (2022). A closer inspection of the top panel of figure 2b shows that the interfacial corrugations are not always evenly distributed along the rim. 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 [START_REF] Stone | Relaxation and breakup of an initially extended drop in an otherwise quiescent fluid[END_REF]), see figure 3c. We note that the interfacial dynamics closely resembles experimental observations by [START_REF] Zhang | Wavelength selection in the crown splash[END_REF], [START_REF] Josserand | Droplet splashing on a thin liquid film[END_REF] and [START_REF] Deegan | Complexities of splashing[END_REF], as well as the dynamics of retracting sheets reported by [START_REF] Wang | Unsteady sheet fragmentation: droplet sizes and speeds[END_REF] 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 idea 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, 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 [START_REF] Asaki | Effect of an insoluble surfactant on capillary oscillations of bubbles in water: Observation of a maximum in the damping[END_REF], 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 higher-σ radially diverging to the lower-σ converging regions. This slows down the interfacial dynamics, and the formation of instabilities Surfactant-induced Marangoni stresses only delay the development of the ligament from the rim; these eventually end-pinch to form droplets, as shown in figure 3f,g. These predictions agree well with the experimental observations of [START_REF] Che | Impact of droplets on liquid films in the presence of surfactant[END_REF]. Figure 4 shows the selected instability wavelength in terms of β s and a 0 . Indeed, a close inspection indicates that the surfactant does not 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 prior to break up developed in the presence of surfactants, as observed in figure 3e. This indicates that Marangoni stresses promote the escape from end-pinching, as previously demonstrated by Constante-Amores et al. (2020a, 2022) Two-dimensional projections, in the x -z plane (y = 4D), 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 ejected 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). Figure 5c shows the presence of large τ peaks at both sides of the rim due to the accumulation of Γ , and suddenly drops at the base of the rim, suggesting a fluid transport from the rim to the sheet. In addition, a lager τ 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, shows 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. 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 shown 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.

t = 5 t = 10 t = 15 t = 20 Surf actant -f ree (a) (b) (c) (d) βs = 0.1 (e) (f) (g) (h) βs = 0.5 (i) (j) (k) (l) βs = 0.7 (m) (n) (o) (p)
Figure 6c shows that the increase of β 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.

Figure 7 shows 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 shows the late stages of a surfactant-laden ligament before the end-pinching mechanism. We proceed to calculate the local Ohnesorge number for the ligament (i.e. Oh = µ l / √ ρ l σ s R, here R stands for the average radius of the ligament) to explain the escape phenomenon observed in figure 8a. The resulting local Oh ∼ 0.07 (e.g., low Oh-regime), and subsequently the local flow dynamics affecting the ligament (i.e., surfactant-driven escape from end-pinching) can be explained by the mechanism proposed by Constante-Amores et al. (2020a) and [START_REF] Kamat | Surfactantdriven escape from endpinching during contraction of nearly inviscid filaments[END_REF]. Figures 8b,c and8d,e show Γ and τ , and the tangential velocity components u t for the framed panels in figure 8a, respectively. Figure 8b,d show that, near the neck (marked with a blue diamond symbol) τ 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. (2020a), two stagnation points sandwiched between the neck boundaries are needed to induce endpinching, 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 surfactantfree 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.

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. In the second mode we observe the initial ejection of filaments that are rapidly recombined into a single filament, see figure 10(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 9, ligaments merge due to the presence of local cusps arising from a nonuniform distribution of mass per unit length (in agreement with [START_REF] Gordillo | On the cusps bordering liquid sheets[END_REF] and [START_REF] Wang | Unsteady sheet fragmentation: droplet sizes and speeds[END_REF]), 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 9a. 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 [START_REF] Wang | Unsteady sheet fragmentation: droplet sizes and speeds[END_REF] for more details). Under the conditions of figure 9, 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 9e), which succumbs to end-pinching resulting in drop formation (not shown). 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. Figure 10 shows the temporal evolution of two protrusions that eventually become ligaments, and merge prior to pinchoff. Figure 10a shows 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. 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 10f shows a threedimensional reconstruction of the |τ | profile with respect to the arc-length, s, across a plane cutting the rim seen in 10a 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 10 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 [START_REF] Che | Impact of droplets on liquid films in the presence of surfactant[END_REF].

Next, we turn our attention to the size distribution of the droplets ejected from impact. Figure 11a shows the probability density function (P DF ), at t = 20, scaled by the initial droplet diameter. It is important to note that as β s increases, a higher droplet size is favoured. Figures 11b-d and 11e-g show the droplet volume, average surfactant concentration, and average velocity for each droplet at t = 15 and t = 20, respectively. 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 11b,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 11e). These small droplets correspond to satellite droplets, which are characterised by V, < |u| > and < Γ >-→ 0, in good agreement with Kovalchuk & Simmons (2018) and [START_REF] Che | Impact of droplets on liquid films in the presence of surfactant[END_REF]. 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 11(c,f) show 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 [START_REF] Blanchard | Concentration of bacteria in jet drops from bursting bubbles[END_REF] and Constante-Amores et al. (2020b) on bursting bubbles. Figure 11c highlights that < Γ > increases with decreasing β s as the surfactant has been highly convected towards the rim. Lastly, figures 11(d,g), which show < |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 < |u| > for the surfactant-laden droplets.

Lastly, we plot in figure 12 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 12a 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 12b, 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 spanwise 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 endpinching, driven by surfactant-induced Marangoni stresses, resulting in longer ligaments. We identify three modes of droplet ejection, and show that surfactant-induced Marangoni stresses result in the bridging of the spanwise surfactant concentration gradients between consecutive 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). , and (772) 3 , respectively. As shown in figure 13, 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.

Figure 1 :

 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 772 3 ; (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 = 186).

Figure 2 :

 2 Figure 2: Wavelength selection in the crown splash regime. Panel (a) shows 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) 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 shown, highlighting the predicted wavelength between crests.

Figure 3 :

 3 Figure 3: Effect of β s on the drop impact dynamics for insoluble surfactants. Spatiotemporal evolution of the three-dimensional interface shape for surfactant-free cases, (a)-(d), and surfactant-laden cases for β 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, W e = 800 and h = 0.2. For the surfactant-laden cases, P e s = 100 and Γ = Γ ∞ /2, the colour indicates the value of Γ , and the legend is shown in panel (e).

Figure 4 :

 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 :

 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 shown. 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. All parameters remain unchanged from figure 3.

Figure 6 :

 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. All parameters remain unchanged from figure 3.

Figure 7 :

 7 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 indicated in the right panel (the axial distance has been normalised as a function of the sheet thickness for each case). All parameters remain unchanged from figure 3.

Figure 8 :

 8 Figure 8: (a) Spatio-temporal evolution of a surfactant-laden ligament and its escape 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.49 and 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 z-direction has been normalised with respect to the value at x = 0 for each 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 :

 9 Figure 9: 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 10 :

 10 Figure 10: 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 11 :

 11 Figure 11: 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. 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. All other parameters remain unchanged from figure 3.

Figure 12 :

 12 Figure 12: 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 13 :

 13 Figure 13: 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

† Email address for correspondence: crc15@ic.ac.uk

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.

Appendix

The code has been benchmarked against: i) simulations of [START_REF] Notz | Dynamics and breakup of a contracting liquid filament[END_REF] on the recoiling of surfactant-free liquid threads, ii) the Rayleigh-Plateau instability model for liquid jets of [START_REF] Plateau | Experimental and Theoretical Statics of Liquids Subject to Molecular Forces Only[END_REF], iii) the scaling laws of pinching off threads by [START_REF] Eggers | Universal pinching of 3d axisymmetric free-surface flow[END_REF] and [START_REF] Lister | Capillary breakup of a viscous thread surrounded by another viscous fluid[END_REF], and iv) the work of [START_REF] Blanchette | Partial coalescence of drops at liquid interfaces[END_REF] on droplet impact at vanishing velocity (droplet coalescence on a flat surface). These validations are found elsewhere, i.e. Constante-Amores et al. (2021, 2020a).

Figure 13 shows the temporal evolution of the kinetic energy, E k , and the conservation of liquid volume for the surfactant-free case (Re = 1000 and W e = 800) for different mesh