3.1. Linear solution ($i=1$)

At ${O}(\epsilon )$, $\mathscr {N}_i(r,T_0,T_2)=\mathscr {M}_i(r,T_0,T_2)=0$ which reduces (3.8a) to a homogeneous equation along with initial conditions obtained from (3.9) for $i=1$, i.e.

(3.10a,b)\begin{equation} p_1^{(\,j)}(0,0) = 0, \quad D_0p_1^{(\,j)}(0,0) ={-}\delta_{qj}. \end{equation}

This leads to

(3.11)\begin{equation} p_1^{(q)}(T_0,T_2) = (\mu_{1q}(T_2)\cos T_0+\nu_{1q}(T_2)\sin T_0)\\ \end{equation}

and

(3.12)\begin{equation} \left.\begin{array}{c@{}} a_1^{(q)}(T_0,T_2) = (\mu_{1q}(T_2)\sin T_0-\nu_{1q}(T_2)\cos T_0), \\ p_1^{(\,j)}(T_0,T_2) = 0, \quad a_1^{(\,j)}(T_0,T_2) = 0\quad \forall j\neq q, \end{array}\right\} \end{equation}

where $\mu _{1q}(T_2)$ and $\nu _{1q}(T_2)$ in 3.11 and 3.12 are functions of the slow time scale $T_2$ and satisfy initial conditions $\mu _{1q}(0)=0$, $\nu _{1q}(0)=-1$, obtained from (3.10). Hence at ${O}(\epsilon )$, $\phi _1$ and $\eta _1$ are obtained as

(3.13)\begin{equation} \phi_1(r,z,T_0,T_2) = (\mu_{1q}(T_2)\cos T_0+\nu_{1q}(T_2)\sin T_0)J_0(r)\exp (z) \end{equation}

and

(3.14)\begin{equation} \eta_1(r,T_0,T_2) = (\mu_{1q}(T_2)\sin T_0-\nu_{1q}(T_2)\cos T_0)J_0(r). \end{equation}

Note that if we were to terminate our calculation at this order, then $\mu _{1q}$ and $\nu _{1q}$ would be treated as constants being $\mu _{1q}=0,\nu _{1q}=-1$ and substituting these in (3.13) and (3.14), we obtain the standing wave that is predicted at linear order, i.e.

(3.15a,b)\begin{equation} \phi(r,z,T_0) ={-}\epsilon\sin(T_0)J_{0}(r)\exp(z),\quad \eta(r,T_0) = \epsilon\cos(T_0)J_{0}(r). \end{equation}

Equations (3.15) predict that the interface $\eta$ evolves temporally as a standing wave with unit frequency (in non-dimensional sense). Expectedly, at all time, there is only the primary mode and no transfer of energy to modes other than this mode (index $q$) is predicted at this order. Consequently, there is no jet seen in the time evolution of the standing wave. At the axis of symmetry, the initial perturbation represented by $\eta$ in (3.15) has an amplitude $\epsilon$ at $T_0=0$ and after one time period at $T_0=2{\rm \pi}$, once again rises to the same amplitude at this location. While we have seen earlier that the radial velocity at the node ($r = r^{*}$) is non-zero and radially inward already at this order, this velocity does not affect the time evolution of the standing wave. To theoretically predict the radial inward movement of nodes which leads to the jet seen in figure 4, our calculation must be extended to nonlinear order. We expect to obtain amplitude equations governing $\mu _{1q}(T_2)$ and $\nu _{1q}(T_2)$ appearing in 3.13 and 3.14 by going up to nonlinear order, and this is presented next.

3.2. Nonlinear solution ($i=2,3$)

On a longer time scale $T_2$, the amplitude of the primary mode ( $j=q$) gets modulated, this being dictated by nonlinear equations governing $\mu _{1q}(T_2)$ and $\nu _{1q}(T_2)$. To determine these equations, it becomes necessary to obtain expressions up to ${O}(\epsilon ^3)$ for $\mathscr {M}_2$ and $\mathscr {N}_2$ as well as $\mathscr {M}_3$ and $\mathscr {N}_3$), where resonant forcing of the primary mode is encountered and eliminated. The nonlinear calculation proceeds in an analogous manner to the linear one outlined above, although the algebra, expectedly, becomes increasingly tedious. We outline only the important steps here.

