# Focusing of concentric free-surface waves

Lohit Kayal1  Vatsal Sanjay2  Nikhil Yewale1  Anil Kumar1  Ratul Dasgupta1
<EMAIL_ADDRESS>1Chemical Engineering, Indian Institute of
Technology, Bombay, India 400 076,2Physics of Fluids Group, Max Planck Center
for Complex Fluid Dynamics, Department of Science and Technology, and J. M.
Burgers Centre for Fluid Dynamics, University of Twente, P. O. Box 217, 7500
AE Enschede, The Netherlands

###### Abstract

Gravito-capillary waves at free-surfaces are ubiquitous in several natural and
industrial processes involving quiescent liquid pools bounded by cylindrical
walls. These waves emanate from the relaxation of initial interface
distortions, which often take the form of a cavity (depression) centred on the
symmetry axis of the container. These surface waves reflect from the container
walls leading to a radially inward propagating wave-train converging
(focussing) onto the symmetry axis. Under the inviscid approximation and for
sufficiently shallow cavities, the relaxation is well-described by the
linearised potential-flow equations. Naturally, adding viscosity to such a
system introduces viscous dissipation that enervates energy and dampens the
oscillations at the symmetry axis. However, for viscous liquids and deeper
cavities, these equations are qualitatively inaccurate. In this study, we
elucidate a modal approach to study the initial-value problem for concentric
gravito-capillary waves generated on a free-surface for inviscid as well as
viscous liquids. For a sufficiently deep cavity, the inward focusing of waves
results in large interfacial oscillations at the axis, necessitating a second-
order nonlinear theory. We demonstrate that this theory effectively models the
interfacial behavior and highlights the crucial role of nonlinearity near the
symmetry axis. Contrary to expectations, the addition of slight viscosity
further intensifies the oscillations at the symmetry axis. This finding
underscores the limitations of the potential flow model and suggests avenues
for more accurate modelling of such complex free-surface flows.

###### keywords:

Surface waves, nonlinear waves, jet formation, wave focussing, Cauchy-Poisson
problem, viscous waves

##  1 Introduction to wave focussing

Focusing of moderate amplitude, progressive surface waves can often in turn
produce unexpectedly large waves. At oceanic scales, spatial wave focusing,
where large amplitude waves form persistently in specific regions (Torres et
al., 2022; Chavarria et al., 2018), can produce waves powerful enough to
damage or capsize ships. A famous example is the Aghulas current region
(Britannica, The Editors of Encyclopaedia, 2024) known for giant waves and
shipping accidents (Mallory, 1974; Smith, 1976). The role of current generated
refractive focusing of waves leading to the birth of such giant waves,
specifically in the Agulhas, was anticipated by Peregrine (1976) (also see
fig. $8$ in Dysthe et al. (2008) and section $2$ in White & Fornberg (1998)).
Refractive focusing of surface waves (Peregrine, 1986) has also been exploited
to design ‘lenses’ i.e., submerged structures in a water basin which focus
incoming divergent, circular waves (see fig. $1$a in Stamnes et al. (1983)),
these being motivated from wave generation of power (McIver, 1985; Murashige &
Kinoshita, 1992).

In addition to spatial focusing, spatio-temporal focusing also occurs (Dysthe
et al., 2008), where large wave amplitudes manifest at specific locations in
space, albeit briefly. Spatio-temporal focussing has obvious relevance not
only towards understanding, for example, rogue (freak) waves in the ocean
(Charlie Wood, 2020) but also to our current study (next section). The
physical mechanisms underlying spatio-temporal focussing have been
distinguished further into linear and nonlinear dispersive focusing (section
$4.2,4.3$, Dysthe et al. (2008)). Linear dispersive focussing of progressive
waves relies on constructive interference exploiting the dispersive nature of
surface gravity waves and is particularly simple to understand in the deep
water limit. For uni-directional wave packets in deep water, generated from a
wave-maker oscillating harmonically at frequency $\Omega$ at one end of a
sufficiently long wave flume, the energy propagation velocity (group velocity)
of the packet is $c_{g}=\frac{g}{2\Omega}$ where $g$ is acceleration due to
gravity. If the wavemaker frequency varies linearly from $\Omega_{1}$ to
$\Omega_{2}$ ($\Omega_{1}>\Omega_{2}$) following
$\frac{d\Omega}{dt}=-\frac{g}{2x_{f}}$ within the time interval
$\left[t_{1},t_{2}\right]$, Longuet-Higgins (1974) showed that the energy of
each wave packet emitted during this period will converge at $x=x_{f}$
simultaneously at $t=t_{f}$ (see Brown & Jensen (2001)). This focussing of
wave energy thus causes a momentary but significant increase in energy density
at $x_{f}$ manifested as a transient, large amplitude wave at that location
around time $t_{f}$. This technique has been discussed in Davis & Zarnick
(1964) and its variants have been employed extensively to generate breaking
waves in the laboratory in a predictable manner in two (Rapp & Melville, 1990)
and three dimensions (Wu & Nepf, 2002; Johannessen & Swan, 2001; McAllister et
al., 2022) as well as in other related contexts such as generation of
parasitic capillary on large amplitude waves (Xu & Perlin, 2023).

On the other hand, in non-linear dispersive focusing, the modulational
instability (Benjamin & Feir, 1967) of a uniform, finite-amplitude wavetrain
(Stokes wave) plays a crucial role. This instability can cause the wave train
to split into groups, where focusing within a group can produce a wave
significantly larger than the others (Zakharov et al., 2006). For further
details on nonlinear focusing, we refer readers to the review by Onorato et
al. (2013).

### Spatio-temporal focussing at gravito-capillary scales

Following this brief introduction to large-scale focusing, we now focus on
length scales where gravitational and capillary restoring forces are nearly
equivalent. Our study aims to achieve an analytical understanding of wave
focusing at these shorter scales. Below, we illustrate two examples where such
small-scale focusing can be readily observed.

Stuhlman Jr (1932) investigated the formation of drops from collapsing bubbles
with diameters under $0.12$ cm in water-air interfaces and $0.15$ cm in
benzene-air interfaces. He hypothesised that these drops emerged from
Worthington jets created by the collapse of the bubble cavity. However,
contemporary research identifies this as just one of two mechanisms
responsible for drop generation (Villermaux et al., 2022). The first high-
speed ( $\approx 6000$ frames per second) images of jet formation were
reported by MacIntyre (1968, 1972) (see original experiments by Kientzler et
al. (1954)). Interestingly, these studies demonstrated that the surface
ripples are created by the retraction of the circular rim of the relaxing
bubble cavity. These ripples travel towards the cavity base before the jet
emerges. In the words of MacIntyre (1972) (see abstract) “..an irrotational
solitary capillary ripple precedes the main toroidal rim transporting mass
along the surface at about $90\%$ of its phase velocity. The convergence of
this flow creates opposed jets…”. The seminal work by Duchemin et al. (2002)
of collapsing bubbles (much smaller than their capillary length scale) at a
gas-liquid interface was able to resolve this focussing process, via direct
numerical simulations (DNS) of the axisymmetric Navier-Stokes equations
without gravity. Figure 1 depicts the generation of an axisymmetric, wave-
train focussing towards the base of the bubble cavity (also the symmetry axis)
for two different Ohnesorge numbers $(Oh)$ and at a fixed Bond number $(Bo)$.
The Bond number $Bo\equiv\frac{\rho^{L}g\hat{R}_{b}^{2}}{T}$ determines the
bubble shape, and Ohnesorge number
$Oh\equiv\frac{\mu^{L}}{\sqrt{\rho^{L}T\hat{R}_{b}}}$ accounts for the ratio
of viscous to capillary forces. Here $\rho^{L},\,\mu^{L},\,T,\hat{R}_{b}$ are
the lower fluid density, lower fluid viscosity, coefficient of surface tension
and equivalent radius of the bubble respectively. We refer the readers to
Deike (2022); Sanjay (2022); Gordillo & Blanco-Rodríguez (2023) for recent
advances on study of bubble collapse and jet formation mechanisms.

Another example of axisymmetric focussing of surface waves was highlighted in
the study by Longuet-Higgins (1990), where several interesting observations
were noted. Longuet-Higgins (1990) studied the inverted conical shaped ‘impact
cavities’ seen in experiments and simulations (Oguz & Prosperetti, 1990) of a
liquid droplet falling on a liquid pool. The author compared these cavities to
an exact solution to the potential flow equations without surface tension or
gravity (Longuet-Higgins, 1983), where the free-surface (gas-liquid interface)
took the form of a cone at all time. The apex of this cone (i.e. the impact
cavity) is often seen to contain a bulge (see fig. $2$a in Longuet-Higgins
(1990)) and the formation of this was attributed to (we quote, section $6$
first paragraph in Longuet-Higgins (1990)) “a ripple on the surface of the
cone converging towards the axis of symmetry”, thus highlighting the role of
wave focussing once again. Longuet-Higgins (1990) insightfully remarked that
this convergence process would be similar to the radially inward propagation
of a circular ripple on a water surface. The interface shape could thus be
approximated as being due to the linear superposition of an initial, localised
wave packet (generated by distorting an initiallly flat surface) whose
Fourier-Bessel representation $F(k)$ ($k$ being the wavenumber) slowly varies
on a time-scale $\bar{t}$ (i.e. slow compared to the wave packet propagation
time-scale $t$). Longuet-Higgins (1990) thus posits that the shape of the
perturbed interface $\eta(r,t,\bar{t})$ may be represented as

$\displaystyle\eta(r,t,\bar{t})=\int_{\Delta
k}F(k,\bar{t}){\mathrm{J}}_{0}(kr)\exp\left(I\sigma(k)t\right)kdk,$ (1)

