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. Reference Lai, Eggers and Deike2018).

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. (Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Séon2018) or figure 4a in Blanco-Rodríguez & Gordillo (Reference Blanco-Rodríguez and Gordillo2021)), 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. Reference Duchemin, Popinet, Josserand and Zaleski2002; Gañán-Calvo Reference Gañán-Calvo2017; Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2018, Reference Gordillo and Rodríguez-Rodríguez2019; Gañán-Calvo & López-Herrera Reference Gañán-Calvo and López-Herrera2021). 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. (Reference Basak, Farsoiya and Dasgupta2021) 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. Reference Mathur, DasGupta, Selvi, John, Kulkarni and Govindarajan2007) 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. (Reference Basak, Farsoiya and Dasgupta2021) 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. (Reference Basak, Farsoiya and Dasgupta2021), 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. Reference Basak, Farsoiya and Dasgupta2021) $\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. (Reference Basak, Farsoiya and Dasgupta2021)) 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. Reference Basak, Farsoiya and Dasgupta2021) 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. Reference Zeff, Kleber, Fineberg and Lathrop2000; Duchemin et al. Reference Duchemin, Popinet, Josserand and Zaleski2002; Gañán-Calvo Reference Gañán-Calvo2017; Lai et al. Reference Lai, Eggers and Deike2018). 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 Reference Gekle and Gordillo2010; Lai et al. Reference Lai, Eggers and Deike2018), 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 (Reference Keller and Miksis1983) 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 (Reference Keller and Miksis1983), 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 Reference Keller and Miksis1983).

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 Reference Debnath1994; Kang & Cho Reference Kang and Cho2019) 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 (Reference Lamb1924, art. 238). Applying the arguments of Lamb (Reference Lamb1924) and Debnath (Reference Debnath1994) 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 Reference Debnath1994) 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 Reference Miles1968; Debnath Reference Debnath1994), 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. (Reference Gordillo, Onuki and Tagawa2020). 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. (Reference Gordillo, Onuki and Tagawa2020). 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. (Reference Antkowiak, Bremond, Le Dizes and Villermaux2007) and other similar experiments (Bergmann et al. Reference Bergmann, De Jong, Choimet, Van Der Meer and Lohse2008; Kiyama et al. Reference Kiyama, Tagawa, Ando and Kameda2016). The theory in Gordillo et al. (Reference Gordillo, Onuki and Tagawa2020) 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. (Reference Gordillo, Onuki and Tagawa2020)) 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. (Reference Gordillo, Onuki and Tagawa2020). 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. (Reference Gordillo, Onuki and Tagawa2020) 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. (Reference Gordillo, Onuki and Tagawa2020), 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. (Reference Gordillo, Onuki and Tagawa2020). 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. (Reference Gordillo, Onuki and Tagawa2020) 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 (Reference Gordillo, Onuki and Tagawa2020) 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. (Reference Basak, Farsoiya and Dasgupta2021) (see § 4 in the supplementary material in Basak et al. (Reference Basak, Farsoiya and Dasgupta2021)), where a qualitative analogy of the flow structure within the jet was found with that of Gekle et al. (Reference Gekle, Gordillo, van der Meer and Lohse2009). The recent discussion (Gañán-Calvo Reference Gañán-Calvo2017, Reference Gañán-Calvo2018b; Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2018, Reference Gordillo and Rodríguez-Rodríguez2019; Gañán-Calvo & López-Herrera Reference Gañán-Calvo and López-Herrera2021) 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 (Reference Gordillo and Rodríguez-Rodríguez2019)) or the main steep wavefront (using the same terminology as in Gañán-Calvo & López-Herrera (Reference Gañán-Calvo and López-Herrera2021); 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 (Reference Gordillo and Rodríguez-Rodríguez2019) 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. (Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Séon2018) or figure 4b in Gordillo & Rodríguez-Rodríguez (Reference Gordillo and Rodríguez-Rodríguez2019)). 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. (Reference Duchemin, Popinet, Josserand and Zaleski2002) and the subsequent numerical and experimental phase boundary on the ${Oh}$ and ${Bo}$ plot reported by Walls, Henaux & Bird (Reference Walls, Henaux and Bird2015), 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. (Reference Walls, Henaux and Bird2015)):

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. Reference Duchemin, Popinet, Josserand and Zaleski2002; Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2019).

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. (Reference Walls, Henaux and Bird2015) (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. Reference Berny, Deike, Séon and Popinet2020) 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. Reference Puthenveettil, Saha, Krishnan and Hopfinger2018). 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. (Reference Walls, Henaux and Bird2015)). 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 Reference Berney2022), which generates the static bubble shape solving equations analogous to Toba (Reference Toba1959).

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 (Reference Boulton-Stone and Blake1993) (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 Reference Huh and Scriven1971). 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. (Reference Shao, Wilson, Saylor and Bostwick2021). As shown in figure 9(a), we choose the $15$th axisymmetric eigenmode obtained by following the calculation in Shao et al. (Reference Shao, Wilson, Saylor and Bostwick2021). 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 Reference Day, Hinch and Lister1998).