For $i=2$, $\mathscr {M}_2(r,T_0,T_2)\neq 0$ and $\mathscr {N}_2(r,T_0,T_2)\neq 0$, and expressions for these are provided in the Appendix. As a result, (3.8a) governing $p_2^{(i)}(T_0,T_2)$ is an inhomogeneous equation. To prevent very lengthy calculation for solving this equation, we introduce an approximation as follows. While solving (3.8a) for $p_2^{(\,j)}(T_2)\;(\,j=0,1,2,3,\ldots )$, constants of integration appear which, strictly speaking, are not constants but functions of the slow time scale $T_2$. These govern the modulation of the amplitude of the secondary modes. To determine the equations governing this modulation, one needs to go to higher order (presumably beyond third order). Our purpose here is to obtain a nonlinear approximation which contains enough physics to model the inception of the jet seen earlier. Thus, to prevent very lengthy, higher-order calculations, we ignore the temporal variation of these functions treating them instead as constants to be determined from initial conditions. Our analytical solution thus contains the modulation of the amplitude of the linear solution (i.e. the primary mode with index $q$), but neglects amplitude modulation of the secondary modes generated in the spectrum, via nonlinearity. The justification for this will be obtained a posteriori, when we compare our analytical predictions with numerical simulations. After lengthy algebra, the ${O}(\epsilon ^2)$ corrections are found to be

(3.16)\begin{align} \phi_2(r,z,T_0,T_2) &= \left[\frac{\nu_{1q}(T_2)^2-\mu_{1q}(T_2)^2}{2}\sin (2T_0)+\mu_{1q}(T_2)\nu_{1q}(T_2)\cos (2T_0)\right]J_0^2(l_q)\nonumber\\ &\quad +\sum_{n = 1}^{\infty}[\xi_{n,q}^{(2)}\sin(\omega_{n,q}T_0) + \xi_{n,q}^{(3)}(\mu_{1q}(T_2),\nu_{1q}(T_2))\cos(2T_0) \nonumber\\ &\quad +\xi_{n,q}^{(4)}(\mu_{1q}(T_2),\nu_{1q}(T_2))\sin(2T_0)] \exp(\alpha_{n,q}z)J_0(\alpha_{n,q}r) \end{align}

and

(3.17)\begin{align} &\eta_2(r,T_0,T_2) = \sum_{n = 1}^{\infty}[\zeta_{n,q}^{(1)}\cos(\omega_{n,q}T_0) + \zeta_{n,q}^{(3)}(\mu_{1q}(T_2),\nu_{1q}(T_2))\cos(2T_0) \nonumber\\ &\quad +\zeta_{n,q}^{(4)}(\mu_{1q}(T_2),\nu_{1q}(T_2))\sin(2T_0) +\zeta_{n,q}^{(5)}(\mu_{1q}(T_2),\nu_{1q}(T_2))]J_0(\alpha_{n,q}r). \end{align}

In (3.16), $\xi _{n,q}^{(1)}=0,\xi _{n,q}^{(2)}$ are constants (part of the complementary function for (3.8a), these being assumed to be constants instead of being functions of $T_2$, as explained earlier). Expressions for $\xi _{n,q}^{(2)}$, $\xi _{n,q}^{(3)}(\mu _{1q}(T_2),\nu _{1q}(T_2))$ and $\xi _{n,q}^{(4)}(\mu _{1q}(T_2),\nu _{1q}(T_2))$ are provided in the Appendix. Similarly, in (3.17), $\zeta _{n,q}^{(1)}$ and $\zeta _{n,q}^{(2)}=0$ are constants while expressions for $\zeta _{n,q}^{(3)}(\mu _{1q}(T_2),\nu _{1q}(T_2))$, $\zeta _{n,q}^{(4)}(\mu _{1q}(T_2),\nu _{1q}(T_2))$ and $\zeta _{n,q}^{(5)}(\mu _{1q}(T_2),\nu _{1q}(T_2))$ are provided in the Appendix. Note that these expressions depend on the unknowns $\mu _{1q}(T_2)$ and $\nu _{1q}(T_2)$. Equations for these two unknowns are obtained at ${O}(\epsilon ^3)$ via elimination of resonant forcing terms for the primary mode.

