Fragmentation of colliding liquid rims

Abstract We present direct numerical simulations of the splashing process between two cylindrical liquid rims. This belongs to a class of impact and collision problems with a wide range of applications in science and engineering, and motivated here by splashing of breaking ocean waves. Interfacial perturbations with a truncated white noise frequency profile are introduced to the rims before their collision, whose subsequent morphological development is simulated by solving the two-phase incompressible Navier–Stokes equation with the adaptive mesh refinement technique, within the Basilisk software environment. We first derive analytical solutions predicting the unsteady interfacial and velocity profiles of the expanding sheet forming between the two rims, and develop scaling laws for the evolution of the lamella rim under capillary deceleration. We then analyse the formation and growth of transverse ligaments ejected from the lamella rims, which we find to originate from the initial corrugated geometry of the perturbed rim surface. Novel scaling models are proposed for predicting the decay of the ligament number density due to the ongoing ligament merging phenomenon, and found to agree well with the numerical results presented here. The role of the mechanism in breaking waves is discussed further and necessary next steps in the problem are identified.


Introduction
Liquid atomisation is a class of challenging multiphase problems (Obenauf & Sojka 2021) featuring a large separation of scales and various interacting physical mechanisms, which is of significance to numerous fields of application including meteorology (Villermaux & Bossa 2009), ink-jet printing (Castrejón-Pita et al. 2015, 2021;Lohse 2022), internal combustion engines (Yarin 2006), and pharmaceutical manufacturing (Mehta et al. 2017).Within the context of air-sea interactions, liquid fragmentation is primarily associated with wave breaking events, and gives rise to ocean sprays.These spray drops are then transported within the atmospheric boundary layer while exchanging with the latter moisture, momentum and heat during their lifetime; thus leaving their impact on both global and regional climates (Lhuissier & Villermaux 2012;Veron 2015;Deike 2022).Atomisation involves topological changes of the liquid bulk driven by external forces, typically followed by formation of corrugated ligaments subject to capillary breakup, and ends with a number of fragments featuring a broad size distribution, the knowledge of which is crucial for various areas of applications listed above (Villermaux 2007).
Among various types of liquid atomisation problems, the impact of liquid droplets has received much scholarly attention for its ubiquitous presence, rich dynamics and vast range of applications (Yarin 2006;Villermaux 2007;Cheng et al. 2022) since the pioneering experimental study of Worthington (1877).The impact process features a competition between inertial and capillary forces, which together with the characteristics of the impacting object and the surrounding gas phase shapes the final outcome of the original droplet: bouncing, spreading or splashing.Empowered by the rapid pace of sensor developments and increasing computational capacities, past works have elucidated considerable details about the ephemeral kinematic and morphological development of drops impacting with various types of surfaces (Cheng et al. 2022), including liquid films (Thoraval et al. 2013), deep pools (Agbaglah et al. 2015;Wang et al. 2023), smooth solid surfaces (Wildeman et al. 2016;Cimpeanu & Papageorgiou 2018), rough solid surfaces with friction (García-Geijo et al. 2021), or an identical droplet (He et al. 2019(He et al. , 2022)).Some recent works have also probed the dynamic properties of drop impact including the distribution of impact force and stresses, providing an alternative approach to investigate impact dynamics at early times (Cheng et al. 2022).Specifically, there have been a series of recent experimental works studying the high-speed impact of droplets with a surface of comparable size as a canonical unsteady fragmentation problem (Wang et al. 2018;Wang & Bourouiba 2017, 2018, 2021, 2022), providing valuable insights into various aspects of the impact process in the limit of large W e, including the self-similar evolution of the liquid-phase thickness and velocity profile (Wang & Bourouiba 2017), the dynamics of the rim bordering the expanding liquid sheet (Wang et al. 2018), the growth of liquid ligaments and the detachment of liquid drops from their tips (Wang & Bourouiba 2018, 2021), and the partition of mass, momentum and energy during the entire collision process (Wang & Bourouiba 2022).
However, the impact-induced fragmentation of liquid bulks featuring non-spherical initial shapes is also attested, which has received considerably less attention and remains less understood compared to drop impact problems (Liu et al. 2021).Among these is the collision of liquid rims, which has seen some recent investigations experimentally (Néel et al. 2020) and numerically (Agbaglah 2021) and is also the focus of the present work.In the two works cited above, an initially intact liquid film is perforated to form small holes, which then expand under surface tension and develop bordering rims travelling at the Taylor-Culick velocity (Taylor 1959;Culick 1960).Neighbouring film holes merge with one another when their bordering rims collide, oscillate and break up into small fragments, a process commonly observed during the rupture of films, which may be induced by rapid radial expansion of liquid shells (Vledouts et al. 2016) or the inflation of a liquid drop interacting with a surrounding airflow (Jackiw & Ashgriz 2022;Tang et al. 2023;Ling & Mahmood 2023).Agbaglah (2021) placed the two holes immediately adjacent to each other so that the two liquid rims collide at low impact Weber number W e; and the fused liquid bridge is found to pinch off under oscillation and form only a few small droplets.Néel et al. (2020) were able to investigate the rim collision phenomena by varying the initial distance between the two perforation sites, thus varying the impact W e value within the range of 50 ⩽ W e ⩽ 200.Apart from the primary capillary breakup mechanism of fused liquid bridges as discussed by Agbaglah (2021), Néel et al. (2020) also identified a critical W e value of 66 beyond which the decelerating transverse lamella rims develop elongating ligaments under the Rayleigh-Taylor (RT) instability, which then produce many secondary fine drops under capillary instabilities featuring a skewed size distribution.Additionally, there have been some early theoretical analyses on the capillary-driven coalescence of two liquid cylinders (Hopper 1993a,b;Eggers et al. 1999); although these are conducted at the creeping flow limit with negligible liquid bulk velocity, and thus may not be directly applicable to the current problem featuring finite collision speeds.
Apart from its presence during the rupture of thin films, the collision of liquid rims is also of specific interest due to its strong resemblance to the secondary wave splashing phenomena observed in ocean wave breaking events (Kiger & Duncan 2012; Mostert et al. 2022; Erinin et al. 2023a).Namely, two consecutive well-defined splashing phases have been identified (Kiger & Duncan 2012) during the lifetime of a deep-water plunging breaker.The first is the 'forward splashing' mechanism occurring right upon the reconnection of the overturning wave front and the sea surface below, which according to Mostert et al. (2022) only produces small amounts of droplets.After this, as shown in fig.1d, the wave bulk catches up with the decelerated splash-up generated from the initial impact at a relative speed u r close to the phase speed of the unbroken wave U g .At the indentation region between these two structures, the shape of the wave bulk and the splash-up can be approximated as two cylinders, whose cross sections feature radii of curvature r b and r s which are typically different.The rapid closure of the indentation region leads to the 'secondary splashing' phenomenon, which is characterised by a wall of vertically projected small droplets along the transverse direction as reproduced in figs.1a-1c, accompanied by air entrainment within the former indentation region (Kiger & Duncan 2012).Under experimental conditions, fragments generated via this mechanism comprise about one third of the total amount of droplets produced over the entire wavebreaking process (Erinin et al. 2023a); and Mostert et al. (2022) found that this splashing mechanism produces many fragments and can be curbed by strong surface tension (small Bo values).While Wang et al. (2016) noted the connection between the corrugated surface of the splash-up and the vortical structures beneath the wave surface, to our knowledge no existing study has analysed the physical mechanism governing the formation of these fragments, and their contribution to the overall droplet distribution associated with wave breaking remains unknown (Andreas et al. 1995;Kiger & Duncan 2012;Veron 2015).Furthermore, this splashing mechanism is found to produce many fragments close to the minimum grid size of Mostert et al. (2022); together with the highly transient nature of wave breaking and the presence of other fragmentation mechanisms, this indicates considerable difficulty in investigating the secondary splashing phenomena directly within the context of wave breaking.While we do not yet reproduce the fragment statistics seen in the studies above, the present work serves as a first step towards understanding the more complex wave splashing phenomena by retaining the major generation mechanism of splash fragments while leaving out many complicating factors, including size difference between the wave bulk and the splash-up evolving with time, and internal turbulent flow within the liquid phase; as a similar approach taken by Gao et al. (2021) reveals the connection between the bubble size distributions of destabilising air cylinders and air cavities entrained by plunging breakers.
In this work, we conduct a comprehensive investigation of the collision between two liquid cylinders with identical size, covering the entire deformation and fragmentation period.The direct comparison and establishment of connections between the rim collision results and the statistics of the secondary wave splashing phenomenon are left for future work, together with the role of gravity and the difference between the sizes of the wave bulk (r b ) and the splash-up (r s ) which complicate the early-time rim dynamics.Twophase numerical simulations are conducted to derive detailed flow field information during this highly transient collision process.Our study is structured as follows.We first present in §2.1 the problem configuration and the parameter space of the current work, and then introduce the numerical method in §2.2.After providing an overview of the rim collision phenomena in §3, we quantitatively analyse the development of each part of the expanding liquid bulk successively following a spatial order, namely the kinematics of the spreading liquid sheet ( §4.1) and its bordering rim ( §4.2), the growth and merge of transverse ligaments topping the rim ( §5), and the statistics of fragments shed from the ligaments ( §6).We conclude the study in §7 with some remarks on future work.

