1. 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 (Chahine Reference Chahine1977; Serra et al. Reference Serra, Colina, Fernández-Pradas, Sevilla and Morenza2004; Duocastella et al. Reference Duocastella, Colina, Fernández-Pradas, Serra and Morenza2007; Thoroddsen et al. Reference Thoroddsen, Takehara, Etoh and Ohl2009; Pain et al. Reference Pain, Hui Terence Goh, Klaseboer, Ohl and Cheong Khoo2012; Supponen et al. Reference Supponen, Obreschkow, Tinguely, Kobel, Dorsaz and Farhat2016; Jalaal et al. Reference Jalaal, Klein Schaarsberg, Visser and Lohse2019a,Reference Jalaal, Li, Klein Schaarsberg, Qin and Lohseb; Serra & Piqué Reference Serra and Piqué2019; Bempedelis et al. Reference Bempedelis, Zhou, Andersson and Ventikos2021). Figure 1 shows a typical sequence of events observed in a laboratory experiment related to the LIFT process: due to optical/thermal 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 approximately 
$100\ \mathrm {\mu }\textrm {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.
Figure 1. A pulsed laser is focused in water, close to the free surface. At 
$t = 0 \ \mathrm {\mu }\textrm {s}$, a highly pressurized bubble is generated, emitting a shock wave. The bubble then expands aspherically into an egg shape, with the tip pointing upwards (
$t = 60 \ \mathrm {\mu }\textrm {s}$). This leads to an asymmetric collapse, where the bubble is pierced by an inner jet (
$t = 110 \ \mathrm {\mu }\textrm {s}$), which then breaks it into a toroid (
$t = 120 \ \mathrm {\mu }\textrm {s}$). At the same time, a jet is protruding upwards. The toroidal bubble further breaks into two toroids (
$t = 130 \ \mathrm {\mu }\textrm {s}$), and expands again while travelling downwards (
$t = 250 \ \mathrm {\mu }\textrm {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 laboratory, University of Twente.)
Several examples of jets associated with cavitation bubbles can be found in the literature. Obreschkow et al. (Reference Obreschkow, Kobel, Dorsaz, De Bosset, Nicollier and Farhat2006) and Robert et al. (Reference Robert, Lettry, Farhat, Monkewitz and Avellan2007) 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 planar. Instead of an axisymmetric crown, the authors observed 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. (Reference Tagawa, Oudalov, Visser, Peters, van der Meer, Sun, Prosperetti and Lohse2012) studied jet formation inside a capillary tube due to laser-induced cavitation. The authors highlighted the compressibility effects and argued that the jet is probably due to the impingement of the shock wave on the initially curved meniscus. Zeng et al. (Reference Zeng, Gonzalez-Avila, Voorde and Ohl2018) studied cavitation inside a droplet, leading to a corrugated free surface due to the Rayleigh–Taylor instability. The authors attributed the subsequent jetting to pressure impulses that focus the flow on suitably curved interfaces. Kiyama et al. (Reference Kiyama, Tagawa, Ando and Kameda2016) performed tube impact tests as in Antkowiak et al. (Reference Antkowiak, Bremond, Le Dizès and Villermaux2007), but for stronger impact. They observed cavitation in the tube and then the formation of crowns around the central jet. They argued that the interaction of expansion and compression waves with the tube wall and the curved interface results in 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, Zhang & Li (Reference Han, Zhang and Li2014) and Liu et al. (Reference Liu, Cui, Ren and Zhang2017), who both used an incompressible boundary element method (BEM) to study similar configurations. Owing 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. (Reference Han, Zhang and Li2014), 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 exist. Koukouvinis et al. (Reference Koukouvinis, Gavaises, Supponen and Farhat2016) 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 did 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. (Reference Li, Zhang, Wang, Li and Liu2019) 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. (Reference Liu, Cui, Ren and Zhang2017)] 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. (Reference Koukouvinis, Gavaises, Supponen and Farhat2016) and Li et al. (Reference Li, Zhang, Wang, Li and Liu2019), 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-interface initial distance on the dynamics.
The paper is organized as follows. In § 2, we present the governing equations, the numerical method, the numerical set-up 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
2.1. Set-up
The numerical set-up 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}$ higher than the surrounding pressure 
$P_\infty$. The bottom boundary has the conditions of a rigid wall, i.e. 
$\boldsymbol {u}=0$, and a zero pressure gradient in the normal direction (
$\partial p/\partial z = 0$). The top and right boundaries have outflow conditions where we impose the pressure as 
$p = P_\infty$, zero normal velocity gradients and vanishing shear stresses (top, 
$\partial v_r/\partial z = 0$, 
$\partial v_z/\partial z = 0$; right, 
$\partial v_z/\partial r = 0$, 
$\partial v_r/\partial r = 0$). Gravity is neglected (
$g = 0$) since we are studying millimetric bubbles with a short lifetime where buoyancy effects are small (Blake & Gibson Reference Blake and Gibson1981), as opposed to UNDEX (underwater explosion) problems where the bubble sizes are much larger and buoyant forces become dominant (Zhang et al. Reference Zhang, Cui, Cui and Wang2015). Hence, we use the constant atmospheric pressure 
$P_\infty$ 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, in both 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 Reference Fuster and Popinet2018).
Figure 2. The numerical set-up 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 
$\rho _1/\rho _2 = 1000$. The dimensions shown are not to scale.
2.2. Governing equations
The governing equations of this two-phase flow (subscript 
$i$ denoting the different phases, 1 in the liquid and 2 in the gas), in their conservative form, are written as
\begin{gather} \frac{\partial\rho_i}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho_i\boldsymbol{u}_i) = 0, \end{gather}
\begin{gather} \frac{\partial\rho_i\boldsymbol{u}_i}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho_i\boldsymbol{u}_i\boldsymbol{u}_i) = \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\mathsf{T}}_i + \sigma\kappa\delta_s\boldsymbol{n}, \end{gather}
which reflect conservation of mass and momentum, respectively. In the equations above, 
$\rho$ is the density and 
$\boldsymbol {u}$ is the velocity vector. The last term in (2.2) deals with the capillary forces and is written in a discrete form, where 
$\delta _s$ is a characteristic function defined as 
$1$ in cells containing an interface and 
$0$ otherwise; 
$\sigma$ is the surface tension coefficient; 
$\kappa$ is the interface curvature; and 
$\boldsymbol {n}$ is the normal to the interface. The stress tensor is
\begin{equation} \boldsymbol{\mathsf{T}}_i ={-}(p_i + \tfrac{2}{3}\mu_i\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_i)\boldsymbol{\mathsf{I}} + \mu_i(\boldsymbol{\nabla}\boldsymbol{u}_i + \boldsymbol{\nabla}\boldsymbol{u}_i{}^\textrm{T}), \end{equation}
where 
$p$ is the pressure and 
$\mu$ is the shear viscosity. Note that, following Stokes’ hypothesis, the bulk viscosity is neglected (Stokes Reference Stokes1845). In the present work, mass transfer and thermal diffusion effects are neglected. The total energy equation is then written as
\begin{equation} \frac{\partial}{\partial t}[\rho_i(e_i + \tfrac{1}{2}\boldsymbol{u}_i^{2})] + \boldsymbol{\nabla}\boldsymbol{\cdot}[\rho_i(e_i + \tfrac{1}{2}\boldsymbol{u}_i^{2})\boldsymbol{u}_i] = \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\mathsf{T}}_i\boldsymbol{\cdot}\boldsymbol{u}_i), \end{equation}
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 Reference Saurel and Abgrall1999),
\begin{equation} \rho_i e_i = \frac{p_i + \varGamma_i\varPi_i}{\varGamma_i - 1}, \end{equation}from which the speed of sound follows as
\begin{equation} c_i = \left[\varGamma_i\left(\frac{p_i + \varPi_i}{\rho_i}\right)\right]^{1/2}. \end{equation}
These relations reproduce the behaviour of the gas, assumed to be diatomic and perfect, by taking 
$\varPi _2=0$ and 
$\varGamma _2=\gamma = 1.4$, with 
$\gamma$ the ratio of the specific heats. For the liquid, the choice 
$\varGamma _1=5.5$ and 
$\varPi _1= 492.115$ MPa is empirically found to reproduce the properties of water (Johnsen & Colonius Reference Johnsen and Colonius2006).
To non-dimensionalize the above governing equations, we use the initial pressure difference 
$\Delta p_0 = P_{b,0} - P_\infty$, the initial radius of the bubble 
$R_0$ and the liquid density 
$\rho _1$. Hence, we find an inertial time scale 
$\tau = R_0(\rho _1/\Delta p_0)^{1/2}$ and a characteristic velocity 
$U = R_0/\tau = (\Delta p_0/\rho _1)^{1/2}$. Times will therefore be multiples of 
$\tau$ and velocities multiples of 
$U$. Using these scales, we obtain the following dimensionless groups:
\begin{equation} {Re} = \frac{R_0 (\rho_1\Delta p_0)^{1/2}}{\mu_1}, \quad We = \frac{\Delta p_0R_0}{\sigma}, \quad Ma = \left(\frac{\Delta p_0}{\rho_1c_1^{2}}\right)^{1/2}, \end{equation}which are the Reynolds, Weber and Mach numbers, respectively. The additional dimensionless groups of the problem are
\begin{equation} PR = \frac{P_{b,0}}{P_\infty}, \quad \chi = \frac{H}{R_0}, \end{equation}which are the pressure ratio (or compression ratio) and the dimensionless initial distance of the bubble to the fluids’ interface, respectively.
2.3. Numerical scheme
We use a finite-volume compressible solver for multiphase flows, introduced by Fuster & Popinet (Reference Fuster and Popinet2018). The interface has a sharp representation, being tracked by a geometric VoF method (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999), 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:
\begin{equation} \frac{\partial c}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(c\boldsymbol{u}) = c\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}, \end{equation}
which advances the interface due to the local velocity field. The fluids’ properties such as density and viscosity are defined as arithmetic means 
$\{\rho ,\mu \}=c\{\rho _1,\mu _1\} + (1-c)\{\rho _2,\mu _2\}$. This scheme has the advantage of being fully conservative due to the simultaneous advection of the conserved quantities (e.g. density, momentum and total energy) 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, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992). The chosen platform for implementation is the free software program basilisk (Popinet Reference Popinet2015), a successor of gerris (Popinet Reference Popinet2003), which 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.
BEMs, 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 Reference Wang and Blake2010, Reference Wang and Blake2011). The advantage of the present scheme is that it includes viscosity, allowing its effect on the dynamics to be studied. 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 BEMs, 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 BEMs.
3. Phenomenology
3.1. 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, Hilgenfeldt & Lohse Reference Brenner, Hilgenfeldt and Lohse2002; Beig, Aboulhasanzadeh & Johnsen Reference Beig, Aboulhasanzadeh and Johnsen2018). 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. Reference Duocastella, Patrascioiu, Fernández-Pradas, Morenza and Serra2010; Jalaal et al. Reference Jalaal, Klein Schaarsberg, Visser and Lohse2019a) and numerical (Li et al. Reference Li, Zhang, Wang, Li and Liu2019; Bempedelis et al. Reference Bempedelis, Zhou, Andersson and Ventikos2021), shows that the eventual crown formation after the central jet has 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 come to any conclusions on the origin of crown formation. Therefore, to settle this question, we turn to numerical simulations, where one can easily vary the control parameters over a large span, in particular including low pressure ratios (
$PR$) 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 set-up, 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. Reference Koukouvinis, Gavaises, Supponen and Farhat2016). 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\ \textrm{mm}$, a pressure ratio of 
$2180$ was found to reproduce the experimental maximum equivalent radius 
$5.2\ \textrm {mm}$. Note that we do not aim to present one-to-one comparisons with the experiments where 
$PR \sim 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. (Reference Li, Zhang, Wang, Li and Liu2019).
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. At 
$\tilde {t} = t/\tau = 0$, a spherical bubble is initialized below the free surface. In experiments, when the laser is focused, plasma forms, emitting a shock wave (Vogel, Busch & Parlitz Reference Vogel, Busch and Parlitz1996). 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 Reference Blake and Gibson1987; Ando, Liu & Ohl Reference Ando, Liu and Ohl2012). The bubbles in this cloud can burst at the free surface (Boulton-Stone & Blake Reference Boulton-Stone and Blake1993; Duchemin et al. Reference Duchemin, Popinet, Josserand and Zaleski2002; Deike et al. Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Séon2018), ejecting small droplets which would ruin the print quality in LIFT (Jalaal et al. Reference Jalaal, Li, Klein Schaarsberg, Qin and Lohse2019b). 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.
Figure 3. A numerical sequence of events due to a pressurized bubble in the vicinity of the free surface. For each instant, the respective snapshot is divided into two halves: on the left, the magnitude of the vertical velocity in the liquid is shown as a colour code; whereas on the right, the magnitude of the radial velocity is presented. The red colour means positive velocity, in the direction of the defined arrows; while blue stands for negative velocities; and stagnation points are shown in white. For this simulation, 
${Re}\rightarrow \infty$, 
$We = 1000$, 
$Ma = 0.05$, 
$PR = 20$ and 
$\chi = 2$.
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, 
$\tilde {t} = 7.6$, the bubble has an egg-like shape, with the tip pointing upwards. At 
$\tilde {t} = 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 (Reference Blake and Cerone1982) and Blake (Reference Blake1988). 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 at 
$\tilde {t} = 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. At 
$\tilde {t} = 19.2$, the inner jet impacts the bottom wall of the bubble and breaks it into a toroidal structure. At 
$\tilde {t} = 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 \ \mathrm {\mu }\textrm {s}$ in figure 1). At 
$\tilde {t} = 32$, we observe the formation of a crown, and the onset of a Rayleigh–Plateau instability at the tip of the central jet, which will eventually lead to pinch-off of a droplet.
3.2. Crown formation
We now focus on the details of the crown formation. To that end we plot, in figure 4(a), a typical example of the normalized volume of the bubble 
$V/V_0$ as a function of time 
$\tilde {t} = t/\tau$. Figure 4(b) shows the free surface at 
$\tilde {t} = 20$, when no crown has formed yet. We thus focus on a region delimited by the black dashed box and record the maximum of 
$-\kappa$, with 
$\kappa (r,z) = f_{rr}/(1+f_r^{2})^{3/2}$ being the in-plane curvature of the interface 
$z = f(r)$, positive in the way the latter is defined. Therefore, 
$\max (-\kappa )$ 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 4(c) at 
$\tilde {t} = 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 4(a) along with the bubble volume. Figures 4(b) and 4(c), marking pre- and post-crown formation states, respectively, are extracted at times shown with thin dashed lines in figure 4(a). 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.
Figure 4. (a) Relative bubble volume 
$V/V_0$ (blue) and in-plane curvature 
$\kappa$ (orange) as functions of time. The thick dashed line marks the instant of crown formation, where the value of 
$\max (-\kappa )$ suddenly increases. (b) Free surface at 
$\tilde {t} = 20$. The dashed black box shows the location where the curvature 
$\kappa$ is recorded, from which 
$\max (-\kappa )$ is determined. (c) Free surface at 
$\tilde {t} = 27$, when the crown has formed. The two instants 
$\tilde {t} = 20$ and 
$\tilde {t} = 27$ of (b) and (c) are marked with thin dashed lines in (a). 
$\mathcal{R}$ and 
$\mathcal{Z}$ are the dimensionless coordinates. For this simulation, 
${Re} = 2000$, 
$We = 1000$, 
$Ma = 0.05$, 
$PR = 20$ and 
$\chi = 2$.
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. Reference Supponen, Kobel, Obreschkow and Farhat2015, Reference Supponen, Obreschkow, Kobel, Tinguely, Dorsaz and Farhat2017; Beig et al. Reference Beig, Aboulhasanzadeh and Johnsen2018). 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 
$PR$, 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 a negligible effect in the process of crown formation.
Figure 5. (a) Relative bubble volume 
$V/V_0$ as a function of time. (b) Instantaneous pressure signals. The red curve shows the maximum pressure whereas the blue curve records the minimum. (c–e) Free surfaces and bubble shapes at 
$\tilde {t} = 19.8$, 
$\tilde {t} = 24.3$ and 
$\tilde {t} = 31$, respectively. In (d), pressure isobars are plotted as well. The solid lines with arrows help in marking the respective instants in time. The colour code is for the pressure field. For this simulation, 
${Re} = 2000$, 
$We = 1000$, 
$Ma = 0.05$, 
$PR = 20$ and 
$\chi = 2$.
A feature that 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 5(d), 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 5(e) shows the pressure field at 
$\tilde {t} = 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 available at https://doi.org/10.1017/jfm.2021.676.
Figure 6 shows snapshots of the free surface at times 
$\tilde {t}\in [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 three dimensions, the curvature comprises two parts: an in-plane curvature plotted in figure 4(a) and an out-of plane curvature 
$\kappa = 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 (2.2). At 
$\tilde {t} = 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 5(c). As the bubble collapses, its internal pressure increases, and so does the gradient of pressure normal to the interface, which peaks around 
$\tilde {t} =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 why two in-plane micro-jets instead of an axisymmetric crown were observed in Obreschkow et al. (Reference Obreschkow, Kobel, Dorsaz, De Bosset, Nicollier and Farhat2006) and Robert et al. (Reference Robert, Lettry, Farhat, Monkewitz and Avellan2007), since the axisymmetry of such an annulus is broken in that configuration.
Figure 6. Shapes of the free surface at 
$\tilde {t} = 19.5,20.5,21.5,22.5,23.5,24.5,25.5,27.0$, from (a) to (h), respectively. The interface is coloured by the value of the normal gradient of pressure. In each snapshot, a blue line indicates the location of zero curvature. For this simulation, 
${Re} = 2000$, 
$We = 1000$, 
$Ma = 0.05$, 
$PR = 20$ and 
$\chi = 2$.
At high 
$We$, when inertia is much stronger than the capillary forces, we also observe the onset of a ‘secondary crown’. An example is shown in figure 7(a). Figure 7(b) tracks the topological changes in the respective regions of interest, similarly to the dashed boxes in figures 4(b) and 4(c). 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.
Figure 7. (a) Free surface at 
$\tilde {t} = 40$. (b) The bubble volume as a function of time (blue curve) and the value of 
$\max (-\kappa )$ in the region where the crown is formed (solid orange curve). The dashed orange line records it in the region where the ‘secondary crown’ is formed. For this simulation, 
${Re}\rightarrow \infty$, 
$We = 4000$, 
$Ma = 0.05$, 
$PR = 10$ and 
$\chi = 2$.
4. 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 § 2.2, we varied the Weber number 
$We$, the Reynolds number 
${Re}$, the initial pressure ratio 
$PR$ and the dimensionless bubble distance to the free surface 
$\chi$. 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}\sim O(10^{4})$, so one can fairly use an inviscid flow approximation (Blake & Gibson Reference Blake and Gibson1981) if one's aim is to compare one-to-one with the experiments. Furthermore, the pressure ratio 
$PR$ increases with the laser energy. Finally, one can vary 
$\chi$ by focusing the laser at different depths. In this parametric study, we perform all simulations for a fixed Mach number 
$Ma = 0.05$ in the weakly compressible surrounding liquid.
4.1. Dependence on the Weber number 
$We$
Figure 8(a) shows the effect of the Weber number on the jet and crown formation. The larger the 
$We$ (the smaller the 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 
$We$, the jets thin similarly as they extend upwards. In addition, as 
$We$ increases, the crown gets more pronounced, i.e. it increasingly resembles a growing thin jet instead of a bulge of liquid. Figure 8(b) shows the instantaneous radial and axial velocities of the crown tips for different Weber numbers. The larger the 
$We$, the larger the crown velocity. At small 
$We$, 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 8(c) depicts the trajectories followed by the crown tips. For relatively low 
$We$, we see again the crown's tendency to rejoin the central jet (arrows pointing towards the 
$z$-axis). However, for large 
$We$, the crown stays separate and eventually pinches off into a detached torus, as is suggested by figure 7(a). An axisymmetric model cannot evidently simulate the secondary breakup of the crown into droplets by the Rayleigh–Plateau instability. A full three-dimensional 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, which is shown in figure 8(a). Figure 8(d) shows the instantaneous dimensionless velocity profiles 
$\mathcal {V}_{jet} = \dot {z}(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 small 
$\tilde {t}$, 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 
$We$ collapse. The velocity then reaches a maximum value before deceleration. The higher the 
$We$, the larger this maximum velocity. The jets eventually reach a constant speed where they rise due to the liquid's acquired inertia. At large 
$We$, where inertial effects dominate, the jet velocity asymptotically reaches a plateau, i.e. the velocity becomes independent of the capillary effects.
Figure 8. (a) Free surface shape at 
$\tilde {t} = 40$ for different Weber numbers 
$We$. (b) Instantaneous radial (dashed) and axial (solid) velocities of the crown's tip for different 
$We$. (c) Trajectories followed by the crown's tip for different 
$We$. The dots mark the beginning of the time interval whereas the arrows mark its end as well as the direction of motion. (d) Instantaneous central jet velocity measured from its tip, for different 
$We$. For these simulations, 
${Re}\rightarrow \infty$, 
$Ma = 0.05$, 
$PR = 20$ and 
$\chi = 2$.
The curves in figure 8(d) also feature a second peak (at 
$\tilde {t}\sim 9$ for this set of parameters), but with a lower amplitude. To explain the origin of this kink, in figure 9(a) we plot the axial tip (
$r = 0$) velocities of both the free surface and the upper wall of the bubble. At 
$\tilde {t}\sim 7.6$, the bubble has reached its maximum expansion where its upper wall stops its upward motion (
$\mathcal {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. (Reference Robinson, Blake, Kodama, Shima and Tomita2001) and Liu et al. (Reference Liu, Cui, Ren and Zhang2017). Physically, as the bubble's tip reverses its motion, a stagnation point develops, associated with a higher static pressure, shown in figure 9(b). 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 8(d), 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 9. (a) The instantaneous jet velocity (solid line) measured from the motion of its tip; and the velocity of the upper tip of the bubble (dashed line), in absolute value. (b) Interfaces at 
$\tilde {t} = 8$. The colour code is for the pressure field. For this simulation, 
${Re}\rightarrow \infty$, 
$We = 600$, 
$Ma = 0.05$, 
$PR = 20$ and 
$\chi = 2$.
4.2. Dependence on the Reynolds number ${Re}$
Figure 10(a) 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 8(a) for a case of high surface tension (e.g. 
$We = 600$) where capillary forces eventually flatten the regions of high curvature and merge the crown with the jet. Similar to figure 8(a), a varying 
${Re}$ does not affect the maximum thickness of the jets. Figure 10(b) has the same traits as figure 8(d). 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}\rightarrow \infty$. 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.
Figure 10. (a) Free surface shape at 
$\tilde {t} = 50$ for three different Reynolds numbers 
${Re}$. (b) Instantaneous central jet velocity measured from its tip, for different 
${Re}$. For these simulations, 
$We\rightarrow \infty$, 
$Ma = 0.05$, 
$PR = 20$ and 
$\chi = 2$.
4.3. Dependence on the pressure ratio $PR$
Figure 11(a) shows the effect of the pressure ratio on the dynamics. In contrast to the Weber and Reynolds numbers, increasing 
$PR$ leads to thicker jets. One must think of 
$PR$ as the amount of energy initially stored inside the bubble. Hence, the larger the 
$PR$, the larger this energy, and the more the bubble is going to expand. Given the same bubble position (i.e. the same value of 
$\chi$), 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.
Figure 11. (a) The free interface for different values of the pressure ratio 
$PR$ shown at analogous instants in time. (b) Egg-shape dependence on the pressure ratio 
$PR$. (c) Location of the bubble's centroid as a function of time, for the different pressure ratios, up until the respective ends of the first oscillation cycles, where the bubble starts its second expansion. (d) Position of the upper tip of the bubble as a function of time, for the different pressure ratios, up until the respective moments of inner jet impacts leading to the bubble's metamorphosis into a toroidal shape. For these simulations, 
${Re}\rightarrow \infty$, 
$We = 1000$, 
$Ma = 0.05$ and 
$\chi = 2$.
An important quantity in the analysis of cavitation bubbles is the Kelvin impulse defined as
\begin{equation} \boldsymbol{I} = \rho\int_S\phi\boldsymbol{n}\,\textrm{d}S, \end{equation}
where 
$S$ is the surface of the bubble, 
$\phi$ is the velocity potential and 
$\boldsymbol {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 
$\boldsymbol{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 (Reference Blake and Cerone1982), Blake & Gibson (Reference Blake and Gibson1987) and Blake (Reference Blake1988). Supponen et al. (Reference Supponen, Obreschkow, Tinguely, Kobel, Dorsaz and Farhat2016) derived the expression of the Kelvin impulse for the asymmetric growth and collapse of a bubble in the vicinity of a free surface,
\begin{equation} \boldsymbol{I} = 0.934R_m^{3}\sqrt{\rho_l(P_\infty - P_v)}\left(\frac{H}{R_m}\right)^{{-}2}\boldsymbol{n}\simeq 0.934R_m^{3}\sqrt{\rho_lP_\infty}\left(\frac{H}{R_m}\right)^{{-}2}\boldsymbol{n}, \end{equation}
where 
$R_m$ is the maximum radius achieved by the bubble and 
$P_v$ is the vapour pressure. With increasing 
$PR$, figure 11(b) shows that further inflation is achieved by the bubble, and thus higher values of 
$R_m$. By inspection of (4.2), one readily sees that 
$\boldsymbol{I}$ increases with increasing 
$R_m$, given that all other parameters are kept the same. This means that, with increasing 
$PR$, the bubble travels further downwards, away from the free surface. Figure 11(c) 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 
$PR$. This is why one does not observe ‘secondary crowns’ with higher pressure ratios, as previously mentioned. Figure 11(d) 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 
$PR$.
4.4. Dependence on the dimensionless bubble–interface distance $\chi$
Figure 12(a) shows the effect of the dimensionless bubble distance from the free surface, 
$\chi$, on the dynamics of both jet and crown. As 
$\chi$ 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 
$\chi$, when only a bump occurs on the free surface (e.g. for 
$\chi = 4.0$). In order to compensate for the lack of proximity to the free surface (too large 
$\chi$), one can increase the pressure ratio 
$PR$, and thus achieve a further expansion of the bubble, then again leading to an egg shape. Figure 12(b) shows that 
$\chi$ 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 8(d) and 10(b). As the distance to the free surface gets smaller, the central jet velocity becomes larger. On the other hand, from the inset figure 12(c), one sees that, as 
$\chi$ 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_\infty$. As it inflates, the bubble pressure 
$P_b$ decreases until the point where 
$P_b = P_\infty$. 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 
$\chi$ 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 
$\chi$.
Figure 12. (a) Interfaces at 
$\tilde {t} = 35$ for three different dimensionless bubble–interface distances, 
$\chi$. (b) Instantaneous jet velocity for different values of 
$\chi$. (c) The inset shows the maximum volume reached by the expanding bubble as a function of 
$\chi$. For these simulations, 
${Re}\rightarrow \infty$, 
$We = 1000$, 
$Ma = 0.05$ and 
$PR = 20$.
5. 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 
$We$, the Reynolds number 
${Re}$, the pressure ratio 
$PR$ and the dimensionless bubble distance from the free surface 
$\chi$. For a given 
$\chi$, the central jet's thickness seems to be affected only by the pressure ratio, which 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 
$\chi$ 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 LIFT in complex fluids, i.e. viscoplastic or viscoelastic fluids, similar to what have been experimentally studied by Unger et al. (Reference Unger, Gruene, Koch, Koch and Chichkov2011), Turkoz et al. (Reference Turkoz, Perazzo, Kim, Stone and Arnold2018) and Jalaal et al. (Reference Jalaal, Klein Schaarsberg, Visser and Lohse2019a). 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.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2021.676.
Acknowledgements
The authors would like to thank D. Fuster for stimulating discussions, and D. Kemper for providing the experimental images for this paper.
Funding
The research leading to these results has received funding from the European Union's Horizon 2020 Research and Innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 813766. This work was carried out on the national e-infrastructure of SURFsara, a subsidiary of the SURF cooperation, the collaborative ICT organization for Dutch education and research.
Declaration of interest
The authors report no conflict of interest.
Appendix A. Numerical details
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 (NRBCs) are not implemented in the current version of the solver, this coarsening serves as a numerical dissipation of the emitted or reflected shock/pressure waves that are mostly generated at locations where the grid is at maximum refinement. Therefore, the highly resolved and localized 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 § 2.1, the size of the numerical domain is chosen large enough to mitigate the lack of NRBCs. Figure 13 shows a couple of simulations depicting the effect of the numerical domain size where ‘
$1{\times }$ domain’ refers to the actual domain size employed for this paper's simulations, and ‘
$2{\times }$ domain’ twice that size. Figure 13(a) shows the shape of the interface at 
$\tilde {t} = 30$. The two 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 13(b) shows the instantaneous velocity 
$\dot {z}(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 9(a), 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.
Figure 13. (a) Shape of the interface at 
$\tilde {t} = 30$ for two different domain sizes. (b) Instantaneous velocities of the free surface jet, and of the upper tip of the bubble, the latter shown in absolute value, for two different domain sizes. For both simulations, 
${Re} = 2000$, 
$We = 1000$, 
$Ma = 0.05$, 
$PR = 20$ and 
$\chi = 2$.
The maximum resolution used is 
$\sim 40$ cells per initial bubble radius, and then progressively coarsened as previously discussed. The time step is governed by an acoustic Courant–Friedrichs–Lewy (CFL) condition based on the maximum speed of sound in the least compressible fluid, i.e. the liquid,
\begin{equation} C = \max\{c(r,z)\}\frac{\Delta t}{\Delta x} = \max\left\{\left[\varGamma_1\left(\frac{p_1(r,z) + \varPi_1}{\rho_1}\right)\right]^{1/2}\right\}\frac{\Delta t}{\Delta x} = 0.5. \end{equation}This condition guarantees that the pressure waves, which travel at the speed of sound, are well captured.
Open access