3.3. Amplitude equation(s)

To obtain the equations governing $\mu _{1q}(T_2)$ and $\nu _{1q}(T_2)$, it is necessary to carry out the calculation till ${O}(\epsilon ^3)$ and seek terms which resonate with the primary mode. As the natural frequency on the left-hand side of (3.18a) is unity for $j=q$ ($\omega _{q,q}^2=1$), we look for terms on the right-hand side in (3.18a) which are proportional to $\cos (T_0)$ or $\sin (T_0)$. Writing (3.8a) for $i=3$ and $j=q$,

(3.18)\begin{equation} (D_0^2 + 1)p_3^{(q)}(T_0,T_2) = \frac{2}{l_q^2 J_0^2(l_q)}\int_{0}^{l_q} {\rm d}r\, r J_0(r) \{D_0\mathscr{N}_3 + \mathscr{M}_3 \}. \end{equation}

Expressions for $\mathscr {M}_3(r,T_0,T_2),\mathscr {N}_3(r,T_0,T_2)$ necessary to evaluate the right-hand side of (3.18) are provided in the Appendix. Note in particular that $\mathscr {M}_3$ and $\mathscr {N}_3$ both contain terms like $D_2\eta _1$ and $D_2\phi _1$ which lead to time derivatives of $\mu _{1q}(T_2)$ and $\nu _{1q}(T_2)$. Elimination of the coefficients of $\cos (T_0)$ and $\sin (T_0)$ appearing on the right-hand side of (3.18) lead us to two nonlinear, coupled ordinary differential equations governing evolution of $\mu _{1q}(T_2)$ and $\nu _{1q}(T_2)$ over the slow time scale $T_2$. The algebra is again quite lengthy and we provide only the amplitude equations here. These are of the form

(3.19a,b)\begin{equation} \left.\begin{array}{c} D_2\mu_{1q}=r_1\mu_{1q}+r_2\nu_{1q},\quad D_2\nu_{1q}=r_3\mu_{1q}-r_1\nu_{1q},\\ \quad \mu_{1q}(0) = 0,\ \nu_{1q}(0) ={-}1 \end{array}\right\}. \end{equation}

Here, $r_{1}$, $r_{2}$ and $r_{3}$ in (3.19a) and (3.19b) are nonlinear functions of $\mu _{1q}$ and $\nu _{1q}$ whose expressions are provided in the Appendix. Equations (3.19a) and (3.19b) are solved numerically in MATLAB and together with expressions at linear and quadratic order provided earlier in (3.13), (3.14), (3.16) and (3.17), we obtain a second-order solution to the primary CP problem for our initial condition. This analytical prediction is free from any fitting parameters and is tested against numerical simulations of the incompressible Euler equation with gravity, in the next section.

3.4. Comparison of theory with numerical simulations

In this section, the results from numerical simulations are compared against the theoretical predictions made earlier. To do this, we truncate the infinite series in (3.16) and (3.17) at $n=70$. This number is twice the index of the primary mode, which in the present study has been chosen to be $q=35$. Figure 12 in Appendix B presents the Hankel transform of the interface for $\epsilon =0.5$ from our ${O}(\epsilon ^2)$ theory, comparing with the Hankel transform of the interface obtained from numerical simulations. This is done at $\hat {t}=0.483$ s and the interface at this instant of time is depicted in figure 6(d). It is clearly seen that at this instant of time, when the interface shows more than $37\,\%$ overshoot at the symmetry axis (figure 6d), there is a perceptible amount of potential energy around the second harmonic of the primary mode, i.e. $q=70$, see inset of figure 12. However, our theoretical model predicts very little potential energy beyond this and consequently, we have truncated our expansion at $n=70$ in expressions (3.16) and (3.17).

Figure 6.