Problem description
The geometrical configuration for the rim collision problem is shown in fig.2a, where two infinitely long cylindrical liquid rims with diameter d 0 , density ρ l and viscosity µ l are aligned along the x axis, and set to travel along the z axis with uniform velocities of opposite signs and the same magnitude U 0 .The liquid cylinders are surrounded by an inert gas phase with density ρ g and viscosity µ g , and the liquid-gas interface is characterised by a surface tension coefficient σ.It is noted that gravitational effects have been neglected in the current setup; and differing from the configurations of Néel et al. (2020) and Agbaglah (2021), there is no interstitial film connecting the two approaching cylinders.Consequently, four non-dimensional controlling parameters can be written for this problem: where W e and Oh are respectively the Weber and Ohnesorge number comparing inertial and viscous effects to capillary forces, and ρ * and µ * are respectively the density and viscosity ratios of the liquid and gas phase.In this work, Oh is set as 0.01 in most of the simulations, whereas its influence on the fragment statistics is briefly discussed in §A.ρ * and µ * are set as 830 and 55, respectively, which are typical values for the air-water system (Pairetti et al. 2018).
We set the width D of the cubic simulation domain as 10d 0 to allow enough space for the morphological evolution of the coalesced liquid structure.Utilising the symmetry of the splashing phenomena about the xz-plane, we only model the merging of the upper halves of the two liquid cylinders to save computational resources.A symmetric boundary condition is therefore applied at the bottom, and an outflow boundary condition is imposed on the top boundary so that fragments produced from the collision can leave the domain from there at late time; while the other boundary conditions are set as periodic.
To investigate the sensitivity of the fragmentation process to the initial conditions (Liu & Bothe 2016;Berny et al. 2022), and also taking into account the surface corrugation of the splash-up in wave breaking events (Kiger & Duncan 2012) which still has not been quantified according to our knowledge; we introduce random transverse perturbation within a certain wavelength range on the cross-sectional radius of the two cylindrical rims (Pal et al. 2021), in the form of filtered white noise signals characterised by the following two parameters, where ε 0 is the non-dimensionalised characteristic amplitude of perturbation, and N max defines the highest wavenumber among the spectrum of the perturbation signal.The filtered white noise signal is the default type of initial interface perturbation we impose on the rims, as Zhang et al. (2010) did for analysing the linear stability of the crown splash.Occasionally we also impose single-wavelength sinusoidal perturbations, or a superposition of sinusoidal perturbations with wavelengths λ = D/8, D/16 and D/32 for comparison, which will be explicitly denoted by 'Sing.' and 'Sup.' when reported.
In the case of single-wavelength perturbations, N max corresponds to the perturbation wavenumber.
As discussed above, W e, ε 0 and N max constitute the parameter space for the present study.Among these, W e is varied between 60 and 280 where the coalesced liquid bulk expands vertically to form a lamella, and ε 0 is set as 0.02, 0.04 and 0.06, within the limit of small radial perturbations.N max varies between 15 and 80, whose influence will be discussed in detail in §5.2.

Numerical method
The open-source scientific computation toolbox Basilisk (Popinet 2019) is used in this work to solve the two-phase nonlinear, incompressible, variable-density Navier-Stokes equations.A second-order accurate discretisation is applied in both space and time, and a geometric volume-of-fluid (VOF) method in a momentum-conserving formulation is used to maintain a sharp representation of the liquid-gas interface while minimising the parasitic currents induced by surface tension (Popinet 2018;Tang et al. 2021).Capillary effects are modelled as source terms in the Navier-Stokes equations using an adaptation of Brackbill's method (Brackbill et al. 1992;Popinet 2009), which calculates the interface curvature by taking the finite-difference discretisation of the derivatives of interface height functions (Popinet 2009).The octree-based adaptive mesh refinement (AMR) scheme based on the estimation of local discretisation errors of the VOF function f and flow velocity u is adopted so as to reduce the computational cost at high resolution levels L, which is defined using the minimum grid size, 3) The results presented in the main body of this work are obtained from three-dimensional simulations at L = 10, at which the late-time evolution of interface profile and liquidphase energetics have reached grid independence.The numerical convergence of fragment statistics is also established for fragments with diameter larger than 4∆ 10 , which is discussed in detail in §A.Results from some two-dimensional simulations are also presented for comparison with three-dimensional rim dynamics ( §4.2), and to investigate lamella foot formation at very early time ( §B).
Since the present study focuses on the liquid-phase morphological development after the rims begin to coalesce, the distance between the symmetric axes of the two cylindrical rims are set at initialisation as ∆z c = 0.95d 0 , so that the slightly overlapped rims form a line of contact.The interfacial perturbations on the rims are introduced as follows.Firstly, discrete white noise signals with unit variance σ 2 are produced using the random number generator provided in Basilisk.Fig. 2b shows the ensemble-averaged power density spectra of the white noise signals generated in Basilisk at different number of realisations N samp .It is observed that the power density P (f ) is close to the theoretical value of σ 2 at all frequencies f , matching the requirement of frequency independence for white noise signals.Next, we apply a low-pass filter on these signals so that only the lowest N max wavelengths are preserved, and the filtered signal is normalised so that its standard deviation becomes ϵ 0 = 0.5εd 0 .The normalised signals η are then mapped onto the transverse radius profile R(x) of the liquid rims, in the form of R(x) = R 0 + η(x).The two liquid cylinders being positioned close to each other ensures that they coalesce immediately when the simulation starts; and there is not sufficient time for capillary effects to smooth out the perturbations, or amplify them to trigger the Rayleigh-Plateau (RP) instability (Pal et al. 2021) so that the cylinders break up prematurely.

Overview of rim splashing
Here we qualitatively describe the rim splashing process as observed in the simulations before analysing the detailed dynamics at each stage in the following sections.Fig. 3 presents the isometric view of the splashing phenomenon at W e = 200, ε = 0.06 and N max = 25, whereas fig. 4 shows the side view for a few different (W e, ε) configurations with N max = 25.Tiny air bubbles are entrained within the indentation space between the two rims following the impact (Thoraval et al. 2013;Josserand et al. 2016;Erinin et al. 2023a), which in our case have no known effects on the subsequent development of liquid-phase flow fields.The rims coalesce rapidly, causing a dramatic increase in the local liquid-phase pressure (Néel et al. 2020).A very large pressure gradient arises at the surface 'indentation' between the two cylinders since the air-phase pressure remains unchanged, leading to the vertical acceleration of fluids near the surface (Longuet-Higgins 2001).The magnitude of the vertical velocity near the contact line can be a few times larger than the initial horizontal velocity U 0 , causing the contact line to advance upwards rapidly with fluid particles nearby converging to it, while the high-pressure region follows it closely rather than remaining on the axis of symmetry (Philippi et al. 2016).A thin transverse liquid lamella (Néel et al. 2020) appears as shown in figs.3b, and the location of the lamella foot agrees well with Eq. (B 12) developed in §B using potential flow theory (Riboux & Gordillo 2014).The liquid lamella is also observed in the experimental study of Goswami & Hardalupas (2023) where two neighbouring drops impact with a solid surface simultaneously and expand into contact.Within the parameter space of the present study, this transverse lamella continues to expand vertically under the pressure difference between the two phases and consumes the two impacting cylinders, evolving into a thin film aligned with the xy-plane; in the meantime it is continuously decelerated by capillary force, forming a thick rim at its upper border.
Without the initial interface perturbation, the liquid rim would remain intact and mostly smooth during the entire simulation period up to t/τ cap = 2.73, where τ cap ≡ ρ l d 3 0 /8σ is the capillary time scale corresponding to the initial rim diameter.When perturbation is introduced, we observe the amplification of transverse perturbation waves on the rim as shown in fig. 4. Depending on the value of W e, these waves will either slowly increase in amplitude and reach nonlinear development, whose growth rate increases with the initial perturbation amplitude ε 0 , as shown in the first and last row of fig.4; or generate slender ligaments which continue to elongate on the top of 'cusp' structures  4, the number density of transverse ligaments (or the characteristic perturbation wave number when no such ligaments form) decreases over time rather than remaining a constant; a phenomenon which we will analyse in more details in §5.In the meantime, lamella expansion gradually slows down as the vertical position of its bordering rim reaches saturation, and begins to slightly retract at low W e values as can be seen in figs.4b and 4c.As is the paradigm of droplet formation in many previous fragmentation studies (Villermaux 2007), here the transverse ligaments act as the direct and only source of fine drops as the latter intermittently detach from the ligament tips, whose statistics will be presented in §6.Under the restoring effects of capillary force, most of the fragments undergo prolate-oblate shape oscillations after their detachment, and a small portion of them may cross path and merge to form larger droplets during their flight (Tang et al. 2023).Note that in our simulations all fragments keep moving upwards before crossing over the top boundary and leaving the simulation domain; while in wave splashing scenarios without wind forcing they will ultimately fall back to the sea surface under gravity and be destroyed (Mostert et al. 2022).

