Jet from a very large, axisymmetric, surface-gravity wave

Abstract We demonstrate that gravity acting alone at large length scales (compared to the capillary length) can produce a jet from a sufficiently steep, axisymmetric surface deformation imposed on a quiescent, deep pool of liquid. Mechanistically, the jet owes it origin to the focusing of a concentric, surface wave towards the axis of symmetry, quite analogous to such focusing of capillary waves and resultant jet formation observed during bubble collapse at small scales. A weakly nonlinear theory based on the method of multiple scales in the potential flow limit is presented for a modal (single-mode) initial condition representing the solution to the primary Cauchy–Poisson problem. A pair of novel, coupled, amplitude equations are derived governing the modulation of the primary mode. For moderate values of the perturbation parameter $\epsilon$ (a measure of the initial perturbation steepness), our second-order theory captures the overshoot (incipient jet) at the axis of symmetry quite well, demonstrating good agreement with numerical simulation of the incompressible, Euler equation with gravity (Popinet 2014, Basilisk. http://basilisk.fr) and no surface tension. We demonstrate that the underlying wave focusing mechanism may be understood in terms of radially inward motion of nodal points of a linearised, axisymmetric, standing wave. This explanation rationalises the ubiquitous observation of such jets accompanying cavity collapse phenomena, spanning length scales from microns to several metres. Expectedly, our theory becomes inaccurate as $\epsilon$ approaches unity. In this strongly nonlinear regime, slender jets form with surface accelerations exceeding gravity by more than an order of magnitude. In this inertial regime, we compare the jets in our simulations with the inertial, self-similar, analytical solution by Longuet-Higgins (J. Fluid Mech., 1983, vol. 127, pp. 103–121) and find qualitative agreement with the same. This analysis demonstrates, from first principles, an example of a jet created purely under gravity from a smooth initial perturbation and provides support to the analytical model of Longuet-Higgins (J. Fluid Mech., 1983, vol. 127, pp. 103–121).


Introduction
The dynamics of waves and jets at geophysical scales has been of recurrent interest.The study by Miles (1968), for example, refers to an insightful quote from Van Dorn (1968): '. . . the concentric circulate ridges that surround the crater Orientale at lat. 20 • S and long.95 • W on the moon may have been initiated as gravity waves on a viscous liquid under the impact of a meteorite'.Miles (1968) subsequently presents an analytical solution to the viscous, linear, Cauchy-Poisson problem motivated by the need to validate this hypothesis.The inviscid, irrotational, linear, Cauchy-Poisson initial value problem (IVP) and its solution (Cauchy 1816;Poisson 1818) governing linear wave evolution on a pool of deep liquid was reported more than one hundred and fifty years before the viscous extension to the same by Miles (1968) for a localised surface perturbation.Our topic of interest here is the gravity-induced collapse of a large cavity and a resultant (Worthington-like) jet which may be ejected during such collapse, at scales large enough to be geophysically relevant.As will be seen below, the Cauchy-Poisson (CP hereafter) linear solution provides a useful reference, against which the formation of this jet may be compared.Throughout this study, when we characterise lengths to be large, this is relative to both the capillary length scale (≈2.7 mm) as well as the amplitude of the initial perturbation.Relevant non-dimensional numbers are quoted appropriately.
Jets can often be seen associated with waves at geophysical scales.For example, numerical models of asteroid impact often depict the generation of a thin ejecta sheet shooting upwards due to impact; see Range et al. (2022, figure 1(a)).This sheet may be treated as being similar to a Worthington jet (in two dimensions) and can be modelled as an infinitely tall spike of cross-sectional area A 0 .Within the CP linear framework, one can ask what surface waves may be produced due to this ejecta sheet?As explained by Lamb (1924) (sections 238-241) for the infinite depth case and also by Pidduck (1912), the air-water interface represented by the variable η(x, t) (x is the direction of wave propagation and t is time) becomes wavy due to such an initial disturbance and evolves self-similarly under acceleration due to gravity g.The evolution is predicted by the CP solution as πx η(x, t) where the functional form of f (•) in (1.1) is provided by Lamb (1924).In contrast to the two-dimensional description above, our interest here lies in axisymmetric, surface gravity waves and jets produced therein, in cylindrical geometry.In cylindrical coordinates, for an initial surface perturbation η(r, t = 0) (with φ(r, ẑ = 0, t = 0) = 0, where r is the radial coordinate, η represents the perturbation of the free surface and φ is the perturbation velocity potential) on a radially unbounded pool of infinite depth, the classical solution to the linearised CP problem predicts (Debnath 1994) η(r, t) = ∞ 0 dk kJ 0 (kr)H[ η(r, 0)] cos( gkt). (1.2) Here, H[•] represents the zeroth-order Hankel transform, J 0 (•) is the zeroth-order Bessel function of the first kind and k represents wavenumber.When the initial surface perturbation, η(r, t = 0), has a finite width (i.e. has a characteristic length scale), the CP integral in (1.2) needs to be evaluated numerically and, unlike (1.1), does not lead to self-similar evolution of the interface.We are specifically interested here in initial conditions (i.e.η(r, 0)) which lead to jets, yet are simple enough to permit analytical Here r is the radial coordinate, â0 is measure of the initial hump amplitude and d, a measure of its width.A jet (red curve) is formed due to the relaxation of the cavity, rising sharply upwards at the symmetry axis, r = 0. (b) A tiny bubble (blue) of radius rc = 8.57 × 10 −3 cm (Bond number Bo ≡ ρgr 2 c /T = 0.001, where T is surface tension, ρ is the water density and g is gravity) at the air-water interface whose collapse produces a jet (red curve) (Gordillo & Rodríguez-Rodríguez 2019;Duchemin et al. 2002).In (b), η ≡ η/r c and r ≡ r/r c (variables with hats are dimensional) where rc is the radius of the bubble.Due to the length scale disparity between (a) and (b), the cavity relaxation is dictated nearly entirely by gravity in (a) and almost entirely by surface tension in (b).The dimensional width of the cavity in (a) (blue curve) is ≈60 cm compared with the bubble diameter (blue curve) in (b) which is 2r c ≈ 1.714 × 10 −2 cm.The (non-dimensional) time corresponding to the jet (red curve) in (a) is t = (t g/ d)/2π = 0.717.In (b), the non-dimensional time corresponding to the jet (red curve) is t = (t T/ρr 3 c )/2π = 0.0963.Note that the simulation in (a) is started from a crest (inset of figure 2a) unlike the trough-like initial shape of the bubble in (b).Thus, a difference of approximately 0.5 is expected in the non-dimensional time of jet formation between (a) and (b).
progress into the nonlinear regime.For such initial conditions, we aim to go beyond the linear regime dictated by (1.2) and solve the weakly nonlinear, axisymmetric CP initial-value problem to understand what aspects of jet formation are contained in such an analytical solution.
1.1.Jets from initial surface deformations Figure 1 depicts two such jets obtained from different initial conditions in axisymmetric geometry.Figure 1(a) depicts a cavity at an air-water interface (blue curve) which has been generated starting from an axisymmetric hump (crest, see inset of figure 2a), by solving the incompressible Euler equations with gravity (and zero surface tension) using the open-source code Basilisk (Popinet 2014).The cavity relaxes producing a Worthington-like jet, indicated by the red curve in this figure.The initial hump is given by the formula η(r, 0) = â0 exp(−r 2 / d2 )[1 − (r/ d) 2 ], â0 , d > 0 with characteristic height â0 and width d (see inset of figure 2a). Figure 1(b), however, also depicts a Worthington jet (red) but is obtained from the relaxation of a nearly spherical cavity (blue), representative of a bubble bursting at an air-water interface (flat line in blue) for low Bond number with zero viscosity.One notes the generation of jets in both situations despite the nearly three orders of magnitude length scale separation and different initial conditions.These jets share qualitatively similar features: the one in figure 1(a) rises significantly beyond its maximum initial amplitude (unity in non-dimensional scale).Similarly, for the jet in figure 1 paragraph, the droplet ejection velocity from such a jet can be up to twenty times the capillary velocity based on the initial bubble radius.
The initial condition which is used to generate the jet in figure 1(a) was proposed by Miles (1968) and represents a volume conserving, free-surface deformation.The linear dimensions of this perturbation (see figure 3 captions for parameter values) was chosen in the simulation to be far greater than the capillary-gravity length (2.7 mm) for air-water.The cavity is thus expected to relax dominantly under the influence of gravity with surface tension playing a relatively insignificant role, particularly in the early stage of relaxation, thereby justifying our neglect of surface tension in the simulation.When the non-dimensional aspect ratio of the cavity β ≡ â0 / d0 1, the collapse is a linear process governed by the CP solution in (1.2).This is seen from the excellent agreement in figure 2(a,b) between numerical simulations and linearised prediction from solving the integral in (1.2) numerically.In contrast, when β is increased to O(1), the cavity collapse once again generates a bump earlier, but which now evolves into a thin, sharply shooting jet.This is shown in figure 3(b).As noted earlier, this jet rises significantly higher (≈20 cm s, see figure 3b) compared with the initial hump height â0 = 14.6 cm.In this case, the temporal evolution is distinctly nonlinear as seen from the rather poor match observed between the numerical simulations and the CP linear solution (1.2) (indicated as 'Lin' in the figure legend).

Literature survey: initial value problems (IVPs)
These numerical results, presented in figure 3, emphasise the need for solving the nonlinear CP problem (i.e.predictions beyond (1.2)) for a prescribed surface deformation η(r, 0) and zero surface potential φ(r, ẑ = 0, t = 0) = 0, particularly when one seeks a first-principles understanding of the generation mechanism of the jet seen in figure 3(b).The theory literature on the solution to such IVPs particularly in axisymmetric coordinates is quite sparse with the bulk of the rich analytical work being focused on understanding of time-periodic, nonlinear, standing (Strutt 1915) or travelling wave solutions (Stokes 1847(Stokes , 1880)).Commencing from these seminal studies (Stokes 1847;Strutt 1915), a rich literature has since developed on finite amplitude, time-periodic surface waves, both in two-dimensional coordinates (Penney et al. 1952;Taylor 1953;Tadjbakhsh & Keller 1960;Fultz 1962;Schwartz & Whitney 1981;Schwartz & Fenton 1982;Dalzell 1999) as well as cylindrical, axisymmetric coordinates (Mack 1962;Tsai & Yue 1987).Additionally, several studies of the stability of these finite amplitude solutions (Benjamin & Feir 1967;Mercer & Roberts 1992) have also been reported.
Here our focus, however, is not on these aforementioned time-periodic solutions but rather on those initial surface deformations which generate a sharply shooting jet, viz.one which rises significantly higher than its parent perturbation, similar to the one seen in figure 3(b).An additional requirement is that the prescribed initial surface deformation η(r, 0) should be simple enough to permit analysis in the nonlinear regime, at least perturbatively.The initial condition presented by Miles (1968) and discussed earlier in figures 2 and 3 excites a continuum of modes (radially) initially.Extending the surface response due to such an initial perturbation beyond the CP linear regime described by (1.2) is technically demanding; this is mainly due to the necessity of accounting for interactions of a continuum of modes, interacting quadratically in a radial, axisymmetric geometry.In further analysis, we opt instead to study an alternative and apparently simpler initial condition in a radially confined geometry.The confinement assumption is mainly for analytical ease as it causes the radial part of the spectrum to be discrete instead of continuous.
In the initial condition that we study here, the free surface of a pool of very large depth (compared with the wavelength of the surface perturbation) is deformed at time t = 0 as a single (radial) eigenmode to the cylindrical, axisymmetric, Laplacian operator, i.e. the zeroth-order Bessel function.This surface deformation thus has the form η(r, 0) = â0 J 0 (l q (r/ R0 )), with â0 being the initial perturbation amplitude, r the radial coordinate, R0 being the domain radius and l q (q = 1, 2, . ..) a root of the Bessel function J 1 (•), this being necessary to satisfy no-penetration at the domain boundary.The length scales are chosen appropriately: the domain radius R0 as well as the width of the initial perturbation around the symmetry axis (≈2π R0 l −1 q ) are chosen to be quite large, i.e. 6 m and ≈0.4 m, respectively, compared with the capillary-gravity length scale of 2.7 mm for the air-water interface (see figure 4).At sufficiently large â0 (>18 times the air-water capillary length, see table 1), numerical simulation of the incompressible Euler equation with gravity (and no surface tension) with this initial surface deformation reveals the formation of a sharply shooting jet at the symmetry axis (r = 0) in figure 4. We notice the qualitative similarity of this jet to the one seen earlier in figure 3(c), i.e. both rise far beyond their respective initial, maximum perturbation height â0 .We will understand the mechanism of generation of this jet from first principles here by solving the corresponding IVP.To place our current study in perspective, we briefly summarise the literature on initial-value problems below.The width of the initial hump is ≈40 cm while the width of the radial domain is R0 = 600 cm, only a part of which is shown here.Note that we depict this as a large amplitude perturbation as its non-dimensional steepness ≡ a 0 l q / R0 = 1.5 > 1.The resultant jet rises >60 cm at its maximum, far exceeding the initial amplitude â0 = 8.13 cm (i.e.approximately seven times the initial perturbation amplitude).Note the formation of droplets from the tip of the jet close to t = 0.5 s (non-dimensional t = 6.28, see (3.1a-c) for temporal scale).The jet profile at t = 0.58 s (t = 7.81) is shown in the inset with vertical and horizontal axes non-dimensional as per (3.1a-c).Parameters correspond to Case 5 in table 1.For this simulation, = 1.5, l q = l 35 = 110.74.The drops visible at the tip of the jet in the inset have been manually removed in the main figure.
Case There are two classes of initial conditions for the solution to the linear CP problem.These are (a) waves are generated from fluid at rest and a specified initial surface deformation (termed the primary CP problem recently by Tyvand, Mulstad & Bestehorn (2021a) and of current interest; (b) waves and fluid motion are generated from a flat surface due to an initial surface impulse (see expressions 3.2.10 and 3.2.12 in Debnath 1994).This is the secondary CP problem (Tyvand et al. 2021a).
The primary and the secondary CP problems in two-dimensional coordinates were investigated numerically using the potential flow equations by Saffman & Yuen (1979).For the latter case, the authors applied a sinusoidal, in space and time, pressure force for a certain duration on a flat interface, and reported the subsequent evolution of the free surface after this forcing was turned off.For the primary CP problem, they used the the finite amplitude solution of Price & Penney (1952) as their initial surface deformation.In both cases, they investigated waves of highest amplitude, obtaining good agreement with the well-known experiments of Taylor (1953).In two-dimensional Cartesian coordinates, the primary CP problem was subsequently numerically investigated by Longuet-Higgins (2001) and Longuet-Higgins & Dommermuth (2001b).The authors imposed a prescribed surface deformation in the shape of a circular trough, flanked by circular crests of larger radii.A trough as it collapsed in time was numerically found to generate a jet rising upwards, see Longuet-Higgins & Dommermuth (2001b, figure 7).Notably, peak accelerations exceeding 10g were observed at the base of their cavity prior to the formation of the jet.In a follow-up work, the secondary CP problem was also investigated numerically by Longuet-Higgins & Dommermuth (2001a).These authors numerically solved the IVP employing an initial condition, which corresponds to a standing wave solution to the linearised problem.The horizontal (û(x, ẑ, t)) and vertical ( v(x, ẑ, t)) components of fluid motion in their simulation were prescribed in two-dimensional coordinates, initially as û(x, ẑ, 0) = C exp(kẑ) cos(kx), v(x, ẑ, 0) = C exp(kẑ) sin(kx) with the surface being flat initially ( η(x, 0) = 0).For C 1, a linear standing wave was observed.However, at larger C = O(1), the authors reported formation of a jet-like structure at x = 0 during the downward motion of the crest (see Longuet-Higgins & Dommermuth 2001a, figure 9b).Recently, the nonlinear, secondary CP problem in two-dimensional Cartesian coordinates has also been solved analytically, employing small time expansion (Tyvand et al. 2021a;Tyvand, Mulstad & Bestehorn 2021b).These authors demonstrate significant differences between the linear and the nonlinear CP solution, with increasing amplitude of the initial pressure impulse.However, a clear jet-like structure, analogous to what was presented in the numerical solution by Longuet-Higgins & Dommermuth (2001a), is not discernable in their figure 6(c) (Tyvand et al. 2021a).Of note are also several experimental studies which have investigated the emergence of jets in set-ups which closely resemble the secondary CP problem (i.e. via application of surface impulse), mostly at capillarity dominated length scales.These include the studies by Antkowiak et al. (2007), Bergmann et al. (2008), Gordillo, Onuki & Tagawa (2020) of the tubular jet as well as several versions of the so-called Pokrowski's experiment (Lavrentiev & Chabat 1980).
The modal initial condition of interest to us in this study, i.e. η(r, 0) = â0 J 0 (l q (r/ R0 )), corresponds to the primary CP problem described above (Tyvand et al. 2021a), albeit in cylindrical, axisymmetric coordinates.This initial condition was first studied by Farsoiya, Mayya & Dasgupta (2017) to analytically solve the viscous, linear, IVP at length scales (approximately a few centimetres) chosen such that surface tension as well as gravity were equally important.It was observed in the direct numerical simulations (DNS) reported by Farsoiya et al. (2017) that by systematically increasing the perturbation amplitude â0 , a capillarity-gravity dominated jet emerged at the symmetry axis rising significantly higher than â0 .The viscous, linear theory presented by Farsoiya et al. (2017) was unable to describe the formation of this jet or even account for its inception.In a subsequent study (Basak, Farsoiya & Dasgupta 2021), the weakly nonlinear solution to the IVP within an inviscid framework and accurate up to O( 2 ) ( ≡ â0 l q / R0 ) corresponding to the same initial condition as that of Farsoiya et al. (2017) was developed.Here too the length scales of interest were similar to that of Farsoiya et al. (2017) with gravity and surface tension forces being equally strong and both forces were accounted for in the theory.This nonlinear theory (Basak et al. 2021) was a significant improvement over the linear model of Farsoiya et al. (2017) and was able to predict the inception of the jet, comparing well with inviscid simulations of Euler's equations with gravity and surface tension (Basak et al. 2021).The physical mechanism underlying this jet formation was however not presented by Basak et al. (2021).Additionally, for a given choice of l q , surface tension and gravity, the analytical theory of Basak et al. (2021) becomes singular for certain values of the domain size R0 , this being related to triadic internal resonance.It was demonstrated (Basak et al. 2021) that these singularities owed their origin to the presence of both surface tension as well as gravity in the theoretical model and represented the cylindrical counterparts of such singularities, better known in the context of Wilton ripples (Wilton 1915) in two-dimensional coordinates.Now, it is already known from experimental as well as simulation studies on bubble bursting at millimetric (Tagawa et al. 2012;Deike et al. 2018;Gordillo & Rodríguez-Rodríguez 2019) and micron length scales (Lee et al. 2011) that jets accompany such cavity collapse quite routinely at these scales, where gravity is insignificant compared with surface tension during the collapse.Motivated partly by this observation, simulations of the incompressible, Euler's equation with only surface tension (no gravity) with the aforementioned initial condition of Farsoiya et al. (2017) have been reported recently by Kayal, Basak & Dasgupta (2022).It was demonstrated clearly in these simulations that surface tension alone at sufficiently small scales can produce a jet similar to what was observed at much larger capillary-gravity length scales by Farsoiya et al. (2017) and Basak et al. (2021).Concomitantly, the weakly nonlinear solution to the IVP for this initial condition and taking into account only surface tension has also been developed by Kayal et al. (2022).Comparison of this analytical theory against numerical simulations demonstrated very good agreement (Kayal et al. 2022), additionally also shedding light into the physical mechanism at work driven by gradients of curvature.It was proven (Kayal et al. 2022) that to describe the inception of the jet in the surface tension driven case, one needs to account for nonlinear effects of curvature (the gradient of which drives the flow).To capture this accurately, it was shown that the analytical theory needs to be at least third-order accurate, i.e.O( 3).This third-order accurate solution was reported by Kayal et al. (2022) and has been found to be free from the aforementioned internal resonance related singularities by Basak et al. (2021), thereby also alleviating a practical deficiency of the capillary-gravity model studied earlier by Basak et al. (2021).
Our present study treats the converse case of Kayal et al. (2022) considering large amplitude surface waves (we typically analyse cases where the initial wave steepness > 0.5) under the influence of gravity (and no surface tension, i.e. infinite Bond number).To be consistent, we choose our length scales of interest to be typically in the metre range (radial domain size is 6 m), far greater than the air-water capillary length and this represents a scale up of nearly two orders of magnitude in length, as compared to all our earlier studies (Farsoiya et al. 2017;Basak et al. 2021;Kayal et al. 2022).Due to the length scale of our initial perturbation, far exceeding the capillary scale, we expect the relaxation of the perturbation to be dominated by gravity with surface tension playing a negligible  6 m in height, via focusing of waves at the symmetry axis of a cylindrical pool.Note that in their case, the focusing was achieved not directly via cavity collapse but via wavemakers at the periphery of a very large cylindrical tank (25 m diameter) which generated wave components arranged to be in phase at the symmetry axis, leading to a spike wave there.McAllister et al. (2022) also demonstrated that the shape of the spike wave (with the exception of a small region around the symmetry axis) could be accounted for by linear, dispersive focusing.Given our numerical observation of purely capillarity driven jets (sans gravity) for the surface deformation η(r, t = 0) = â0 J 0 (l q (r/ R0 )) at millimetric length scales in Kayal et al. (2022), we ask here if gravity acting alone, at typical length scales of tens of metres, may also produce a jet similar to the one previously seen by us (Kayal et al. 2022), employing the same initial condition.We also note here that the experimental study by McAllister et al. (2022) demonstrated recently that similar waves, purely due to gravity, may be generated via focusing.We will show here, theoretically and computationally, that the answer to our question is in the affirmative and also explain the physical mechanism at work.Analogous to the purely surface tension driven, nonlinear theory developed by us in Kayal et al. (2022), the purely gravity driven, second-order, nonlinear theory developed here is devoid of singularities, arising from internal resonance.This theory developed via multiple scale analysis is then compared with numerical simulations obtaining very good agreement.

Physical mechanism of jet formation: wave focusing
Figure 5, illustrated further in Appendix C, depicts the radially inward motion of the crests, highlighted by horizontal arrows.To understand physically this radially inward motion seen in figure 5 (Appendix C), we shall employ mass conservation via the kinematic boundary condition.Recall that the solution to the linear CP problem provided in (1.2) for the modal surface deformation, i.e. η(r, t = 0) = −â 0 J 0 (l q (r/ R0 )) and φ(r, ẑ = 0, 0 where ûr and ûz are the radial and axial components of perturbation velocity.Here a negative sign is deliberately used in the prescription for η(r, t = 0) to generate an initial trough (instead of a hump) around r = 0 resembling locally a cavity around the symmetry axis.The linearised kinematic boundary condition, i.e. (∂ η/∂t) = ûz (ẑ = 0, r, t), implies that the temporal evolution of the interface (at linear order) is determined only by the vertical component of perturbation velocity ûz at the undisturbed liquid level ẑ = 0. Thus for this initial condition, the zeros of the initial surface deformation, viz. the solution to J 0 (l q (r/ R0 )) = 0, behave as nodes (at linear order).
Let us denote the radial location of the first zero of J 0 (l q (r/ R0 )) = 0 as r * .Observe from (2.1c) that ûr (r = r * , ẑ = 0, 0 < t < T0 /4) / = 0, being a small negative value in the first quarter of the oscillation cycle with time period T0 , i.e. the early stages of cavity collapse.Thus, linear theory predicts a radially inward, horizontal velocity at the location r = r * , as seen from the negative sign of (2.1c) for r = r * , ẑ = 0 and t < T0 /4.This radially inward velocity affects the solution only at nonlinear order leading to radial inward motion of the first node, as may be seen via consideration of the nonlinear term in the kinematic boundary condition.One can now correlate this radially inward motion of the node, with the radially inward motion highlighted earlier in figure 5 (also see Appendix C).Notably, our allusion to the nonlinear term in the kinematic boundary condition implies that the physical mechanism of inward radial motion of the nodes may be understood purely from a kinematic viewpoint (using only mass conservation) without appealing to forces or pressure gradients, which drive the flow.This thus provides insight into why one should generically expect jets from such large amplitude monochromatic perturbations, independent of whether one is investigating the phenomena at capillarity driven length scales (Kayal et al. 2022) or far larger, gravity driven scales as is our present case.
Our analytical argument above complements the well accepted viewpoint that these jets are a consequence of flow focusing and thus appear in a wide range of situations.Examples include large cavity collapse (Ghabache, Séon & Antkowiak 2014), impact-driven liquid jets (Antkowiak et al. 2007), needle-less injection (Tagawa et al. 2012) or collapse in granular jets (Lohse et al. 2004).

Weakly nonlinear theory: multiple scale approach
In this section, we develop an O( 2 ) weakly nonlinear theory employing the method of multiple scales and ≡ â0 l q / R0 as perturbation parameter to predict from first principles how the jet develops.Note that in our earlier studies (Basak et al. 2021;Kayal et al. 2022), we have used the closely related method of strained coordinates (Lindstedt-Poincare technique).Our usage of multiple scales here is motivated by the fact that at sufficiently early time, the effect of nonlinearity should be small and the interface should behave as a linear standing wave.At a longer time window, however, the effect of nonlinearity becomes apparent including the production of free and bound mode components absent initially.We derive here the amplitude equations governing the modulation of the primary mode as this was not obtained earlier in the strained coordinate method (Basak et al. 2021;Kayal et al. 2022).As our study is restricted to jets on a water pool only, for pure water, the corresponding viscous-gravity length scale (ν 2 /g) 1/3 1 mm is quite small compared to the length scales of interest to us here.Consequently, we resort to the inviscid, irrotational, potential flow approximation as is standard practise in analysing large-scale surface gravity waves (McAllister et al. 2022).These potential flow equations are non-dimensionalised using the following scales: φ, (3.1a-c) which yields the following set of equations written in axisymmetric, cylindrical coordinates.Here (non-dimensional) φ represents the disturbance velocity potential and (non-dimensional) η represents the disturbed free surface.The governing equations and boundary conditions are The number l q (q ∈ Z + ) in the initial condition η(r, t = 0) = â0 J 0 (l q (r/ R0 )) satisfies J 1 (l q ) = 0.The chosen integer value of q in ≡ l q (â 0 / R0 ) is related to the number of extremas of J 0 within the radial domain R0 .The numerical value of R0 l −1 q provides a rough measure of the wavelength of the initial perturbation while nonlinearity is controlled by the magnitude of the perturbation amplitude â0 relative to R0 l −1 q .To minimise the effect of the confining radial boundary R0 on the cavity collapse process, we require R0 l −1 q R0 and this is ensured by choosing l q 1.In this study, we choose q = 35 implying l 35 = 110.74.The variables φ, η and t are expanded in a power series in ( 1) for finite l q .Employing multiple scale analysis, we replace the temporal dependencies in all dependent variables as T n ≡ n t (up to second order).Thus, we have up to second order, where T 0 ≡ t and T 2 ≡ 2 t.After performing a Taylor series expansion of (3.2b) and (3.2c) about the unperturbed interface z = 0 and thereafter substituting expansions (3.3) and (3.4) in (3.2a-j), we extract the following set of linear equations governing φ i (r, z, T 0 , T 2 ) and η i (r, T 0 , T 2 ) at the ith order for all i = 1, 2, 3, . ... Hence, at every O( i ) and using the short-hand symbol D n ≡ ∂/∂T n , we have and Here, δ 1i is the Kronecker delta while M i (r, T 0 , T 2 ) and N i (r, T 0 , T 2 ) are nonlinear terms involving products of lower-order solutions of φ i and η i with M 1 (r, T 0 , T 2 ) = N 1 (r, T 0 , T 2 ) = 0. Due to this at the linear order, we have a homogeneous set of equations.
The expressions for M 2 (r, T 0 , T 2 ), N 2 (r, T 0 , T 2 ) as well as M 3 (r, T 0 , T 2 ), N 3 (r, T 0 , T 2 ) (necessary for eliminating resonant forcing) are provided in the Appendix.As all φ i satisfy the Laplace equation (3.5a), we expand φ i and η i at every order in a Dini series (Mack 1962) in the following form (note that q is a given fixed integer used to indicate the primary Bessel mode that is excited initially): and By construction, each term in the expansion in (3.6) satisfies the Laplace equation (3.5a), while each term in (3.7) satisfies the mass conservation equation (3.5g).In addition, (3.6) and (3.7) together respect the finiteness condition as well as the no-penetration and the free-edge conditions (3.5d,e, f ).Our task thus reduces to ensuring that (3.5b,c) are satisfied, and these in turn determine equations governing p ( j) i (T 0 , T 2 ) and a ( j) i (T 0 , T 2 ) in the expansions (3.6) and (3.7).These equations have to be solved subject to initial conditions (3.5h,i,j).Substituting (3.6) and (3.7) into (3.5b,c),taking inner products with J 0 (α n,q )r dr and using the orthogonality relations for Bessel functions (see the supplementary material of Basak et al. (2021) for the relevant orthogonality relations), we obtain the following equations at each O( i ): Here, ω 2 j,q ≡ α j,q is the non-dimensional form of the dispersion law for pure gravity waves on a free surface.Equation (3.8a) is a second-order, partial differential equation which is solved subject to the following initial conditions (these are obtained by substituting (3.6) and (3.7) into initial conditions (3.5h,i,j)): The expression for a ( j) i (T 0 , T 2 ) may then be obtained from (3.8b), knowing the expression for p ( j) i (T 0 , T 2 ) from the solution to (3.8a).
(3.15a,b) Equations (3.15) predict that the interface η 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 η in (3.15) has an amplitude at T 0 = 0 and after one time period at T 0 = 2π, 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 μ 1q (T 2 ) and ν 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 μ 1q (T 2 ) and ν 1q (T 2 ).To determine these equations, it becomes necessary to obtain expressions up to O( 3 ) for M 2 and N 2 as well as M 3 and 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, M 2 (r, T 0 , T 2 ) / = 0 and N 2 (r, T 0 , T 2 ) / = 0, and expressions for these are provided in the Appendix.As a result, (3.8a) governing p (i) 2 (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 ( j) 2 (T 2 ) ( j = 0, 1, 2, 3, . ..), 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( 2 ) corrections are found to be 5) n,q (μ 1q (T 2 ), ν 1q (T 2 ))]J 0 (α n,q r). (3.17) n,q 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 ξ (2)  n,q , ξ (3) n,q (μ 1q (T 2 ), ν 1q (T 2 )) and ξ ( 4) n,q (μ 1q (T 2 ), ν 1q (T 2 )) are provided in the Appendix.Similarly, in (3.17), ζ (1) n,q and ζ (2) n,q = 0 are constants while expressions for n,q (μ 1q (T 2 ), ν 1q (T 2 )) are provided in the Appendix.Note that these expressions depend on the unknowns μ 1q (T 2 ) and ν 1q (T 2 ).Equations for these two unknowns are obtained at O( 3 ) via elimination of resonant forcing terms for the primary mode.

Amplitude equation(s)
To obtain the equations governing μ 1q (T 2 ) and ν 1q (T 2 ), it is necessary to carry out the calculation till O( 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 (ω 2 q,q = 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, (D 2 0 + 1)p Expressions for M 3 (r, T 0 , T 2 ), 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 M 3 and N 3 both contain terms like D 2 η 1 and D 2 φ 1 which lead to time derivatives of μ 1q (T 2 ) and ν 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 μ 1q (T 2 ) and ν 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 Here, r 1 , r 2 and r 3 in (3.19a) and (3.19b) are nonlinear functions of μ 1q and ν 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.

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 = 0.5 from our O( 2 ) theory, comparing with the Hankel transform of the interface obtained from numerical simulations.This is done at 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).The numerical simulations have been carried out using the open source code Basilisk (Popinet 2014) employing a cylindrical domain of radius R0 = 600 cm and depth Ĥ = 300 cm.As the chosen value of q = 35, this implies that the (approximate) wavelength of the initial perturbation is far smaller than Ĥ = 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 R0 at an angle of π/2 at all times.
Figure 6 shows the comparison of linear and weakly nonlinear theory with the simulation for = 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( 2) theory.This overshoot is the precursor 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 η(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.
to the slender jet that develops at the axis of symmetry, seen for far larger values of = 1.5, please refer to figure 4. Figure 7 shows the inward movement of the first zero crossing of the interface η (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 to 0.9 (Case 3 in table 1).It is seen that although the O( 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 ∼ 0.9 represents a simulation which lies outside the range of utility of the O( 2 ) theory for describing these jets.While the theory is able to describe the inception (up to = 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 > 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. (2021) (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 .This figure demonstrates that the weakly nonlinear theory presented here may be used reliably up to ∼ 0.6 beyond which it becomes quite inaccurate.Figure 9(b) demonstrates that the displacement of the tip of the jet, viz.η(r = 0, t), displays a parabolic dependence in time, in agreement with the earlier observations of Ghabache et al. (2014).
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. 2002;Gordillo & Rodríguez-Rodríguez 2019).Predictions from our theoretical model vis-a-vis numerical Figure 8. Interface for = 0.9, t = 0.47 s (non-dimensional t = 6.33) before the maximum height is reached.The O( 2 ) theory is unable to describe the slender jet seen in simulations although it manages to get the overshoot correct qualitatively. which translates in dimensional variables to ηovershoot ∼ â3 0 / R2 0 (Ghabache et al. 2014).Our weakly nonlinear model estimates the overshoot with an error of <20 % up to ∼ 0.6.Up to this = 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 2 + 0.19 + 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. (2014, figure 2).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.

Strongly nonlinear regime, inertial regime
In the previous section, we have studied simulations for < 1 and the weakly nonlinear theory was able to describe the inception of the jet up to ≈ 0.5.We now move to the strongly nonlinear regime ( > 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 (1983) in axisymmetric, extensional flow.It will be seen that this analytical model of Longuet-Higgins (1983) provides a reasonably good approximation to the jets that arise in our simulations, when > 1.2.
Interestingly, this exact solution is singular at t = 0 where it predicts a divergent velocity (ċ ∼ t−1/3 ) and a divergent acceleration (c ∼ t−4/3 ) at the tip of the free surface shaped as a hyperboloid, viz. at ẑ = c(t).In a time window around this initial 'jerky start' to motion, the half-angle θ between the asymptotes to the hyperbola are predicted to evolve from its initial value θ( 0 where t0 is time shift that will be necessary for comparison with our simulations, but is zero for Longuet-Higgins (1983).
To compare our strongly nonlinear jets ( > 1) with the solution of Longuet-Higgins (1983), 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.η(r = 0, t) versus time for various > 1.We fit parabolas (in red and green) of the form αt 2 + O(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 ηe (0, t) = αt 2 + O(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 ηl (0, t)).The value of 2α from each of these datasets provides an approximate measure of the acceleration of the interface viz.(∂ 2 η/∂t 2 )(r = 0, 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 η(r = 0, 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 η(r ≈ 1.8, t) and is plotted on the left vertical axes of figures 10(e) and 10(g).
Note from these figures that around the time of jet formation, the value of 2α (red parabolas) climbs from 5236 cm s −2 (i.e.5×g) at = 1.2 (figure 10a) to approximately 21 684 cm s −2 (≈20×g) at = 1.6 (figure 10g).Also note that the green parabolas in these figures predict a value of 2α ranging from −965 cm s −2 for figure 10(a) to ≈ − 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 = 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 > 1, there exists a time window when gravity may be neglected in our simulations, enabling comparison with the solution of Longuet-Higgins (1983).To compare the angle θ obtained from our simulations with the prediction from (3.20), we set t0 to be the instant of time when θ is closest to the initial angle in the model of Longuet-Higgins (1983), i.e. ≈54.5 • .This definition is also consistent with recent experimental measurements of large-scale jets  by McAllister et al. (2022), who too compared their jet angle with the model of Longuet-Higgins (1983).Note that the model of Longuet-Higgins (1983) 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 ( ∼ O(1)) at the instant of such peak acceleration is significantly larger than a typical gravity based velocity scale √ gl, l being the displacement of the interface at r = 0.These justify the relevance of the model of Longuet-Higgins (1983) 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 t = t0 .Here, U is estimated from the (red) parabolic fits at t0 in 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 (1962), Dalzell (1967), Tsai & Yue (1987) 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 (1999) who studied quadratic interactions for two surface gravity modes in deep water (Cartesian geometry) and, particularly, also investigated the standing wave limit; Dalzell (1999, (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.

Conclusions
We have shown in this study, numerically as well as from a first-principles theory, that gravity driven focusing of a nonlinear surface wave towards the axis of symmetry can produce jets quite analogous to what also occurs at small scales for pure capillary waves.Our theory developed using multiple scale analysis leads to novel amplitude equations governing the modulation of the primary mode, excited initially.The solution to these equations are found to be able to describe the formation of the jet up to ∼ 0.5, capturing the overshoot qualitatively, albeit quantitative differences with numerical simulations are also detected.As is increased to nearly O(1), the weakly nonlinear theory becomes expectedly inaccurate, while simulations in this regime show slender jets shooting up with intense accelerations.Here, we find qualitative agreement of the time evolution of these jets around the symmetry axis, with a self-similar solution due to Longuet-Higgins (1983).Remarkably, this solution was derived by Longuet-Higgins (1983) for the purely inertial regime with zero gravity and we find evidence of such a time window in our simulations.In marked contrast to the pure surface tension driven jets at small scales studied recently by Kayal et al. (2022), where surface tension is a part of the self-similar evolution via Keller-Miksis scales (Keller & Miksis 1983), we find here that in the pure gravity driven case and for > 1, there is a time window when gravitational acceleration at the jet tip becomes small compared to the pressure gradient driven, inertial acceleration of the fluid.During this, the interface around the jet tip evolves self-similarly showing qualitative agreement with the exact solution provided by Longuet-Higgins (1983).Our study discusses reasons for why such an agreement may be expected and also explains the physical mechanism for jet formation based on mass conservation, thus independent of dynamics of the problem.These, in turn, provide insights into why these jets may be generically observed accompanying cavity collapse phenomena across scales, spanning several orders of magnitude (Lee et al. 2011;McAllister et al. 2022).
In conclusion, we note a remark by Longuet-Higgins (2001) (page 496, first paragraph) who insightfully observed that in axisymmetric geometry, the flow focusing responsible for jet formation would be very intense.We hypothesise that accelerations significantly more than those reported here may be momentarily generated in this geometry by increasing the value of further.Finally, we note that the theory developed here may be viewed as a weakly nonlinear solution to the primary CP problem in axisymmetric coordinates, which provides insights into development of sharply shooting jets in a radial confined geometry.Our study has several implications for such large-scale jets generated, for example, in the ocean as well as in other geophysical contexts such as impact craters (Lherm et al. 2022).The plot presents the Hankel transform of the analytical solution, i.e. η 1 + 2 η 2 for = 0.5 (Case 2 in table 1) at t = 0.483 s when the interface reaches nearly its maximum height (see figure 6d).For comparison, the Hankel transform of the interface obtained from numerical simulation of Euler's equation at this time instant is also presented.It is seen that there is no significant potential energy in modes beyond n = 70, thereby justifying truncation of the infinite series beyond this n.The inset is a blowup of the region highlighted in green and shows gravitational potential energy present in modes with n ≤ 70. where dr rJ ν 1 (α m 1 ,q r)J ν 2 (α m 2 ,q r)J ν 3 (α m 3 ,q r) . . ., {m 1 , m 2 , . ..} ∈ Z + .(A15)

Appendix B: Hankel transformation
Figure 12 shows the Hankel transformation of the interface η(r, t).We have used the definition of Basak et al. (2021) and it is defined as follows and is used in figure 12:

Figure 1 .
Figure 1.Jet formation due to cavity relaxation at length scales separated by approximately three orders of magnitude.(a) A cavity (blue) at an air-water interface generated numerically (Popinet 2014) from an initial hump (see inset of figure 2a) of the form η ≡ η(r, 0)/â 0 = exp(−r 2 )[1 − r 2 ] with r ≡ r/ d, non-dimensional parameters β ≡ â0 / d ≡ 0.584 and Bond number Bo = ∞ (zero surface tension).Here r is the radial coordinate, â0 is measure of the initial hump amplitude and d, a measure of its width.A jet (red curve) is formed due to the relaxation of the cavity, rising sharply upwards at the symmetry axis, r = 0. (b) A tiny bubble (blue) of radius rc = 8.57 × 10 −3 cm (Bond number Bo ≡ ρgr 2 c /T = 0.001, where T is surface tension, ρ is the water density and g is gravity) at the air-water interface whose collapse produces a jet (red curve)(Gordillo & Rodríguez-Rodríguez 2019;Duchemin et al. 2002).In (b), η ≡ η/r c and r ≡ r/r c (variables with hats are dimensional) where rc is the radius of the bubble.Due to the length scale disparity between (a) and (b), the cavity relaxation is dictated nearly entirely by gravity in (a) and almost entirely by surface tension in (b).The dimensional width of the cavity in (a) (blue curve) is ≈60 cm compared with the bubble diameter (blue curve) in (b) which is 2r c ≈ 1.714 × 10 −2 cm.The (non-dimensional) time corresponding to the jet (red curve)

Figure 2 .Figure 3 .
Figure 2. Linear evolution of a cavity from a hump of the initial form η ≡ η(r, 0)/â 0 = exp(−r 2 )[1 − r 2 ], r ≡ r/ d.In CGS units, â0 = 0.54, d = 25, β ≡ â0 / d = 0.0216 1.(a) t = 0.32 s when a cavity is formed from the initial perturbation (shown in inset).Lin is the linear prediction obtained from expression (1.2).(b) Downward motion of a bump formed at the symmetry axis.Note the very good agreement between simulations and the linear CP prediction.The physical parameters are chosen corresponding to air and water.(a) t = t g/ d/2π = 0.314 and (b) t = t g/ d/2π = 0.717.

Figure 4 .
Figure4.Development of a jet starting from a surface deformation in the shape of single Bessel mode at t = 0, i.e. η(r, 0) = â0 J 0 (l 35 (r/ R0 )).The width of the initial hump is ≈40 cm while the width of the radial domain is R0 = 600 cm, only a part of which is shown here.Note that we depict this as a large amplitude perturbation as its non-dimensional steepness ≡ a 0 l q / R0 = 1.5 > 1.The resultant jet rises >60 cm at its maximum, far exceeding the initial amplitude â0 = 8.13 cm (i.e.approximately seven times the initial perturbation amplitude).Note the formation of droplets from the tip of the jet close to t = 0.5 s (non-dimensional t = 6.28, see (3.1a-c) for temporal scale).The jet profile at t = 0.58 s (t = 7.81) is shown in the inset with vertical and horizontal axes non-dimensional as per (3.1a-c).Parameters correspond to Case 5 in table 1.For this simulation, = 1.5, l q = l 35 = 110.74.The drops visible at the tip of the jet in the inset have been manually removed in the main figure.
role.The motivation to search for jets at such scales comes from the rogue wave literature where a very interesting recent experimental observation of McAllister et al. (2022) reported a 'spike wave'.McAllister et al. (2022) generate a spike wave which rises up to

Figure 5 .
Figure 5. Temporal evolution and development of a jet at the symmetry axis.This is from simulations with = 1.5, l q = l 35 = 110.74(Case 5 in table 1).The arrows indicate the local direction of motion of the interface at that time instant.The inset is a non-dimensional representation of the main figure.
η(r, 0) = J 0 (r), ∂η ∂t (r, 0) = 0, φ(r, z, 0) = 0. (3.2h,i, j) Equation (3.2a) is the Laplace equation in cylindrical axisymmetric coordinates, (3.2b,c) are the constant pressure condition (from Bernoulli's equation) and the kinematic boundary condition, respectively, both at the free surface.Equation (3.2d) is the finiteness of the velocity potential at infinite depth (deep water limit is assumed for simplicity).Equation (3.2e, f ) are the no-penetration and free-edge boundary conditions at the outer radial boundary r = R0 .Equation (3.2g) is the overall mass conservation while (3.2h,i,j) are initial conditions corresponding to the primary CP problem.

Figure 6
Figure 6.(a-d) Comparison of the interface profile η/ at various instants with linear and weakly nonlinear, analytical solutions.Solid red, O( 2 ) weakly nonlinear theory; solid brown, O( ) linear theory and numerical simulation (Sim, solid blue line) for = 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 (≈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) t = 0 (non-dimensional t = 0), (b) t = 0.266 s (t = 3.58), (c) t = 0.337 s (t = 5.07) and (d) t = 0.483 s (t = 6.5).

Figure 7 .
Figure7.Radial inward motion of the first zero crossing between figures 6(b) and 6(c).Note from figure6that 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 figure6.As noted earlier, the inward motion of the first zero crossing of η(r = r * , t) may be interpreted as radial inward focusing of the two humps shown earlier in figure5.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 9 .
Figure 9. (a) Comparison of the maximum height ( ηmax /â 0 ) at the symmetry axis, as a function of increasing .Note that η overshoot ≡ ηmax /â 0 − 1 = O(2) which translates in dimensional variables to ηovershoot ∼ â3 0 / R2 0(Ghabache et al. 2014).Our weakly nonlinear model estimates the overshoot with an error of <20 % up to ∼ 0.6.Up to this = 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 2 + 0.19 + 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 ofGhabache et al. (2014, figure  2).

Figure 12 .
Figure12.Justification of the neglect of terms in the infinite series beyond n = 70, i.e. 2q in (3.16) and (3.17).The plot presents the Hankel transform of the analytical solution, i.e. η 1 + 2 η 2 for = 0.5 (Case 2 in table 1) at t = 0.483 s when the interface reaches nearly its maximum height (see figure6d).For comparison, the Hankel transform of the interface obtained from numerical simulation of Euler's equation at this time instant is also presented.It is seen that there is no significant potential energy in modes beyond n = 70, thereby justifying truncation of the infinite series beyond this n.The inset is a blowup of the region highlighted in green and shows gravitational potential energy present in modes with n ≤ 70.
Radial inward focusing of surface waves is depicted by inward arrows in figure13(a-d).

Figure 13 .Figure 14 Figure 15 Figure 16 Figure 17 Figure 18 .
Figure 13.Radial inward focusing corresponding to the simulation presented in figure 5. (a-d) Same simulation as figure 5, providing a closer view of the inward focusing of the concentric, large, amplitude surface-gravity wave towards the symmetry axis.The inward motion of the crests (in blue) leads to the eventual emergence of the jet at t ≈ 0.389 s.(a) t = 0.285, (b) t = 0.35 s, (c) t = 0.38 s and (d) t = 0.39 s.