(a–d) Comparison of the interface profile ${\eta }/{\epsilon }$ at various instants with linear and weakly nonlinear, analytical solutions. Solid red, ${O}(\epsilon ^2)$ weakly nonlinear theory; solid brown, ${O}(\epsilon )$ linear theory and numerical simulation (Sim, solid blue line) for $\epsilon =0.5$ and $l_{35}=110.74$ (Case $2$ in table 1). (b) Inset depicts a close-up view of the first zero crossing of the perturbed interface. Note the qualitative differences between the simulation and linear theory particularly in (c) and (d). The shaded region in blue in (c) depicts the wave crest whose inward focusing results in an overshoot. Significant overshoot ($\approx$37 %) is seen in the numerical simulation, i.e. above the dashed pink line in (d), at $r=0$. Unlike linear theory, the weakly nonlinear theory predicts this overshoot nearly correctly as seen in (d). (a) $\hat {t}=0$ (non-dimensional $t=0$), (b) $\hat {t}=0.266$ s ($t=3.58$), (c) $\hat {t}=0.337$ s ($t=5.07$) and (d) $\hat {t}=0.483$ s ($t=6.5$).

The numerical simulations have been carried out using the open source code Basilisk (Popinet Reference Popinet2014) employing a cylindrical domain of radius $\hat {R}_0=600$ cm and depth $\hat {H}=300$ cm. As the chosen value of $q=35$, this implies that the (approximate) wavelength of the initial perturbation is far smaller than $\hat {H}=300$, justifying the deep water approximation that we make in our theory. The parameters of various simulations are summarised in table 1. In all simulations, we have employed free-edge boundary conditions at the wall implying that the interface intersects the solid boundary at $\hat {R}_0$ at an angle of ${\rm \pi} /2$ at all times.

Figure 6 shows the comparison of linear and weakly nonlinear theory with the simulation for $\epsilon =0.5$. It can be seen that the weakly nonlinear theory captures the simulation profiles quite well. Note in particular figure 6(d), showing clearly that the overshoot seen in the simulation at the symmetry axis over and above unity (pink dashed line in the figure) is predicted well by the ${O}(\epsilon ^2)$ theory. This overshoot is the precursor to the slender jet that develops at the axis of symmetry, seen for far larger values of $\epsilon =1.5$, please refer to figure 4. Figure 7 shows the inward movement of the first zero crossing of the interface $\eta$ (around $r=-2.4$ in figure 6b) from panels (b) to (c). Note the good agreement in figure 7 between nonlinear theory and simulation. During this time window, the interface around the axis of symmetry is shaped like a cavity which collapses subsequently, leading to the overshoot seen in figure 6(d). In figure 8, we test the limit of our theory, extending the magnitude of $\epsilon$ to $0.9$ (Case $3$ in table 1). It is seen that although the ${O}(\epsilon ^2)$ theory shows the overshoot, it significantly underdescribes the jet being unable to capture its slenderness well. For reference, the initial condition is also provided in this figure in green. This value of $\epsilon \sim 0.9$ represents a simulation which lies outside the range of utility of the ${O}(\epsilon ^2)$ theory for describing these jets. While the theory is able to describe the inception (up to $\epsilon =0.5$), it is unsuitable in this parametric regime where nonlinearity becomes very strong. This regime will be analysed in detail in the next sub-section. We note that the inability of the current multiple scale analysis to capture the overshoot accurately for $\epsilon > 0.6$ is unrelated to the approximation discussed above in (3.16). To check this, we have compared our current solution with the stretched coordinate solution reported earlier by Basak et al. (Reference Basak, Farsoiya and Dasgupta2021) (setting surface tension equal to zero in the latter) and found that predictions from the two methods are indistinguishable. Figure 9(a) compares the overshoot seen at the symmetry axis ($r=0$) as a function of $\epsilon$. This figure demonstrates that the weakly nonlinear theory presented here may be used reliably up to $\epsilon \sim 0.6$ beyond which it becomes quite inaccurate. Figure 9(b) demonstrates that the displacement of the tip of the jet, viz. $\hat {\eta }(\hat {r}=0,\hat {t})$, displays a parabolic dependence in time, in agreement with the earlier observations of Ghabache et al. (Reference Ghabache, Séon and Antkowiak2014).

Figure 7.