Liquid lamella expansion
In this section, we study the dynamics of the liquid lamellae consisting of the expanding sheet ( §4.1) and its bordering rim ( §4.2).To simplify the problem which features random initial perturbation, here we do not consider the formation of liquid ligaments.Instead, we average the transverse cross section of the coalesced liquid bulk along the x axis following Wang & Bourouiba (2017) so that the lamella expansion process can be described in a quasi-one-dimensional manner, characterised by one-dimensional velocity and thickness profiles u y (y, t) and h(y, t).

Liquid sheet kinematics
The unsteady evolution of the expanding sheet profile is of both fundamental and practical importance for impact problems, especially for predicting their maximum spread radius (Yarin 2006;Josserand & Thoroddsen 2016;Wang & Bourouiba 2017), and the first step towards modelling the profile evolution of expanding sheets is to understand the fluid motion within them.
We first derive the velocity profile within liquid lamella sheets in Cartesian coordinates for W e ≫ 1 and Oh ≪ 1.At early times, the vertically expanding sheet is bounded at its lower end by the line of collision between the two cylinders, which we call the lamella foot.The trajectory of this foot is given by y n = 2 U 0 t/R 0 , per the analysis given in §B and following Gordillo et al. (2019).At later times, this lamella foot becomes increasingly indistinct as the colliding cylinders merge.
We now proceed to discuss the kinematics of the fluid within the vertical lamella sheet itself, subject to these considerations, following the general analytical strategy of Wang & Bourouiba (2017).After neglecting viscous, compressibility and capillary effects, the quasi-one-dimensional momentum equation in the vertical direction can be written as which may be presented alternatively in the Lagrangian form along the characteristics where D/Dt is the total derivative.(4.2) suggests that the velocity of a fluid particle remains unchanged within the liquid sheet as it travels vertically upwards.We can then integrate along the characteristics and obtain the motion of fluid particles within a Lagrangian frame following it, where u y and ξ are the initial vertical velocity and position of fluid particles.Following Yarin & Weiss (1995), Wang & Bourouiba (2017) assume that the initial velocity is proportional to the initial position, namely u y = βξ, similar to the velocity field of a stagnation-point flow.Consequently, (4.3) may be rewritten in the Eulerian reference frame as where 1/β corresponds to the time needed to set up the post-collisional velocity profile within the liquid bulk, which according to Wang & Bourouiba (2017) is at the same order of the collision timescale d/U 0 .As sheet expansion further progresses so that t ≫ 1/β, (4.4) asymptotes to We measure the vertical velocity profile within the expanding sheet from our numerical simulations at different times for W e = 120, and first plot them in the inset of fig.5a.Consistent with our assumption, the vertical velocity u y is observed to scale linearly with the vertical position y, with the slope decreasing over time.In the main plot we rescale the horizontal axis by U 0 t, after which the velocity profiles at different times all collapse onto a single straight line y = x for y ⩽ 2U 0 t, agreeing with the theoretical model (4.5).
For W e = 120, the collision timescale is d/U 0 ≈ 0.52τ cap , while fig.5a indicates that the velocity profile at t/τ cap = 0.73 is already very close to the analytical solution (4.5); which suggests that in fact the initialisation timescale 1/β ≪ d/U 0 , and model (4.5) is applicable to the early collision stage where t ≈ d/U 0 .Note that as time elapses, the vertical velocity deviates from Eq. (4.5) at progressively smaller values of y/U 0 t, where the fluid parcels move away from the sheet towards the bordering rim.Indeed, in fig.5b we also plot the vertical velocity normalised by y/t at different times, W e values and perturbation waveforms.All results presented therein feature a range where the normalised velocity u y t/y = 1, suggesting that Eq. (4.5) remains valid at different initial configurations within our current parameter space.
Having established the liquid velocity profile within the lamella sheet, we can further solve for its thickness profile utilising the continuity equation, which can be written as follows for a thin sheet expanding in the y direction, combined with Eq. (4.5), this yields This can be solved by separation of variables in the form of h(y, t) = f (t)g(y) to obtain the evolution of sheet thickness h, which is written in a non-dimensional formulation as note that this result differs from the axisymmetric configuration where hU 2 0 t 2 /R 3 0 evolves self-similarly without an explicit time dependence (Wang & Bourouiba 2017;Gordillo et al. 2019).
We plot the lamella thickness profiles for W e = 120 in fig.6a, where it can be seen that the 'bulge' centred around y = 0 gradually flattens into an extended thin liquid sheet, pushing the bordering rim further along the vertical direction.The bordering rim is connected to the sheet via a neck, reminiscent of capillary waves upstream of inviscid liquid rims receding at the Taylor-Culick velocity (Savva & Bush 2009).Fig. 6b further shows the profiles rescaled by Eq. 4.8.The 'bulk' region where y/R 0 ⩽ 2U 0 t/R 0 initially retains its cylindrical shape during the initialisation period t ∼ 1/β, and only comes to agreement with Eq. (4.8) when the lamella foot disappears and the bulk can no longer be decisively told apart from the lamella.The non-dimensionalised lamella profile in fig.6b is found to be well described by the following functional form f (x), (4.9) Lastly, we note that although the liquid lamella expands vertically to form a thin film, the latter does not suffer from spontaneous perforation during our simulation period, although this is observed in fig.14 of Vledouts et al. (2016) and marks the onset of bag film fragmentation in droplet aerobreakup problems (Tang et al. 2023;Ling & Mahmood 2023), where the liquid films are subject to radial accelerations.Neither does the lamella film experience destabilisation under shear force arising from its interaction with the surrounding gas phase (Villermaux & Bossa 2011;Riboux & Gordillo 2015), which may arise at larger W e values.We therefore do not take special measures to artificially stabilise (Liu & Bothe 2016) or perforate (Chirco et al. 2022) the expanding lamellae in this study.