where ${\mathrm{J}}_{0}$ is the Bessel function, $r$ is the radial coordinate
and the spectrum of the surface perturbation $F(k,\bar{t})$ evolves slowly on
a time-scale $\bar{t}$, $I\equiv\sqrt{-1}$ and $\sigma(k)$ satisfies the
dispersion relation for capillary waves (see eqn. $6.2$ in Longuet-Higgins
(1990)). Note that if the slow variation of $F(k,\bar{t})$ over $\bar{t}$ is
supressed, eqn. 1 represents the solution to the linearised Cauchy-Poisson
problem with an initial surface distortion whose Hankel transform is $F(k)$.
Longuet-Higgins (1990) however did not report any systematic comparison of
available experimental or simulational data (Oguz & Prosperetti, 1990) with
eqn. 1 although the author anticipated that nonlinearity could become
important during the convergence; see last para in page $405$ of Longuet-
Higgins (1990).

Our current study is partly motivated by the aforementioned observations of
Longuet-Higgins (1990) and Duchemin et al. (2002) and aims at obtaining an
analytical description of spatio-temporal wave focussing at these short
scales. We seek an initial, localised surface distortion which produces a
wave-train, and whose radial convergence may be studied analytically, at least
in the potential flow limit. We refer the reader to the review by Eggers et
al. (2024) where this limit corresponding to Ohnesorge $Oh=0$ is discussed. In
the next section we present a localised initial surface distortion which is
expressable as a linear superposition of eigenmodes. It will be seen that this
distortion generates a surface wave-train which focuses towards the symmetry
axis of the container. We emphasize that the wave-trains or the solitory
ripple seen in Kientzler et al. (1954) and Longuet-Higgins (1990)
respectively, have different physical origins compared to the ones we study
here. However, following Longuet-Higgins (1990) we intuitively expect there
are aspects of their convergence which do not sensitively depend on how these
are generated in the first place.

We further develop an inviscid nonlinear theory for the focusing of a
concentric wave-train resulting from the aforementioned a priori imposed free-
surface deformation. The theory developed from first principles here has no
fitting parameters and helps delineate those aspects of focusing which may be
accounted for by linear theory compared to nonlinear features. In a series of
earlier theoretical and computational studies from our group (Farsoiya et al.,
2017; Basak et al., 2021; Kayal et al., 2022; Kayal & Dasgupta, 2023), we have
solved the initial-value problem corresponding to delocalised, initial
interface distortions in the form of a single eigenmode
$({\mathrm{J}}_{0}(kr))$ at gravity dominated large-scales (Kayal & Dasgupta
(2023)), gravito-capillary intermediate scales (Farsoiya et al., 2017; Basak
et al., 2021) and capillarity dominated small-scales (Kayal et al. (2022))
(also see the recent study in Dhote et al. (2024) for a delocalised initial
perturbation on a sessile bubble). In contrast to these studies, we study here
a localised, multimodal initial excitation. Apart from the obvious advantage
of easier experimental realisation of this initial condition (see Ghabache et
al. (2014b) for experimental realisation at gravity dominated scales), this
initial condition has the additional advantage that already at linear order, a
radially propagating concentric wave-train is obtained and one can ask how
does this converge at the axis of symmetry? In contrast, for a single mode
initial condition as in Basak et al. (2021); Kayal et al. (2022), at linear
order one obtains only a standing wave and it is necessary to proceed to
quadratic order and beyond to generate other modes via non-linear interactions
whose superposition can generate a focussing wave-train.

The manuscript is structured as follows: § 2 illustrates the time evolution of
a relaxing cavity and introduces the analytical equations for wave evolution.
§ 3 compares these analytical results with direct numerical simulations (DNS).
Finally, the paper culminates with discussions and outlook in § 4.

(a)

(b)

Figure 1: An example of capillary wave focussing obtained from direct
numerical simulations (DNS) conducted using the open-source code Basilisk
Popinet & collaborators (2013–2024). The initial cavity shape (inset in first
figure of upper and lower panels) is obtained by solving the Young-Laplace
equation with gravity to determine the shape of a static bubble at the free
surface (without its cap). In CGS units, initial bubble radius $0.075$,
surface tension $T=72$, gravity $g=981$, density $\rho^{\text{L}}=1.0$ and
$\rho^{\text{U}}=0.001$ for upper and lower fluid. Upper panel (blue)
simulations are conducted using zero viscosity for both gas (above) and liquid
(below). (Red, lower panel) simulations have dynamic viscosity
$\mu^{\text{U}}=0.0001$ and $\mu^{\text{L}}=0.01$. Axes are non-
dimensionalised using initial bubble radius. Time is non-dimensionalised using
the capillary time-scale
$t=\frac{\hat{t}}{\sqrt{\frac{\rho\hat{R}_{b}^{3}}{T}}}$. For the upper panel
$Bo\equiv\frac{\rho^{\text{L}}g\hat{R}_{b}^{2}}{T}=0.076$ and
Oh$=\frac{\mu^{\text{L}}}{\sqrt{\rho^{\text{L}}T\hat{R}_{b}}}=0$. For the
lower panel Bo$=0.076$ and Oh$=0.0043$.

##  2 Time evolution of a relaxing cavity

Figure 2: A (not to scale) cross-sectional representation of the initial
interface distortion $\hat{\eta}(\hat{r},0)$ shaped as a cavity of half-width
$\hat{b}$ and depth $\hat{a}_{0}$ in a cylinder of radius $\hat{R}$ filled
with liquid (in blue). The functional form chosen for $\hat{\eta}(\hat{r},0)$
was first proposed by Miles (1968) and represents a volume preserving
distortion. The red dotted line indicates the unperturbed level of the free-
surface of the liquid pool. The gas-liquid surface tension is $T$. Liquid
density and viscosity are $\rho$ and $\mu$ respectively, $g$ is gravity.

As shown in fig. 2, the system consists of a cylindrical container of radius
$\hat{R}$ filled with quiescent liquid (indicated in blue). As we do not model
the upper fluid in our theory, here onwards the superscript L is dropped from
the variables representing fluid properties. For simplicity of analytical
calculation, the cylinder is assumed to be infinitely deep and the gas-liquid
density ratio is fixed at $0.001$ to model air-water configuration. In our
theoretical calculations, we approximate the gas-liquid interface as a free-
surface and neglect any motion in the gas phase (although, it is modeled in
our DNS). Some of the relevant length scales are the gravito-capillary length
$l_{c}\equiv\sqrt{T/\rho g}\approx 2.7\,$\mathrm{mm}$$ and the visco-capillary
length scale $l_{\mu}\equiv\mu^{2}/\rho T\approx
0.01\,$\mathrm{\SIUnitSymbolMicro m}$$. For our chosen sizes of the initial
interface perturbation ($0.26-3.9\,$\mathrm{mm}$$), these length scales
justify the inclusion of both capillarity as well as gravity in the
theoretical calculation while neglecting viscosity at the leading order.
However, we stress that viscosity is known to have a non-monotonic effect on
wave focussing in a collapsing bubble, as demonstrated by Ghabache et al.
(2014a). Their fig. $3$ shows that the jet velocity during bubble bursting
varies non-monotonically with increasing viscosity. Thus, the fastest jets
occur not in an inviscid system but at an ‘optimal’ viscosity. In what
follows, we employ potential flow equations in our theory and do not treat the
boundary layers expected to be generated at the air-water interface and the
cylinder walls (Mei & Liu, 1973). We will address the inclusion of viscous
effects later in the study.

Before delving into the theoretical formulation, it is instructive to discuss
the phenomenology of the problem. Fig. 3, panels (a)-(i) depict the interface
at various time instants as obtained from DNS. These are obtained by solving
the inviscid, axisymmetric, and incompressible Euler’s equations with surface-
tension and gravity in cylindrical coordinates (Basilisk, Popinet &
collaborators, 2013–2024) (script file is available as supplementary material
(Kayal, 2024)). The images in fig. 3 are obtained by generating the surface of
revolution of axisymmetric DNS data. As shown in panel (a), the interface is
initially distorted in the shape of an axisymmetric, stationary, and localised
perturbation. As this cavity relaxes, waves are generated which travel outward
reflecting off the wall (between panels (e) and (f)). This produces a wave-
train which focusses at the symmetry axis of the container ($r=0$). One notes
the formation of a small dimple-like structure at the symmetry axis in panel
(h). In § 3, we will demonstrate that neither the dimple nor other interface
features around the symmetry axis can be explained by the linear theory.

###  2.1 Governing equations: potential flow

We now turn to the theoretical analysis of the phenomenology illustrated in
fig. 3. In the base state, we consider a quiescent pool of liquid with density
$\rho$ and surface tension $T$ contained in a cylinder of radius $\hat{R}$.
For analytical simplicity, we assume this pool is infinitely deep compared to
the wavelength of the excited interface waves. For further simplicity, we
assume that the solid-liquid contact angle at the cylinder wall is always
fixed at $\pi/2$ and the contact line is free to move ($\partial_{n}v_{t}=0$).
This is the simplest contact line condition which allows for reflection of
waves at the boundary without complicating the analytical treatment of the
problem. The variables $\hat{\eta}\left(\hat{r},\hat{t}\right)$ are used to
represent the axisymmetric perturbed interface (see fig. 1) and
$\hat{\phi}\left(\hat{r},\hat{z},\hat{t}\right)$ is the disturbance velocity
potential; $\hat{r}$ and $\hat{z}$ being the radial and axial coordinates in
cylindrical geometry respectively. Variables with the dimensions of length
(e.g. $\hat{r},\hat{z},\hat{\eta}$) and time ($\hat{t}$) are scaled using
length and time-scales $L\equiv\hat{R}$ and
$T_{0}\equiv\sqrt{\frac{\hat{R}}{g}}$, respectively. The velocity potential
$\hat{\phi}$ is non-dimensionalised using the scale $L^{2}/T_{0}$. Under the
potential flow approximation, the nondimensional governing equations and
boundary conditions governing perturbed quantities are,