Radial inward motion of the first zero crossing between figures 6(b) and 6(c). Note from figure 6 that within this time window, the interface is shaped as a cavity around $r=0$ and the collapse of this cavity occurs nonlinearly, leading to the overshoot seen in figure 6. As noted earlier, the inward motion of the first zero crossing of $\eta (r=r^{*},t)$ may be interpreted as radial inward focusing of the two humps shown earlier in figure 5. As seen, linear theory (pink symbols) does not allow any radial displacement of the zero crossings. However, the weakly nonlinear theory is able to describe the inward motion of the nodes correctly qualitatively, although some quantitative differences persist between theory and simulations. The inset depicts the main figure with dimensional axes.

Figure 8.

Interface for $\epsilon =0.9$, $\hat {t}=0.47$ s (non-dimensional $t=6.33$) before the maximum height is reached. The ${O}(\epsilon ^2)$ theory is unable to describe the slender jet seen in simulations although it manages to get the overshoot correct qualitatively. The line in green represents the initial condition.

Figure 9.

(a) Comparison of the maximum height ($\hat {\eta }_{max}/\hat {a}_0$) at the symmetry axis, as a function of increasing $\epsilon$. Note that $\eta _{overshoot}\equiv {\hat {\eta }_{max}}/{\hat {a}_0} - 1 = {O}(\epsilon ^2)$ which translates in dimensional variables to $\hat {\eta }_{overshoot} \sim {\hat {a}_0^3}/{\hat {R}_0^2}$ (Ghabache et al. Reference Ghabache, Séon and Antkowiak2014). Our weakly nonlinear model estimates the overshoot with an error of <20 % up to $\epsilon \sim 0.6$. Up to this $\epsilon =0.6$, the shape of the jet in theory and simulation agree reasonably well. Beyond this, the weakly nonlinear model becomes quite inaccurate. Note that the analytical model displays overshoot earlier in time compared with the simulation. The thin blue line indicates the linear solution which does not have any overshoot. The green curve has the equation $0.32\epsilon ^2 + 0.19\epsilon + 1.1$ and is obtained by fitting a parabola to the analytical model (diamonds). (b) Variation of non-dimensional amplitude with time at the axis of symmetry. Dots represent simulation data while solid curves are parabolic fits to the same, in qualitative agreement with observations of Ghabache et al. (Reference Ghabache, Séon and Antkowiak2014, figure 2).

We have argued earlier that the development of the jet requires radially inward focusing of surface gravity waves at the axis of symmetry, much akin to the phenomena of bubble collapse albeit at far smaller capillary scales, where radially inward focusing of capillary waves are implicated in the formation of a slender jet (Duchemin et al. Reference Duchemin, Popinet, Josserand and Zaleski2002; Gordillo & Rodríguez-Rodríguez Reference Gordillo and Rodríguez-Rodríguez2019). Predictions from our theoretical model vis-a-vis numerical simulations presented in figure 6, where an inward motion of the first nodes is visible, see panel (b) inset for the weakly nonlinear theory as well as simulations, demonstrate that this radially inward flow focusing at large length scales may also be interpreted as the reason for the generation of jets at these large length scales. Alternatively, this process may also be interpreted as constructive interference at the symmetry axis from a large number of (free) modes in the spectrum. In our study, these free modes are produced via quadratic interactions of the monochromatic initial condition.

3.5. Strongly nonlinear regime, inertial regime

In the previous section, we have studied simulations for $\epsilon < 1$ and the weakly nonlinear theory was able to describe the inception of the jet up to $\epsilon \approx 0.5$. We now move to the strongly nonlinear regime ($\epsilon >1$) where, due to stronger nonlinearity, the theory becomes inadequate and a slender jet emerges in the numerical simulations (figure 4). Here, qualitative understanding may be obtained by referring to the analytical solution on hyperboloidal jets by Longuet-Higgins (Reference Longuet-Higgins1983) in axisymmetric, extensional flow. It will be seen that this analytical model of Longuet-Higgins (Reference Longuet-Higgins1983) provides a reasonably good approximation to the jets that arise in our simulations, when $\epsilon >1.2$.