Bordering rim evolution
As is discussed in §4.1, while the expanding lamella sheet abides by the velocity profile (4.5) and the self-similar thickness profile (4.8), these two models break down for the bordering rim, which demands separate scaling laws to describe its kinemamics.Similar rim structures are ubiquitous in impact problems and act as the crucial link between the expanding sheet and shedding droplets (Wang et al. 2018), and understanding their motion lays the foundation for further theoretical analysis of ligament merging, which we will perform in §5.2.
In a quasi-one-dimensional framework, the kinematics of the bordering rim can be characterised by two parameters, namely its average vertical position y rim and diameter b rim .We first present their evolution at different W e values in figs.7a and 7b, respectively; where time and length are scaled using U 0 /R 0 and R 0 .It is seen that both y rim and b rim first increase with time and eventually saturate; and y rim and b rim show different dependence on W e, as the former increases and the latter decreases with W e, although these W e-dependencies have become very subtle by W e = 200.Further, the evolution of y rim for W e = 200, ε = 0.06 and W e = 200, ε = 0.04 in fig.7a is virtually the same, which shows that the dynamic behaviour of the rim is largely independent of the initial perturbation amplitude ε.
We now seek to further compare our numerical results obtained in figs.7a and 7b with available theoretical models.Following Gordillo et al. (2019), the following mass-and momentum-conservation equations can be proposed in the non-dimensionalised form to describe the dynamic behaviour of the advancing lamella rim in the inviscid limit, ) where the cylinder radius R 0 , initial impact velocity U 0 and their quotient R 0 /U 0 are chosen as the reference scales for length, velocity and time.The numerical solutions of the ordinary differential equation (ODE) system (4.10)-(4.12) at various W e values, with initial conditions as defined in §B, are presented in figs.7a and 7b respectively as transparent solid lines.Since most of our measurements from three-dimensional simulations are taken when U 0 t/R 0 > 1, we also present results of two-dimensional simulations at W e = 80 and 160 using solid dots, which extend to much earlier times.An excellent agreement between the predictions of (4.10)-(4.12)and the two-dimensional results is observed.Three-dimensional measurements of y rim and b rim are found to saturate earlier and become smaller than their two-dimensional counterparts and theoretical predictions as time elapses.These observations can be explained by taking into account the formation of ligaments on the lamella rim.This mechanism arises due to highly nonlinear transverse perturbations imposed on the liquid rim, which is not present in two-dimensional simulations.At sufficiently large times, the growth of ligaments causes a substantial mass flux away from the liquid rim, which becomes more prominent at higher W e values as is shown in fig. 4. Other factors may also play a role, e.g., the initial overlapping between the two rims in three-dimensional simulations.
In figs.7c and 7d we present the evolution of y rim and b rim again, where time is now non-dimensionalised using the capillary timescale τ cap .The growth of both y rim and b rim in three-dimensional simulations is consistent with a power law of t/τ cap up to t/τ cap ≈ 1.4, after which deviation from this power law is observed.It is also found that prefactors of √ W e (for y rim ) and W e −1/4 (for b rim ) can collapse the evolution data reasonably well, especially for b rim as a single master curve is clearly observed in fig.7d, leading to the following scaling arguments, which we will use for further theoretical analysis of the ligament merging phenomenon in §5.2, while a complete determination for (4.13) valid for very late time remains for future work.The scaling of W e in (4.13) agree with the experimental results of Villermaux & Bossa (2011) and Wang et al. (2018) for drops impacting a small surface; although the dependence on time is different, as their models aim at describing the entire sheet expansion-retraction motion.Nonetheless, similar empirical √ t scalings have been proposed by Mongruel et al. (2009), Thoroddsen et al. (2012) and Visser et al. (2015) for quantifying the radial position of the expanding lamellae in the inertial regime, indicating that this simple form of time dependence can still describe the rim kinematics reasonably well between its formation and the onset of the retraction motion.Note that (4.9) implies the evolution of the liquid momentum carried by the rim p rim ≡ πρ l Db 2 rim ẏrim is independent of the values of W e, even though its average vertical velocity ẏrim and volume Ω rim ≡ πb 2 rim D do depend on W e. This contrasts with the axisymmetric results of Wang & Bourouiba (2022) where the rim volume remains independent of W e due to their axisymmetric configuration.

Formation and growth
In the scenario considered by Wang & Bourouiba (2021), liquid ligaments grow slowly out of the corrugated bordering rim along its azimuthal direction, which they ascribed to a combination of local geometry, pulling effects of inertial force associated with rim deceleration and the global liquid-phase mass conservation.For this study, and given our perturbation profile, we find the transverse ligaments form very early for W e ⩾ 120, nearly at the same time when the lamella is born out of the indentation region between the two cylinders, as presented in the simulation snapshots of fig. 8.A closer look at fig. 8b reveals that the ligaments are produced preferentially from the concave regions along the two perturbed cylinders, suggesting that the ligament formation process is closely associated with the initial rim perturbation waveform.Indeed, a previous investigation by Gordillo et al. (2020) suggests that the ligaments originate from the non-uniform initial distribution of normal interface velocity, which is in turn determined by the upstream liquid velocity and the curved initial cavity profile at the moment of impact.Note also that the nascent ligaments do not always project exactly vertically along the vertical yaxis; rather, they can display complex twisting and surface oscillation motions, and may grow in an oblique direction.At the same time, the liquid sheet beneath its bordering rim also features a wavy surface not perfectly aligned with the xy-plane.These are most likely due to the difference in the initial interface perturbations imposed on the two colliding rims, which gives rise to oscillatory motions under capillary effects.
We now move on to discuss the subsequent growth of the height of these ligaments after their formation, which is shown in fig.9a for rim collision at W e = 120 and 160, where we present the height evolution of three individual ligaments at each W e. The shedding of the first fine drop from these ligaments is characterised by a kink around t = 0.4τ cap .While this initial pinch-off may happen at even earlier times as W e increases, as shown in fig.7 of Wang & Bourouiba (2018); this is still much later than the onset of the micro-splashing phenomena investigated by Thoroddsen et al. (2012), which happens at t/τ cap = O(10 −4 ) (see e.g., their fig.5b).Before this first pinch-off event, the height h lig of different ligaments increases following a similar trend, and the growth rate at W e = 160 is higher than that at W e = 120.Given that the initial phase of growth of these ligaments is governed primarily by inertio-capillary effects, we compare it in fig.9a with the selfsimilar power law of h ∝ t 2/3 , proposed by Lai et al. (2018) in their investigation of inertio-capillary-dominated collapse of small surface bubbles.It is found that the height increase of the majority of transverse ligaments (except Ligament 1 at W e = 120) agrees better with the linear growth model.This is most likely because the formation of fast jets observed by Lai et al. (2018) is preceded by focusing of interfacial capillary waves at the bottom of the bubble cavity while the liquid bulk remains quiescent; whereas here the bulk velocity plays a vital role in driving the closure of rim indentation, and may thus modify the rate of jet growth.In addition, a recent study by Gordillo & Blanco-Rodríguez (2023) suggests that the exponent of power laws dictating the evolution of jet radius and speed is dependent on the initial geometry of the collapsing air cavity, which may also account for the difference between our results and those of Lai et al. (2018) since the concave regions on our perturbed cylindrical surfaces do not feature a uniform radius of curvature.
To better understand the fluid motion within the growing ligaments, we measure the liquid-phase vertical velocity u y from the bottom of the expanding sheet up to the tip of a single ligament, and plot it in fig.9b as a function of the vertical coordinate y.Two linear scaling regimes are observed; the first one is well described by u y = y/t, which corresponds back to the velocity profile (4.5) we established for the expanding sheet.As the fluid particles move higher up and away from the sheet, its velocity first increases and then abruptly decreases as it enters the neck and rim region respectively; a similar abrupt deceleration is also noted by Wang & Bourouiba (2022) for drop collision problems.Interestingly, the combination of the neck and rim causes a constant decrease in the vertical velocity of approximately 0.7U 0 for different W e values, times and perturbation waveforms.This pattern is not explicable from the scaling model (4.13) since according to it, the fluid velocity should be equal to the average rim velocity ẏrim after deceleration, which is always one half of the sheet velocity y rim /t before deceleration.The implication is that the fluid velocity at the ligament root is always faster than the average rim velocity, which Wang & Bourouiba (2021) ascribed to the additional acceleration due to the interface curvature at the rim-ligament junction.After this constant offset at the bordering rim which is most likely a viscous effect (Ghabache et al. 2014), u y grows linearly again with the vertical position y with the same slope as the sheet region.This second linear scaling regime most likely corresponds to the ballistic region (although in the present study gravity is not included) in Worthington jets identified by Gekle & Gordillo (2010), where the liquid particles travel at constant speed upwards before being slowed down once again at the bulb by capillary effects.
While Wang & Bourouiba (2022) were able to demonstrate theoretically that energy dissipation at the bordering rim is responsible for the local rapid deceleration of fluid particles, which we also observed in fig.9b, the detailed nonlinear dissipation mechanism is not captured by their one-dimensional model.As discussed in §1, there is still a dissipation deficit of 15% not covered by their calculations occurring at early-time.To help elucidate this discrepancy, we present contour plots in fig. 10 for a simulation case at W e = 200, showing the distribution of the liquid-phase viscous dissipation rate ϵ d ≡ µ l (∂u i /∂x j )(∂u j /∂x i ) within the centre-plane x = 0 for t/τ cap ⩽ 0.36, covering the early deformation period t/τ cap ⩽ 0.2 where Wang & Bourouiba (2022) observed the dissipation deficit (see their fig.19).It is found that when the lamella is born at very early time and its bordering rim has not yet fully developed (fig.10a), there is an extremely high concentration of ϵ d located at the lamella foot, agreeing with the simulation results of Wildeman et al. (2016) for drops impacting a smooth surface (see their fig.4) and Fudge et al. (2023) for drops impacting a liquid pool (see their fig.8b).As time elapses, the bordering rim takes shape with its neck featuring relatively high concentration of ϵ d ; whereas the dissipation at the lamella foot weakens and eventually becomes negligible by t/τ cap = 0.364 (fig.10d), matching the saturation trend shown in fig.19 of Wang & Bourouiba (2022) for the dissipation deficit.This decay pattern of liquid-phase dissipation may be explained as follows.The liquid particles feeding the lamella at early time mostly comes from a very narrow boundary straddling the corners on either side of the lamella foot (Riboux & Gordillo 2014); and they generate vorticity (Batchelor 2000;Li et al. 2018) and experience capillary deceleration while traversing the highly-curved free surface, leading to very large values of viscous dissipation.Since this early-time deceleration is geometry-induced, Fudge et al. (2023) further hypothesised that the magnitude of the velocity gradient remains largely unchanged at different flow configurations so that the early-time dissipation rate is proportional to the liquid viscosity µ l ; although this is not directly verified in the present work.As the liquid sheet extends further upwards, the interface curvature at the lamella foot decreases and capillary deceleration becomes much weaker, hence the decrease in the dissipation rate ϵ d .Eventually, as the lamella foot is no longer discernable from the flattened liquid bulk and the capillary effects become negligible, the inviscid sheet velocity profile (4.5) will be established.
These observations of the liquid-phase dissipation evolution therefore suggest that the unidentified dissipation deficit found by Wang & Bourouiba (2022) is most likely linked with the early-time lamella formation which also features strong viscous effects; and according to Wildeman et al. (2016) and Ó Náraigh & Mairal (2023), this kind of dissipation is a type of 'general head loss' imposed by the deformation mode of the impacting object and independent of the detailed impact parameters.The head loss in impact problems dissipates a fixed fraction of the initial kinetic energy via recirculating flows (Wildeman et al. 2016;Ó Náraigh & Mairal 2023), which matches the observations of Wang & Bourouiba (2022).It is also noted that besides the early-time lamella foot and the late-time rim neck, the flow field elsewhere within the coalesced liquid bulk shown in 10 is nearly inviscid, supporting our derivation of the centerline velocity profile within the liquid sheet (4.5) based on inviscid flow assumptions.While we have not fully explored the influence of Oh on the lamella expansion process, it can be expected that in the inviscid limit where Oh ≪ 1, its influence will be confined to the narrow 'viscous' regions identified in this section, and thus will not significantly affect the deformation of the droplet bulk, or the ligament dynamics and fragmentation mechanisms to be discussed below.

