Crown formation from a cavitating bubble close to a free surface

A rapidly growing bubble close to a free surface induces jetting: a central jet protruding outwards and a crown surrounding it at later stages. While the formation mechanism of the central jet is known and documented, that of the crown remains unsettled. We perform axisymmetric simulations of the problem using the free software program basilisk, where a finite-volume compressible solver has been implemented, that uses a geometric Volume-of-Fluid method (VoF) for the tracking of the interface. We show that the mechanism of crown formation is a combination of a pressure distortion over the curved interface, inducing flow focusing, and of a flow reversal, caused by the second expansion of the toroidal bubble that drives the crown. The work culminates in a parametric study with the Weber number, the Reynolds number, the pressure ratio and the dimensionless bubble distance to the free surface as control parameters. Their effects on both the central jet and the crown are explored. For high Weber numbers, we observe the formation of weaker"secondary crowns", highly correlated with the third oscillation cycle of the bubble.


Introduction
Free surface jetting is a widely observed phenomenon in nature with several practical applications. Of particular relevance to this paper is its occurrence in the laser-induced forward transfer (LIFT) process, a digital printing technique where parts of a donor film are transferred to a receiving substrate (Serra et al. 2004;Duocastella et al. 2007;Serra & Piqué 2019;Supponen et al. 2016;Jalaal et al. 2019a,b;Chahine 1977;Thoroddsen et al. 2009;Pain et al. 2012;Bempedelis et al. 2021). Figure 1 shows a typical sequence of events observed in a laboratory experiment related to the LIFT process: due to optical/thermal Figure 1. A pulsed laser is focused in water, close to the free surface. At t = 0 µs, a highly pressurized bubble is generated, emitting a shock wave. The bubble then aspherically expands into an egg-shape, with the tip pointing upwards (t = 60µs). This leads to an asymmetric collapse, where the bubble is pierced by an inner jet (t = 110 µs) which then breaks it into a toroid (t = 120 µs). At the same time, a jet is protruding upwards. The toroidal bubble further breaks into two toroids (t = 130 µs), and expands again while travelling downwards (t = 250 µs). Then the formation of an axisymmetric crown around the central jet is observed, see the arrow in the last snapshot. The experimental snapshots were provided by Dave Kemper (Physics of Fluids lab, University of Twente).
breakdown via focusing of a pulsed laser, a bubble nucleates, expands and collapses while a jet is ejected upwards from the free surface. In addition, an unexpected and rather remarkable phenomenon, namely the formation of an axisymmetric crown around the rim of the jet, is observed about 100 µs after the ejection of the central jet. While it is clear that the central jet is caused by the rapidly expanding bubble, the mechanism responsible for the crown observed in this as well as in other circumstances is still controversial.
Several examples of jets associated with cavitation bubbles can be found in the literature. Obreschkow et al. (2006) and Robert et al. (2007) studied laser-induced bubbles inside a cylindrical water jet close to the free surface, which may be considered as a variation of the experiment shown in figure 1 in which the free surface is initially plane. Instead of an axisymmetric crown, the authors observe two in-plane micro jets. They speculated that this might be the effect of a shock wave impinging on the curved water interface. Tagawa et al. (2012) studied the jet formation inside a capillary tube due to laser-induced cavitation. The authors highlighted the compressibility effects and argued that the jet is likely due to the impingement of the shock wave on the initially curved meniscus. Zeng et al. (2018) studied cavitation inside a droplet, leading to a corrugated free surface due to the Rayleigh-Taylor instability. The authors attribute the subsequent jetting to pressure impulses that focus the flow on suitably curved interfaces. Kiyama et al. (2016) performed tube impact tests as in Antkowiak et al. (2007), but for stronger impact. They observed cavitation in the tube and then the formation of crowns around the central jet. They argue that the interaction of expansion and compression waves with the tube wall and the curved interface results in the crown formation.
In LIFT, the focus on liquid compressibility as a determining factor in the formation of the secondary jet, i.e., the crown, is put into question by the work of Han et al. (2014) and Liu et al. (2017), who both used an incompressible Boundary Element Method (BEM) to study similar configurations. Due to the challenges that the BEM encounters in simulating the bubble after it breaks into a toroid, the former authors simply removed the toroidal bubble. The latter authors studied the interaction of two, vertically adjacent, oscillating bubbles close to a free surface but, unlike Han et al. (2014), they maintained the bubble even after it became a toroid and they observed a trace of a crown.
Other numerical studies dispensing with the assumption of liquid incompressibility Figure 2. The numerical setup for the study of a cavitation bubble in the vicinity of a free surface. The blue depicts the heavier fluid (e.g. water) while grey represents the lighter one (e.g. air), respectively referred to as fluids 1 and 2 and having an initial density ratio of ρ1/ρ2 = 1000. The shown dimensions are not to scale.
exist. Koukouvinis et al. (2016) performed Volume-of-Fluid (VoF) simulations with the aim of replicating their own experimental results. They achieved an overall good agreement, both qualitatively and quantitatively. Although their numerical results exhibit a crown, they do not discuss it nor make any comment on the mechanism of its formation or on how it depends on the control parameters. Li et al. (2019) also used a VoF method to simulate an oscillating bubble close to a free surface, also achieving a rather good agreement with their experiments. In addition, they performed a parametric study, varying the initial bubble-to-interface distance, observing bursting behaviour in some instances. On the subject of crown formation, they write "when the toroidal bubble starts to rebound, the pushing effect from the rebounding bubble, similar to the effect from the initial expanding bubble, tends to induce the upward deformation of the free surface again, which is the cause of the secondary spike". They add "since [Liu et al. (2017)] ignored the compressibility of the water, acoustic emissions were absent in their simulation. Therefore, we argue that the acoustic emissions are not important in the formation of the crown spike". Beyond these statements, however, they did not provide any quantitative evidence of the role of compressibility.
In the present work, we perform Direct Numerical Simulations (DNS) of a cavitation bubble in the vicinity of a free surface, focusing in particular, unlike the work of Koukouvinis et al. (2016) and Li et al. (2019), on the formation of the crown, and studying to what extent liquid compressibility contributes to the phenomenon. Furthermore, we perform a parametric study, analysing the effect of surface tension, viscosity, pressure ratio and bubble-to-interface initial distance on the dynamics.
The paper is organised as follows. In §2, we present the governing equations, the numerical method, the numerical setup, and the relevant dimensionless numbers. In §3, we reiterate in more detail the above sketched dynamics of the phenomenon, focusing on the origin of crown formation. In §4, we present the parametric study. Finally, we summarize our findings in §5.
2. Geometry, governing equations, control parameters, and numerical scheme The numerical setup is shown in figure 2. A bubble of initial radius R 0 is placed in a liquid at a distance H between its centre and the liquid's interface. The bubble has an initial internal pressure P b,0 larger than the surrounding pressure P ∞ . The bottom boundary has the conditions of a rigid wall, i.e. u = 0, and a zero pressure gradient in the normal direction (∂p/∂z = 0). The top and right boundaries have outflow conditions where we impose the pressure as p = P ∞ , zero normal velocity gradients and vanishing shear stresses (top: ∂v r /∂z = 0, ∂v z /∂z = 0, right: ∂v z /∂r = 0, ∂v r /∂r = 0). Gravity is neglected (g = 0) since we are studying millimetric bubbles with a short lifetime where buoyancy effects are small (Blake & Gibson 1981), as opposed to UNDEX problems where the bubble sizes are much larger and buoyant forces become dominant (Zhang et al. 2015). Hence, we use the constant atmospheric pressure P ∞ at the right boundary instead of the hydrostatic pressure. All the simulations are performed in axisymmetric configuration. The cylindrical domain is 50 times the initial radius of the bubble, both in radius and height. Its relatively large size is chosen so as to diminish the effect of boundaries on the problem, especially the reflection of pressure waves and their interaction with the physical process (Fuster & Popinet 2018).