Interestingly, this exact solution is singular at $\hat {t}=0$ where it predicts a divergent velocity ($\dot {c}\sim \hat {t}^{-1/3}$) and a divergent acceleration ($\ddot {c}\sim \hat {t}^{-4/3}$) at the tip of the free surface shaped as a hyperboloid, viz. at $\hat {z}=c(\hat {t})$. In a time window around this initial ‘jerky start’ to motion, the half-angle $\theta$ between the asymptotes to the hyperbola are predicted to evolve from its initial value $\theta (0) = \tan ^{-1}{\sqrt {2}} \approx 54.5^{\circ }$ as (Longuet-Higgins Reference Longuet-Higgins1983)

(3.20)\begin{equation} \tan\theta(\hat{t}) \sim \sqrt{2}\left[1 + \frac{U}{L}\left(\frac{\hat{t} - \hat{t}_0}{\sqrt{3}}\right)\right]^{{-}3/2}, \end{equation}

where $\hat {t}_0$ is time shift that will be necessary for comparison with our simulations, but is zero for Longuet-Higgins (Reference Longuet-Higgins1983).

To compare our strongly nonlinear jets ($\epsilon >1$) with the solution of Longuet-Higgins (Reference Longuet-Higgins1983), we first extract the time window where gravity may be neglected in our simulations. Figure 10(a,c,e,g) depict the interface displacement at the symmetry axis, viz. $\hat {\eta }(\hat {r}=0,\hat {t})$ versus time for various $\epsilon > 1$. We fit parabolas (in red and green) of the form $\alpha \hat {t}^2 + {O}(\hat {t})$ to each of these datasets to estimate the average acceleration in this time window, numerically. This is done at early time (red parabolas indicating the early time window of jet inception indicated by the formula $\hat {\eta }_e(0,\hat {t}) = \alpha \hat {t}^{2} + {O}(\hat {t})$) and at a later time (green parabolas indicating a late time window when a well-developed jet has already formed and indicated by a similar formula $\hat {\eta }_l(0,\hat {t})$). The value of $2\alpha$ from each of these datasets provides an approximate measure of the acceleration of the interface viz. $({\partial ^2\hat {\eta }}/{\partial \hat {t}^2})(\hat {r}=0,\hat {t})$ at the symmetry axis. The jerkiness in the data for figures 10(e) and 10(g) is due to the numerical noise associated with the tendency to eject small droplets at the jet tip, which typically tear off and are subsequently excluded in our estimation for $\hat {\eta }(\hat {r}=0,\hat {t})$. For these two figures, as a slightly better estimate of the jet acceleration, we have also estimated the acceleration as above but by tracking the interface displacement at a radial location slightly off the symmetry axis where drop formation does not occur. This is depicted as $\hat {\eta }(\hat {r}\approx 1.8,\hat {t})$ and is plotted on the left vertical axes of figures 10(e) and 10(g).

Figure 10.

(a,c,e,g) Interface height at the symmetry axis, i.e. $\hat {\eta }(\hat {r}=0,\hat {t})$ versus time measured from simulations at various $\epsilon > 1$. Acceleration at the symmetry axis ($\hat {r}=0$) or slightly off it ($\hat {r}\approx 1.8$) is measured by fitting parabolas (see text) to these data and estimating the second derivative. (b,d,f,h) Comparison of the angle $\theta$ (in degrees) enclosed at the jet crest starting from $\hat {t}=\hat {t}_0$ with the prediction by Longuet-Higgins (Reference Longuet-Higgins1983). The variable $\tau \equiv {U(\hat {t} - \hat {t}_0)}/{L}$ is the non-dimensional time. (a,b) $\epsilon =1.2$, (c,d) $\epsilon =1.3$, (e,f) $\epsilon =1.5$ and (g,h) $\epsilon =1.6$.