$\displaystyle\frac{\partial^{2}\phi}{\partial
r^{2}}+\frac{1}{r}\frac{\partial\phi}{\partial
r}+\frac{\partial^{2}\phi}{\partial z^{2}}=0,$ (2a)
$\displaystyle\frac{\partial\eta}{\partial
t}+\left(\frac{\partial\eta}{\partial
r}\right)\left(\frac{\partial\phi}{\partial
r}\right)_{z=\eta}-\left(\frac{\partial\phi}{\partial z}\right)_{z=\eta}=0,$
(2b) $\displaystyle\left(\frac{\partial\phi}{\partial
t}\right)_{z=\eta}+\eta+\frac{1}{2}\left\\{\left(\frac{\partial\phi}{\partial
r}\right)^{2}+\left(\frac{\partial\phi}{\partial
z}\right)^{2}\right\\}_{z=\eta}-\alpha\left\\{\frac{\frac{\partial^{2}\eta}{\partial
r^{2}}}{\left\\{1+\left(\frac{\partial\eta}{\partial
r}\right)^{2}\right\\}^{\frac{3}{2}}}+\frac{1}{r}\frac{\frac{\partial\eta}{\partial
r}}{\left\\{1+\left(\frac{\partial\eta}{\partial
r}\right)^{2}\right\\}^{\frac{1}{2}}}\right\\}=0,$ (2c)
$\displaystyle\int_{0}^{1}r\eta(r,t)dr=0,\quad\left(\frac{\partial\phi}{\partial
r}\right)_{r=1}=0,$ (2d,e)
$\displaystyle\lim_{z\to-\infty}\phi\to\text{finite}$ (2f)
$\displaystyle\eta(r,t=0)=-\varepsilon\left(1-\frac{r^{2}}{b^{2}}\right)\exp\left(-\frac{r^{2}}{b^{2}}\right)=\sum_{m=1}^{N}\eta_{m}(0){\mathrm{J}}_{0}(k_{m}r),\;\;\frac{\partial\phi}{\partial
n}(r,z=\eta(r,0),t=0)=0,$ (2g,h)

where $\varepsilon>0$ and $n$ in eqn. 2h is a distance coordinate measured
normal to the free-surface at $t=0$. The dimensionless parameters are defined
as follows: $\alpha\equiv\dfrac{T}{\rho g\hat{R}^{2}}$, representing the
inverse Bond number (based on the cylinder radius);
$b\equiv\dfrac{\hat{b}}{\hat{R}}$ is the dimensionless measure of cavity
width; and $\varepsilon\equiv\dfrac{\hat{a}_{0}}{\hat{R}}$ is the
dimensionless measure of cavity depth (see Fig. 2 caption for the meaning of
the symbols).

(a) t=0

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 3: Wave focussing observed in DNS from the cavity-shaped interface
distortion at $t=0$ (panel (a)). The figure is to be read left to right and
top to bottom for progression of time. After the waves reflect off the
cylinder wall (between panels (e) and (f); the confining walls are not shown),
they focus inwards towards $r=0$ producing strongly nonlinear oscillations of
increasing amplitude. The arrows indicate the instantaneous direction of wave
motion. The DNS parameters may be read from Case $1$ in table 1.

In cylindrical, axisymmetric coordinates. eqn. 2a is the Laplace equation, 2b
and 2c are the kinematic boundary condition and the Bernoulli equation applied
at the free surface respectively. Eqn. 2d restricts initial interfacial
distortions to those which are volume conserving while 2e enforces no-
penetration at the cylinder wall. Eqn. 2f is the finiteness condition at
infinite depth.

Eqns. 2 g & h represent the initial conditions. We decompose the initial
interface distortion i.e.
$\eta(r,t=0)=-\varepsilon\left(1-\frac{r^{2}}{b^{2}}\right)\exp\left(-\frac{r^{2}}{b^{2}}\right)$
(Miles, 1968), into a sum of Bessel modes as indicated by the second equality
sign in eqn. 2g and ${\mathrm{J}}_{1}(k_{m})=0$ for $m\in\mathbb{Z}^{+}$. The
numerical values of the modal coefficients at $t=0$ i.e.
$\eta_{m}(0)\,(m=1,2,3\ldots)$ in eqn. 2g are determined from the
orthogonality relation between Bessel functions i.e.
$\eta_{m}(0)=\frac{\int_{0}^{1}dr\;r\textrm{J}_{0}(k_{m}r)\eta(r,0)}{\int_{0}^{1}dr\;r\textrm{J}_{0}^{2}(k_{m}r)}$.
A sample representation of the initial condition and its modal coefficients is
presented in fig. 4a and 4b respectively where it is seen that about $17$
modes are excited initially. Subject to these initial and boundary conditions
presented in eqns. 2 a-h, we need to determine the modal amplitudes
$\eta_{m}(t),\,m=1,2,3\ldots$ as a function of time and this is carried out
next.

(a) Cavity shape at $t=0$

(b) Modal coefficients

Figure 4: Panel (a) The gas-liquid interface initially deformed as a cavity of
half-width $b=\dfrac{\hat{b}}{\hat{R}}$ and depth
$\varepsilon\equiv\dfrac{\hat{a_{0}}}{\hat{R}}$. Panel (b) The modal
cofficients $\eta_{m}(0)$ obtained by decomposing the initial distorted
interface. For this initial distortion, $\varepsilon=0.091,b=0.187$. It is
seen that only the first ten modes or so are excited initially. For accuracy,
we consider the energy in the first seventeen modes initially ($m=1,2,3\ldots
17$).

###  2.2 Equations for modal amplitude $\eta_{j}(t)$

In this section we solve the initial, boundary-value problem posed in eqns. 3
a-h. We derive equations governing the time evolution of modal coefficients
$\eta_{j}(t)$ upto quadratic order (i.e. terms which are cubic or higher in
the modal coefficients are neglected). The approach for doing this is
classical and was laid out in Hasselmann (1962) in Cartesian coordinates
although their initial conditions were random functions in contrast to the
deterministic initial distortion posed in eqn. 2g. The procedure below closely
follows the approach of Nayfeh (1987), who derived similar equations (his
eqns. $14$ and $15$) in the context of the Faraday instability (i.e. with
vertical oscillatory forcing) including gravity but not surface tension
(Nayfeh, 1987) in his analysis. In contrast to forced waves being studied by
Nayfeh (1987), we consider free waves in our current study and include both
surface-tension and gravity in the analysis. We first expand $\phi$ and $\eta$
in eqns. 2 as

$\displaystyle\phi(r,z,t)=\sum_{m=1}^{\infty}\phi_{m}(t)\textrm{J}_{0}(k_{m}r)\exp(k_{m}z),\quad\quad\eta(r,t)=\sum_{m=1}^{\infty}\eta_{m}(t)\textrm{J}_{0}(k_{m}r)$
(3a,b)

By construction, each term in the expansion in 3 satisfies the Laplace
equation 1a, eqns. 2 (d) and (e) as well as the finiteness condition 2f.
Taylor expanding eqns. 2b and c about $z=0$ we obtain

$\displaystyle\frac{\partial\eta}{\partial
t}-\left(\frac{\partial\phi}{\partial
z}\right)_{z=0}-\left(\frac{\partial^{2}\phi}{\partial
z^{2}}\right)_{z=0}\eta+\frac{\partial\eta}{\partial
r}\left(\frac{\partial\phi}{\partial r}\right)_{z=0}+\text{H.O.T}=0$ (4a)
$\displaystyle\left(\frac{\partial\phi}{\partial
t}\right)_{z=0}+\eta\left(\frac{\partial^{2}\phi}{\partial t\partial
z}\right)_{z=0}+\eta+\frac{1}{2}\left\\{\left(\frac{\partial\phi}{\partial
r}\right)^{2}+\left(\frac{\partial\phi}{\partial z}\right)^{2}\right\\}_{z=0}$
$\displaystyle-\alpha\left\\{\frac{\frac{\partial^{2}\eta}{\partial
r^{2}}}{\left\\{1+\left(\frac{\partial\eta}{\partial
r}\right)^{2}\right\\}^{\frac{3}{2}}}+\frac{1}{r}\frac{\frac{\partial\eta}{\partial
r}}{\left\\{1+\left(\frac{\partial\eta}{\partial
r}\right)^{2}\right\\}^{\frac{1}{2}}}\right\\}+\text{H.O.T}=0$ (4b)

where H.O.T represents higher order terms. Substituting expansions 3a & b into
4a,b and using orthogonality relations between Bessel functions we obtain for
$n,p,m\in\mathbb{Z}^{+}$

$\displaystyle\frac{d\eta_{n}}{dt}-k_{n}\phi_{n}(t)+\sum_{m,p}\left(D_{npm}-k_{m}^{2}C_{npm}\right)\phi_{m}(t)\eta_{p}(t)=0$
(5a) $\displaystyle\frac{d\phi_{n}}{dt}+(1+\alpha
k_{n}^{2})\eta_{n}(t)+\sum_{m,p}k_{m}C_{npm}\left(\frac{d\phi_{m}}{dt}\right)\eta_{p}(t)+\frac{1}{2}\sum_{m,p}\left(D_{npm}+k_{m}k_{p}C_{npm}\right)\phi_{m}(t)\phi_{p}(t)=0$
(5b) $\displaystyle n=1,2,3\ldots$

The nonlinear interaction coefficients $C_{npm}$ and $D_{npm}$ in eqn. 5 are
related as (Nayfeh, 1987):

$\displaystyle
D_{npm}=\frac{1}{2}\left(k_{p}^{2}+k_{m}^{2}-k_{n}^{2}\right)C_{npm}$ (6)

