Dimples, jets and self-similarity in nonlinear capillary waves

Abstract

Numerical studies of dimple and jet formation from a collapsing cavity often model the initial cavity shape as a truncated sphere, mimicking a bursting bubble. In this study, we present a minimal model containing only nonlinear inertial and capillary forces, which produces dimples and jets from a collapsing capillary wave trough. The trough in our simulation develops from a smooth initial perturbation, chosen to be an eigenmode to the linearised ${O}(\epsilon )$ problem ($\epsilon$ is the non-dimensional amplitude). We explain the physical mechanism of dimple formation and demonstrate that, for moderate $\epsilon$, the sharp dimple seen in simulations is well captured by the weakly nonlinear ${O}(\epsilon ^3)$ theory developed here. For $\epsilon \gg 1$ the regime is strongly nonlinear, spreading surface energy into many modes, and the precursor dimple now develops into a sharply rising jet. Here, simulations reveal a novel localised window (in space and time) where the jet evolves self-similarly following inviscid (Keller & Miksis, SIAM J. Appl. Maths, vol. 43, issue 2, 1983, pp. 268–277) scales. We develop an analogy of this regime to a self-similar solution of the first kind, for linearised capillary waves. Our first-principles study demonstrates that, at sufficiently small scales, dimples and jets can form from radial inward focusing of capillary waves, and the formation of this may be described by a relatively simple model employing (nonlinear) inertial and capillary effects. Viscosity and gravity can, however, significantly influence the focusing process, either intensifying the singularity or weakening it (Walls et al., Phys. Rev. E, vol. 92, issue 2, 2015, 021002; Gordillo & Rodríguez-Rodríguez, J. Fluid Mech., vol. 867, 2019, pp. 556–571). This leads, in particular, to critical values of Ohnesorge and Bond numbers, which cannot be obtained from our minimal model.

1. Introduction

Nonlinear surface waves, ranging in wavelengths from metres to micrometres, remain a bountiful source of curious observations. A phenomenon of sustained attention has been jet formation during cavity collapse (Kientzler et al. 1954). Owing to its relevance to applications such as the generation of sea spray from bursting oceanic bubbles (Blanchard 1963) or the spread of aroma (Séon & Liger-Belair 2017), research on this has been sustained (Blanchard 1963; MacIntyre 1972) and intense (Zeff et al. 2000; Duchemin et al. 2002; Bolaños-Jiménez et al. 2008; Gordillo 2008; Gañán-Calvo 2017; Lai, Eggers & Deike 2018; Gordillo & Rodríguez-Rodríguez 2019), continuing actively into the present (Ji, Yang & Feng 2021). The event appears quite ubiquitously and we refer the reader to table 1 in Basak, Farsoiya & Dasgupta (2021) where a comprehensive literature has been summarised.

An estimate of the span of lengths where these jets are observed is instructive. Consider the ‘spike wave’ rising up to $6$ m in height (see figure 1 in McAllister et al. (2022)) at the centre of a circular wave tank. At these large scales, the jet dynamics is governed by fluid inertia and gravity (McAllister et al. 2022). At the lower extreme in length scales, surface tension dictates the dynamics of an ejecting ethanol jet rising to about $100\ \mathrm {\mu }{\rm m}$ due to the collapse of a $35\ \mathrm {\mu }{\rm m}$ bubble (see figure 1, top row, in Lee et al. (2011) and supplementary movie 1 in that paper). In this study, we focus on small-scale jets where the dominant restoring force is surface tension, with gravity playing a negligible role. In contrast to prior numerical studies, which obtain the jet by deforming the initial interface in the shape of a truncated spherical cavity (Duchemin et al. 2002; Lai et al. 2018), the initial deformation here is chosen to be an eigenmode to the linearised problem, viz. the zeroth-order Bessel function (figure 1). This choice of the initial condition is crucial from a theoretical perspective as we justify below.

Figure 1. Three-dimensional rendition of an axisymmetric simulation (generated via rotation of the pink curve about the symmetry axis in all plots). (a) Cross-section of the interface at $\hat {t}=0$. (b) Trough with precursor dimple at $\hat {t}=0.7947\,\hat {T}$. (c) Rising jet at $\hat {t}=1.0526\,\hat {T}$. Here, $\epsilon \approx 2.3$, $l_{15}=47.9$, the full simulation domain is not shown and $\hat {T}$ is the linear time period of the eigenmode that has been excited.

From a modal expansion viewpoint, a truncated cavity with relatively sharp corners and overhang (see figure 3, lower panel, in Duchemin et al. (2002)), while being a faithful representation of an incipient bursting bubble, is a multi-valued, highly nonlinear surface deformation, which does not lend itself to eigenfunction expansion (in the radial coordinate). In contrast, the smooth perturbation that we study here excites a single eigenmode initially. Further, unlike a truncated nearly spherical bubble cavity (at small Bond number; see Duchemin et al. (2002) or figure 2 in Gordillo & Rodríguez-Rodríguez (2019)), no precursor capillary waves are generated in our simulations. It will be seen that, similar to the truncated cavity (Duchemin et al. 2002), this initial perturbation produces a dimple at moderate amplitude, giving way to a jet at larger amplitudes. This initial deformation further justifies its choice by yielding to a first-principles (weakly) nonlinear analysis employing eigenfunction expansion. The advantage of the theory, devoid of any fitting parameters, is that it describes the onset of dimple formation from a capillary wave trough. In addition, it also provides a simple physical picture explaining the birth of the dimple due to wave focusing.

In previous studies from our group, this initial condition has been explored in two regimes, viz. the ‘inviscid, gravity–capillarity’ regime in Basak et al. (2021), where an ${O}(\epsilon ^2)$ theory was presented, and the ‘viscous, linear, gravity–capillarity regime’ in Farsoiya, Mayya & Dasgupta (2017), where the ${O}(\epsilon )$ linear problem was studied. In particular, it was shown in Basak et al. (2021) that, while the ${O}(\epsilon ^2)$ theory can describe the jet qualitatively, it does not capture the onset of sharp dimple formation seen in our numerical simulations. An additional difficulty with the capillarity–gravity calculation in Basak et al. (2021) is the presence of singularities in the theory due to second harmonic resonance.