Note from these figures that around the time of jet formation, the value of $2\alpha$ (red parabolas) climbs from $5236$ cm s$^{-2}$ (i.e. $5{\times }g$) at $\epsilon =1.2$ (figure 10a) to approximately 21 684 cm s$^{-2}$ (${\approx } 20{\times }g$) at $\epsilon =1.6$ (figure 10g). Also note that the green parabolas in these figures predict a value of $2\alpha$ ranging from $-$965 cm s$^{-2}$ for figure 10(a) to ${\approx } -524$ cm s$^{-2}$ for figure 10(g) (indicating that the jet is still experiencing upward acceleration and has not relaxed yet to $g=-980$). One must note that these numbers only comprise rough estimates for acceleration, this being due to estimating the second derivative from noisy data introducing significant uncertainties. We also caution the reader that for the simulations for $\epsilon =1.2- 1.6$ in figures 20–23 (Appendix E), the detailed shape of the jet in a small region around the symmetry axis ($r=0$) remains grid dependent at the current resolution of $2048^2$, further affecting the accuracy of the asymptotes and the jet tip acceleration reported here. Subject to these uncertainties, the qualitative conclusion inferred from the above analysis is that for $\epsilon > 1$, there exists a time window when gravity may be neglected in our simulations, enabling comparison with the solution of Longuet-Higgins (Reference Longuet-Higgins1983). To compare the angle $\theta$ obtained from our simulations with the prediction from (3.20), we set $\hat {t}_0$ to be the instant of time when $\theta$ is closest to the initial angle in the model of Longuet-Higgins (Reference Longuet-Higgins1983), i.e. ${\approx } 54.5^{\circ }$. This definition is also consistent with recent experimental measurements of large-scale jets by McAllister et al. (Reference McAllister, Draycott, Davey, Yang, Adcock, Liao and van den Bremer2022), who too compared their jet angle with the model of Longuet-Higgins (Reference Longuet-Higgins1983). Note that the model of Longuet-Higgins (Reference Longuet-Higgins1983) implies not only an infinite initial acceleration but also infinite initial velocity at the tip of the hyperboloid. We have checked that the velocity of the jet in our simulation ($\epsilon \sim O(1)$) at the instant of such peak acceleration is significantly larger than a typical gravity based velocity scale $\sqrt {gl}$, $l$ being the displacement of the interface at $\hat {r}=0$. These justify the relevance of the model of Longuet-Higgins (Reference Longuet-Higgins1983) for our case and the comparison with this model, which we describe next.

For comparing our simulations with the prediction in (3.20), we choose $U$ and $L$ to be equal to the jet velocity and the jet width, respectively, in the simulation, both values obtained at $\hat {t}=\hat {t}_0$. Here, $U$ is estimated from the (red) parabolic fits at $\hat {t}_0$ in figure 10(a,c,e,g) and $L$ is a measure of the width of the jet at the same time instant. The angle $\theta$ is measured in our simulations by visually fitting straight lines to the jet tip as depicted in figure 11 (also see Appendix D). The inset to this figure shows a hyperboloid of two sheets with parameter values indicated in the caption to this figure. Note that the angle $2\theta$ is the angle made by the two straight lines which represent the asymptote to the hyperboloid at large distance. Hence, while generating these fits, we have excluded a small region at the jet tip where the slope of the interface is nearly zero or has a droplet-like bulge, whose dimensions are close to the capillary length, see e.g. figure 16(a) in Appendix D. This is consistent with Longuet-Higgins (Reference Longuet-Higgins1983, figure $3$). Figure 10(b,d,f,h) present the temporal evolution of the angle $\theta$ with time $\tau \equiv {U(\hat {t} - \hat {t}_0)}/{L}$. The L-H in these figure legends indicates the prediction in (3.20) by Longuet-Higgins (Reference Longuet-Higgins1983). It is observed that for $\epsilon \geq 1.3$, our jets evolve in the qualitative manner predicted by (3.20), the agreement getting better with increasing $\epsilon$. Despite this apparent good agreement, we must note that our angle estimation has a degree of subjectivity associated with our visual choice of asymptotes (see sample asymptotes provided in Appendix D). Consequently, significant subjectivity persists in the degree of agreement that may be inferred between the jets in our simulations with those of Longuet-Higgins (Reference Longuet-Higgins1983). This is partly due to the formation and ejection of droplet-like features from the jet tip in our simulations (without surface tension), which render accurate estimation of jet acceleration challenging. In addition, the jet of Longuet-Higgins (Reference Longuet-Higgins1983) is shaped as a hyperboloid at all times, in contrast to our jet which arises as a superposition of Bessel modes. This implies that the model of Longuet-Higgins (Reference Longuet-Higgins1983) devoid of gravity and capillarity is, at best, only a local approximation to our jet. This local region of applicability where asymptote angles may be defined has been performed visually in the present case, as illustrated in Appendix D. Importantly, we have neglected droplets at the jet tip, when their size is close to the air–water capillary length for consistency with the exclusion of surface tension in the simulations. Subject to these varied approximations, one may thus conclude only a rough qualitative agreement of our jets at large $\epsilon$, with the theory of Longuet-Higgins (Reference Longuet-Higgins1983).