and
$C_{npm}=\dfrac{\int_{0}^{1}r\textrm{J}_{0}(k_{n}r)\textrm{J}_{0}(k_{p}r)\textrm{J}_{0}(k_{m}r)dr}{\int_{0}^{1}r\textrm{J}^{2}_{0}(k_{n}r)dr}$.
For the benefit of the reader, the detailed proof of 6 is provided in Appendix
A. Retaining self-consistently up to quadratic order terms, eqns. 5 a and b
may be combined into a second order equation for $\eta_{n}$ alone. This is:

$\displaystyle\dfrac{d^{2}\eta_{n}}{dt^{2}}+\omega_{n}^{2}\eta_{n}+k_{n}\sum_{m,p}\left[2+\frac{k_{p}^{2}-k_{m}^{2}-k_{n}^{2}}{2k_{m}k_{n}}\right]C_{npm}\left(\dfrac{d^{2}\eta_{m}}{dt^{2}}\right)\eta_{p}$
(7) $\displaystyle+$
$\displaystyle\dfrac{1}{2}k_{n}\sum_{m,p}\left[1+\dfrac{k_{p}^{2}+k_{m}^{2}-k_{n}^{2}}{2k_{m}k_{p}}+\dfrac{k_{p}^{2}-k_{m}^{2}-k_{n}^{2}}{k_{m}k_{n}}\right]C_{npm}\left(\dfrac{d\eta_{m}}{dt}\right)\left(\frac{d\eta_{p}}{dt}\right)=0$

Note that $\omega_{n}$ is the linear oscillation frequency of the
$n^{\text{th}}$ mode, viz. $\omega_{n}\equiv\sqrt{k_{n}\left(1+\alpha
k_{n}^{2}\right)}$. We solve the coupled ordinary differential eqns. 7
numerically subject to the initial conditions discussed earlier for
$n=1,2,3..\ldots 34$ (i.e. twice the initial number of excited modes, see fig.
4b) using ‘DifferentialEquations.jl’, an open-source package by Rackauckas et
al. (2017) and collaborators. The ‘DifferentialEquations.jl‘ automatically
chooses an ODE solver based on stiffness detection algorithms as described by
Rackauckas & Nie (2019). The Julia script file can be found in Kayal (2024).
We note that while numerically solving eqn. 7, we compute
$\frac{d^{2}\eta_{m}}{dt^{2}}$ in the third term of the equation (the
nonlinear term) via the linear estimate, viz,
$\frac{d^{2}\eta_{m}}{dt^{2}}=-\omega_{m}^{2}\eta_{m}$. Interestingly, the
solution to eqn. 7 shows instability at large times for sufficiently large
modal amplitudes, in qualitative agreement with the stability analysis of Zhu
et al. (2003) of axisymmetric standing waves (see their fig. $12$). In further
analysis, we restrict ourselves to numerical solutions to eqn. 7 which do not
contain this instability within the time period of our interest in this study
i.e. the time taken of focussing of the wave-train.

As benchmarking of our numerical solution procedure, we first solve eqns. 7
employing the single mode initial surface distortion that was studied in Basak
et al. (2021) i.e. in our current notation
$\eta(r,t=0)=\varepsilon{\mathrm{J}}_{0}(l_{5}\;r),\;\varepsilon>0$ where
$l_{5}=16.4706$ is the fifth non-trivial root of the Bessel function
${\mathrm{J}}_{1}$. For this initial condition, the second-order accurate
solution is expectedly of the form

$\displaystyle\eta(r,t)=\varepsilon\eta_{1}(r,t)+\varepsilon^{2}\eta_{2}(r,t)$
(8)

where explicit expressions for $\eta_{1}$ and $\eta_{2}$ were provided in
Basak et al. (2021) (we note the slight difference in non-dimensionalisation
of length between the current study and the one by Basak et al. (2021)
involving a factor of $l_{q}$). Fig. 5, demonstrates a comparison between the
prediction of eqn. 8 (indicated in the figure as ‘B21’ for Basak et al.
(2021)), the solution obtained from solving eqn. 7 with the same initial
condition (labelled in the figure as ‘Analytical’) and the numerical
simulation from Basilisk (depicted as ‘Simulation’). Fig. 5 demonstrates good
agreement between the three, thereby providing confidence on our numerical
procedure for solving eqns. 7.

Figure 5: Benchmarking of our solution procedure for solving the coupled
O.D.E’s in eqns. 7 against inviscid DNS (indicated as ‘Simulation’ in the
legend of panel (a)) and analytical predictions by Basak et al. (2021),
indicated as ‘B21’. For DNS, the dimensionless parameters are
$\varepsilon\equiv\dfrac{a_{0}}{\hat{R}}=\frac{0.5}{16.4706}=0.03,\alpha=0.004$
and $Oh=0$. Note that the initial condition here has a crest around $r=0$, see
inset of panel (a).

##  3 Comparison of direct numerical simulations (DNS) with theory

In this section, we compare results from our direct numerical simulations with
the theory discussed in § 2. Before this, it is instructive to rationalize the
reflection process and estimate its duration. To do this, we observe that the
Fourier-Bessel spectrum of the initial interface distortion prominently
features a mode at $m=4$ (see Fig. 4b). A rough estimate of the time required
for the energy associated with any mode excited in the initial spectrum to
complete a return trip (from $\hat{r}=0$ to the wall and back) can be derived
from linear theory. When this return time is estimated for the dominant mode
in the initial spectrum, we expect the numerical value to roughly coincide
with the generation time of the largest amplitude oscillation at $\hat{r}=0$
during the focusing process. This is illustrated in Fig. 6, where the time
signal from tracking the interface at $\hat{r}=0$ is presented (case 2 in
Table 1). Note that this figure uses dimensional variables, denoted with hats.
After the outward travelling waves move away, the interface at $\hat{r}=0$
remains relatively quiescent, as indicated by the nearly flat time signal
around $\hat{t}=0.2$ s. As a result of reflection, the energy associated with
every mode of wavenumber $k$ present initially focusses back to $\hat{r}=0$,
this return trip is carried out with its group-velocity
$\hat{c}_{g}=\frac{g+3\left(T/\rho\right)k^{2}}{2\sqrt{gk+Tk^{3}/\rho}}$. In
fig. 4b, the dominant mode is $k_{d}=\frac{l_{4}}{\hat{R}}$ and the largest
oscillation at $\hat{r}=0$ during the focussing process is seen to be
generated at $\hat{t}_{\text{peak}}=0.384$ s from fig. 6. Using the linear
estimate $\hat{t}_{\text{peak}}\approx\frac{2\hat{R}}{\hat{c}_{gd}}$ where
$\hat{c}_{gd}$ is the group-velocity of the dominant mode, we obtain the value
$0.403$ s which is reasonably close to the observed
$\hat{t}_{\text{peak}}=0.384$ s.

Figure 6: Time signal of the interface at $\hat{r}=0$. The green line indicates approximately the time window when focussing takes place at $\hat{r}=0$ $\mathrm{Case}$ | $\varepsilon\equiv\frac{\hat{a}_{0}}{\hat{R}}$ | $Oh\equiv\frac{\mu}{\sqrt{\rho T\hat{b}}}$ | $\hat{a}_{0}$ | $\nu$
---|---|---|---|---
1 | 0.061 | 0 | 0.26 | 0
2 | 0.091 | 0 | 0.39 | 0
3 | 0.091 | $1.17\times 10^{-5}$ | 0.39 | $8.9\times 10^{-5}$
4 | 0.091 | $1.17\times 10^{-4}$ | 0.39 | $8.9\times 10^{-4}$
5 | 0.091 | $1.17\times 10^{-3}$ | 0.39 | $8.9\times 10^{-3}$
6 | 0.091 | $1.17\times 10^{-2}$ | 0.39 | $8.9\times 10^{-2}$
7 | 0.006 | 0 | 0.026 | 0
8 | 0.006 | $1.17\times 10^{-5}$ | 0.026 | $8.9\times 10^{-5}$
9 | 0.006 | $1.17\times 10^{-4}$ | 0.026 | $8.9\times 10^{-4}$
10 | 0.006 | $1.17\times 10^{-3}$ | 0.026 | $8.9\times 10^{-3}$
11 | 0.006 | $3.7\times 10^{-3}$ | 0.026 | $2.81\times 10^{-2}$
12 | 0.006 | $1.17\times 10^{-2}$ | 0.026 | $8.9\times 10^{-2}$
Table 1: All dimensional lengths are indicated with a hat. Values are quoted
in CGS units. In all of the cases we have used $\hat{R}=4.282$ cm,
$\hat{b}=0.8$ cm, $T=72$ dyne$/$cm, $g=-981$ cm$/\text{s}^{2}$, $\rho=1$
gm$/\textrm{cm}^{3}$. These imply dimensionless values
$b\equiv\frac{\hat{b}}{\hat{R}}=0.187$, $\alpha\equiv\frac{T}{\rho
g\hat{R}^{2}}=0.004$