Our current nonlinear minimal model based only on capillary and inertial forces surmounts the deficiency of the second-order theory of Basak et al. (2021), obviates its singularities (Kochurin et al. 2020) and allows a mechanistic interpretation of dimple formation. We present a minimal model that demonstrates dimple and jet formation at sufficiently small scales, retaining only nonlinear curvature and inertial contributions but no contributions from gravitational or viscous forces. The absence of viscosity in our model is particularly important, as in recent literature there has been significant debate on the physical mechanism of jet formation and the role of viscosity in bubble bursting (Gañán-Calvo 2017, 2018a; Gordillo & Rodríguez-Rodríguez 2018, 2019; Gañán-Calvo & López-Herrera 2021).

While it is known now that the fastest jets from bubble bursting occur at finite viscosity (Duchemin et al. 2002; Deike et al. 2018), our study demonstrates that one can also obtain jets from inviscid wave focusing. Viscosity or gravity, both of which are neglected in our study, however, can significantly influence the jetting process, leading to regions on the Bond number versus Ohnesorge number plane (${Bo}$–${Oh}$) where jetting or droplet ejection is not expected. Owing to the neglect of gravity and viscosity in our study, these phase boundaries, however, cannot be predicted from our model. The nonlinear capillary waves and the accompanying flow generated in our simulations arise purely due to initial interfacial curvature, requiring accurate estimation of the same via the nonlinear theory. In what follows, gravity effects are neglected by restricting attention to wavelengths much smaller than the air–water capillary length scale (${\approx }2.7$ mm). Similarly, viscous effects are ignored, as the capillary-viscous length scale for air–water, $\rho \nu ^2/T \approx 0.01\ \mathrm {\mu }{\rm m}$ (Duchemin et al. 2002), is much smaller than wavelengths of interest.

Figure 1 represents three snapshots (at different times) obtained from axisymmetric simulations of the inviscid, incompressible Euler's equation with surface tension using the open-source code Basilisk (Popinet 2009, 2014). The undisturbed interface is flat, with quiescent water in a cylindrical container of radius $\hat {R}_0$. Figure 1(a) represents the initial interface excitation $\hat {\eta }(\hat {r},0) = \hat {a}_0 {\mathrm {J}}_{0}(l_q({\hat {r}}/{\hat {R}_0}))$, with ${\mathrm {J}}_0(\,{\cdot }\,)$ being the zeroth-order Bessel function and $q \in \mathbb {Z}^{+}$ indicates the number of zero crossings of ${\mathrm {J}}_0(\,{\cdot }\,)$ within the computational domain. The blue surface in the background is the axisymmetric rendition of the interface (pink curve). Evolving from the initial perturbation, the interface develops a trough (figure 1b), and as this trough rises, a dimple forms at its base. This dimple subsequently develops into a jet (figure 1c), rising significantly more than the initial perturbation at $\hat {r}=0$.

For simplicity, in our theoretical analysis, the container is assumed to be deep, implying that the longest capillary wave ($\hat {\lambda }_{{max}}$) in our domain satisfies $\hat {\lambda }_{{max}} \ll \hat {H}$, where $\hat {H}$ is the undisturbed water depth (dimensional variables have a hat). The simulations are described by two non-dimensional parameters, viz. $\epsilon \equiv \hat {a}_0({l_q}/{\hat {R}_0})$ and a number $l_q$, with $q$ being the index of the primary Bessel mode that we excite initially. These provide a non-dimensional measure of the initial perturbation amplitude $\hat {a}_0$ and wavelength $\hat {\lambda }$ (see figure 1), respectively. The numbers $l_q$ satisfy ${\mathrm {J}}_1(l_q)=0$ (Weisstein 2020), necessary for respecting the no-penetration condition at $\hat {r}=\hat {R}_0$ (see next section). A rough estimate of the width of the crest in figure 1, viz. $\hat {\lambda } \approx {2{\rm \pi} \hat {R}_0}/{l_q}$, is obtained from the asymptotic expression of ${\mathrm {J}}_{0}(\,{\cdot }\,)$ at large radius (Abramowitz & Stegun 1972). We may treat $\hat {\lambda }$ as a measure of the initial perturbation wavelength, noting that this is an approximation, as the zeros of ${\mathrm {J}}_{0}(\,{\cdot }\,)$ are not exactly equally spaced. In order to minimise container wall effects on dimple and jet formation, we require ${\hat {R}_0}/{\hat {\lambda }} \approx {l_q}/{2{\rm \pi} } \gg 1$, which can be ensured by choosing $q$ to be sufficiently large. In this study, we have chosen $q \geqslant 15$ to stay consistent with this requirement.

1.1. Physical mechanism of dimple formation

Figure 2 explains the physical mechanism of dimple formation. Plotted in figure 2(a) is $\hat {a}_0 {\mathrm {J}}_{0}(l_q({\hat {r}}/{\hat {R}_0}))$ for large and small amplitudes, viz. $\epsilon = \hat {a}_0({l_q}/{\hat {R}_0}) \approx 2.7$ (brown curve) and $\epsilon \approx 0.27$ (orange curve), with $l_{15}\approx 47.9$ for both cases. Figure 2(b) shows the initial curvature profile $\hat {\kappa }(\hat {r},0)$ for the two Bessel modes in figure 2(a) using the same colour legend. With $\hat {\eta }(\hat {r},0) = \hat {a}_0{\mathrm {J}}_{0}(l_q({\hat {r}}/{\hat {R}_0}))$ and $\epsilon \ll 1$, the linearised approximation to curvature is $-{\partial ^2\hat {\eta }}/{\partial \hat {r}^2} - ({1}/{\hat {r}})({\partial \hat {\eta }}/{\partial \hat {r}})\propto \hat {a}_0{\mathrm {J}}_{0}(l_q({\hat {r}}/{\hat {R}_0}))$ (using the Bessel equation). Thus for a small-amplitude Bessel mode (figure 2a, orange curve with $\epsilon =0.27$), its curvature profile (figure 2b, orange curve) is also a Bessel function, implying that the curvature is zero at the location(s) where the interface has zero crossings with the undisturbed liquid level. Cyan dots in figure 2(a) show the first zero crossing, and the corresponding dots in figure 2(b) verify that the curvature at this radial location is (nearly) zero at small $\epsilon$. As the flow due to this initial perturbation is driven by interfacial curvature, the above implies that, if the interface is initialised as a small-amplitude Bessel mode (e.g. $\epsilon =0.27$), we expect a standing wave whose shape at all subsequent time remains proportional to the initial Bessel mode. In particular, there are well-defined nodes in the linear description (cyan dots in figure 2a) where the velocity of the interface is zero at all time within a linear description.