Ligament merging phenomenon
As is noted in §3, the number of transverse ligaments on the rim will decrease over time as they merge with their neighbours.This merging process is observed only when the initial perturbation is not monochromatic such that the nascent ligaments are not equidistantly spaced, and it turns out crucial in maintaining the growth of ligaments and the continuation of fragment shedding.Figs.11a-11c show the development of ligaments formed out of monochromatic perturbation waveforms, and it is observed that after two rounds of pinching-off events at their tips, the ligaments can no longer sustain their own growth and are re-absorbed back into the underlying liquid rim; whereas ligaments formed out of filtered white noise perturbations at the same values of W e, ε and N max merge with their neighbours and survive end-pinching, continuing to grow and shed fragments until the end of simulation, as observed in figs.11d-11f.
The ligament merging process is shown in more details in the snapshots of fig.12, where the merging ligaments are observed to be located on rim 'cusp' structures extruding from the liquid sheet beneath, an indication of the non-uniform incoming mass distribution along the bordering rim (Gordillo et al. 2014;Wang & Bourouiba 2018).When ligament merging occurs, the roots of two neighbouring ligaments approach each other; liquid is then drawn upwards from the underlying cusp in between the two ligaments, causing them to coalesce into a thicker and more corrugated ligament; while the length of the fused ligament does not differ significantly from those of its parents.Ligament merging is thus capable of delaying the next end-pinching event since the fused ligament takes longer to be stretched, thus allowing incoming fluid from the rim to sustain their growth and evade re-absorption into the rim when the depleted ligament becomes too short, as observed in fig.11.Note that the monotonically decreasing trend of ligament numbers observed by us and Wang & Bourouiba (2021) differs from the early experimental work of Thoroddsen & Sakakibara (1998) on drop impact, where they also observed splitting of liquid fingers aside from their merging behaviour, so that the total finger number remains approximately constant.
To the knowledge of the authors, a scaling law for the ligament numbers accounting for their merging dynamics is not yet available.The early work of Marmanis & Thoroddsen (1996) found that the number of liquid fingers at the maximum spread radius for high-speed drop impacts scales with Re 3/4 (Re being the impact Reynolds number), without accounting for their evolution over time; whereas the recent W e 3/8 scaling model proposed by Wang & Bourouiba (2021) is based on the understanding that the ligaments are formed from a subset of rim corrugations arising from a combined Rayleigh-Taylor (RT) and Rayleigh-Plateau (RP) instability; a physical mechanism which is most likely not yet active within our parameter space as we find the formation of ligaments more closely linked with the initial interface perturbation geometry.The early theoretical analysis of Yarin & Weiss (1995) predicted that perturbed free rims will spontaneously develop 'cusp' structures due to nonlinearity, where two neighbouring rim sections impinge and give rise to free jets.This physical picture closely resembles our present observations, although Yarin & Weiss (1995) did not proceed to develop scaling models for the splashing fragments.We therefore seek to derive a new ligament number scaling model accounting for the merging dynamics in this section, and compare it with the simulation results.
As the first step towards quantifying ligament dynamics, we seek to determine the evolution of the average ligament width w lig over time by inspecting the evolution of fragment diameter d frag , which can be easily reconstructed using the droplet-tracking algorithm of Chan et al. (2021).The connection between these two quantities is established in fig.13a, where we plot the ratio between d frag and w lig at the instant of pinch-off.The ligament diameters are measured at the cross section corresponding to one half of the total ligament length.Most of the measured data are found to scatter between 1 and 2, centred around 1.4 which is close to the average value of 1.5 as found by Wang & Bourouiba (2018).These are below the theoretical value of 1.89 as predicted by the RP instability (Gordillo & Gekle 2010), indicating that end-pinching is indeed the dominant fragment production mechanism for the present study, and that the diameter of fragments d frag remains in proportion to the width of their parent ligaments w lig .The inset of fig.13b shows the diameters of the fragments d frag versus their time of formation, where the data for each configuration of (W e, ε) have been averaged across three individual realisations to reduce the range of scatter.It is found that larger fragments are generally produced at later times and smaller W e values.The scatter in data increases over time, which is most likely due to the increase in the diameter difference and the surface corrugation of remaining ligaments as they merge with one other.The main plot suggests that the evolution of individual fragment diameter in fig.13b roughly scales with W e −3/4 t/τ cap , while noting that our fitted prefactor W e −3/4 is not conclusive and remains to be further validated at higher W e values.Since the fragment size remains in proportional to the width of their parent ligament according to fig.13a, it can be inferred that (5.1) Gordillo et al. (2014) and Wang & Bourouiba (2018) proposed the following drift velocity of ligaments u drift on top of a liquid cusp to characterise their migration, where u s (y rim ) = y rim /t is the velocity at which liquids enters the rim from the expanding liquid sheet, calculated from the self-similar expanding sheet velocity profile (4.5); and u rim is the vertical rim velocity determined by differentiating the scaling law for y rim in (4.13).It is noted that both u s (y rim ) and u rim scale as √ W e • t −1/2 , and α is the angle between the local rim and the horizontal plane as shown in fig.14.Overall, (5.2) suggests that the migration of ligaments is driven by the tangential projection of the net incoming fluid velocity along the corrugated rim.
As u drift drives the ligament migration and causes the ligament number density N lig to decrease, the average transverse ligament spacing 1/N lig on the rim consequently suggesting that the rim slope α remains unchanged over time and depends only upon the impact Weber number.Taking into account that tan α ≈ ε rim N lig , α remaining constant also indicates that as the average ligament spacing 1/N lig becomes larger over time owing to the ongoing merging of ligaments, the rim corrugation ε rim also increases proportionally to maintain a constant local slope.By further incorporating (5.3) and (5.2), the following model predicting the evolution of the ligament number density N lig can be derived, (5.6) In the main plots of fig.15 we show the decay of the ligament number density N lig at different values of N max and W e. Fig. 15a shows that while increasing N max leads to the formation of more ligaments at early time, N lig appears to reach saturation and does not increase proportionally with N max when t/τ cap = 0.18; and it is likely that viscous damping effects are particularly strong when the lamella is ejected at the early impact stage, which may smooth out short-wavelength components in the initial perturbation spectra (see our fig.10 and relevant discussions).Larger N max also leads to faster decay of N lig as expected, since the average distance between neighbouring ligaments becomes smaller; and for all values of N max , fig.15a suggests that the evolution of N lig becomes largely similar for t/τ cap > 1.5, where the remaining ligaments take much longer to migrate and merge.Fig. 15b shows that the decay of N lig does not appear to depend strongly on W e and thus lends some support to the prediction of (5.6) that it is W eindependent, although ensemble averaging would be needed to ascertain this due to the randomness in the initial perturbation waveform.The insets of fig.15 show the evolution of N lig again in log-log axes, which is also compared with model (5.6).The inset of fig.15a indicates that the measured results collapse well and agree with (5.6) for N max ⩾ 35 throughout the period of measurement, whereas at N max = 25 the decay of N lig is initially slower, and only matches the prediction of t −1/2 for t/τ cap ⩾ 1.The measurement of N lig extends to earlier times in fig.15b, and its inset suggests that the evolution of N lig at different W e is initially slower and matches model (5.6)only after t/τ cap ⩾ 0.2.This slower decay at early times observed in both insets arises most likely because the rim slope α takes a finite period of time to develop before reaching the steady-state value given by (5.5).Overall, these results indicate that model (5.6) offers a good working description of the ligament merging phenomenon occurring within our parameter space.However, the approximations we make for its derivation restricts its application to the scenario when N lig is large and the slope α of the corrugated rim has reached the steady-state value predicted by (5.5).