In the collage of images in figs. 7 and 8, we present the shape of the
interface as a function of time for case $1$ and $2$ in table 1 respectively,
comparing this to linear and nonlinear theoretical predictions. The only
difference between these two figures is in the value of $\varepsilon$, all
other dimensionless numbers remaining the same. Here linear theory implies
solution to eqn. 7 without the nonlinear terms. Note that this is equivalent
to superposition of the form
$\eta(r,t)=\sum_{m=1}^{17}{\mathrm{J}}_{0}(k_{m}r)\cos(\omega_{m}t)$ where
$\omega_{m}(k_{m})$ satisfies the gravito-capillary dispersion relation for
deeep water. In fig. 7, the transition from outward propagating waves to
inward propagating ones occur between panels (c) and (d). For panels (a), (b)
and (c) it is evident that linear theory represents the outgoing waves
accurately. However as focussing commences from panel (d) onwards, we notice
significant differences between linear theory and (inviscid) DNS.
Interestingly, second order theory seems to predict the shape of the interface
around $r=0$ quite well. Fig. 8 shows a more intense scenario than Fig. 7,
featuring a larger $\varepsilon=0.091$. The transition from predominantly
linear to nonlinear behavior occurs between panels (c) and (d), representing
outgoing and incoming waves, respectively. Notably, sharp dimple-like
structures emerge around $r=0$, as seen in panel (h), which are well described
by nonlinear theory. Additionally, the tendency to form jets, as seen in the
final panel, is noteworthy, although the nonlinear theory is only
qualitatively accurate in this context. We refer the reader to the
accompanying Movie $\\#1$ ($\epsilon=0.061$) and $\\#2$ ($\epsilon=0.091$),
see additional supplementary material which visualises these.

Figure 7: Waves generated from the cavity shaped interface distortion at $t=0$
(inset of panel (a)). We compare the interface shape as a function of time as
predicted by linear theory (L, solid blue line), second-order nonlinear theory
(N, solid green line) and (inviscid) DNS (Sim, red symbols). The waves
reflect-off the cylinder wall at $r=1$ (not shown) and focus back towards
$r=0$ generating oscillations of increasing amplitude. This corresponds to
case $1$ of table 1 with $\varepsilon=0.061$. To highlight the difference
between linear and nonlinear predictions, the figures have been plotted upto
$r=0.5$ instead of the entire radial domain up to $r=1$. The arrows depict the
instantaneous direction of motion of the waves. Figure 8: The same as fig. 7,
but for $\varepsilon=0.091$ corresponding to case $2$ in table 1. Note the
good qualitative agreement between nonlinear theory and (inviscid) DNS but not
linear theory, in capturing the dimple in panel (g). Also note the large
amplitude oscillations at $r=0$ with a tendendency to generate narrow jet-like
structure (panel (h)).

###  3.1 Role of nonlinearity at $r=0$

Figs. 7 and 8 show that although the linear solution is a reasonable model for
the interface evolution before reflection, it shows deviation from the fully
nonlinear simulation at the axis of symmetry during radial convergence of the
wave-train. In this subsection, we explain the apparent significance of
nonlinearity around the symmetry axis. To do this, we revisit results for the
single-mode interface distortion described by
$\eta(r,0)=\varepsilon\textrm{J}_{0}(k_{5}r)$ (where $\varepsilon>0$
corresponds to an initial crest at $r=0$ and $q=5$ is the primary mode), as
studied in Basak et al. (2021). For this initial condition, the expression for
$\eta(r,t)$ was analytically derived up to $\mathcal{O}(\varepsilon^{2})$ in
Basak et al. (2021) as:

$\eta(r,t)=\underbrace{\varepsilon\textrm{J}_{0}(k_{5}\;r)\cos(\omega_{5}t)}_{\text{Primary
mode}}\;+\;\varepsilon^{2}\sum_{j=1}^{\infty}\left[\underbrace{\zeta_{1}^{(j)}\cos(\omega_{j}t)}_{\text{Free
modes}}+\overbrace{\zeta_{2}^{(j)}\cos(2\omega_{5}t)+\zeta_{3}^{(j)}}^{\text{Bound
modes}}\right]\textrm{J}_{0}(k_{j}r),$ (9)

where $\zeta_{1}^{(j)}+\zeta_{2}^{(j)}+\zeta_{3}^{(j)}=0,\forall
j\in\mathbb{Z}^{+}$ to ensure that the initial condition is satisfied. Note
that that expression 9 has been suitably modified from Basak et al. (2021) to
make this compatible with the length and time scales in the present analysis.
Here $\varepsilon=\frac{\hat{a}_{0}}{\hat{R}}$, frequency
$\omega_{j}=\sqrt{k_{j}(1+\alpha k_{j}^{2})}$ and expressions for
$\zeta_{1}^{(j)},\,\zeta_{2}^{(j)}$ and $\zeta_{3}^{(j)}$ are provided in the
Appendix of (Basak et al., 2021). As highlighted in eqn. 9, the expression for
$\eta(r,t)$ comprises of three qualitatively different parts. The first term
on the right hand side of eqn. 9 represents the primary mode which is excited
at $t=0$. This mode has wavenumber $k_{5}$ and oscillates harmonically with
frequency $\omega_{5}$. Self interaction of the primary mode generates other
modes; the so-called “free modes” in eqn. 9 arise at
$\mathcal{O}(\varepsilon^{2})$ and their temporal frequency satisfy the
dispersion relation i.e. modes with wavenumber $k_{j}$ have frequency
$\omega_{j}$. The third kind viz. the “bound modes” which too appear at
$\mathcal{O}(\varepsilon^{2})$ vary temporally with the frequency
$2\omega_{5}\left(i.e.\;\omega_{5}+\omega_{5}\right)$ or
$0\left(i.e.\;\omega_{5}-\omega_{5}\right)$ independent of their wavenumber
$k_{j}$.

In Fig. 9, the interface from inviscid DNS with the initial condition
$\eta(r,0)=\varepsilon\textrm{J}_{0}(k_{5}r),\;\varepsilon=0.03>0$ is shown at
an instant when it forms a dimple-like protrusion at $r=0$. This is
represented by the curve with red dots, labelled as ‘Simulation’. In this
figure, we also plot the formula from Basak et al. (2021), excluding the bound
components (labelled as ‘Primary + Free’), i.e., setting
$\zeta_{2}^{(j)}=\zeta_{3}^{(j)}=0$ in equation (9). It is evident that this
approximation does not capture the dimple, which is otherwise predicted by the
full nonlinear expression (indicated as ‘Nonlinear’ in the figure caption and
referring to eqn. 9).

Figure 9: Various approximations for describing the dimple produced from a
single mode initial perturbation

The above exercise can also be carried out when the initial interface
deformation takes the shape of a cavity. For this initial condition,
$\eta(r,0)=\sum_{n=1}^{\infty}\eta_{n}(t)\textrm{J}_{0}(k_{n}r)$, as
previously shown in fig. 4. From the numerical solution to eqns. 7, the
temporal frequency spectrum at $r=0$ is obtained. We track the time series
generated by $\eta_{n}(t)$ and eliminate the frequencies $2\omega_{n}$ and $0$
from its Fourier spectrum, thus removing the bound wave components for all
initially excited modes. Fig. 10 demonstrates that after the removal of these
bound modes, the interface (labeled ‘Primary + Free’) fails to capture the
dimple shape. In contrast, the full numerical solution to eqns. 7 faithfully
reproduces the dimple.

Figure 10: Shape of a dimple for a cavity with $\epsilon=0.091$

###  3.2 Viscous effects: comparison with linear theory

In this section, we analyze viscous effects for the chosen initial condition.
Using cylindrical coordinates, Miles (1968) solved the problem of free-surface
waves on a viscous liquid in a linear regime within a radially unbounded
domain for multi-modal initial excitation. Farsoiya et al. (2017) extended
this theory to internal waves, considering viscosity and density due to both
upper and lower fluids, focusing on a single-mode initial perturbation. Due to
the availability of superposition in the linear regime, the results of
Farsoiya et al. (2017) are easily extended to multi-modal initial excitations
as well. In Cartesian geometry, the single mode initial excitation case was
first explicitly studied by Prosperetti (1976) treating free-surface waves and
by Prosperetti (1981) treating internal waves. In the Laplace domain and in
cylindrical axisymmetric coordinates, the solution to the evolution of a
single mode was shown in Farsoiya et al. (2017) to be given by,

$\tilde{\eta}_{m}(s)=\hat{\eta}_{m}(0)\dfrac{s+\left(4\tilde{k}_{m}^{2}\nu-\frac{4\tilde{k}_{m}^{3}\nu}{\tilde{k}_{m}+\sqrt{\tilde{k}_{m}^{2}+s/\nu}}\right)}{s^{2}+\left(4\tilde{k}_{m}^{2}\nu-\frac{4\tilde{k}_{m}^{3}\nu}{\tilde{k}_{m}+\sqrt{\tilde{k}_{m}^{2}+s/\nu}}\right)s+\hat{\omega}_{m}^{2}},\quad\hat{\omega}_{m}^{2}\equiv
g\tilde{k}_{m}+T\tilde{k}_{m}^{3}/\rho,\quad\tilde{k}_{m}\equiv\frac{k_{m}}{\hat{R}}$
(10)

Employing linear superposition, the corresponding (dimensional) expression for
the interface evolution in the time domain for the current multi-modal case is

$\hat{\eta}(\hat{r},\hat{t})=\sum_{m=1}^{17}\hat{\eta}_{m}(\hat{t}){\mathrm{J}}_{0}\left(k_{m}\frac{\hat{r}}{\hat{R}}\right),\quad\hat{\eta}_{m}(\hat{t})\equiv\mathcal{L}^{-1}\left[\tilde{\eta}_{m}(s)\right]$
(11)

Here $\mathcal{L}^{-1}$ is the inverse Laplace operator. We stress that
expression 10 accounts for dissipation in the bulk liquid and boundary layer,
as demonstrated by Prosperetti (1976) in Cartesian coordinates and by Farsoiya
et al. (2017) in cylindrical coordinates. Equation 11 is compared with DNS for
two different values of $\varepsilon$ and $Oh$ in figs. 11h and 12, where
inverse Laplace transforms were performed using the Cohen method by Henri
Cohen & Zagier (2000) which is a default method in mpmath (2023), a free
Python library for arbitrary-precision floating-point arithmetic. See Kayal
(2024) for the code. Fig. 11h benchmarks the theory at a relatively small
$\varepsilon=0.006$, where linear viscous theory is expected to be accurate.
Excellent agreement with linear viscous theory is observed in Fig. 11h.
Conversely, fig. 12 shows a clear distinction between linear and nonlinear
predictions.