Figure 2. (a) Plot of $\hat {\eta }(\hat {r},0) = \hat {a}_0 {\mathrm {J}}_{0}(l_q({\hat {r}}/{\hat {R}_0}))$ with $\epsilon \approx 2.7$ (brown curve) and $\epsilon \approx 0.27$ (orange curve). (b) Curvature profile $\hat {\kappa }(\hat {r},0)$ of the initial interface using the same colour legend as panel (a). (c) Radial inward movement (yellow arrow) of the interface zero crossing in the numerical simulation for $\epsilon = 2.7$, at $\hat {t}=0$ (purple curve) and $\hat {t}=30\Delta \hat {t}$ (brown curve); the blue arrows show the velocity field. (d–g) Radial inward focusing of the wave (indicated in blue) leading to dimple formation.

The description changes qualitatively as $\epsilon$ is increased towards unity and exceeds it. For $\epsilon > 1$ ($\epsilon =2.7$), the curvature profile indicated by the brown curve in figure 2(b) appears quite distinct from the corresponding shape of the interface in figure 2(a) (brown curve), the latter being just a large-amplitude Bessel function. Notice that, unlike the small-$\epsilon$ case described in the previous paragraph, the interface in figure 2(a) has finite positive curvature at its zero crossing (brown curve in figure 2b). A finite curvature at the zero crossing implies, in general, that a finite velocity will arise subsequently at this location (which is equivalent to saying that the zero crossings of a nonlinear interface do not behave as nodes). The local velocity field obtained from numerical simulations, around the zero crossing of the interface for the case $\epsilon = 2.7$, is depicted in figure 2(c) at $\hat {t} = \Delta \hat {t}$ and $\hat {t} = 30\Delta \hat {t}$ (where $\Delta \hat {t}$ is the computational time step). It is seen from that figure that the effect of this velocity field is to move the interface locally vertically downwards in such a manner that its zero crossing moves radially inwards – the yellow arrow in figure 2(c) highlights the inward movement of the zero crossings. This radially inward motion persists further, eventually producing a trough (figure 2d), which is significantly narrower than what would result had the initial perturbation evolved as a linear standing wave.

Figures 2(d)–(g) depict the subsequent inward motion of the zero crossing(s) and the two capillary humps (shaded in blue) in time. Note that figure 2(d) depicts the trough that forms from our large-amplitude Bessel mode excited at time $t=0$ (in brown in figure 2a). This trough is saddled by two capillary humps on either side, as seen in figure 2(d), and the radially inward motion of these two capillary humps is strongly reminiscent of capillary wave focusing in bubble bursting (Lai et al. 2018).

An important difference to note is that, unlike the bubble bursting case where a train of precursor capillary waves are typically observed especially at low Ohnesorge numbers (see figure 2, upper panel, in Deike et al. (2018) or figure 4a in Blanco-Rodríguez & Gordillo (2021)), here only two parent waves (shaded in blue) are seen in figure 2(d). Their inward focusing produces a region of large curvature at the symmetry axis, which subsequently triggers the dimple (shaded pink region in figures 2e–g). Analogous focusing of capillary waves is implicated in the formation of dimples and jets from a bursting bubble (Duchemin et al. 2002; Gañán-Calvo 2017; Gordillo & Rodríguez-Rodríguez 2018, 2019; Gañán-Calvo & López-Herrera 2021). In the following section, we develop a weakly nonlinear perturbative approach to mathematically describe the formation of this dimple. The motivation for doing this lies in the observation that the dimple is already seen at moderate values of $\epsilon \approx 0.8$ (figure 3), rendering its formation accessible to first-principles theory.

Figure 3. (a–d) Interface shape at different time instants from simulations and from second- and third-order theory for $\epsilon = 0.8$ and $l_{15} = 47.9$. The green arrows indicate the instantaneous direction of interface motion. The black dotted curve in panel (d) indicates the wave trough at an earlier instant $\hat {t}=0.45\,\hat {T}$, before the emergence of the dimple (the bottom of the trough is not seen within the vertical scale of the figure and extends further downwards). (e) Hankel transform of the interface following (5.3) in Basak et al. (2021) at the instant of dimple formation in panel (d). The vertical axis in panel (e) is a measure of the surface energy in higher-order modes ($j=15$ is the primary mode here but not shown in the figure). The ${O}(\epsilon ^3)$ theory (red symbols) describes the surface energy nearly up to $j=45 = 3\times 15$, but the ${O}(\epsilon ^2)$ theory (blue symbols) does not, and predicts zero surface energy in modes beyond $j=30=2\times 15$.

1.2. Weakly nonlinear regime: $\epsilon < 1$

We use $\epsilon \equiv {\hat {a}_0l_q}/{\hat {R}_0} < 1$ to solve the potential flow equations with cylindrical, axisymmetric surface tension terms (Mathur et al. 2007) perturbatively up to ${O}(\epsilon ^3)$. Equations (1.1a–i) are non-dimensionalised using the capillary scales, $l_q/\hat {R}_0$ and $({T l_q^3}/{\rho \hat {R}_0^3})^{1/2}$, respectively, with surface tension $T$, density $\rho$ and container radius $\hat {R}_0$. Equation (1.1a) is the Laplace equation governing the perturbation velocity potential $\phi$, (1.1b,c) are the kinematic boundary condition and the Bernoulli equation at the interface, respectively. Equations (1.1d,e) are the no-penetration and free-edge boundary conditions at the wall, respectively, ensuring that the contact angle is ${\rm \pi} /2$. Equation (1.1f) ensures that the perturbation conserves volume, while (1.1g–i) represent initial conditions:

(1.1a)$$\begin{gather} \nabla^2\phi = 0, \end{gather}$$

(1.1b)$$\begin{gather}\frac{\partial \eta}{\partial t} + \frac{\partial \eta}{\partial r}\left(\frac{\partial\phi}{\partial r}\right)_{z = \eta} = \left(\frac{\partial\phi}{\partial z}\right)_{z = \eta}, \end{gather}$$

(1.1c)$$\begin{gather}\left(\frac{\partial\phi}{\partial t}\right)_{z = \eta} \! + \frac{1}{2} \left[\left(\dfrac{\partial\phi}{\partial r}\right)^2 + \left(\dfrac{\partial\phi}{\partial z}\right)^2\right]_{z=\eta} \! - \left[\dfrac{\dfrac{\partial^2\eta}{\partial r^2}}{\left\{1 + \left(\dfrac{\partial\eta}{\partial r}\right)^2\right\}^{3/2}} + \dfrac{1}{r}\dfrac{\dfrac{\partial\eta}{\partial r}}{\left\{1 + \left(\dfrac{\partial\eta}{\partial r}\right)^2\right\}^{1/2}}\right] \nonumber\\ =0 , \end{gather}$$

(1.1d–i)\begin{gather} \left.\begin{gathered} \left(\frac{\partial\phi}{\partial r}\right)_{r=l_q}=0,\quad \left(\frac{\partial\eta}{\partial r}\right)_{r=l_q} = 0, \quad \int_{0}^{l_q} {\rm d}r\,r\eta(r,t) = 0,\\ \eta(r,0) = \epsilon \mathrm{J_0}(r), \quad \frac{\partial\eta}{\partial t}(r,0) =0,\quad \phi(r,z,0) = 0. \end{gathered}\right\} \end{gather}

These equations contain the dimensionless group $\epsilon$ and parameter $l_q$. The analytical procedure is similar to the capillarity–gravity case in Basak et al. (2021) with the important difference that we obtain corrections up to cubic order here. The necessity of going up to ${O}(\epsilon ^3)$ is seen from the exact curvature expressions in (1.1c). For $\eta (r,0) \propto \epsilon {\mathrm {J}}_{0}(r)$, the first nonlinear contribution from curvature appears at ${O}(\epsilon ^3)$.