Droplet generation and characteristics
Since the expanding sheet remains intact during our simulation period, fragments are formed solely through the breakup of liquid ligaments originating from the bordering rim.The majority of fragments are produced via the end-pinching mechanism, a process that coincides with the ligament merging phenomena and is already visible in the snapshots presented in figs.4 and 12. Namely, the ligament tips are decelerated by capillary force and produce surface corrugations, which in turn creates a local pressure gradient within the ligament neck that drains the liquid towards the tip and triggers pinch-off.After this, the enlarged ligament tip detaches as a primary drop (Gordillo & Gekle 2010), whose diameter is proportional to the parent ligament width as shown in fig.13a.Occasionally, satellite drops are from the small amount of remnant liquid within the neck before it can be fully reabsorbed into the ligament after end-pinching, as also shown in fig.5b of Wang & Bourouiba (2018).In contrast with the primary drops, the production of these satellite drops is sensitive to initial liquid-phase velocity perturbations, thus introducing randomness to the jet breakup process.However, they do not affect the size and velocity of the subsequent primary drops ejected, as shown by Berny et al. (2022) in the instance of jet-droplet production by bubble-bursting.Possibly due to the difference in their generation mechanisms, our ligaments can grow much longer than those Wang & Bourouiba (2018) observed in their experiments; and at late times some particularly long ligaments may break up due to the RP instability and shed multiple primary fragments at a time, as can be seen in the rightmost ligament in fig.4e and fig.2f of Liu et al. (2022).This appears inconsistent with the conclusion of Wang & Bourouiba (2021) that a ligament can only pinch off to produce one primary drop at a time under the chaotic dripping regime.We first analyse the time evolution of various fragment properties during the rim collision process.Fig. 16a shows the evolution of the total number density of primary fragments N frag , where it is observed to generally increase over time towards saturation, despite infrequent decreases due to coalescence of primary fragments with another fragment or a neighbouring ligament.It is also noted that N frag increases with both W e and ε, as larger W e and ε encourages the growth and subsequent pinch-off of liquid ligaments.
The decrease of drop production rate ∆N frag /∆t can be explained as follows.The ongoing ligament merging phenomenon causes the ligament number density to decrease, hence fewer fragments can detach at the same time.In the meanwhile, merged ligaments become more corrugated and thicker, thus end-pinching events occur at larger ligament widths as time elapses.Consequently, the time interval between successive end-pinching events also becomes longer, since the necking timescale increases with w i according to the experimental results of Wang & Bourouiba (2018).Fig. 16b shows the evolution of the fragment ejection speed u ej , the speed of the drop at the moment it detaches from its parent ligament.The fastest drops found in our simulations feature u ej slightly over 3U 0 , which is comparable to the fragment ejection velocity in the prompt splashing phenomena (Burzynski et al. 2020).Regardless of the detailed geometrical features of the parent ligaments, when scaled with the initial collision speed U 0 , the decay of u ej over time does not significantly depend on either W e or ε 0 .We further compare the measured u ej values with the predictions of both the rim dynamic equations (4.10)-(4.12)and the scaling model (4.13).Both capture the early-time evolution of u ej up to t/τ cap ≈ 0.5, after which the decrease of u ej slightly steepens and shows a better agreement with the predictions of (4.10)-(4.12).This is most likely due to the late-time capillary deceleration not well-represented by (4.13).We note that while the rim velocity predicted by (4.10)-(4.12) is closer to u ej , they tend to over-predict the threedimensional simulation results, as already observed in fig.7; suggesting the existence of a gap between the actual rim velocity and the fragment ejection velocity.This agrees with the earlier results of Liu et al. (2022), where their fig.5 also shows fragment speeds remaining higher than the rim velocity.Similar measurements of the fragment ejection velocity have also been reported by Thoroddsen et al. (2012) for micro-splashing in drop impact problems at much larger values of W e, which were well explained by the theory of Riboux & Gordillo (2015).However, in Riboux & Gordillo (2015) the lamella rim fragments under Rayleigh-Taylor and capillary instabilities, and therefore the velocity and size of ejected droplets are directly determined by the rim.In Thoroddsen et al. (2012) fragmentation also occurs shortly after the emergence of the lamella sheet, based on which Riboux & Gordillo (2015) modelled the droplet ejection velocity using flow field information at the lamella foot.These differ from our scenario where the ejection of droplets are governed by the end-pinching of ligaments erupting on the lamella rim (Wang & Bourouiba 2018).Development of theoretical models capable of predicting fragment ejection speed in our case thus requires detailed knowledge of ligament growth and merging dynamics, which is out of the scope of the current work.
Since rim collision is a transient process where both the total number and the individual size of fragments increase over time, as shown in Fig. 13, it is of interest to determine how the fragment size and velocity distribution functions evolve with time.We first show the fragment size distributions n(r/R 0 ) in fig.17a at different times t c /τ cap , which are computed by sampling over a time window of t c − 0.45τ cap ⩽ t ⩽ t c + 0.45τ cap .To ensure statistical convergence of the data, three individual realisations are computed for each initial configuration (W e, ε, N max ); and all results presented in figs.17 and 18 have been averaged across these ensembles.It can be seen that initially fragments with r ⩽ 0.2R 0 are produced at early time.While the number density n within this range remains largely unchanged over time and does not appear to depend strongly on r, this should be treated with caution due to a lack of fully established grid independence of fragment statistics at small sizes, as discussed in §A.As time elapses, the number of larger fragments increases and causes the falling tail of the distribution to move further to the right, indicating that larger drops are fewer and produced later in time.This can be attributed to the ongoing ligament merging process, as it increases the thickness of individual ligaments.
We now seek to compare our numerical results with the experimental data and the theoretical model of Néel et al. (2020).It is noted therein that the variation in the fragment size distribution for the rim collision problem arises from two sources, namely those of the transverse ligament size and of the fragments produced from the breakup of a single ligament.A linear superposition of these two effects yields the following size distribution function, We find that setting the coefficient χ to 5 leads to an excellent match between our simulation results with (W e, ε, N max ) = (200, 0.06, 25) and the re-normalised data of Néel et al. (2020) for r ⩾ 4∆ 10 .Fig. 13b of Néel et al. (2020) suggests a χ value of approximately 25 for their controlled rim production setup, which is larger than our fitted value by a factor of 5.This might be because rim fragmentation has completed in the experiments, and the larger fragments produced at later times increases the value of d.Our numerical results differ from Néel et al. (2020) for r ⩽ 4∆ 10 , where their model (6.2) exhibits a fall-off not found in our data.This more uniform portion of fragment size distribution for small r values is also found in fig.23b of Lhuissier & Villermaux (2012), where they attributed it to the transverse impact between adjacent rim ligaments.The differences between their experimental and our numerical configurations may also play an important role, as their expanding rims feature toroidal shapes and therefore coalesce within a finite period of time.Last but not least, the two liquid rims are connected by an interstitial thin film in the configuration of Néel et al. (2020), whereas in our case the rims come into contact directly.
The dependence of the time-averaged fragment size distribution on the controlling parameters W e, ε and N max is further shown in figs.17c and 17d.It can be seen that the number density n of small fragments with r ⩽ 4∆ 10 continues to increase with W e, ε and N max , consistent with our observations in fig.16a.More specifically, here the sensitivity of n to ε and N max , as shown in fig.17d, further supports our understanding that the differences between the numerical and experimental initial conditions causes the difference between the corresponding results.The number density of intermediate fragments with 4∆ 10 ⩽ r ⩽ 0.2R 0 also increases with W e in fig.17c; but different from that of smaller fragments, it appears to asymptote to (r/R 0 ) −1/2 for W e ⩾ 240 or N max ⩾ 60.The tail of the distribution functions consisting of even larger fragments with r ⩾ 0.2R 0 remains approximately independent of the controlling parameters, and its decaying trend agrees well with a power law of (r/R 0 ) −7/2 .While this finding differs from that of Néel et al. (2020), where m and n decrease with increasing W e and cause the slope of the tail to decrease correspondingly, this may again be due to differences in detail of the initial conditions.Overall, figs.17c and 17d suggest that as W e increases, the fragment size distribution within r ⩾ 4∆ asymptotes to a regime independent of the controlling parameter, and well described by a power-law decay with a break in slope; which may originate from the insensitivity of ligament merging and breakup phenomena to the initial perturbation configurations.It is noted that a similar transition between two power law regimes has been observed in the droplet size distributions associated with wave breaking by Mostert et al. (2022); Erinin et al. (2023a), although their distributions feature steeper slopes, with the power law transitioning from r −2 to r −6 .
Here we show that the power-law scaling (r/R 0 ) −1/2 we observed in figs.17c and 17d at small fragment sizes can be derived from the ligament merging dynamics previously established in this work.The ligament merging timescale is defined as the ratio between the average ligament spacing L 0 /N lig and the drift velocity u drift , which can be evaluated based on Eq. (5.2), The rate of droplet shedding is controlled by the ligament necking process, thus ∆t shed = t neck as given in Eq. (6.1) (Wang & Bourouiba 2018).Consequently, the total number of fragments shed from all ligaments between two consecutive merging events can be estimated as where w is the average width of ligaments.
As we observed in figs.13a and 13b, w ∝ r ∝ t/τ cap , whereas Eq. (5.6) suggests N lig ∝ (t/τ cap ) −1/2 .Substituting these two scalings into Eq.(6.5) leads to These derivations suggest that the (r/R 0 ) −1/2 scaling at small fragment sizes can be explained as a competition between ligament merging and end-pinching, whereas the (r/R 0 ) −7/2 scaling found at larger fragment size ranges still awaits further analysis.It is likely that this steeper dependence on r originates from the presence of corrugations on individual ligaments, or the variation of ligament width across the bordering rim (Néel et al. 2020); neither of which has been considered in the derivations above.
Lastly, we discuss the fragment velocity statistics.Figs.18a and 18b present the ensemble-averaged vertical (fig.18a) and in-plane (fig.18b) components u y and u xz of the fragment ejection velocity, which are plotted as functions of the fragment radius r.
Since most of the initial liquid momentum is deflected in the vertical direction over the collision process, u y remains a few times to a decade larger than u xz .Fig. 18a shows that the u y values for fragments within the size range of 0.1 ⩽ r/R 0 ⩽ 0.4 collapse reasonably well when scaled by the initial rim velocity U 0 , and is well predicted by a power-law decay model u y ∝ (r/R 0 ) −1 .The velocity distribution at larger fragment sizes deviate from this power-law scaling, most likely due to a combination of late-time effects including amplified ligament corrugations, rim deceleration and the onset of RP instability on the ligaments.The in-plane velocity component values u xz measured in fig.18b are more scattered compared with fig.18a, but the distribution can still be roughly described by the same power-law decay model u xz ∝ (r/R 0 ) −1 .These observations further corroborate the fragment size distribution model (6.6) we proposed, since the power-law exponent of −1 for u y can be derived based on our observations in figs.16b and Eq.(5.1), As the only source of the liquid in-plane motion is the transverse drifting of ligaments, the in-plane velocity component u xz can be estimated by the drifting velocity u drift .Combining Eqs.(5.2) and (5.6) similarly leads to The good agreement between the velocity scaling models derived above and our numerical results once again highlights the importance of the rim ligament merging phenomenon in determining the size and velocity statistics of splashing fragments.Figs.18c and 18d show the number density of fragments as a function of their travelling speed u frag .Fig. 18c suggests that while most of the fragments produced at early time feature a skewed velocity distribution peaking at u frag ≈ 2.8U 0 , the maximum fragment speed decreases and more fragments travelling at lower speeds are recorded as time elapses, and the distribution function has developed a plateau by t/τ cap = 1.36.This suggests that as the ligaments continue to grow in length and break up, the droplets produced come to span uniformly across a large range of travelling speeds, with fewer drops featuring very slow or particularly fast speeds.Fig. 18d shows the velocity distribution at different W e values averaged over the entire collision process.It is found that as W e increases the velocity distribution becomes broader; and for W e beyond 240, the right tail of the velocity distribution appears to reach a W e-independent regime, similar to our observation in fig.17c for the size distribution of fragments; and the fastest speeds recorded is around u frag ≈ 3.2U 0 .While a direct comparison with breaking wave statistics (Mostert et al. 2022;Erinin et al. 2023a) is out of the scope of the current work, it is noted that the shapes of velocity distribution functions obtained here in figs.18c and 18d differ from their counterparts in wave breaking phenomena (see, e.g.fig.16d of Mostert et al. (2022) and 12 of Erinin et al. (2023a)), as the latter are skewed and narrower than our distributions.This may be due to the presence of gravity in wave breaking, which may arrest the ligament fragmentation process and define a short timescale for completing the splashing phenomenon, thereby reducing the total number of fragments.These splash fragments are themselves also decelerated by gravity.In the velocity distribution, this would appear as a narrowing of the distribution, with lower velocities at the peak and the higher velocities represented by a long tail.Other fragmentation mechanisms in the wave-breaking phenomena may also alter the shape of velocity distribution, e.g., the bursting of surface bubbles (Lhuissier & Villermaux 2012;Berny et al. 2021), whose contribution to the production of droplets is known to be significant after the wave splashing phase.