To further investigate the impact of viscosity, fig. 13 presents the
interfacial velocity at $r=0$ from DNS for various $Oh$ values. The most
notable observation is that the highest peak velocity during wave focusing
occurs in the viscous simulation rather than the inviscid one. This non-
monotonous behavior as a function of the Ohnesorge number is a well-known
phenomenon in other contexts (Ghabache et al., 2014a), indicating a
significant effect of viscosity. In our analysis of converging waves, we
attribute the observed non-monotonic behavior to viscous dissipation within
the boundary layer at the gas-liquid interface. Even as the Ohnesorge number
approaches zero ($Oh=0^{+}$), this boundary layer remains significant, similar
to the dissipative anomaly seen in fully developed turbulence (Prandtl, 1904;
Onsager, 1949; Dubrulle, 2019; Eggers, 2018) and recently explored in contexts
such as sheet retraction (Sanjay et al., 2022) and drop impact (Sanjay et al.,
2023) interfacial flows. Consequently, this non-zero viscous dissipation
intensifies the focusing of capillary waves, thereby increasing the velocity
at the center ($r=0$). To validate this hypothesis, in the next section, we
next employ the viscous potential flow approach, which accounts for bulk
viscous dissipation but neglects dissipation in the gas-liquid boundary layer,
to model the converging waves.

###  3.3 Viscous potential flow

To further elucidate viscous effects, we incorporate viscosity into the
nonlinear equations using the viscous potential flow model (Joseph, 2006).
Unlike the linear case discussed previously, this method does not account for
the boundary layer formed at the free surface, since it does not enforce the
zero shear stress boundary condition (Moore, 1963). As is well-known, in this
approach the normal stress boundary condition (eqn. 2 c) is modified to
incorporate the effect of bulk viscous damping to obtain

$\displaystyle\left(\frac{\partial\phi}{\partial
t}\right)_{z=\eta}+\eta+2\;b\;Oh\;\sqrt{\alpha}\;\left(\frac{\partial^{2}\phi}{\partial
z^{2}}\right)_{z=\eta}+\frac{1}{2}\left\\{\left(\frac{\partial\phi}{\partial
r}\right)^{2}+\left(\frac{\partial\phi}{\partial
z}\right)^{2}\right\\}_{z=\eta}$
$\displaystyle-\alpha\left\\{\frac{\frac{\partial^{2}\eta}{\partial
r^{2}}}{\left\\{1+\left(\frac{\partial\eta}{\partial
r}\right)^{2}\right\\}^{\frac{3}{2}}}+\frac{1}{r}\frac{\frac{\partial\eta}{\partial
r}}{\left\\{1+\left(\frac{\partial\eta}{\partial
r}\right)^{2}\right\\}^{\frac{1}{2}}}\right\\}=0.$ (12)

We follow the same strategy as the inviscid case and obtain a modified
differential equation for $\eta_{n}$, i.e. the viscous counterpart of eqn. 7
leading to

$\displaystyle\dfrac{d^{2}\eta_{n}}{dt^{2}}+\omega_{n}^{2}\eta_{n}+2\;b\;Oh\;\sqrt{\alpha}\;k_{n}^{2}\frac{d\eta_{n}}{dt}+2b\;Oh\;\sqrt{\alpha}\;k_{n}\sum_{m,p}k_{m}^{2}C_{npm}\frac{d\eta_{m}}{dt}\eta_{p}$
(13)
$\displaystyle+k_{n}\sum_{m,p}\left[2+\frac{k_{p}^{2}-k_{m}^{2}-k_{n}^{2}}{2k_{m}k_{n}}\right]C_{npm}\left(\dfrac{d^{2}\eta_{m}}{dt^{2}}\right)\eta_{p}$
$\displaystyle+$
$\displaystyle\dfrac{1}{2}k_{n}\sum_{m,p}\left[1+\dfrac{k_{p}^{2}+k_{m}^{2}-k_{n}^{2}}{2k_{m}k_{p}}+\dfrac{k_{p}^{2}-k_{m}^{2}-k_{n}^{2}}{k_{m}k_{n}}\right]C_{npm}\left(\dfrac{d\eta_{m}}{dt}\right)\left(\frac{d\eta_{p}}{dt}\right)=0$

In fig. 14, we compare the nonlinear analytical inviscid solution (referred to
as ‘Inviscid’ in the legend), the viscous potential flow (VPF) solution for
$\varepsilon=0.091$ and the viscous DNS (referred to as ‘Simulation’) for Case
$4$ in table 1. It is seen that the VPF solution, is indistinguishable from
the inviscid one in the limit of $Oh=0^{+}$, highlighting the importance of
resolving the viscous boundary layer in theory.

To further quantify the comparison between these cases, in fig. 15a shows the
velocity at the axis of symmetry within a shallow cavity during focusing. The
linear viscous theory, which accounts for the boundary layer at the free
surface, describes the change in $v_{z}$ with Ohnesorge number slightly better
than the VPF model. Fig. 15b presents results for a deeper cavity where non-
linearity plays a significant role, and the non-monotonic behavior observed in
fig. 13 as a function of $Oh$ is evident. The VPF model fails to capture this
non-monotonic behavior, highlighting the importance of resolving the boundary
layer at the gas-liquid interface, as discussed in § 3.2. We propose
developing a nonlinear-viscous theory superior to the VPF model to explain the
observations in Fig. 13 and Fig. 15b.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 11: Viscous DNS (indicated as ‘Simulation’ with red dots in the legend
to panel (a)

) with $\varepsilon=0.006$ and $Oh=1.17\times 10^{-3}$ corresponding to case
$10$ in table 1. One notes the excellent agreement with linear, viscous theory
(green dots, ’Linear’, eqn. 11 in text) with hardly any nonlinear
contribution.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 12: Viscous DNS (indicated as ‘Simulation’ with red dots in the legend
to panel (a)) with $\varepsilon=0.091$ and $Oh=1.17\times 10^{-4}$
corresponding to case $4$ in table 1. In contrast to fig. 11h, increasing the
value of $\varepsilon$ and a corresponding reduction in viscosity, has a
dramatic effect in the simulations. We note that viscous linear theory is no
longer adequate particularly during the focussing process in panels (f)-(h).
In panel (h), we also provide a comparison of the interface at this time-
instant, for the inviscid numerical simulation ($Oh=0$) with the same
$\varepsilon$. It is seen that the viscous simulation has a crest which at
this instant is taller than the one obtained from the inviscid simulation.

Figure 13: Velocity at the interface at $\hat{r}=0$ for different values of
$Oh$ and fixed $\varepsilon=0.091$. Note that the viscous DNS for
$Oh=1.17\times 10^{-4}$ (solid blue line) produces the largest velocity peak
around $t\approx 5.7$. Note in particular that the inviscid signal ($Oh=0$,
red solid line) has a peak which is shorter by a factor of half compared to
the signal corresponding to $Oh=1.17\times 10^{-4}$. We refer the reader to
Appendix B where the grid convergence for this (and other) simulations are
provided

.

(a)

(b)

(c)

(d)

Figure 14: Comparison of the viscous potential solution (VPF, black dotted
line), inviscid solution (green solid line) and DNS (red dots) at
$\varepsilon=0.091$ and $Oh=1.17\times 10^{-4}$, case $4$ in table 1.

(a) Shallow cavity with $\varepsilon=0.006$

(b) Deep cavity with $\varepsilon=0.091$

Figure 15: Comparison of the maximum velocity at $r=0$ after reflection for
different Ohnesorge number for a shallow cavity, cases $7-12$ in table 1
(panel a) and for a deep cavity, case $2-5$ (panel b). Panel (a): ‘+’ symbols
represent DNS with finite viscosity. Black dotted symbols represent DNS with
zero viscosity. Red symbols represent the linear viscous solution obtained by
numerical inversion of eqn. 10. Green symbols indicate viscous potential flow
(VPF) approximation obtained from solving eqn. 13. At the $Oh=0$ limit, VPF
(green dashed line) and linear viscous theory (red dashed line) coincide with
the linear inviscid theory (blue dashed line). Panel (b): Symbols have the
same meaning as panel (a), the only difference is that we have employed
nonlinear inviscid theory in this case. Note that non-monotonicity in the
velocity at $r=0$ as a function of $Oh$. The viscous potential flow
approximation (VPF) despite being nonlinear is unable to describe this non-
monotonicity, presumably because of its inability to resolve the boundary
layer at the free-surface. In both of the panels, the dotted black line
represents the velocity of inviscid DNS.

##  4 Conclusion & outlook

In this study, we have discussed the dynamics of a localised free-surface
perturbation in a cylindrical pool of liquid, which generates a train of
waves. These waves, upon reflecting from the container walls, converge back
towards the axis of symmetry, leading to progressively increasing free-surface
oscillations at the center. Using the potential flow approximation, we derived
a set of ordinary differential equations governing the evolution of modal
amplitudes up to second order in the amplitude.

For shallow cavities, linear theory suffices to explain the wave evolution.
However, as the cavity depth increases, the limitations of linear theory
become evident, particularly in predicting the focusing effects at $r=0$. Our
findings demonstrate that linear dispersive focusing alone is inadequate to
describe the intricate dimple shape forming at the axis of symmetry for deeper
cavities. A nonlinear theory that accounts for energy transfer among various
modes is essential for accurately modeling the cavity shape evolution. The
role of bound modes is critical in capturing the emergence of the dimple.