To eliminate resonant forcing of primary mode ($q$), we have $\tau = t[1 + \epsilon ^2 \varOmega _2 + \cdots ]$. Following Basak et al. (2021), we expand both $\phi (r,z,\tau ) = \sum _{i=1}^{3}\epsilon ^{i}\phi _i(r,z,\tau )$ and $\eta (r,\tau ) = \sum _{i=1}^{3}\epsilon ^{i}\eta _i(r,\tau )$. Taylor-expanding all the $z$-dependent terms about $z$=$0$, we obtain coupled linear partial differential equations at every order ${O}(\epsilon ^{i})$, with nonlinearity appearing at $i \geqslant 2$. Using the Dini series (Basak et al. 2021) $\phi _i(r,z,\tau ) = \sum _{j = 0}^{\infty } p_i^{(j)}(\tau )\exp (\alpha _{j,q} z)\mathrm {J_0}(\alpha _{j,q} r)$ and $\eta _i(r,z,\tau ) = \sum _{j = 0}^{\infty } a_i^{(j)}(\tau )\mathrm {J_0}(\alpha _{j,q} r)$ with $\alpha _{j,q} \equiv l_j/l_q$, the task of determining $\phi _i$ and $\eta _i$ becomes that of determining $p_i^{(j)}(\tau )$ and $a_i^{(j)}(\tau )$. These satisfy simple harmonic oscillator equations with inhomogeneities at ${O}(\epsilon ^2)$ and ${O}(\epsilon ^3)$, which may be solved subject to the initial conditions in (1.1g–i). The final expressions for $\eta _1(r,\tau )$, $\eta _2(r,\tau )$ and $\eta _3(r,\tau )$ are provided below (the expressions for the modal coefficients are lengthy and so are provided in the supplementary material available at https://doi.org/10.1017/jfm.2022.854):

(1.2a)\begin{gather}\eta_1(r,\tau) = \cos(\tau)\mathrm{J_0}(r),\end{gather}

(1.2b)\begin{gather}\eta_2(r,\tau) = \frac{1}{2}\sum_{j = 1}^{\infty} [\zeta_{j,q}^{(1)}\cos(\omega_{j,q}\tau) + \zeta_{j,q}^{(2)}\cos(2\tau) + \zeta_{j,q}^{(3)}]\mathrm{J_0}(\alpha_{j,q}r),\end{gather}

(1.2c)\begin{gather}\eta_{3}(r, \tau)=\left[\vphantom{\sum_{m=1, m \neq q}^{\infty}}{\mu}^{(q)}\cos (\tau)+{\kappa}^{(q)} \cos (3 \tau) \right.\nonumber\\ \left.+\sum_{m=1, m \neq q}^{\infty}\lbrace{\gamma}_{m}^{(q)} \cos [(\omega_{m, q}+1)\tau]+{\chi}_{m}^{(q)} \cos [(\omega_{m, q}-1) \tau]\rbrace\right] \mathrm{J}_{0}(r)\nonumber\\ +\sum_{j=1, j \neq q}^{\infty}\left[\vphantom{\sum_{m=1, m \neq q}^{\infty}}{\mu}^{(j)} \cos (\tau)+{\kappa}^{(j)} \cos (3 \tau)+ {\nu}^{(j)} \cos (\omega_{j, q} \tau)\right. \nonumber\\ \left. +\sum_{m=1, m \neq q}^{\infty}\lbrace{\gamma}_{m}^{(j)} \cos [(\omega_{m, q}+1) \tau]+{\chi}_{m}^{(j)} \cos [(\omega_{m, q}-1) \tau]\rbrace\right] \mathrm{J}_{0}(\alpha_{j, q} r).\end{gather}

The expressions for $\eta _2(r,\tau )$ and $\eta _3(r,\tau )$ clearly show the excitation of modes not present initially in the spectrum, implying that surface energy initially injected into a single mode (index $q$) gets redistributed among other modes, and this is crucial for resolving dimple formation. Note the presence of bound and free wave components in the expressions for $\eta _{2}(r,\tau )$ and $\eta _{3}(r,\tau )$ in (1.2b,c) ($\omega _{j,q}^2 \equiv \alpha _{j, q}^3$).

Figure 3 compares the shape of the interface at different time instants for moderate $\epsilon =0.8$ between theory and simulations. It is clear that the ${O}(\epsilon ^3)$ theory does significantly better than its second-order counterpart, particularly in capturing the instantaneous shape and location of the dimple (figure 3d). Note that the width of this dimple is substantially smaller than the width of the collapsing wave trough (black dotted curve in figure 3d). Figure 3(e) presents the Hankel transform of the interface (following (5.3) in Basak et al. (2021)) as obtained from simulations and from theory at ${O}(\epsilon ^2)$ and ${O}(\epsilon ^3)$. Note that, for this simulation, the $15$th mode ($l_{15}$) is the primary mode and we expect the ${O}(\epsilon ^3)$ theory to describe well the surface energy in modes up to $l=45 = 3\times 15$. This is indeed seen to be so, as seen from figure 3(e), where it is also clear that the ${O}(\epsilon ^2)$ theory does not capture the surface energy in modes beyond twice the primary mode, i.e. beyond modes with index $30 = 2 \times 15$.

Figure 4 presents the vertical velocity of the dimple at its formation from simulations (red plus signs) and from ${O}(\epsilon ^3)$ theory (blue diamonds). Good agreement between the two can be seen up to $\epsilon = 1.1$. Beyond $\epsilon > 1.1$, sharply rising jets are seen in simulations and the weakly nonlinear theory systematically underestimates its velocity. A better estimate may be obtained by a two-term order-of-magnitude balance, in (1.1c). Approximating the dimple velocity to be predominantly vertical, we have $v_z^2 \sim ({T}/{\rho })({1}/{R_1} + {1}/{R_2})$, where $R_1^{-1}$ and $R_2^{-1}$ are estimated from simulations. This is shown in figure 4(a) (black circles), where it is seen that, for $\epsilon \gg 1$, this approximation does better than the ${O}(\epsilon ^3)$ theory. The limitation of the weakly nonlinear approach as $\epsilon \gg 1$ is also seen clearly in figure 4(b), where we track the interface in time at $r=0$ for $\epsilon =1.8$. It is seen (insets) that the dimple becomes very narrow with increasing $\epsilon$. Such sharp dimples involves energy transfer to modes with indices $\gg 3q$ in the spectrum (Basak et al. 2021) and fall into the strongly nonlinear regime that we discuss next.

Figure 4. (a) Dimple velocity as a function of $\epsilon$ for $l_{15}$; the inset shows a power-law fit to the simulation data (log–log), with the exponent being ${\approx }1.6$. ‘Approx’ indicates the two-term order-of-magnitude balance discussed in the text. (b) Time signal for $\epsilon = 1.8$ and $l_{15} = 47.9$. At $\epsilon \gg 1$, the weakly nonlinear theory becomes inadequate.

1.3. Strongly nonlinear regime ($\epsilon \gg 1$): self-similar evolution

As $\epsilon \gg 1$, we approach the strongly nonlinear regime where the curvature at the symmetry axis, as a consequence of wave focusing, increases sharply (singularity). This is reflected in the progressive shortening of the radial extent of the dimple as $\epsilon$ is increased; see figure 4(b) inset (also see figure 2f). In the context of jets generated from bubble bursting and Faraday waves, such local singularities have been studied and self-similar solutions obtained in their vicinity (Zeff et al. 2000; Duchemin et al. 2002; Gañán-Calvo 2017; Lai et al. 2018). We anticipate an analogous singularity in local curvature as $\epsilon$ is increased, generating sharply shooting jets.

Resorting to dimensional analysis, the functional dependence of $\hat {\eta }$ on the physical parameters of the problem is expressed as $\hat {\eta } = \hat {f}(\hat {r},\hat {t};\hat {a}_0,\hat {R}_0,{T}/{\rho },l_q)$. This may be rewritten non-dimensionally as

(1.2)\begin{equation} \eta \equiv \frac{\hat{\eta}l_q}{\hat{R}_0} = f\left(\frac{l_q}{\hat{R}_0}\hat r, \sqrt{\frac{Tl_q^3}{\rho\hat{R}_0^3}} \, \hat t;\epsilon,l_q\right) = f(r,t;\epsilon,l_q). \end{equation}

We hypothesise that, for sufficiently large $\epsilon$ and $l_q$, the local curvature at the axis of symmetry becomes singular due to wave focusing. Consequently, the temporal evolution of the jet that emerges from the dimple must be independent of $\hat {R}_0$ and $l_q$ (locally in space). This can occur only if, in a narrow region around the jet (where the local singularity is effectively felt), $\eta$ depends only on the combination ${r^3}/{t^2}$, i.e. is independent of the ratio $l_q/\hat {R}_0$. We thus have the ansatz that $\tilde {\eta }(\tilde {r}) \equiv ({\eta - z_b})/{(t - t_0)^{2/3}}= g(\tilde {r};\epsilon,l_q)$ (following Gekle & Gordillo 2010; Lai et al. 2018), where $\tilde {r} \equiv {r}/{(t - t_0)^{2/3}}$, $t_0$ is the instant of dimple formation and $z_b$ is the (scaled) location of the first minimum nearest to the symmetry axis. The scaled data are depicted in figure 5(c,d) for $\epsilon =2.3$ and $l_{15}$ and $l_{24}$. Figure 5(a,b) depict the unscaled data at different times while figure 5(c,d) show the scaled data. The excellent collapse in a region around the jet is apparent.

Figure 5. (a,b) Unscaled data and (c,d) scaled data for (a,c) $\epsilon = 2.3,\ l_{15} = 47.9$, and (b,d) $\epsilon = 2.3,\ l_{24} = 74.18$. The linear time period of the standing wave in units of $t$ is $2{\rm \pi}$. The inset in panel (b) shows the instant of dimple formation $t_0$.

The similarity scales above are the Keller & Miksis (1983) scales noted in the context of the singularity at pinch-off for a fluid wedge. Figure 5 shows that, while the temporal evolution of the jet happens self-similarly, the initial and boundary parameters are not forgotten as the scaled interface shape $\tilde {\eta } = f(\tilde {r})$ depends on the parameters $(\epsilon,l_q)$. This is analogous to Keller & Miksis (1983), where the self-similar wedge shape is not universal, but depends on the initial wedge angle. Interestingly, the localised self-similar behaviour observed here in the strongly nonlinear regime may also be found in the linearised regime for capillary waves (Keller & Miksis 1983).

For a radially unbounded domain (i.e. no confining walls) and with only an initial interface deformation (zero surface impulse), the solution to the axisymmetric linearised Cauchy–Poisson problem (Debnath 1994; Kang & Cho 2019) is $\hat {\eta }(\hat {r},\hat {t}) =\int _{0}^{\infty }k\,\mathrm {J}_0(k\hat {r})\bar {\eta }_0(k)\cos (\omega \hat {t})\,\textrm {d}k$, with $\bar {\eta }_0(k)$ being the Hankel transform of $\hat {\eta }(\hat {r},0)$. This solution lacks adequate length scales and, unless the initial condition introduces these, the resultant waves evolve self-similarly. A classic example is for pure gravity waves in two dimensions. This solution (due to Cauchy and Poisson) is discussed by Lamb (1924, art. 238). Applying the arguments of Lamb (1924) and Debnath (1994) to pure capillary waves on a radial pool, we considered a localised initial perturbation, viz. $\hat {\eta }(\hat {r},0)=({\hat {V}_0}/{2{\rm \pi} \hat {r}})\delta (\hat {r})$, where $\hat {V}_0^{1/3}$ is the only length scale in this initial condition. Dimensional reasoning predicts that, instead of $\hat {\eta }$ being a function of $\hat {r}$ and $\hat {t}$ individually, it should be possible to express the scaled interface in terms of a single scaled variable, i.e. ${\hat {\eta } T^{\prime } \hat {t}^2}/{\hat {r}\hat {V}_0} = \psi ({\hat {r}}/((T^{\prime })^{1/3} \hat {t}^{2/3}))$, where $T^{\prime } \equiv {T}/{\rho }$. The functional form of $\psi (\,{\cdot }\,)$ may be obtained by solving the aforementioned integral numerically. Alternatively, the asymptotic form for $\psi (\,{\cdot }\,)$ is deduced by applying the stationary phase technique ($\hat {r},\hat {t}\rightarrow \infty$ with ${\hat {r}}/{\hat {t}}$ fixed) (Debnath 1994) to evaluate the integral (see supplementary material). This predicts $\tilde {\tilde {\eta }} \equiv {\hat {\eta }T^{\prime }\hat {t}^2}/{\hat {V}_0\hat {r}}\sim ({2\sqrt {2}}/{9{\rm \pi} }) \sin (\frac {4}{27} \tilde {\tilde {r}}^{3})$, with $\tilde {\tilde {r}}\equiv {\hat {r}}/((T')^{1/3}\hat {t}^{2/3})$, and we note the Keller–Miksis scale in this.

Figure 6 compares the stationary phase solution (Asy) with the numerical solution to this problem at different time instants, showing self-similar behaviour. Addition of further length scales into the initial condition breaks this self-similarity. Consider linear capillary waves generated from a volume-conserving perturbation of finite width (Miles 1968; Debnath 1994), i.e. $\hat {\eta }(\hat {r},0) =\hat {d}(1-{\hat {r}^2}/{\hat {a}^2})\exp (-{\hat {r}^2}/{\hat {a}^2})$. This initial condition has two length scales, viz. $\hat {d}$ and $\hat {a}$, and dimensional analysis predicts three groups, viz. ${\rm \pi} _1\equiv {\hat {\eta }(T')^3\hat {t}^6}/{\hat {d}\hat {a}^4\hat {r}^5}$, ${\rm \pi} _2\equiv {\hat {r}}/((T')^{1/3}\hat {t}^{2/3})$ and ${\rm \pi} _3\equiv {\hat {a}^2\hat {r}^4}/((T')^2\hat {t}^4)$, implying from the stationary phase, ${\rm \pi} _1 \sim ({8\sqrt {2}}/{729})\exp (-\frac {4}{81}{\rm \pi} _3)\sin (\frac {4}{27}{\rm \pi} _2^3)$, which is not a self-similar result (see supplementary material).

Figure 6. Plot of $\hat {\eta }(\hat {r},0) =\hat {V}_0 {\delta (\hat {r})}/{2{\rm \pi} \hat {r}}$. The interface evolves self-similarly at all times. ‘Asy’ denotes the asymptotic solution from the stationary phase for large $\hat {r}$ and $\hat {t}$.

1.4. Comparison with results from the literature

In this section, we compare velocity predictions at the instant of dimple formation from our weakly nonlinear theory, and simulations in the strongly nonlinear regime with the inviscid theory presented in Gordillo et al. (2020). In this study, the authors examine a closely related problem of a liquid with a free surface contained in a narrow cylindrical tube, where a jet is triggered at the symmetry axis by focusing a laser pulse in the liquid bulk. This leads to the creation of a vapour bubble, which induces a uniform velocity in the liquid bulk below the free surface; see figure 1(c) in Gordillo et al. (2020). Also closely related to this are similar experiments in the same study where an impulse is imparted to the bottom of the tube, the ‘Pokrovski experiment’ in Antkowiak et al. (2007) and other similar experiments (Bergmann et al. 2008; Kiyama et al. 2016). The theory in Gordillo et al. (2020) adopts a potential flow approach similar to our present case and predicts that in a narrow region around the axis of symmetry, the velocity at the free surface (due to the uniform velocity induced in the bulk) is nearly radial (the initial meniscus is spherical for Gordillo et al. (2020)) and given by $\hat {v}_r = \hat {v}_n\cos (\hat {k}\psi )$ for $\psi \ll 1$, where $\psi$ is the angle measured from the centre of the circle; see figure 1(c) in Gordillo et al. (2020). As this formula is valid for small $\psi$, we compare this velocity with the vertical component of velocity $\hat {v}_z$ obtained from our numerical simulations (for large $\epsilon$) and with our weakly nonlinear theory (at smaller $\epsilon$) at the instant of dimple formation.

Note that in Gordillo et al. (2020) the parameters $\hat {v}_n$ and $\hat {k}$ in the above formula are functions of the contact angle that the meniscus makes in their case with the tube wall and are obtained from the numerical solution to the Laplace equation. For our present study, the contact angle is ${\rm \pi} /2$ and, as discussed in our introduction, we are primarily interested in the regime where the contact angle effects are minimal. This is ensured by suitably choosing $l_q$ for our initial perturbation, which assures that the wavelength $\hat {\lambda }$ of the initial perturbation (see figure 1a) is much smaller than the cylinder radius (see discussion below figure 1). Hence while comparing with the theoretical prediction for the velocity $\hat {v}_r$ with Gordillo et al. (2020), we treat $\hat {v}_n$ and $\hat {k}$ purely as fitting parameters that are influenced not by the actual contact angle but by some equivalent local contact angle. It is useful to recall that in our case the width of the dimple is far smaller than the cylinder radius (figure 2d) and it is thus expected to be relatively independent of the contact angle at the cylinder wall.

Figure 7(a) depicts a circle fitted to match the local shape of the dimple at the instant when it first emerges in our numerical simulations ($\epsilon =2.3$). The angle $\psi$ is defined with respect to the centre of this circle, as also done in Gordillo et al. (2020). Figure 7(b) (comparison with weakly nonlinear theory), figure 7(c) (comparison with strongly nonlinear simulation at $\epsilon =2.3$) and figure 7(d) (comparison with strongly nonlinear simulation at $\epsilon =2.7$) show comparisons with the formula of Gordillo et al. (2020) for specified values of the fitting parameters, depicting in general good agreement.

Figure 7. (a) A circle of radius $\hat {r}_c$ fitted locally to the dimple (blue symbols) arising for $\epsilon = 2.3$ in our numerical simulations at the instant of its formation. The angle $\psi$ is measured from the centre of this circle. (b) Comparison of the radial velocity in Gordillo, Onuki & Tagawa (2020) with the vertical velocity prediction of the dimple from our weakly nonlinear theory for $\epsilon =0.8$, with $\hat {v}_n=86.42,\ \hat {k}=2.7$. (c) Comparison with simulation for $\epsilon =2.3$, with $\hat {v}_n=489,\ \hat {k}=2.2$. (d) Comparison with simulation for $\epsilon =2.7$, with $\hat {v}_n=516.5,\ \hat {k}=3.2$.

1.5. Critical analysis of our model and results

In this section, we critically analyse the strengths and weaknesses of our theoretical and numerical model. A comparison of the vertical component of the instantaneous fluid acceleration was reported earlier for the ${O}(\epsilon ^2)$ inviscid irrotational theory with both gravity and surface tension in Basak et al. (2021) (see § 4 in the supplementary material in Basak et al. (2021)), where a qualitative analogy of the flow structure within the jet was found with that of Gekle et al. (2009). The recent discussion (Gañán-Calvo 2017, 2018b; Gordillo & Rodríguez-Rodríguez 2018, 2019; Gañán-Calvo & López-Herrera 2021) on the physical mechanism of jet formation in the context of bubble bursting has extensively debated the wavelength of the wave that dominates the focusing process. Stated in terms of time, a part of the debate relates to whether the ejection of the jet is dominated by the arrival at the symmetry axis of the smallest non-damped ripple (see discussion below figure 4 in Gordillo & Rodríguez-Rodríguez (2019)) or the main steep wavefront (using the same terminology as in Gañán-Calvo & López-Herrera (2021); see their discussion on p. 4, paragraph 1). Note that the steep wavefront is much longer than the ripples and thus they have different phase speeds.

Also under debate has been the role of viscosity in setting the critical Ohnesorge number at which jetting is observed to peak. In this context, our results provide strong evidence that a nonlinear inviscid irrotational model with only capillarity has the essential ingredients for generating a dimple at small $\epsilon$ and a jet (with tip ejection of droplets) at larger values of this parameter (with fixed $l_q$). This conclusion aligns with the suggestion of Gordillo & Rodríguez-Rodríguez (2019) that the jet is generated by an inertial mechanism (but where viscosity can significantly influence the process; see discussion below). Notably, our initial condition does not generate a precursor capillary wave train (see figure 2d–g) allowing us to identify the wave whose focusing generates the dimple and the jet in our simulations. We note that precursor ripples are typically absent in a bubble only when the viscosity is sufficiently high (see figure 3 in Deike et al. (2018) or figure 4b in Gordillo & Rodríguez-Rodríguez (2019)). As expected, our minimal model has limitations, particularly due to the neglect of viscosity (but also gravity, see below). From the results of Duchemin et al. (2002) and the subsequent numerical and experimental phase boundary on the ${Oh}$ and ${Bo}$ plot reported by Walls, Henaux & Bird (2015), the following important points emerge (for a bubble of equivalent radius $\hat {R}_c$, we define Bond number ${Bo} \equiv {\rho g \hat {R}_c^2}/{T}$ and Ohnesorge number ${Oh} \equiv {\mu }/{\rho T \hat {R}_c}$ Walls et al. (2015)):

1. (i) No jet drops are produced from sufficiently large bubbles (${Bo} > 3$) even in the inviscid (${Oh}=0$) limit.

2. (ii) At any finite viscosity, sufficiently small bubbles (e.g. those satisfying ${Oh} > 0.037$ and ${Bo}\approx 10^{-3}$) do not eject droplets and presumably no jets as well.

3. (iii) Even at ${Bo}\rightarrow 0$, the fastest jets are produced at finite and not zero viscosity, i.e. there is a critical Ohnesorge number at which the jet velocity is maximum (Duchemin et al. 2002; Gordillo & Rodríguez-Rodríguez 2019).

In analogy with bubble bursting simulations, we may also define a Bond number based on the wavelength of our initial perturbation, viz. ${Bo}_{B}\equiv {\rho g \hat {\lambda }^2}/{T}$, and an Ohnesorge number ${Oh}_{B}\equiv {\mu }/{\rho T \hat {\lambda }}$ (the subscript ‘$B$’ for Bessel). An approximate expression relating $\hat {\lambda }$ to the parameter $l_q$ and the radius of the container $\hat {R}_0$ has been discussed in the introduction to this study. Our current results without viscosity or gravity thus correspond to the limit ${Oh}_{B}\rightarrow 0,\ {Bo}_{B}\rightarrow 0$. As we argue below (see next paragraph), $\hat {\lambda }$ in our simulations behaves like the bubble radius $\hat {R}_c$, and we thus expect proportionality between the two sets of non-dimensional numbers, viz. ${Bo}_{B}\propto {Bo}$ and ${Oh}_{B}\propto {Oh}$. Thus our results are consistent with those of Walls et al. (2015) (their figure 4), where, on their Bond number versus Ohnesorge number plane, the ${Oh}\rightarrow 0,\ {Bo}\rightarrow 0$ limit lies within the jetting (jet drop) region. Additionally, a quantitative explanation of why sufficiently large bubbles (the first point above) do not produce drops (or even jets) may also be obtained from our theory, as we explain next.

Refer to figures 8(a) and 8(b), where we compare the cavity that is generated from the initial condition in our numerical simulations to a bubble cavity generated by solving the Young–Laplace equations (Berny et al. 2020) without the bubble cap. It is visually apparent that the wavelength $\hat {\lambda }$ for the cavity generated from our perturbation is qualitatively like the effective diameter $2\hat {R}_c$ of a bubble cavity (based on the volume of an equivalent sphere), while the amplitude $\hat {a}_0$ for our cavity (generated from a Bessel mode) is akin to the submergence depth $Z_c$ of the bubble (Puthenveettil et al. 2018). Consequently, the ratio $\hat {Z}_c/2\hat {R}_c$ for a submerged static bubble behaves like the steepness parameter $\epsilon \equiv {\hat {a}_0l_q}/{\hat {R}_0}$ for our initial perturbation. For a static bubble, $\hat {Z}_c/\hat {R}_c$ is uniquely determined by the Bond number and, as ${Bo}$ is increased, $\hat {Z}_c/\hat {R}_c$ decreases monotonically (the static bubble becomes larger, deforms radially and protrudes out of the liquid more and more as ${Bo}$ increases; see figure 5 in Walls et al. (2015)). This implies that, with increasing Bond number, the equivalent $\epsilon$ for a bubble tends to zero, which is the linear limit in our theory, where jetting is not observed. It may thus be anticipated that, at large Bond number, a relatively large and radially extended bubble generates a spectrum of capillary waves (after bursting) whose subsequent time evolution is well described by linear superposition. As this is a linear description, there is no dimple or jet formed from such a bubble, thereby explaining why large bubbles fail to produce jet drops.

Figure 8. (a) Cavity produced in our simulations starting with a Bessel mode perturbation, with $\epsilon =2.7$ and $l_{15}$. Here, $\hat {\lambda }$ and $\hat {a}_0$ provide a measure of the initial perturbation wavelength and amplitude, respectively, and the resultant cavity dimensions are approximately proportional to these. (b) A bubble cavity for ${Bo}=0.001$ generated using the Basilisk sandbox code (Berney 2022), which generates the static bubble shape solving equations analogous to Toba (1959).

Our theory cannot, however, rationalise the last of the two aforementioned points because these relate to viscosity. We note here that, while the divergence of curvature is sufficient in our model to trigger the dimple, in reality, viscous dissipation also diverges at the symmetry axis as wave focusing proceeds in time. This has been demonstrated clearly by Boulton-Stone & Blake (1993) (their figure 11c,d). As our theory does not model this phenomenon, it cannot describe the critical Ohnesorge number(s) where jetting peaks or ceases. As a first step towards a viscous description employing our initial condition, we have also conducted viscous simulations presented in the Appendix. Note that, in a confined simulation geometry such as ours, it is inconsistent to have a viscous simulation with a moving contact line due to the associated singularity (Huh & Scriven 1971). To alleviate this issue, we have conducted simulations employing eigenmodes which have been obtained from pinning the interface at the cylinder wall. For a partially liquid-filled circular cylinder with a flat interface pinned at the cylinder wall, the radial eigenmodes have been recently reported in Shao et al. (2021). As shown in figure 9(a), we choose the $15$th axisymmetric eigenmode obtained by following the calculation in Shao et al. (2021). This eigenmode is provided as an initial surface perturbation with a sufficiently large perturbation amplitude to our viscous Navier–Stokes solver. It is seen from figure 9(a–d) that, in qualitative agreement with the inviscid simulations presented earlier employing Bessel modes, the pinned viscous simulation also shows a clear dimple and subsequent jet. The theoretical analysis of this situation using a weakly nonlinear viscous theory is underway and will be reported subsequently.

In conclusion, we have employed a single Bessel mode to produce dimples in the weakly nonlinear regime and jets in the strongly nonlinear regime. The choice of initial condition is crucial, as it allows us to theoretically explain wave focusing and radial inward motion of the capillary humps in terms of (nonlinear) generation of new eigenmodes. Our first-principles calculation establishes that a minimal nonlinear inviscid irrotational model with capillary forces only can demonstrate dimple and jet formation. Both viscosity and gravity, however, retain strong influences, especially when $0.001 < {Oh} < 1$ and $0 < {Bo} < 1$, a range that is very important for bursting bubbles at a liquid–gas free surface. The inviscid self-similarity observed in our simulations retains memory of boundary and initial parameters, in contrast to that observed in capillary pinch-off of an axisymmetric drop, where memory of initial parameters is lost (Day, Hinch & Lister 1998).