Conclusions
We have investigated the collision and subsequent fragmentation of perturbed liquid rims, focusing on the range of 120 ⩽ W e ⩽ 280 that allows for the growth and merging of transverse ligaments and production of fine drops from such ligaments via the pinch-off mechanism.We look into different parts of the post-collisional liquid bulk as it evolves over time, and our key findings are summarised as below: Firstly, following the quasi-one-dimensional approach of Wang & Bourouiba (2017), we derive analytical solutions of the liquid velocity and free surface profiles for the vertically expanding lamella sheet, which are shown to be in good agreement with the numerical results.Capillary effects have been neglected in this model for sheet evolution, but prove significant for the dynamics of the bordering rim.We then compared the growth of its vertical position and thickness with the model of Gordillo et al. (2019), and develop scaling laws collapsing the data.
Secondly, we analyse the behaviour of transverse ligaments on top of the lamella rim.The ligaments produce fragments primarily via the end-pinching mechanism, and when the initial perturbation waveform is polychromatic, they migrate on the rim and merge with each other to form thicker and more corrugated ligaments, thus preventing their absorption into the rim and sustaining the fragmentation process.A novel scaling model is derived for predicting the evolution of ligament number density based on the migration speed model of Wang & Bourouiba (2018).
Lastly, we present the size and velocity statistics associated with the rim collision phenomenon.An excellent agreement between our fragment size distribution and the experimental results of Néel et al. (2020) is found within the range of grid convergence (r ⩾ 4∆ 10 ).The fragment size distribution becomes insensitive to the initial configurations when W e or N max further increases, which can be described using a power law with a break in slope.A theoretical model is proposed predicting the power-law distribution observed for r ⩽ 0.2R 0 .Over time, the fragment speed u frag develops a largely uniform spread over the range of 0 ⩽ u frag ⩽ 3.2U 0 as their parent ligaments continue to grow vertically and decelerate to form slower drops.
The implications of the present work are manifold.Firstly, it sheds new light on the fluid physics involved in a liquid impact problem that has not received much attention prior to the recent works of Néel et al. (2020) and Agbaglah (2021).Furthermore, the results we obtained are also of reference value for ongoing research works on spherical drop impacts, especially the influence of initial perturbations on the fragmentation process, which may also be present during the early-time prompt splashing phenomena observed in previous experimental studies (Burzynski et al. 2020;Wang et al. 2023).Lastly, this work serves as a stepping stone towards understanding the secondary splashing phenomenon observed in wave breaking events and the associated fragment statistics (Mostert et al. 2022;Erinin et al. 2023a), and provides the basis for investigating the influence of other physical mechanisms not covered in the present work, e.g., viscosity, gravity and air-phase turbulence effects.