A notable observation is the significant influence of viscosity on the
focusing process. Contrary to expectations, the maximum velocity at the axis
of symmetry is higher for a slightly viscous fluid than for an inviscid one.
This non-monotonic behavior with respect to the Ohnesorge number $Oh$ is not
captured by either the linear viscous model (Farsoiya et al., 2017;
Prosperetti, 1976) or the nonlinear viscous potential flow (VPF) model
(Joseph, 2006). The failure of the VPF model, which neglects the boundary
layer effects, underscores the necessity of accounting for the boundary layer
in the limit of $Oh\to 0^{+}$. The VPF model’s convergence to an inviscid
solution in this limit further strengthens the assertion that boundary layers
are crucial for the observed velocity enhancements.

Future research should develop a comprehensive nonlinear viscous theory that
incorporates both boundary layer effects and the nonlinearity of wave
focusing. This theory would provide a more accurate analytical resolution of
discrepancies observed in current models. Additionally, extending this work to
non-Newtonian fluids, such as viscoplastic or viscoelastic liquids (Sanjay et
al., 2021), could reveal new insights and broaden the applicability of our
theoretical framework.

Acknowledgements. RD thanks Prof. B. Sutherland for sharing the study by Smith
(1976).

Funding. We gratefully acknowledge financial support from DST-SERB (Govt. of
India) grants MTR/2019/001240, CRG/2020/003707 and SPR/2021/000536 on topics
related to waves, jet formation, cavity collapse and the viscous Cauchy-
Poisson problem. The Ph.D. tenure of L.K. is supported by the Prime-Minister’s
Research Fellowship (PMRF), Govt. of India and is gratefully acknowledged.

Declaration of Interests. The authors report no conflict of interest.

## Appendix A

We derive the relation between the nonlinear interaction coefficients
$C_{mnq}$ and $D_{mnq}$ discussed in eqn. 6. This relation has been provided
in Nayfeh (1987) and Miles (1976) without proof and the same is presented
here. Following Nayfeh (1987), we represent eqn. 3 a in (semi) basis-
independent notation as

$\displaystyle\phi(\mathbf{x}{,}z,t)=\sum_{m=1}^{\infty}\phi_{m}(t)\Psi_{m}(\mathbf{x}{)}\exp(k_{m}z)$
(14)

where $\mathbf{x}$ is the horizontal position vector and $\Psi_{m}$ satisfies
the equation $\nabla_{H}^{2}\Psi_{m}+k_{m}^{2}\Psi=0$ as a consequence of
$\phi$ satisfying the Laplace eqn; note that
$\nabla^{2}=\nabla_{H}^{2}+\frac{\partial^{2}}{\partial z^{2}}$. We assume
that $\Psi_{m}(\mathbf{x})$ follow the orthogonality rule $\int\int
dS\;\psi_{m}(\mathbf{x})\psi_{q}(\mathbf{x})=\delta_{mq}S$ where $\delta_{mq}$
is the Kronecker delta. Using Stokes theorem to relate an area integral (over
s) in two dimensions to the line integral, we have for a vector field
$\mathbf{F}(\mathbf{x})$

$\displaystyle\int\int ds\;\bm{\nabla}_{H}\cdot\mathbf{F}=\int
dl\;\left(\mathbf{F}\cdot\mathbf{n}\right).$ (15)

Choosing $\mathbf{F}=\psi_{q}\psi_{m}\bm{\nabla}_{H}\psi_{n}$, eqn. 15 leads
to

$\displaystyle\int\int\left[\psi_{q}\left(\bm{\nabla}_{H}\psi_{m}\cdot\bm{\nabla}_{H}\psi_{n}\right)+\psi_{m}\left(\bm{\nabla}_{H}\psi_{q}\cdot\bm{\nabla}_{H}\psi_{n}\right)+\psi_{q}\psi_{m}\nabla_{H}^{2}\psi_{n}\right]\;ds=0,$
(16)

the right hand side following from the no-penetration condition at the wall.
Following the same notation as Nayfeh (1987), we define,

$\displaystyle\int\int
ds\;\psi_{m}(\mathbf{x}{)}\psi_{n}(\mathbf{x}{)}\psi_{q}(\mathbf{x}{)}\equiv
S\;C_{mnq},\quad\int\int
ds\;\left(\bm{\nabla}_{H}\psi_{m}(\mathbf{x}{)}\cdot\bm{\nabla}_{H}\psi_{n}(\mathbf{x}{)}\right)\psi_{q}(\mathbf{x}{)}\equiv
S\;D_{mnq}.$ (17)

Note that $D_{nmq}=D_{mnq}$. Using 17, eqn. 16 may be written compactly as

$\displaystyle D_{mnq}+D_{qnm}-k_{n}^{2}C_{mnq}=0.$ (18)

Replacing $m\rightarrow n,n\rightarrow q,q\rightarrow m$ in equation 18, we
obtain

$\displaystyle D_{nqm}+D_{mqn}-k_{q}^{2}C_{nqm}=0.$ (19)

which may be rewritten as

$\displaystyle D_{qnm}+D_{qmn}-k_{q}^{2}C_{nqm}=0$ (20)

Using 20 in 18, we obtain

$\displaystyle
D_{mnq}=k_{n}^{2}C_{mnq}-\left(k_{q}^{2}C_{nqm}-D_{qmn}\right).$ (21)

Replacing once again $m\rightarrow q,n\rightarrow m,q\rightarrow n$ in
equation 18, we obtain

$\displaystyle D_{qmn}+D_{nmq}-k_{m}^{2}C_{qmn}=0.$ (22)

Combining 21 and 22 and the fact that $D_{mnq}=D_{nmq}$, we obtain

$\displaystyle
D_{nmq}=\frac{1}{2}\left(k_{n}^{2}+k_{m}^{2}-k_{q}^{2}\right)C_{nmq}$ (23)

After some manipulation, expression 6 follows from the above expression.

## Appendix B

Figs. 16 and 17 present the grid convergence results at three different grid
resolutions ($512^{2},\,1024^{2}$ and $2048^{2}$) for case $4$ and case $2$
respectively in table 1. Figs. 18 presents the grid convergence results for
fig. 13.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 16: Comparison of interface profile for case $4$ in table 1 at three
different grid resolutions, $512^{2}$ (blue dots), $1024^{2}$ (red dots) and
$2048^{2}$ (green dots)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 17: Comparison of interface profile for case $2$ in table 1 for three
different grid resolutions, $512^{2}$ (blue dots), $1024^{2}$ (red dots) and
$2048^{2}$ (green dots)

(a) $Oh=0$

(b) $Oh=1.17\times 10^{-4}$

Figure 18: Comparison of the vertical velocity for case 2 (panel a) and case 4
(panel b) for three different grid resolutions, $512^{2}$ (red solid line),
$1024^{2}$ (blue solid line) and $2048^{2}$ (green solid line).