Governing equations
The governing equations of this two-phase flow (subscript i denoting the different phases, equal to 1 in the liquid, and to 2 in the gas), in their conservative form, are written as follows: which reflect conservation of mass and momentum, respectively. In the equations above, ρ is the density and u is the velocity vector. The last term in equation (2.2) deals with the capillary forces and is written in a discrete form, where δ s is a characteristic function defined as 1 in cells containing an interface, and otherwise 0; σ is the surface tension coefficient, κ the interface curvature and n the normal to the interface. The stress tensor is where p is the pressure and µ is the shear viscosity. Note that, following Stokes' hypothesis, the bulk viscosity is neglected (Stokes 1845). In the present work, mass transfer and thermal diffusion effects are neglected. The total energy equation is then written as: where e is the internal energy of the fluid. The system is closed by an equation of state (EOS) relating the different thermodynamic properties. We use the stiffened-gas EOS, written in the Mie-Grüneisen form (Saurel & Abgrall 1999): from which the speed of sound follows as: These relations reproduce the behaviour of the gas, assumed to be diatomic and perfect, by taking Π 2 = 0 and Γ 2 = γ = 1.4, with γ the ratio of the specific heats. For the liquid, the choice Γ 1 = 5.5 and Π 1 = 492.115 MPa is empirically found to reproduce the properties of water (Johnsen & Colonius 2006). To non-dimensionalize the above governing equations, we use the initial pressure difference ∆p 0 = P b,0 − P ∞ , the initial radius of the bubble R 0 and the liquid density ρ 1 . Hence, we find an inertial timescale τ = R 0 (ρ 1 /∆p 0 ) 1/2 , and a characteristic velocity U = R 0 /τ = (∆p 0 /ρ 1 ) 1/2 . Times will therefore be multiples of τ and velocities multiples of U . Using these scales, we obtain the following dimensionless groups: which are the Reynolds, Weber and Mach numbers, respectively. The additional dimensionless groups of the problem are: which are the pressure ratio (or compression ratio), and the dimensionless initial distance of the bubble to the fluids' interface, respectively.