Declaration of Interests
The authors report no conflict of interest.4∆ Lmax can be regarded as the lower limit of fragment radius above which the fragment size and velocity distributions can be considered fully grid converged, while the statistics of fragments below this threshold are still grid-dependent and should be treated with caution (Mostert et al. 2022).Figs.19e and 19f present the fragment size and velocity distributions at a few different Oh values at t/τ cap = 1.82.Fig. 19e shows that for r ⩾ 4∆ 10 where grid convergence of fragment statistics is fully established, the fragment size distribution is not sensitive to viscous effects.On the other hand, the velocity of fragments with radii near 4∆ 10 as shown in fig.19f decreases slightly with increasing Oh values, which might be because of more significant viscous dissipation at the lamella feet and rim necks, as discussed in §5.1.
Within the same reference frame, the liquid phase velocity decays to 0 far away from the contact region, thus ∇ϕ = 0, z ≫ 1.
(B 3) The final boundary condition comes from considering the flow condition on the drop surface (but outside the contact area) for 1 ≫ y ⩾ y n (t) and z = z d (y).At very early time immediately after the impact, the rim bulks largely retain their cylindrical shape, as shown in fig.20a.Thus, the drop shape in the vicinity of the stagnation point can be approximated as up to first order in y.At high W e values capillary effects can be neglected, and the unsteady Bernoulli Equation can be applied along the drop surface, For small values of t, unsteady effects dominate and Eq.(B 5) can be further simplified using ϕ(t = 0) = 0 at the free surface, The solution of Eq. (B 1) subject to boundary conditions (B 2), (B 3) and (B 6) thus corresponds to the flow field induced by an infinitely long lamina with width 2y n moving perpendicular to its surface.The following solution written in elliptical coordinates is provided in Art.71, 3°of Lamb (1932): where y = a cosh ξ cos η, z = a sinh ξ sin η.Utilising η = 0 at z = 0, we derive the total vertical velocity u z for y > y n (t), z = 0 within the laboratory frame, Thus C = 1/2, and by integration one finds y n = 2 √ t, or written dimensionally as We conduct two-dimensional numerical simulations at a very high resolution level L max = 15 to investigate the evolution of the interface profile close to the contact region, which we show in fig.21a.It can be observed that for tU 0 /R 0 ⩾ 0.0096, a bulge on the interface profile appears at z = 0 representing the nascent lamella, whereas the local minimum in y denotes the lamella foot.Fig. 21b further compares the evolution of the lamella foot location y n measured from the numerical simulations with the theoretical prediction (B 12), where a good agreement is reached for W e ⩾ 160, indicating that the potential flow analysis employed in this section indeed captures the lamella ejection at very early time.
Solving Eqs.(4.10)-(4.12)requires the initial values of (b rim , y rim , v rim ) at the moment of lamella foot formation t e .(B 12) suggests that y rim (t e ) = 2 √ t e , v rim (t e ) = 1/ √ t e .(B 13) We also expect that the nascent rim thickness b rim (t e ) can be approximated using the lamella height function h[y rim (t e ), t e ].Thus the ejection time t e remains the only unknown in the initial conditions, although our numerical simulations suggest that it is close to zero.To regularize the initial conditions, we choose a well-resolved and sufficiently small time t ′ e = 4.5 × 10 −3 τ cap when solving Eqs.(4.10)-(4.12).We have confirmed that changing values of t ′ e does not have significant influences on the solution of Eqs.(4.10)-(4.12).

Figure 1 :
Figure 1: (a)-(c): Wave splashing observed in previous numerical (a,b) and experimental (c) studies, adapted from Wang et al. (2016) (a), Mostert et al. (2022) (b) and Erinin et al. (2023a) (c), respectively.(d): Sketch showing the ensemble-averaged breaking wave profile taken from Erinin et al. (2023b) after the initial moment of impact in the breaking wave, at the moment of secondary splashing, where dashed lines indicate our simplification of the problem as the collision of two cylindrical rims with radii r b and rs.In this study we consider the basic case rs = r b .
Figure 2: (a): Sketch showing the configuration of the liquid rim collision problem; (b): ensembleaveraged power density spectrum of the white noise signal for generating initial interface perturbations on the cylindrical rims.

Figure 5 :
Figure 5: (a): Liquid sheet velocity profile at W e = 120, scaled according to (4.5); (b): verification of (4.5) at different values of W e, time and perturbation waveforms.'Sing.' indicates that the initial perturbation we impose features a single wavenumber Nmax, while 'Sup.' denotes a combination of sinusoidal perturbations with wavelengths λ = D/8, D/16 and D/32.
Figure 7: (a)(b): The evolution of the vertical position yrim (a) and the rim thickness brim (b) over time, compared with solutions of Eqs.(4.10)-(4.12) at corresponding W e values (solid lines).Early-time measurements from two two-dimensional simulations with W e = 80 and 160 are also included.(c)(d): Results in (a) and (b) rescaled using Eq.(4.13).

Figure 8 :
Figure 8: Snapshots taken from a simulation case at W e = 200, ε = 0.06 and Nmax = 25 showing ligaments generated from the 'indentation' region between two colliding rims.From left to right: t/τcap = 0.045, 0.091 and 0.136.

Figure 11 :
Figure 11: Snapshots at W e = 120, ε = 0.06 and Nmax = 25 showing ligament evolution from monochromatic initial perturbations (upper row) and filtered white-noise perturbations (lower row).Re-absorption of ligaments back into the rim is observed for the monochromatic perturbation case after two cycles of drop shedding.From left to right: t/τcap = 0.2, 0.4 and 0.6.

Figure 12 :
Figure 12: Snapshots taken from a simulation case at W e = 160, ε = 0.04 and Nmax = 25 showing ligaments merging on the corrugated rim bordering the expanding sheet, while shedding fragments via the end-pinching mechanism.

Figure 13 :
Figure 13: (a): Measurement of the ratio between the diameter di of the detaching fragment and the width wi of its originating ligament at different ejection times.(b): The evolution of the fragment diameter d frag of ejected fragments.The results in (b) have been ensemble-averaged across three realisations for each (W e, ε) configuration, and rescaled by W e −0.75 in the main plot.

Figure 14 :
Figure 14: Main plot: Sketch showing the quantities defined in §5.2 for developing the ligament merging model (5.6).Inset: Sketch showing the local geometry of the junction region at the ligament base.

Figure 15 :
Figure 15: Evolution of the ligament number density N lig at different values of Nmax with W e = 120 (a) and different W e with Nmax = 60 (b).The insets compare the evolution of N lig with our model (5.6).

Figure 16 :
Figure16: The evolution of the total number density N frag (a) and the ejection velocity (b) of primary fragments, compared with the rim velocity urim derived from solving Eqs.(4.10)-(4.12)(solid transparent lines) and (4.13) (dash-dotted line).

Figure 17 :
Figure 17: (a): The evolution of time-and ensemble-averaged size distribution function n(r/R0) of all fragments produced by colliding rims at W e = 200, ε = 0.06, and Nmax = 25.(b): The fragment size probability distribution function f (r/R0) compared with the experimental data and model of Néel et al. (2020).(c,d): The influence of W e (c) and ε and Nmax (d) on the fragment size distribution function, where ε = 0.06 and Nmax = 25 for all simulation results presented in (c), and W e = 200 for those presented in (d).
6.2) where K m−n is a modified (m−n)-th order Bessel function of the second kind; and m and n reflect the roughness of the distribution of ligament widths and corrugation amplitudes on individual ligaments.Néel et al. (2020) fit their experimental data at W e = 193 using (6.2) with m = 40, n = 5, which we can reproduce in fig.17b together with our fragment size probability density function (pdf) at W e = 200.Note that the fragment sizes were originally normalised by Néel et al. (2020) using the average fragment diameter d in their fig.15b, which is shown in their fig.13b to saturate at large W e values and remain proportional to their interstitial sheet thickness h.The sheet thickness is in turn related to the pre-collisional rim radius R 0 via the rim collision Weber number, given in their work as W e = 16R 0 /h.This allows us to estimate their average fragment diameter as d

Figure 18 :
Figure 18: (a,b): Ensemble-averaged vertical (a) and in-plane (b) components of fragment ejection velocity uy and uxz calculated at various W e values, with Nmax = 25.(c): Evolution of the fragment velocity distribution over time, obtained at W e = 200, ε = 0.06 and Nmax = 25.(d): The probability distribution functions of fragment velocities at different W e values.

Figure 19 :
Figure 19: The size (a,c) and velocity (b,d) distributions of fragments n(r/R0) and v(r/R0) at W e = 200, ε = 0.06 and Nmax = 25 at t/τcap = 0.23 (a,b) and 1.82 (c,d), binned and averaged across three realisations with the same initial configurations at Lmax = 9, 10 and 11.Fig. 19a also includes statistics produced from a single realisation at Lmax = 12.(e,f): The size (e) and velocity (f) distributions of fragments at W e = 280 and different Oh values at t/τcap = 1.82.
Figure 20: (a): Sketch showing the cross-sectional view of the liquid cylinder collision problem.(b): Sketch showing the boundary conditions (B 2), (B 3) and (B 6) under which we solve the Laplace equation (B 1).

(Figure 21 :
Figure 21: (a): Two-dimensional simulation results at Lmax = 15 for W e = 160 showing the evolution of the contact region.(b): Comparison between simulation results at different values of W e and the theoretical prediction (B 12) at very early time.