Figure 11.

(Inset) A hyperboloid of two sheets symmetric about the vertical $z$ axis with parameters $a=50$, $b=50$ and $c=200$, and having the equation $-{\hat {x}^2}/{a(\hat {t})^2} - {\hat {y}^2}/{b(\hat {t})^2} + {\hat {z}^2}/{c(\hat {t})^2} - 1 =0$. Only the lower sheet ($z < 0$) is depicted here and serves as a reference for measuring the jet angle $2\theta$, seen in simulations. The angle $2\theta$ (Longuet-Higgins Reference Longuet-Higgins1983) is defined as the angle between the two asymptotes to the hyperboloid of Longuet-Higgins (Reference Longuet-Higgins1983). Also note that the small bulge around the jet tip is approximately $3$ mm in radius which is close to the air–water capillary length scale and can eject droplets at later time. Hence, a small region close to the tip of the jet is excluded in our estimation of asymptotes. The main figure depicts the interface $\hat {\eta }$ obtained from numerical simulation with $\epsilon =1.5$, a few time steps after jet inception occurs (i.e. a few time steps after $\hat {t}_0$). Further asymptotes are also provided in Appendix D for $\epsilon =1.2$, 1.3, 1.5 and $1.6$.

We also query the relation between our study and that of earlier ones like Mack (Reference Mack1962) and Tsai & Yue (Reference Tsai and Yue1987), both of which have investigated time periodic solutions in this geometry. We emphasise that as we are solving an IVP, our analysis needs to take into account free modes (coming from the complementary function in the solution to (3.8a)) as well as particular integrals which produce bound components. The presence of the former leads to aperiodicity as the frequencies of the free modes are not rational multiples of each other. The nonlinearly produced free modes, when they have the same temporal phase as the primary mode and frequencies close to the primary mode, can reinforce the primary mode when the latter reaches its maximum at the symmetry axis and this provides a rough physical picture of the origin of the overshoot seen in simulations. However, as we are solving an IVP, there are also bound components in the solution (free modes and the bound components are both necessary to satisfy the initial conditions) and the precise role of these in the overshoot requires further study, and is proposed as future work.

In conclusion, we note that here our motivation has been to study a simple initial condition, for which jet formation occurs at sufficient perturbation amplitude and which permits analytical foray into the nonlinear regime. Due to the solution being specific to the initial condition, all our conclusions need not be generic or applicable to jets in scenarios mentioned in the introduction. However, our analysis yields useful insights towards what is necessary to generate an overshoot, an essential characteristic of a jet. An important attribute of our analytical solution is that we do not impose time periodicity and our initial condition excites a single mode only. In future studies, we would like to compare our present approach to those of Mack (Reference Mack1962), Dalzell (Reference Dalzell1967), Tsai & Yue (Reference Tsai and Yue1987) who also studied surface waves in cylindrical geometry, particularly towards understanding the role of free modes vis-a-vis bound components in jet formation and to elucidate the role of nonlinearity further. We conclude by highlighting the theoretical contribution by Dalzell (Reference Dalzell1999) who studied quadratic interactions for two surface gravity modes in deep water (Cartesian geometry) and, particularly, also investigated the standing wave limit; Dalzell (Reference Dalzell1999, ($32$) and ($33$)) provide time-periodic expressions. An investigation of jets in Cartesian geometry using our IVP approach is currently underway and comparisons of the same, vis-a-vis these earlier insightful nonlinear approaches, will be presented in a forthcoming publication.