## References

  * Basak et al. (2021) Basak, S., Farsoiya, P. K. & Dasgupta, R. 2021 Jetting in finite-amplitude, free, capillary-gravity waves. J. Fluid Mech. 909.
  * Benjamin & Feir (1967) Benjamin, T. B. & Feir, J. E. 1967 The disintegration of wave trains on deep water part 1. theory. J. Fluid Mech. 27 (3), 417–430.
  * Britannica, The Editors of Encyclopaedia (2024) Britannica, The Editors of Encyclopaedia 2024 Aghulas current. https://www.britannica.com/place/Agulhas-Current.
  * Brown & Jensen (2001) Brown, Michael G & Jensen, Atle 2001 Experiments on focusing unidirectional water waves. Journal of Geophysical Research: Oceans 106 (C8), 16917–16928.
  * Charlie Wood (2020) Charlie Wood 2020 The grand unified theory of rogue waves. https://www.quantamagazine.org/the-grand-unified-theory-of-rogue-waves-20200205/.
  * Chavarria et al. (2018) Chavarria, Gerardo Ruiz, Le Gal, Patrice & Le Bars, Michael 2018 Geometrical focusing of surface waves. Physical Review Fluids 3 (9), 094803.
  * Davis & Zarnick (1964) Davis, MC & Zarnick, Ernest E 1964 Testing ship models in transient waves. In 5th Symposium on Naval Hydrodynamics, , vol. 507.
  * Deike (2022) Deike, L. 2022 Mass transfer at the ocean–atmosphere interface: the role of wave breaking, droplets, and bubbles. Annu. Rev. Fluid Mech. 54, 191–224.
  * Dhote et al. (2024) Dhote, Yashika, Kumar, Anil, Kayal, Lohit, Goswami, Partha Sarathi & Dasgupta, Ratul 2024 Standing waves and jets on a sessile, incompressible bubble. Physics of Fluids 36 (1).
  * Dubrulle (2019) Dubrulle, B. 2019 Beyond Kolmogorov cascades 867.
  * Duchemin et al. (2002) Duchemin, Laurent, Popinet, Stéphane, Josserand, Christophe & Zaleski, Stéphane 2002 Jet formation in bubbles bursting at a free surface. Physics of fluids 14 (9), 3000–3008.
  * Dysthe et al. (2008) Dysthe, Kristian, Krogstad, Harald E & Müller, Peter 2008 Oceanic rogue waves. Annu. Rev. Fluid Mech. 40, 287–310.
  * Eggers (2018) Eggers, J. 2018 Role of singularities in hydrodynamics 3 (11), 110503.
  * Eggers et al. (2024) Eggers, J. G., Sprittles, J. & Snoeijer, J. H. 2024 Coalescence dynamics. Annu. Rev. Fluid Mech. .
  * Farsoiya et al. (2017) Farsoiya, Palas Kumar, Mayya, YS & Dasgupta, Ratul 2017 Axisymmetric viscous interfacial oscillations–theory and simulations. Journal of Fluid Mechanics 826, 797–818.
  * Ghabache et al. (2014a) Ghabache, Elisabeth, Antkowiak, Arnaud, Josserand, Christophe & Séon, Thomas 2014a On the physics of fizziness: How bubble bursting controls droplets ejection. Physics of Fluids 26 (12).
  * Ghabache et al. (2014b) Ghabache, Élisabeth, Séon, Thomas & Antkowiak, Arnaud 2014b Liquid jet eruption from hollow relaxation. Journal of fluid mechanics 761, 206–219.
  * Gordillo & Blanco-Rodríguez (2023) Gordillo, J. M. & Blanco-Rodríguez, F. J. 2023 Theory of the jets ejected after the inertial collapse of cavities with applications to bubble bursting jets. Phys. Rev. Fluids 8 (7), 073606\.
  * Hasselmann (1962) Hasselmann, K. 1962 On the non-linear energy transfer in a gravity-wave spectrum part 1. general theory. Journal of Fluid Mechanics 12 (4), 481–500.
  * Henri Cohen & Zagier (2000) Henri Cohen, Fernando Rodriguez Villegas & Zagier, Don 2000 Convergence acceleration of alternating series. Experimental Mathematics 9 (1), 3–12.
  * Johannessen & Swan (2001) Johannessen, Thomas B & Swan, Chris 2001 A laboratory study of the focusing of transient and directionally spread surface water waves. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 457 (2008), 971–1006.
  * Joseph (2006) Joseph, Daniel D 2006 Potential flow of viscous fluids: Historical notes. International Journal of Multiphase Flow 32 (3), 285–310.
  * Kayal (2024) Kayal 2024 Codes for cavity evolution. https://github.com/Lohittitas/concentricWave.
  * Kayal et al. (2022) Kayal, Lohit, Basak, Saswata & Dasgupta, Ratul 2022 Dimples, jets and self-similarity in nonlinear capillary waves. Journal of Fluid Mechanics 951, A26.
  * Kayal & Dasgupta (2023) Kayal, Lohit & Dasgupta, Ratul 2023 Jet from a very large, axisymmetric, surface-gravity wave. Journal of Fluid Mechanics 975, A22.
  * Kientzler et al. (1954) Kientzler, CF, Arons, Arnold B, Blanchard, Duncan C & Woodcock, AHhttps 1954 Photographic investigation of the projection of droplets by bubbles bursting at a water surface. Tellus 6 (1), 1–7.
  * Longuet-Higgins (1974) Longuet-Higgins, MS 1974 Breaking waves in deep or shallow water. In Proc. 10th Conf. on Naval Hydrodynamics, , vol. 597, p. 605. MIT.
  * Longuet-Higgins (1983) Longuet-Higgins, Michael S 1983 Bubbles, breaking waves and hyperbolic jets at a free surface. Journal of Fluid Mechanics 127, 103–121.
  * Longuet-Higgins (1990) Longuet-Higgins, Michael S 1990 An analytic model of sound production by raindrops. Journal of Fluid Mechanics 214, 395–410.
  * MacIntyre (1968) MacIntyre, Ferren 1968 Bubbles. boundary-layer” microtome” for micronthick samples of a liquid surface. The Journal of Physical Chemistry 72 (2), 589–592.
  * MacIntyre (1972) MacIntyre, Ferren 1972 Flow patterns in breaking bubbles. Journal of Geophysical Research 77 (27), 5211–5228.
  * Mallory (1974) Mallory, John Kennaway 1974 Abnormal waves on the south east coast of south africa. The International hydrographic review .
  * McAllister et al. (2022) McAllister, ML, Draycott, S, Davey, T, Yang, Y, Adcock, TAA, Liao, S & van den Bremer, TS 2022 Wave breaking and jet formation on axisymmetric surface gravity waves. Journal of Fluid Mechanics 935.
  * McIver (1985) McIver, Maureen 1985 Diffraction of water waves by a moored, horizontal, flat plate. Journal of Engineering Mathematics 19 (4), 297–319.
  * Mei & Liu (1973) Mei, CC & Liu, LF 1973 The damping of surface gravity waves in a bounded liquid. Journal of Fluid Mechanics 59 (2), 239–256.
  * Miles (1968) Miles, John W 1968 The cauchy–poisson problem for a viscous liquid. Journal of Fluid Mechanics 34 (2), 359–370.
  * Miles (1976) Miles, John W. 1976 Nonlinear surface waves in closed basins. Journal of Fluid Mechanics 75 (3), 419–448.
  * Moore (1963) Moore, D. W. 1963 The boundary layer on a spherical gas bubble. J. Fluid Mech. 16 (2), 161–176.
  * mpmath (2023) mpmath 2023 mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 1.3.0). http://mpmath.org/.
  * Murashige & Kinoshita (1992) Murashige, Sunao & Kinoshita, Takeshi 1992 An ideal ocean wave focusing lens and its shape. Applied ocean research 14 (5), 275–290.
  * Nayfeh (1987) Nayfeh, Ali H 1987 Surface waves in closed basins under parametric and internal resonances. The Physics of fluids 30 (10), 2976–2983.
  * Oguz & Prosperetti (1990) Oguz, Hasan N & Prosperetti, Andrea 1990 Bubble entrainment by the impact of drops on liquid surfaces. Journal of Fluid Mechanics 219, 143–179.
  * Onorato et al. (2013) Onorato, Miguel, Residori, S, Bortolozzo, U, Montina, A & Arecchi, FT 2013 Rogue waves and their generating mechanisms in different physical contexts. Physics Reports 528 (2), 47–89.
  * Onsager (1949) Onsager, L. 1949 Statistical hydrodynamics. Il Nuovo Cimento 6 (2), 279–287.
  * Peregrine (1976) Peregrine, D Howell 1976 Interaction of water waves and currents. Advances in applied mechanics 16, 9–117.
  * Peregrine (1986) Peregrine, D. H. 1986 Approximate descriptions of the focussing of water waves. In Coastal Engineering 1986, pp. 675–685. ASCE.
  * Popinet & collaborators (2013–2024) Popinet, S. & collaborators 2013–2024 Basilisk C: volume of fluid method. http://basilisk.fr (Last accessed: August 23, 2023).
  * Prandtl (1904) Prandtl, L. 1904 Über flüssigkeitsbewegung bei sehr kleiner reibung. Math-Kongr, Heidelberg pp. 484–491.
  * Prosperetti (1976) Prosperetti, A. 1976 Viscous effects on small-amplitude surface waves. Phys. Fluids 19 (2), 195–203.
  * Prosperetti (1981) Prosperetti, Andrea 1981 On the motion of two superposed viscous fluids .
  * Rackauckas & Nie (2019) Rackauckas, Christopher & Nie, Qing 2019 Confederated modular differential equation apis for accelerated algorithm development and benchmarking. Advances in Engineering Software 132, 1–6.
  * Rackauckas et al. (2017) Rackauckas, Christopher, Nie, Qing & collaborators 2017 Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia. Journal of Open Research Software 5 (1), 15.
  * Rapp & Melville (1990) Rapp, Ronald James & Melville, W Kendall 1990 Laboratory measurements of deep-water breaking waves. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 331 (1622), 735–800.
  * Sanjay (2022) Sanjay, V. 2022 Viscous Free-Surface Flows. PhD thesis, University of Twente.
  * Sanjay et al. (2023) Sanjay, Vatsal, Chantelot, Pierre & Lohse, Detlef 2023 When does an impacting drop stop bouncing? Journal of Fluid Mechanics 958, A26.
  * Sanjay et al. (2021) Sanjay, V., Lohse, D. & Jalaal, M. 2021 Bursting bubble in a viscoplastic medium. J. Fluid Mech. 922, A2.
  * Sanjay et al. (2022) Sanjay, V., Sen, U., Kant, P. & Lohse, D. 2022 Taylor-Culick retractions and the influence of the surroundings 948, A14.
  * Smith (1976) Smith, Ronald 1976 Giant waves. Journal of Fluid Mechanics 77 (3), 417–431.
  * Stamnes et al. (1983) Stamnes, Jakob J, Løvhaugen, Odd, Spjelkavik, Bjøsrn, Mei, Chiang C, Lo, Edmond & Yue, Dick KP 1983 Nonlinear focusing of surface waves by a lens–theory and experiment. Journal of Fluid Mechanics 135, 71–94.
  * Stuhlman Jr (1932) Stuhlman Jr, Otto 1932 The mechanics of effervescence. Physics 2 (6), 457–466.
  * Torres et al. (2022) Torres, Theo, Lloyd, Max, Dolan, Sam R & Weinfurtner, Silke 2022 Wave focusing by submerged islands and gravitational analogues. Physical Review Research 4 (3), 033210.
  * Villermaux et al. (2022) Villermaux, E., Wang, X. & Deike, L. 2022 Bubbles spray aerosols: Certitudes and mysteries. PNAS Nexus 1 (5), pgac261.
  * White & Fornberg (1998) White, Benjamin S & Fornberg, Bengt 1998 On the chance of freak waves at sea. Journal of fluid mechanics 355, 113–138.
  * Wu & Nepf (2002) Wu, Chin H & Nepf, HM 2002 Breaking criteria and energy losses for three-dimensional wave breaking. Journal of Geophysical Research: Oceans 107 (C10), 41–1.
  * Xu & Perlin (2023) Xu, Chang & Perlin, Marc 2023 Parasitic waves and micro-breaking on highly nonlinear gravity–capillary waves in a convergent channel. Journal of Fluid Mechanics 962, A46.
  * Zakharov et al. (2006) Zakharov, Vladimir E, Dyachenko, Alexander I & Prokofiev, Alexander O 2006 Freak waves as nonlinear stage of stokes wave modulation instability. European Journal of Mechanics-B/Fluids 25 (5), 677–692.
  * Zhu et al. (2003) Zhu, Qiang, Liu, Yuming & Yue, Dick KP 2003 Three-dimensional instability of standing waves. Journal of Fluid Mechanics 496, 213–242.