Numerical scheme
We use a finite-volume compressible solver for multiphase flows, introduced by Fuster & Popinet (2018). The interface has a sharp representation, being tracked by a geometric Volume-of-Fluid (VoF) method (Scardovelli & Zaleski 1999), where the different phases are labelled by a volume fraction c, equal to 1 in fluid 1, to 0 in fluid 2, and to values in between for the cells containing an interface. This volume fraction obeys the following advection equation: which advances the interface due to the local velocity field. The fluids' properties such as density and viscosity are defined as arithmetic means {ρ, µ} = c{ρ 1 , µ 1 } + (1 − c){ρ 2 , µ 2 }. This scheme has the advantage of being fully conservative due to the simultaneous advection of the conserved quantities (e.g. total energy, density, momentum) with the volume fraction c, thus avoiding any numerical diffusion between the two phases during the advection step. Another advantage of the used scheme is that it takes into account viscous and surface tension effects, the latter being modelled as continuum surface forces (CSF) in the cells containing an interface (Brackbill et al. 1992). The chosen platform for implementation is the free software program basilisk (Popinet 2015), a successor of gerris (Popinet 2003), that has a quad/octree discretization allowing both static and adaptive mesh refinement. For a complete overview of the algorithms involved, the reader should refer to the corresponding publications cited in this subsection. Boundary Element Methods (BEM), which assume an inviscid and an incompressible ambient liquid, are often used for the simulation of bubble phenomena. Some of these methods account for the liquid's weak compressibility, at low range of Mach numbers (Wang & Blake 2010. The advantage of the present scheme is that it includes viscosity, allowing to study its effect on the dynamics. Moreover, it is an all-Mach formulation, taking into account the liquid's compressibility and allowing the capture of travelling waves that can sometimes be crucial to the physical phenomena. In addition, VoF methods naturally handle the breakup of the interfaces (e.g. aspherical collapse of a bubble), while in BEM, special care is needed to reconnect the ruptured interfaces. We should note that these advantages come at a cost in terms of speed and computational power, making the DNS approach much more demanding than the BEM.

Asymmetric bubble collapse and jet formation
If potential instabilities of the spherical shape do not occur or can be neglected, a bubble in an infinite liquid domain, subject to a pressure difference with its surrounding liquid, collapses in a spherically symmetric fashion. During this collapse, the internal pressure increases and reaches its peak when the bubble is at its minimum volume. At this specific instant, the internal energy stored by the bubble also reaches a peak, and the high internal pressure generates a compression wave in the liquid which can evolve into a shock or be a shock already at the bubble interface if the collapse is sufficiently violent, i.e., if the initial pressure difference is large enough (Brenner et al. 2002;Beig et al. 2018). In LIFT, however, the bubble is not in an infinite liquid domain, but close to an air-liquid interface, leading to breakage of the spherical symmetry, and to an only axially symmetric collapsing bubble, resulting in an upward jet. Past work on LIFT, both experimental (Duocastella et al. 2010;Jalaal et al. 2019a) and numerical (Li et al. 2019;Bempedelis et al. 2021), shows that the eventual crown formation after the central jet had formed, is correlated with the second expansion of the bubble and the induced outward flow. However, since the shock wave is emitted at the same time, it is difficult to disentangle these two effects and conclude on the origin of crown formation. Therefore, to settle this question, we recur to numerical simulations, where one can easily vary the control parameters over a large span, in particular including low pressure ratios (P R) so as to weaken the compression shock and check whether a crown still forms.
Typically, in an experiment, the initial pressure inside the bubble is unknown. However, it is a key parameter in any numerical simulation, and a good prediction of its value is a prerequisite for any decent comparison between experiments and numerics. A good prediction of P b,0 is one that reproduces the experimental inflation of the bubble, or in other words one that is enough for the bubble to reach its experimental maximum equivalent radius R m , given its initial radius R 0 . If the bubble was in an infinite medium, P b,0 would be well predicted by the Rayleigh-Plesset equation or by one of its weakly compressible variants. Although the spherical symmetry is broken in our current setup, the use of such equations is still useful and greatly limits the number of trial and error iterations needed for a good prediction of P b,0 (Koukouvinis et al. 2016). The latter authors also point out the great computational cost and difficulties associated with simulating the experimental parameters. In their simulation, given an initial radius of 0.1 mm, a pressure ratio of 2180 was found to reproduce the experimental maximum equivalent radius 5.2 mm. Note that we do not aim to present one-to-one comparisons with the experiments where P R ∼ O(10 3 ). This allows us to make further compromises with the experimental parameters such as assuming a bigger R 0 , and thus pressure ratios that are orders of magnitude smaller, which were found sufficient for our current purposes. This renders the simulations more manageable. Furthermore, for higher pressure ratios, bursting behaviour was observed at the interface, such as in Li et al. (2019). Figure 3 shows a typical sequence of events produced by our numerical simulations. The snapshots show a very good qualitative agreement with their respective experimental counterparts in figure 1. Att = t/τ = 0, a spherical bubble is initialized below the free surface. In experiments, when the laser is focused, plasma forms, emitting a shock wave (Vogel et al. 1996). The latter can be seen in the first snapshot of figure 1. The emitted shock wave hits the flat interface, and a rarefaction wave is reflected, carrying negative pressure. If gaseous impurities exist between the free surface and the laser-induced bubble, and if the reflected pressure is lower than the liquid's vapour pressure, the cavitation of a cloud of tiny bubbles can occur (Blake & Gibson 1987;Ando et al. 2012). The bubbles in this cloud can burst at the free surface (Boulton-Stone & Blake 1993;Duchemin et al. 2002;Deike et al. 2018), ejecting small droplets which would ruin the print quality in LIFT (Jalaal et al. 2019b). This explains the small droplets that can be seen at the tip of the central jet in the subsequent snapshots of figure 1. In our simulations, we do not attempt to model the small secondary bubbles, nor does our numerical model allow spontaneous cavitation. In all cases, by the time the crown forms, this cloud would be long gone, pushed away by the expanding bubble, and therefore having no effect on the crown formation which is our principal focus.
The pressurized spherical bubble expands but soon loses its sphericity as its upper half gets elongated and deforms the free interface while the bottom half remains relatively hemispherical. So at the end of the expansion phase,t = 7.6, the bubble has an egglike shape, with the tip pointing upwards. Att = 16.5, an inner jet develops inside the bubble. An important quantity for the explanation of such a jet, and in the analysis of cavitation bubbles in general, is the Kelvin impulse. For details regarding its derivation and its implications on the translatory motion of the bubble, the reader is referred to Blake & Cerone (1982) and Blake (1988). The predictions of this theory fall in line with the law of Bjerknes: a bubble would move towards a rigid boundary but away from a free surface. In other words, the inner jet that pierces the bubble is directed away from a free surface, as can be seen att = 16.5. Meanwhile, the free surface is rising upwards due to the acquired inertia. This results in a stagnation point associated with a high pressure that further drives both the central and inner jets. Att = 19.2, the inner jet impacts the bottom wall of the bubble and breaks it into a toroidal structure. Att = 23, we see that the torus further breaks into two toroids, due to shearing instabilities. The bubble is nearing its minimum volume and locally, the free surface slightly moves downwards for it has to fill the void created by the collapsing bubble and its downwards migration. The same behaviour is noticed in experiments as well (t = 130 µs in figure 1). Att = 32, we observe the formation of a crown, and the onset of a Rayleigh-Plateau instability at the tip of the central jet, that will eventually lead to pinch-off of a droplet.

Crown formation
We now focus on the details of the crown formation. To that end we plot, in figure 4a, a typical example of the normalized volume of the bubble V /V 0 as function of timet = t/τ . Figure 4b shows the free surface att = 20, when no crown has formed yet. We thus focus on a region delimited by the black dashed box and record the maximum of −κ, with κ(r, z) = f rr /(1 + f 2 r ) 3/2 being the in-plane curvature of the interface z = f (r), positive in the way the latter is defined. Therefore, max(−κ) is zero, owing to the virtually flat lines at both ends of the curved interface, and rapidly grows upon the emergence of an inflection point, i.e. upon the formation of the crown. Such a case is illustrated in figure 4c att = 27. We therefore track this quantity in the region of interest, for the time interval in which the crown remains bounded by the black dashed box, and plot it in figure 4a along with the bubble volume. Figures 4b and 4c, marking pre and post crown formation states respectively, are extracted at times shown with thin dashed lines in figure 4a. The formation of the crown seems to be highly correlated with the second expansion of the bubble and happens almost at the same time. We will support this observation with more information on the pressure field.
In figure 5, we plot the life cycle of the same bubble in panel (a) and the recorded peak pressures in panel (b). The strongest shock wave occurs upon the impact of the inner jet and the breakup of the bubble into a toroidal structure (figure 5c). In general, during aspherical collapse of a bubble close to a free interface, and next to the shock wave emitted at maximum compression, a water-hammer shock is generated when the inner jet hits the wall of the bubble (Supponen et al. 2015(Supponen et al. , 2017Beig et al. 2018). In our simulations, this shock wave is evident, and although part of it is blocked by the bubble itself, it still reaches the free surface and impinges on the curved interface; however, no crown forms (figure 5c). At the moment of crown formation, when the bubble is at its minimum volume (figure 5d), we do not see the propagation of a shock wave, due to the relatively low chosen pressure ratio P R, but rather just an increase in pressure, due to the maximum compression of the bubble. Therefore, the compressibility of the ambient liquid and its ability to sustain and propagate a shock wave seem to have negligible effect in the process of crown formation. A feature we believe is crucial in the formation of the crown is the position of the interface at that moment (figure 5d). We previously mentioned that locally, the free surface dips down to fill in the void created by the collapsing bubble. Therefore, when the upper torus expands again, a flow reversal takes place and the surrounding liquid is pushed outwards due to its near incompressibility, resulting in a reversal of the near interface's curvature and direction of motion. Thus, a crown is formed due to the second expansion of the bubble. In addition, the toroids are under high pressure, owing to their maximum compression, i.e. minimum volume. The distorted pressure field around the curved interface, shown by the focusing of isobars in figure 5d, contributes to the process by focusing the flow on the curved interface (higher pressure gradient normal to the region of high curvature), but is not sufficient per se to drive the crown. It thus helps in picking the curved interface as the preferential location of crown appearance. Figure  5e shows the pressure field att = 29 where the crown has already formed. It must be stated that the other waves that we observe at later stages of the bubble life are due to the collapse of smaller bubbles. Movies of the simulation depicting both the pressure and the velocity fields can be found in the supplementary materials. Figure 6 shows snapshots of the free surface at timest ∈ [19.5−27], where the interface is coloured by the value of the gradient of pressure in the normal direction. The blue line in each of the snapshots indicates the location of zero curvature. In 3D, the curvature comprises two parts: an in-plane curvature plotted in figure 4a, and an out-of plane curvature κ = 1/R. If one looks at the central jet in the snapshots, its in-plane curvature is null, depicted by the straight line, whereas its out-of plane curvature is convex and positive. The central jet then connects to the flat free surface (z = 0), with a concave negative in-plane curvature. Therefore, an annulus exists where the two curvatures cancel out, and where, counter-intuitively, surface tension does not act, as can be inferred from equation 2.2. Att = 19.5, we see alternating gradients of pressure at distinct patches of the interface. This is due to the propagating water hammer shock, illustrated in figure 5c. As the bubble collapses, its internal pressure increases, and so does the gradient of pressure normal to the interface, which peaks aroundt = 24.5 when the bubble reaches its minimum volume. The crown seems to initiate from the zero-curvature annulus which is then shifted upwards. This might be the reason as to why one observes two in-plane micro jets instead of an axisymmetric crown in Obreschkow et al. (2006) and Robert et al. (2007), since the axisymmetry of such an annulus is broken in that configuration.
At high W e, when inertia is much stronger than the capillary forces, we also observe the onset of a "secondary crown". An example is shown in figure 7a. Figure 7b tracks the topological changes in the respective regions of interest, similarly to the dashed boxes in figures 4b and 4c. The solid orange line corresponds to the primary crown whereas the dashed orange line corresponds to its "secondary" counterpart. We clearly see that the latter also occurs due to the third expansion of the bubble. However, since the bubble rebounds become weaker with time, these topological features also become less pronounced. For instance, it takes more time to reverse the curvature of the interface and thus the increasing delay between the minima of the blue curve and the associated jumps in curvature. A movie of the simulation is included in the supplementary material. Also, counter-intuitively, the formation of a "secondary crown" is weakened by higher pressure ratios, since then the bubble tends to migrate faster downwards and away from the free surface. This is why one rarely observes these topologies in experiments. This presents further evidence as to how and why the rebound of the bubble drives the crown.

Parametric study
The dynamic processes, explained above, depend on several physical parameters. Hence, in this section, we perform a parametric study. Out of the five dimensionless numbers that we defined in section 2.2, we varied the Weber number W e, the Reynolds number Re, the initial pressure ratio P R and the dimensionless bubble distance to the free surface χ. In varying the Weber number, we study the effect of surface tension. Experimentally speaking, surface tension effects can be changed by changing the liquid or by adding surfactants, provided they instantaneously spread. The viscosity can also be readily changed by mixing water with other fluids (e.g. glycerol), and thus obtaining different Reynolds numbers Re. For the standard operating parameters of the experiments, Re ∼ O(10 4 ), so one can fairly use an inviscid flow approximation (Blake & Gibson 1981) if one's aim is to compare one-to-one with the experiments. Furthermore, the pressure ratio P R increases with the laser energy. Finally, one can vary χ by focusing the laser at different depths. In this parametric study, we perform all simulations for a fixed Mach number M a = 0.05 in the weakly compressible surrounding liquid.  Figure 8a shows the effect of the Weber number on the jet and crown formation. The larger the W e (smaller surface tension), the faster the central jet. The maximum thickness of the jet does not seem to be affected by surface tension, and, for the different values of W e, the jets thin similarly as they extend upwards. In addition, as W e increases, the crown gets more pronounced, i.e., it increasingly resembles a growing thin jet instead of a bulge of liquid. Figure 8b shows the instantaneous radial and axial velocities of the crown tips for different Weber numbers. The larger the W e, the larger the crown velocity. At small W e, surface tension forces are strong and tend to reconnect the crown with the central jet via capillary waves. This can be seen by the negative radial velocity at later times. Figure 8c depicts the trajectories followed by the crown tips. For relatively low W e, we see again the crown's tendency to rejoin the central jet (arrows pointing towards the z-axis). However, for large W e, the crown stays separate and eventually pinches off into a detached torus, as is suggested by figure 7a. An axisymmetric model cannot evidently simulate the secondary break-up of the crown into droplets by the Rayleigh-Plateau instability. A full 3D simulation, allowing variations in the out-of-plane curvature, is therefore needed to capture this phenomenon. In contrast, the development of the Rayleigh-Plateau instability on the central jet can be captured with the present axisymmetric simulations. At long times, this instability results in a spherical droplet at the tip of the central jet, the inset of which is shown in figure 8a. Figure 8d shows the instantaneous dimensionless velocity profiles V jet =ż(t) of the free surface, at r = 0, for different values of the Weber number, ranging from as low as we could numerically go, to infinity. For smallt, the free surface accelerates due to the rapid expansion of the bubble beneath it. Surface tension seems to have no effect on this first phase as all the curves for different W e collapse. The velocity then reaches a maximum value before deceleration. The higher the W e, the larger this maximum velocity. The jets eventually reach a constant speed where they rise due to the liquid's acquired inertia. At large W e, where inertial effects dominate, the jet velocity asymptotically reaches a plateau, i.e. the velocity becomes independent of the capillary effects. The curves in figure 8d also feature a second peak (att ∼ 9 for this set of parameters), but with a lower amplitude. To explain the origin of this kink, in figure 9a we plot the axial tip (r = 0) velocities of both the free surface and the upper wall of the bubble. At t ∼ 7.6, the bubble has reached its maximum expansion where its upper wall stops its upward motion (V bubble = 0) and the inner jet starts developing. This velocity (dashed curve) is plotted as absolute value, in order to stress this particular moment in time. One clearly sees that the second surge in the free surface velocity happens at the same time. This was also observed in the boundary integral simulations of Robinson et al. (2001) and Liu et al. (2017). Physically, as the bubble's tip reverses its motion, a stagnation point develops, associated with a higher static pressure, shown in figure 9b. This pressure repels the fluid axially in both directions, thus accelerating the free surface upwards, and contributing to the formation of both the central and inner jets. Note that all the key features of figure 8d, whether the peak velocity or the kink, seem to happen at the same instant in time. We thus conclude that, for the range of Weber numbers considered here, surface tension has practically no effect on the growth of the bubble which seems to expand in the same way in all the cases that we have considered. Figure 10a shows the effect of a varying Reynolds number Re on the phenomenon. With increasing Re, the jet becomes faster and so does the crown. This specific parametric study was done in the limit of an infinite Weber number so as to illustrate the competition between inertia and viscosity, with capillarity playing no role. One notices that indeed the jets no longer break into droplets due to the absence of capillary forces. This has its effect on the crown as well which in time remains a distinct topology, separated from the central jet. This is to be contrasted with figure 8a for a case of high surface tension (e.g. W e = 600) where capillary forces eventually flatten the regions of high curvature and merge the crown with the jet. Similar to figure 8a, a varying Re does not affect the maximum thicknesses of the jets. Figure 10b has the same traits as figure 8b. For low values of Re, viscous forces increasingly act against the liquid inertia, thus leading to lower jet velocities. The curves tend to converge asymptotically with increasing Re to the inviscid limit Re → ∞. Note that the respective kinks in the curves also occur at the same time, meaning that viscosity has a negligible effect on the bubble expansion for the range of Re studied here.

Dependence on the Reynolds number Re
4.3. Dependence on the pressure ratio P R Figure 11a shows the effect of the pressure ratio on the dynamics. In contrast to the Weber and Reynolds numbers, increasing P R leads to thicker jets. One must think of P R as the amount of energy initially stored inside the bubble. Hence, the larger the P R, the larger this energy, and the more the bubble is going to expand. Given the same bubble position (i.e. the same value of χ), this leads to a more pronounced egg-shape, elongated further in the direction of the free surface (figure 11b). Therefore, when the collapse phase begins, the central jet is formed as previously explained, and due to the altered initial deformation of the free surface, its thickness changes. Again, it is the first expansion phase of the bubble that dictates the shape of the central spike. In addition, we see that as PR increases, the central jet becomes faster.
An important quantity in the analysis of cavitation bubbles is the Kelvin impulse defined as: where S is the surface of the bubble, φ is the velocity potential and n is the outward normal to the fluid (i.e. into the bubble). This impulse is a measure of the liquid's linear momentum. Considering that the bubble has a source-like behaviour, and by ensuring the global conservation of linear momentum in a control volume encompassing the free surface, the value of I can be found. The Kelvin impulse describes the translatory motion of the bubble in a semi-infinite fluid, and its predictions fall in line with the law of Bjerknes, stating that the oscillating bubble migrates away from the free surface. For details regarding the derivation of the Kelvin impulse, the reader is referred to Blake & Cerone (1982), Blake & Gibson (1987), and Blake (1988). Supponen et al. (2016) derived the expression of the Kelvin impulse for the asymmetric growth and collapse of a bubble in the vicinity of a free surface, where R m is the maximum radius achieved by the bubble, and P v is the vapour pressure. With increasing P R, figure 11b shows that further inflation is achieved by the bubble, and thus higher values of R m . By inspection of equation 4.2, one readily sees that I increases with increasing R m , given that all other parameters are kept the same. This means that with increasing P R, the bubble travels further downwards, away from the free surface. Figure 11c shows the location of the bubble's centroid as a function of time, up until the end of the first oscillation cycle, for the different pressure ratios. During the expansion into the egg-shape, the centroid rises and when the collapse phase begins, the inner jet develops and the bubble migrates downwards, away from the free surface (z = 0). The distance travelled by the bubble clearly increases with increasing P R. This is why one does not observe "secondary crowns" with higher pressure ratios, as previously mentioned. Figure 11d shows the instantaneous position of the upper tip of the bubble, up until the inner jet impacts the lower wall of the then toroidal bubble. Again, one sees the expansion phase, and then the inner jet that pierces the bubble, with a velocity that increases with increasing P R.

4.4.
Dependence on the dimensionless bubble-interface distance χ Figure 12a shows the effect of the dimensionless bubble distance from the free surface, χ, on the dynamics of both jet and crown. As χ increases, the bubble retains its spherical shape as compared to an elongated egg-shape when the bubble is closer to the free surface. This is crucial to the formation of the central jet and, subsequently, the crown, as previously explained. These features will be suppressed for too large χ, when only a bump occurs on the free surface (e.g., for χ = 4.0). In order to compensate for the lack of proximity to the free surface (too large χ), one can increase the pressure ratio P R, and thus achieve a further expansion of the bubble, then again leading to an egg-shape. Figure 12b shows that χ seems to be the most sensitive parameter for which the slightest change creates a large velocity difference. All the curves exhibit the same behaviour as in figures 8d and 10b. As the distance to the free surface gets smaller, the central jet velocity becomes larger. On the other hand, from the inset 12c, one sees that as χ increases, the maximum volume achieved by the bubble expansion slightly increases as well, leading to a delay in the formation of the inner jet and thus the appearance of the kink. A bubble expands if its internal pressure P b is higher than the surrounding pressure P ∞ . As it inflates, the bubble pressure P b decreases until the point where P b = P ∞ . From this point on, the bubble continues its expansion due to the acquired inertia of the surrounding liquid. As the distance from the free surface χ increases, the bubble is submerged in a larger liquid mass, and therefore inertial effects are stronger, leading to further inflation. At maximum expansion, the bubble pressure P b is thus higher for smaller values of χ.

Conclusion
In this paper, we performed numerical simulations of a cavitating bubble close to a free surface. We observe the formation of a central jet protruding upwards while an inner jet pierces the bubble and breaks it into a toroid. We also observe the formation of a crown highly correlated with the second expansion of the bubble. Further examination of the pressure signal leads us to the conclusion that such a topological feature is not the result of shock wave impingement on the curved free surface, but rather, a combination of a pressure distortion over the curved interface, leading to a preferential selection of the location of crown formation, and of a flow reversal, caused by the second expansion of the toroidal bubble that drives this crown. For high Weber numbers, we also observe weaker "secondary crowns", correlated with the third oscillation cycle of the bubble.
We performed a parametric study, varying the Weber number W e, the Reynolds number Re, the pressure ratio P R and the dimensionless bubble distance from the free surface χ. For a given χ, the central jet's thickness seems to be only affected by the pressure ratio that dictates the shape of the expanding bubble, and thus the deformation of the free surface. All the studied parameters have a direct effect on the velocity of both the central jet and the crown. Varying χ seems to have the most effect on the formed topologies (jet and crown) which can be reduced merely to bumps provided that the bubble is far enough from the free surface. This study considered Newtonian fluids only. However, further numerical studies could provide insight into laser-induced forward transfer (LIFT) in complex fluids, i.e. viscoplastic/viscoelastic fluids, similar to what have been experimentally studied by Unger et al. (2011), Turkoz et al. (2018), and Jalaal et al. (2019a. This is most useful for engineering and biomedical applications, and could allow us to shed more light on cavitation in soft matter. Another route of relevant extension of this work is towards liquid pools with finite depths, where the bubble's interaction with the bottom wall becomes important. In this appendix, some of the simulations' numerical details will be presented. A static mesh is employed with maximum refinement around the interfaces and the regions of interest. The mesh is then progressively coarsened as the boundaries are approached. Since Non Reflective Boundary Conditions (NRBC) are not implemented in the current version of the solver, this coarsening serves as a numerical dissipation of the emitted/reflected shock/pressure waves that are mostly generated at locations where the grid is at maximum refinement. Therefore, the highly resolved and localised fronts (pressure peaks) of the waves get flattened via interpolation on the coarser grid. It is imperative that this coarsening be made progressively; otherwise, the waves will be spuriously reflected and will contaminate the regions of interest. In addition, and as mentioned in subsection 2.1, the size of the numerical domain is chosen large enough to mitigate the lack of NRBC. Figure 13 shows a couple of simulations depicting the effect of the numerical domain size where"1x domain" refers to the actual domain size employed for this paper's simulations, and "2x domain" twice that size. Figure 13a shows the shape of the interface att = 30. Both simulations are virtually identical, especially the free surface topology, i.e. the central jet and the crown. There are slight discrepancies at the level of the toroidal bubble which clearly do not interfere with the dynamics of the free surface. Figure 13b shows the instantaneous velocityż(t) of both the free surface and the upper tip of the bubble recorded at r = 0. The curves are similar to those of figure 9a, enclosing the same physics, and identical for both domain sizes. This proves that the currently used size, i.e. 50 times the initial radius of the bubble, added to the progressive mesh coarsening strategy, is safe enough to model the physical process at hand.
The maximum resolution used is ∼ 40 cells per initial bubble radius, and then progressively coarsened as previously discussed. The timestep is governed by an acoustic CFL condition based on the maximum speed of sound in the least compressible fluid, i.e. the liquid, C = max{c(r, z)} ∆t ∆x = max Γ 1 p 1 (r, z) + Π 1 ρ 1 1/2 ∆t ∆x = 0.5. (A 1) This condition guarantees that the pressure waves, that travel at the speed of sound, are well captured.