# A boundary perturbation method to simulate nonlinear deformations of a two-
dimensional bubble

Philippe Guyenne  Department of Mathematical Sciences, University of Delaware,
Newark, DE 19716, USA<EMAIL_ADDRESS>
###### Abstract

Nonlinear deformations of a two-dimensional gas bubble are investigated in the
framework of a Hamiltonian formulation involving surface variables alone. The
Dirichlet–Neumann operator is introduced to accomplish this dimensional
reduction and is expressed via a Taylor series expansion. A recursion formula
is derived to determine explicitly each term in this Taylor series up to an
arbitrary order of nonlinearity. Both analytical and numerical strategies are
proposed to deal with this nonlinear free-boundary problem under forced or
freely oscillating conditions. Simplified models are established in various
approximate regimes, including a Rayleigh–Plesset equation for the time
evolution of a purely circular pulsating bubble, and a second-order Stokes
wave solution for weakly nonlinear shape oscillations that rotate steadily on
the bubble surface. In addition, a numerical scheme is developed to simulate
the full governing equations, by exploiting the efficient and accurate
treatment of the Dirichlet–Neumann operator via the fast Fourier transform.
Extensive tests are conducted to assess the numerical convergence of this
operator as a function of various parameters. The performance of this direct
solver is illustrated by applying it to the simulation of cycles of
compression-dilatation for a purely circular bubble under uniform forcing, and
to the computation of freely evolving shape distortions represented by
steadily rotating waves and time-periodic standing waves. The former solutions
are validated against predictions by the Rayleigh–Plesset model, while the
latter solutions are compared to laboratory measurements in the case of mode-2
standing waves.

###### keywords:

Dirichlet–Neumann operator , gas bubbles , Hamiltonian systems , pseudo-
spectral method , Rayleigh–Plesset equation

††journal: …

##  1 Introduction

Gas bubble dynamics has been the subject of intensive research in recent
decades due to its prominent role in a number of physical phenomena. For
example, bubble formation during gas entrainment into water by breaking waves
produces sound over large distances which contributes to ambient noise in the
ocean. Oscillations of bubble volume may lead to collapse and generate shock
waves with damaging effects on the submerged machinery in cavitation problems.
Large distortions of bubble shape may lead to break-up and contribute to
increasing bubble population in the surrounding fluid. Related phenomena such
as sonoluminescence [3], and applications involving the ultrasonic excitation
of microbubbles for medical imaging [39] or the implosion of deuterium-tritium
capsules for inertial confinement fusion [36], have drawn a lot of attention
in recent years. We are aware of a similar large literature on the dynamics of
fluid droplets, which is closely related to the present subject, but for
convenience we will restrict the following introduction to gas bubbles.

In theoretical studies of cavitation for a single bubble, a common approach
relies on the Rayleigh–Plesset (RP) model which is a second-order nonlinear
ordinary differential equation (ODE) describing the time evolution of the
bubble radius under the influence of surface tension, viscosity and an
external pressure field. In general, the latter is meant to represent a source
of disturbances from the surrounding environment (e.g. internal variations of
the fluid flow, perturbations due to atmospheric pressure or waves at the sea
surface, noise generated by submerged machinery or by other bubbles). The
resulting oscillations of the bubble radius are usually referred to as radial
or volume oscillations. The RP model can be derived from the incompressible
Navier–Stokes equations by assuming spherical symmetry (i.e. angular
invariance) with appropriate boundary conditions. Its simple closed form as
compared to the full Navier–Stokes equations makes it well suited for both
mathematical analysis and numerical simulation. In particular, bifurcation
theory for nonlinear ODEs has been successfully applied to characterizing the
rich behavior of solutions under external forcing. Reviews on this model with
a detailed discussion on its properties, predictions and extensions
(incorporating e.g. compressible or thermal effects) can be found in [16; 47;
49].

For radial pulsations of sufficiently large amplitude, distortions of the
bubble shape (with angular variation) tend to develop, induced by parametric
instability as reported in experimental observations. Investigation of these
shape oscillations has also produced an abundant literature, with earlier
mathematical results focusing on the linear case for small distortions [47] or
the weakly nonlinear case without radial motion [52]. The problem of bubble
dynamics with a nontrivial geometry is significantly more difficult to analyze
than the RP framework, and the vast majority of theoretical studies have
considered three-dimensional axisymmetric bubbles in potential flow, meaning
that their shape is assumed to be invariant with respect to the azimuthal
angle in spherical coordinates. Seeking a weakly nonlinear solution in the
form of an asymptotic series for small- to moderate-amplitude distortions and
employing the method of multiple scales, a hierarchy of inhomogeneous linear
equations can be established to determine each order of approximation as a
function of previous orders. Then, by expanding each order in terms of
Legendre polynomials relative to the inclination angle (owing to the
axisymmetry) and by exploiting their orthogonality property, explicit
expressions can be obtained for the shape modes (i.e. the coefficients in
factor of the Legendre polynomials) at any order in the time-periodic setting.
These perturbation calculations are particularly tedious because the
inhomogeneous equations become increasingly complicated with the level of
approximation, and for this reason, they have usually been restricted to a few
leading orders.

In earlier work, Tsamopoulos and Brown [52] derived such an asymptotic
solution up to second order in wave amplitude for the first four shape modes
in the freely oscillating regime (without radial motion or external forcing).
Their analysis reveals the natural wave frequencies at first order (i.e. in
the linear approximation) as well as effects from nonlinear wave interactions
at higher order. More intriguing phenomena can occur if both volume and shape
modes are coupled together, especially under resonant or near-resonant
conditions, with one or more shape modes being excited. On one hand, radial
oscillations driven by an external acoustic field can trigger shape
distortions via parametric instability. This stability problem for a pulsating
spherical bubble has been extensively investigated at lowest order in the
asymptotic procedure, where the volume mode typically obeys the RP equation
while the distortion modes satisfy a Mathieu-type equation with non-autonomous
coefficients depending on the radial motion [3; 47]. Exponential growth is
predicted for this instability at the linear level, but the inclusion of
higher-order nonlinear corrections has been found to promote its saturation
[22; 31]. On the other hand, second-order interactions among shape modes can
transfer energy to the volume mode, potentially leading to a monopole emission
of sound [40] or an erratic motion of the bubble [2; 11].

Aside from theoretical results based on asymptotic solutions or reduced
models, there is also a large literature on direct numerical simulations of
bubble dynamics, owing to advances in computer power and numerical methods
over the last few decades. Typically solving the Navier–Stokes equations and
using interface reconstruction techniques on Cartesian grids, these studies
have considered a variety of configurations involving multiple bubbles and
complex flow conditions or complex surrounding domains, e.g. [4; 14; 19; 20;
35] to cite a few references. While such computations have produced impressive
results, being able to handle extreme situations such as bubble merging or
break-up, they may not be the best option when it comes to accurately
simulating shape oscillations of a bubble because the computed solution is
particularly prone to numerical diffusion. Moreover, their computational cost
is usually high because they require solving the governing equations over the
entire physical domain. If advective inertial effects are assumed to be small
compared to viscous effects, then the Navier–Stokes equations can be
simplified to those for Stokes flow and boundary integral methods have been
developed for their simulation, e.g. in view of medical applications to
microbubbles for drug delivery in blood flow [21]. However, the flow linearity
associated with this approach makes it inadequate for the description of
nonlinear bubble deformations.

Closer to the present problem, direct numerical simulations based on nonlinear
potential flow theory have also been performed via boundary integral and
conformal mapping methods [38; 43; 48; 50; 53; 58]. These exploit coordinate
transformations that are motivated by complex analysis or by the bubble
geometry, and they only involve dynamical variables representing the bubble
surface. Both methods can resolve strong shape distortions (without change in
topology), with boundary integral techniques being more adaptable to various
flow configurations, while the conformal mapping approach is restricted to the
two-dimensional setting. In particular, the latter method enables the exact
derivation of steady solutions or the time-dependent computation of rotating
and standing waves for nonlinear oscillations on a bubble surface, similar to
Crapper’s solutions for the capillary water wave problem [9; 13; 54]. Using a
series expansion in Legendre polynomials without assuming small amplitudes,
McDougald and Leal [42] solved numerically the nonlinear system of ordinary
differential equations for the radial and shape modes of a three-dimensional
axisymmetric bubble. However, the computational cost was high and these
authors only examined the interaction of a few leading modes.

In this paper, we present a combination of mathematical and numerical results
on the surface dynamics of a gas bubble immersed in a two-dimensional liquid
of infinite extent. Volume and shape oscillations of the bubble are considered
under forced or unforced conditions by a far-field pressure. The starting
point is the Hamiltonian formulation for nonlinear potential flow around a
deformable bubble under the influence of surface tension without gravity, as
proposed by Benjamin [1]. The corresponding conjugate variables turn out to be
natural choices of surface quantities for this free-boundary problem, thus
allowing for dimensionality reduction. In this theoretical framework, our new
contributions include:

  1. 1.

Detailed restatement of this Hamiltonian formulation by introducing the
Dirichlet–Neumann operator (DNO) to express the full governing equations as a
closed system in terms of surface variables alone. In doing so, the connection
between these governing equations and the conserved Hamiltonian is also
clarified. Specifics of this two-dimensional problem were not treated by
Benjamin [1].

  2. 2.

Taylor series representation of the DNO associated with the Laplace problem in
polar coordinates. Each term in this Taylor series is determined explicitly by
a recursion formula involving the bubble deformations relative to a quiescent
circular state.

  3. 3.

Analytical calculation of a second-order Stokes wave solution for steadily
rotating shape oscillations in the weakly nonlinear regime, without volume
change. Asymptotic expressions are obtained for the surface displacement,
velocity potential and angular wave speed.

  4. 4.

Derivation of a two-dimensional inviscid version of the RP equation in the
purely circular geometry. Its independent Hamiltonian structure and linear
dispersion relation for radial pulsations driven by a uniform pressure field
are also established.

  5. 5.

Development of an accurate and efficient numerical scheme for the direct
simulation of nonlinear bubble deformations, coupling shape oscillations and
radial pulsations together. It is applicable to surface distortions of
moderate steepness, nonetheless the computations can be performed up to an
arbitrary order of nonlinearity with an arbitrary number of shape modes in an
automatic manner via the fast Fourier transform.

  6. 6.

Extensive tests to examine the numerical convergence of the DNO. Validation
against predictions by the RP equation is also provided in the purely circular
case under uniform pressure forcing.

  7. 7.

Computations of nonlinear shape oscillations represented by steadily rotating
waves and time-periodic standing waves on the bubble surface in the unforced
regime without volume change. For the latter solutions, comparison with
existing laboratory measurements is presented on the frequency of mode-2
oscillations as a function of the bubble’s maximum aspect ratio.

While the DNO has been a common tool in the formulation, analysis and
simulation of the water wave problem [7; 8; 25; 46; 55], its applications have
focused mostly on cases where the reference geometry is rectangular with
Cartesian coordinates. We are only aware of a few exceptions. De la Llave and
Panayotaros [10] considered nonlinear gravity waves on the surface of a sphere
and derived a series expansion for this operator in terms of spherical
coordinates. Their study was strictly mathematical and did not produce any
numerical result. Guyenne and Părău [30] adopted a similar approach in the
axisymmetric cylindrical setting to compute solitary waves on the surface of a
ferrofluid jet. This formalism has also been applied to scattering problems in
acoustics and electromagnetics with non-Cartesian (e.g. polar or spherical)
parameterizations of the irregular domain [15; 45]. However, these
investigations have been restricted to linear time-harmonic waves in the
presence of a stationary object with a fixed shape. To our knowledge, it is
the first time here that the DNO is introduced to describe this nonlinear
hydrodynamic problem with a moving boundary in polar coordinates. Furthermore,
unlike all previous applications, the DNO in this case contains an additional
component to allow for radial pulsations of the bubble (i.e. volume changes
such as compression or dilatation).

While the problem under consideration is two-dimensional in space, this is not
a limitation of the mathematical formulation and numerical procedure that we
advocate to analyze it. The present study contributes to their development in
this idealized situation before tackling the more general three-dimensional
case in the future. This two-dimensional problem is challenging and
interesting in its own right, but has not been examined much in previous
modeling work as opposed to e.g. the three-dimensional axisymmetric
configuration. Indeed, experimental observations have revealed possible
significant differences in bubble dynamics between the two- and three-
dimensional geometries, e.g. regarding the condensation, collapse and
sonoluminescence of explosive bubbles [12]. Even for the three-dimensional
axisymmetric problem which essentially reduces to a planar analysis, there are
qualitative and quantitative differences as compared to the purely two-
dimensional case. The present results may thus serve as benchmark solutions to
test other mathematical or numerical models in this area. Lastly,
understanding the dynamics of two-dimensional bubbles may be relevant in
applications to microfluidics where the technology is based on small thin
devices (e.g. PDMS chips). Recent research has investigated ways to remove
unwanted air bubbles from microfluidic systems [32] or, on the contrary, to
exploit their presence for various purposes such as micropumps, micromixers,
microvalves and microactuators [37].

The remainder of this paper is organized as follows. In Section 2, we present
the mathematical formulation of this two-dimensional problem on the nonlinear
deformations of a gas bubble, and we elaborate on its Hamiltonian structure by
introducing the DNO to accomplish the reduction to surface variables. In
Section 3, we derive a Taylor series representation for the DNO with two
distinct components to enable both shape distortions and volume variations of
the bubble. In Section 4, we discuss various analytical approximations
including a second-order Stokes wave solution for steadily rotating shape
oscillations and a RP model for purely circular pulsations. In Section 5, we
describe the numerical methods for space discretization and time integration
to solve the full governing equations in Hamiltonian form. In Section 6, we
show numerical tests on the convergence of the DNO as well as direct
computations of rotating and standing waves on the bubble surface. For
standing waves, a comparison with laboratory measurements and other
theoretical results is provided. Compression-dilatation cycles under the
excitation of a far-field pressure are also simulated and validated against
predictions by the RP model. Finally, concluding remarks are given in Section
7.

##  2 Mathematical formulation

###  2.1 Governing equations

We consider the motion of a single gas bubble immersed in a two-dimensional
liquid (e.g. water) spanning a domain ${\cal D}$ of infinite extent. Given the
geometry of this problem, we adopt a polar coordinate system $(r,\theta)$
whose origin coincides with the bubble center such that

${\cal D}=\mathbb{R}^{2}\setminus\\{0\leq\theta<2\pi,0\leq
r<R+\eta(\theta,t)\\}\,,$

where $R>0$ is the radius of the unperturbed circular bubble and
$\eta(\theta,t)$ denotes the surface perturbation relative to this simple
configuration at any time $t$. We disregard any motion of the bubble center
and restrict our attention to its surface dynamics, which implies in
particular that body forces such as gravity or buoyancy are neglected. The
exterior flow is assumed to be incompressible, inviscid and irrotational so
that the fluid velocity is given by ${\bf u}(r,\theta,t)=\nabla\varphi$ where
the velocity potential $\varphi(r,\theta,t)$ satisfies the Laplace equation

$\nabla^{2}\varphi=0\,,\quad\mbox{in}\quad{\cal D}\,,$ (1)

with $\nabla=(\partial_{r},r^{-1}\partial_{\theta})^{\top}$ being the spatial
gradient in polar coordinates. At the free surface

${\cal S}=\\{0\leq\theta<2\pi,r=s(\theta,t)=R+\eta(\theta,t)\\}\,,$

there are two boundary conditions, namely the kinematic condition

$\partial_{t}\eta=\partial_{r}\varphi-s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\varphi)\,,\quad\mbox{on}\quad{\cal
S}\,,$ (2)

and the dynamic (or Bernoulli’s) condition

$\partial_{t}\varphi=-\frac{1}{2}|\nabla\varphi|^{2}+\frac{\sigma}{\rho}\kappa+\frac{1}{\rho}\Delta
p\,,\quad\mbox{on}\quad{\cal S}\,,$ (3)

where

$\kappa=\frac{s^{2}+2(\partial_{\theta}\eta)^{2}-s\partial_{\theta}^{2}\eta}{\big{(}s^{2}+(\partial_{\theta}\eta)^{2}\big{)}^{3/2}}\,,$
(4)

denotes the mean curvature at any point on the free surface. Keep in mind that
$\partial_{t}s=\partial_{t}\eta$, $\partial_{\theta}s=\partial_{\theta}\eta$
and similarly for higher derivatives in (2)–(4). The parameters $\sigma$ and
$\rho$ represent the surface tension and fluid density respectively. Typical
values are $\sigma=75$ dyn cm-1 and $\rho=1$ g cm-3 for water. Because the
bubble interior is not empty, the contribution $\Delta p=p_{\infty}-p_{B}$ in
(3) denotes the difference between the fluid pressure $p_{\infty}$ in the far
field and the bubble pressure $p_{B}$ exerted by the internal gas. The
competition between these two pressure disturbances drives the bubble
deformations. More details on $p_{\infty}$ and $p_{B}$ will be provided in a
subsequent section. Finally, the vanishing condition

$|\nabla\varphi|\to 0\,,\quad\mbox{as}\quad r\to+\infty\,,$ (5)

is imposed in the far field. Note that the boundary conditions in $\theta$ are
naturally periodic in this geometric configuration. The present problem is
nonlocal and nonlinear due to the dependence on ${\cal S}$ and associated
boundary conditions (2)–(3).

An interesting situation arises in the absence of pressure disturbances, i.e.
$\Delta p={\rm constant}$, which may be set to $\Delta p=0$ without loss of
generality through a gauge transformation on the velocity potential $\varphi$
[7]. If so, this system of equations possesses a number of invariants of
motion, notably the energy

$H=\int_{\cal D}\frac{1}{2}|\nabla\varphi|^{2}dA+\int_{\cal
S}\frac{\sigma}{\rho}d\ell=K+P\,,$ (6)

where $dA=rdrd\theta$ is the elementary area over ${\cal D}$ and
$d\ell=\sqrt{s^{2}+(\partial_{\theta}\eta)^{2}}d\theta$ is the elementary
arclength along ${\cal S}$. The first term $K$ in (6) is the kinetic part
associated with fluid motion, while the second term $P$ is the potential part
associated with surface tension. Other invariants include the mean surface
level

$Q=\int_{0}^{2\pi}s\,d\theta\,,$ (7)

the mass (or volume)

$V=\int_{\cal D}dA\,,$ (8)

and the angular momentum (or impulse)

$I=\int_{\cal D}({\bf r}\times\nabla\varphi)\,dA\,,$ (9)

where ${\bf r}=r\,{\bf e}_{r}$ denotes the position vector. It turns out that

$I=\int_{\cal
D}(\partial_{\theta}\varphi)\,rdrd\theta=\int_{s}^{+\infty}\left(\varphi\big{|}_{0}^{2\pi}\right)rdr=0\,,$

by virtue of the periodic boundary conditions in $\theta$. Furthermore,
because the fluid mass as defined in (8) is clearly infinite, it may be
substituted by its complement over the bounded area of the gas bubble, more
specifically

$V=\int_{0}^{2\pi}\int_{0}^{s}rdrd\theta=\int_{0}^{2\pi}\frac{1}{2}s^{2}d\theta\,,$
(10)

which is also conserved over time.

###  2.2 Hamiltonian formulation

As shown by Benjamin [1], Eqs. (1)–(5) with $\Delta p=0$ can be re-expressed
as a canonical Hamiltonian system

$\partial_{t}\begin{pmatrix}\eta\\\ \xi\end{pmatrix}=\begin{pmatrix}0&-1\\\
1&0\end{pmatrix}\begin{pmatrix}\partial_{\eta}H\\\
\partial_{\xi}H\end{pmatrix}\,,$ (11)

in terms of the two conjugate variables $\eta(\theta,t)$ and

$\xi(\theta,t)=\varphi(s(\theta,t),\theta,t)\,,$ (12)

the latter being the trace of the velocity potential evaluated at the free
surface ${\cal S}$. Note the sign difference in the symplectic matrix of (11)
as compared to the standard canonical form (see e.g. the water wave problem
[57]), which is explained by the fact that ${\cal S}$ is the inner boundary of
the fluid domain ${\cal D}$ with respect to the coordinate system. The
Hamiltonian $H$ in (11) coincides with the energy (6). Because Benjamin [1]
did not specifically cover this two-dimensional case in polar coordinates, we
will elaborate further on it here. In doing so, we provide a new perspective
on this problem by introducing the Dirichlet–Neumann operator (DNO)

$G(\eta)\xi=(-1,s^{-1}\partial_{\theta}\eta)^{\top}\cdot\nabla\varphi\big{|}_{r=s}=\sqrt{1+s^{-2}(\partial_{\theta}\eta)^{2}}(\nabla\varphi\cdot{\bf
n})\big{|}_{r=s}\,,$ (13)

which is the singular integral operator that takes Dirichlet data $\xi$ on
${\cal S}$, solves the Laplace equation (1) subject to (5), and returns the
corresponding Neumann data (i.e. the normal fluid velocity on ${\cal S}$). The
outward unit vector ${\bf n}$ normal to ${\cal S}$ is given by

${\bf
n}=\frac{(-1,s^{-1}\partial_{\theta}\eta)^{\top}}{\sqrt{1+s^{-2}(\partial_{\theta}\eta)^{2}}}=\frac{s}{\sqrt{s^{2}+(\partial_{\theta}\eta)^{2}}}(-1,s^{-1}\partial_{\theta}\eta)^{\top}\,,$
(14)

which points in the opposite $r$-direction, consistent with the sign
difference in the canonical form (11) as mentioned above. The choice of
definition (13) for the DNO will become more evident in subsequent
calculations. It is a linear operator in $\xi$ but depends nonlinearly on
$\eta$. Similar to studies on water waves [8; 24; 27; 55], an advantage of
using the DNO is that the dependence on $\eta$ and $\xi$, which are the two
conjugate variables, is made more explicit in the equations of motion (2)–(3)
and in the Hamiltonian (6). Indeed, the variable $\eta$ appears as part of the
domain of integration ${\cal D}$ in the original form (6) of $H$, while the
variable $\xi$ does not even appear explicitly there, and so it is not clear
at this stage what the variational derivatives of $H$ with respect to
$(\eta,\xi)$ mean as indicated by (11).

In the following, we will write (11) more explicitly via the DNO. With the
definitions (12) and (13) at hand, we get the identities

$\partial_{t}\xi=\partial_{t}\varphi+(\partial_{t}\eta)(\partial_{r}\varphi)\big{|}_{r=s}\,,\quad\partial_{\theta}\xi=\partial_{\theta}\varphi+(\partial_{\theta}\eta)(\partial_{r}\varphi)\big{|}_{r=s}\,,$
(15)

by differentiating (12) and applying the chain rule. Together with

$G(\eta)\xi=-\partial_{r}\varphi+s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\varphi)\big{|}_{r=s}\,,$
(16)

from (13), we deduce

$\partial_{\theta}\varphi\big{|}_{r=s}=\partial_{\theta}\xi-(\partial_{\theta}\eta)(\partial_{r}\varphi)\big{|}_{r=s}\,,$
(17)

and

$\displaystyle\partial_{r}\varphi\big{|}_{r=s}$ $\displaystyle=$
$\displaystyle-G(\eta)\xi+s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\varphi)\big{|}_{r=s}=-G(\eta)\xi+s^{-2}(\partial_{\theta}\eta)\big{[}\partial_{\theta}\xi-(\partial_{\theta}\eta)(\partial_{r}\varphi)\big{]}\big{|}_{r=s}\,,$
$\displaystyle=$
$\displaystyle-G(\eta)\xi+s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\xi)-s^{-2}(\partial_{\theta}\eta)^{2}(\partial_{r}\varphi)\big{|}_{r=s}\,,$

which implies

$\partial_{r}\varphi\big{|}_{r=s}=\frac{1}{1+s^{-2}(\partial_{\theta}\eta)^{2}}\big{[}s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\xi)-G(\eta)\xi\big{]}\,.$
(18)

Then substituting (18) back into (15) and (17), we find

$\displaystyle\partial_{\theta}\varphi\big{|}_{r=s}$ $\displaystyle=$
$\displaystyle\partial_{\theta}\xi-\frac{\partial_{\theta}\eta}{1+s^{-2}(\partial_{\theta}\eta)^{2}}\big{[}s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\xi)-G(\eta)\xi\big{]}\,,$
(19) $\displaystyle=$
$\displaystyle\frac{1}{1+s^{-2}(\partial_{\theta}\eta)^{2}}\big{[}\partial_{\theta}\xi+(\partial_{\theta}\eta)G(\eta)\xi\big{]}\,,$

and

$\displaystyle\partial_{t}\varphi\big{|}_{r=s}$ $\displaystyle=$
$\displaystyle\partial_{t}\xi-(\partial_{t}\eta)(\partial_{r}\varphi)\big{|}_{r=s}\,,$
(20) $\displaystyle=$
$\displaystyle\partial_{t}\xi+\frac{G(\eta)\xi}{1+s^{-2}(\partial_{\theta}\eta)^{2}}\big{[}s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\xi)-G(\eta)\xi\big{]}\,,$

where we have used the relation

$\partial_{t}\eta=\partial_{r}\varphi-s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\varphi)\big{|}_{r=s}=-G(\eta)\xi\,,$

according to the kinematic condition (2) and the definition (16) of the DNO.
Combining the squares of (18) and (19), i.e.

$\displaystyle(\partial_{r}\varphi)^{2}\big{|}_{r=s}$ $\displaystyle=$
$\displaystyle\frac{1}{\big{(}1+s^{-2}(\partial_{\theta}\eta)^{2}\big{)}^{2}}\Big{[}s^{-4}(\partial_{\theta}\eta)^{2}(\partial_{\theta}\xi)^{2}-2s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\xi)G(\eta)\xi+\big{(}G(\eta)\xi\big{)}^{2}\Big{]}\,,$
$\displaystyle(\partial_{\theta}\varphi)^{2}\big{|}_{r=s}$ $\displaystyle=$
$\displaystyle\frac{1}{\big{(}1+s^{-2}(\partial_{\theta}\eta)^{2}\big{)}^{2}}\Big{[}(\partial_{\theta}\xi)^{2}+2(\partial_{\theta}\eta)(\partial_{\theta}\xi)G(\eta)\xi+(\partial_{\theta}\eta)^{2}\big{(}G(\eta)\xi\big{)}^{2}\Big{]}\,,$

yields

$\displaystyle(\partial_{r}\varphi)^{2}+s^{-2}(\partial_{\theta}\varphi)^{2}\big{|}_{r=s}$
$\displaystyle=$
$\displaystyle\frac{1}{\big{(}1+s^{-2}(\partial_{\theta}\eta)^{2}\big{)}^{2}}\Big{[}s^{-2}\big{(}1+s^{-2}(\partial_{\theta}\eta)^{2}\big{)}(\partial_{\theta}\xi)^{2}+\big{(}1+s^{-2}(\partial_{\theta}\eta)^{2}\big{)}\big{(}G(\eta)\xi\big{)}^{2}\Big{]}\,,$
$\displaystyle=$
$\displaystyle\frac{1}{1+s^{-2}(\partial_{\theta}\eta)^{2}}\Big{[}s^{-2}(\partial_{\theta}\xi)^{2}+\big{(}G(\eta)\xi\big{)}^{2}\Big{]}\,.$

Finally, putting all these expressions together in (2)–(3) leads to a closed
system of two equations

$\displaystyle\partial_{t}\eta$ $\displaystyle=$ $\displaystyle-G(\eta)\xi\,,$
(21) $\displaystyle\partial_{t}\xi$ $\displaystyle=$
$\displaystyle-\frac{1}{2\big{(}1+s^{-2}(\partial_{\theta}\eta)^{2}\big{)}}\Big{[}s^{-2}(\partial_{\theta}\xi)^{2}+2s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\xi)G(\eta)\xi-\big{(}G(\eta)\xi\big{)}^{2}\Big{]}$
(22) $\displaystyle+\frac{\sigma}{\rho}\kappa\,,$

for $\eta$ and $\xi$, which gives a more explicit form of (11) and is
completely equivalent to the original nonlinear formulation (1)–(5) (with
$\Delta p=0$). Recall that, by definition, the solution of the Laplace
equation (1) subject to the far-field vanishing condition (5) is encoded in
the DNO. We also point out that, in addition to being Hamiltonian, Eqs.
(21)–(22) is a lower-dimensional version of (1)–(5) in terms of surface
variables alone. It is thus more appealing for mathematical analysis and
numerical simulation. In this spirit, a numerical model will be proposed here
to solve (21)–(22). Of course, the contribution $\Delta p$ can be added to
(22) as in (3), but then the Hamiltonian structure would not be preserved in
general.

###  2.3 Variational derivatives of the Hamiltonian

We now apply Benjamin’s formalism [1] together with the DNO introduced above
to show the connection between the variational derivatives of $H$ in (11) and
the lower-dimensional system (21)–(22) that we just derived. For this purpose,
we rewrite the kinetic part of the Hamiltonian (6) as

$K=\int_{\cal D}\frac{1}{2}|\nabla\varphi|^{2}dA=\int_{\cal
D}\frac{1}{2}\big{[}\nabla\cdot(\varphi\nabla\varphi)-\varphi\nabla^{2}\varphi\big{]}dA\,,$

where the last term vanishes due to (1). Then using the divergence theorem, we
obtain

$\displaystyle 2K$ $\displaystyle=$ $\displaystyle\int_{\cal
D}\nabla\cdot(\varphi\nabla\varphi)\,dA=\int_{\cal
S}\varphi(\nabla\varphi\cdot{\bf n})\,d\ell\,,$ $\displaystyle=$
$\displaystyle\int_{\cal S}\varphi(\nabla\varphi\cdot{\bf
n})\sqrt{s^{2}+(\partial_{\theta}\eta)^{2}}\,d\theta=\int_{0}^{2\pi}\xi
G(\eta)\xi\,sd\theta\,,$

by virtue of (12) and (13). Following Benjamin [1], a crucial step is to
introduce the weighted element $d\mu=\sqrt{{\cal G}/{\cal
G}_{11}}d\theta=sd\theta$ along some interval ${\cal C}$, where ${\cal
G}=r^{2}$ and ${\cal G}_{11}=1$ are respectively the determinant and first
diagonal entry of the metric tensor for polar coordinates, evaluated on ${\cal
S}$. The Hamiltonian then reads

$H=K+P=\frac{1}{2}\int_{\cal C}\xi
G(\eta)\xi\,d\mu+\frac{\sigma}{\rho}\int_{0}^{2\pi}\sqrt{s^{2}+(\partial_{\theta}\eta)^{2}}\,d\theta\,,$
(23)

in terms of the DNO and the two conjugate variables $(\eta,\xi)$. Given such a
functional, the variation $\partial_{u}H$ is defined via the Gâteaux
derivative

$\langle\partial_{u}H,v\rangle=\frac{d}{d\lambda}H(u+\lambda
v)\Big{|}_{\lambda=0}\,,$ (24)

with respect to the inner product $\langle f,h\rangle=\int_{\cal C}fh\,d\mu$
for any functions $u$, $v$ in some Hilbert space. From the alternate form (23)
of $H$ and the definition (24) of the variational derivative, we deduce

$\displaystyle\partial_{\xi}H$ $\displaystyle=$
$\displaystyle\partial_{\xi}K=G(\eta)\xi\,,$ (25)
$\displaystyle\partial_{\eta}K$ $\displaystyle=$
$\displaystyle-\frac{1}{2\big{(}1+s^{-2}(\partial_{\theta}\eta)^{2}\big{)}}\Big{[}s^{-2}(\partial_{\theta}\xi)^{2}+2s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\xi)G(\eta)\xi-\big{(}G(\eta)\xi\big{)}^{2}\Big{]}\,.$

Equation (25) for $\partial_{\xi}H$ follows from the self-adjointness of the
DNO [15; 45; 46]. Establishing the other identity for $\partial_{\eta}K$ is a
lengthy calculation but is closely related to that for the water wave problem.
Therefore, we refer the reader to [7] for more details.

Regarding the capillary component, we note that $\partial_{u}H$ can be
expressed alternatively as

$\int_{\cal
C}(\partial_{u}H)v\,d\mu=\langle\partial_{u}H,v\rangle=\int_{0}^{2\pi}(\partial_{u}H)^{*}v\,d\theta\,,$

which implies that $\partial_{u}H=s^{-1}(\partial_{u}H)^{*}$ where
$(\partial_{u}H)^{*}$ is the variational derivative (24) with respect to the
inner product $\langle f,h\rangle=\int_{0}^{2\pi}fh\,d\theta$. Applying this
idea to the potential part of (23), we find

$\frac{\rho}{\sigma}(\partial_{\eta}P)^{*}=\frac{s^{3}+2s(\partial_{\theta}\eta)^{2}-s^{2}\partial_{\theta}^{2}\eta}{\big{(}s^{2}+(\partial_{\theta}\eta)^{2}\big{)}^{3/2}}\,,$

hence $\partial_{\eta}P=s^{-1}(\partial_{\eta}P)^{*}=\sigma\kappa/\rho$ with
$\kappa$ as presented in (4). Combining these results for
$\partial_{\eta}H=\partial_{\eta}K+\partial_{\eta}P$ together with (25) for
$\partial_{\xi}H$ shows the equivalence between (11) and (21)–(22).

###  2.4 Variational derivatives of the momentum

For completeness and reference in a subsequent section, we use the same
formalism to evaluate variational derivatives of the angular momentum (9) with
respect to the surface variables $(\eta,\xi)$. Similar to $H$, we first
rewrite $I$ explicitly in terms of these variables as follows

$I=\int_{\cal D}\partial_{\theta}\varphi\,dA=\int_{\cal
D}\nabla(r^{2}\varphi)\cdot\nabla\theta\,dA=\int_{\cal
D}\nabla\cdot(r^{2}\varphi\nabla\theta)\,dA\,,$

by virtue of the identity

$\nabla\cdot(r^{2}\varphi\nabla\theta)=\nabla(r^{2}\varphi)\cdot\nabla\theta+r^{2}\varphi\nabla^{2}\theta\,,$

together with the fact that $\nabla^{2}\theta=0$. The divergence theorem
implies

$\displaystyle I$ $\displaystyle=$ $\displaystyle\int_{\cal
S}s^{2}\varphi(\nabla\theta\cdot{\bf n})\,d\ell=\int_{\cal
S}s^{2}\varphi(\nabla\theta\cdot{\bf
n})\sqrt{s^{2}+(\partial_{\theta}\eta)^{2}}\,d\theta\,,$ (26) $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}\xi\partial_{\theta}\eta\,sd\theta=\int_{\cal
C}\xi\partial_{\theta}\eta\,d\mu\,,$

via (12) and (14), after recognizing that $\nabla\theta=(0,r^{-1})^{\top}$. It
is then a straightforward calculation to show that

$\partial_{\xi}I=\partial_{\theta}\eta\,,\quad\partial_{\eta}I=-\partial_{\theta}\xi\,,$
(27)

based on the definition (24). For $\partial_{\eta}I$, because Eq. (26)
involves $\partial_{\theta}\eta$ and also depends on $\eta$ through $s$, an
integration by parts in $\theta$ is required to obtain its expression in (27).

Likewise, we easily see that

$\partial_{\xi}Q=0\,,\quad\partial_{\eta}Q=1\,.$ (28)

##  3 Dirichlet–Neumann operator

The question now is how to calculate the DNO given the boundary $\eta$ and
Dirichlet data $\xi$ at any time $t$, in order to solve (21)–(22). The
approach that we advocate here is of boundary perturbation type and is based
on the fact that the DNO is analytic in $\eta$ under certain (relatively mild)
regularity conditions, as first shown by Coifman and Meyer [6]. Rigorous
proofs of this analyticity property can be found in [45] and [15] for the two-
dimensional circular and three-dimensional spherical cases, respectively.

The starting point is the Laplace equation (1) in polar coordinates

$\partial_{r}^{2}\varphi+r^{-1}\partial_{r}\varphi+r^{-2}\partial_{\theta}^{2}\varphi=0\,,$

for which we consider two different elementary solutions

$\varphi_{1}=\ln r\,,\quad\varphi_{2}=r^{-n}e^{{\rm i}n\theta}\,,$ (29)

with $n\in\mathbb{N}$ (non-negative integers). The time dependence is omitted
because the Laplace equation is solved in a frozen domain at each time $t$.
The first solution $\varphi_{1}$ describes bubble motion in the $r$-direction
alone (i.e. purely circular compression or dilatation) while the second
solution $\varphi_{2}$ allows for wave development on the bubble surface (i.e.
shape distortions) in the $\theta$-direction. Both phenomena are physically
relevant, hence the importance to examine both solutions. In particular, the
radial motion is inherent to this physical problem (even in the absence of
waves on the bubble surface) owing to surface tension associated with the
bubble curvature.

Note that both $\varphi_{1}$ and $\varphi_{2}$ satisfy the periodic boundary
conditions in $\theta$ as well as the vanishing condition (5) as
$r\to+\infty$. This explains why the option $\varphi_{2}$ with $n<0$ is ruled
out. There is no singularity at $r=0$ for either $\varphi_{1}$ or
$\varphi_{2}$ because this location is outside of the fluid domain ${\cal D}$.

Next we exploit the analyticity property of the DNO by seeking a Taylor series
representation in $\eta$, namely

$G(\eta)=\sum_{j=0}^{+\infty}G_{j}(\eta)\,,$ (30)

about the reference circular geometry $\eta=0$, i.e. $r=R$ with $R$ being a
constant. By construction, each term $G_{j}$ in (30) is homogeneous of degree
$j$ in $\eta$ and its action as an operator on any basis function $e^{{\rm
i}n\theta}$ (associated with Fourier mode or wavenumber $n$) can be determined
in a recursive manner. This derivation of the DNO relies on the choice of
harmonic function (i.e. solution to the Laplace equation) and, because
$\varphi_{1}$ and $\varphi_{2}$ are two different harmonic functions, we will
present details for each case separately. We will first consider $\varphi_{2}$
which has both $(r,\theta)$-dependences, before inspecting the univariate case
with $\varphi_{1}$. Incidentally, we recognize that $\varphi_{2}$ should be
specified as a real-valued function (corresponding to a velocity potential).
However, its complex-valued form in (29) is more convenient for the purposes
of our derivation and, because the DNO is a linear operator in $\xi$, it is
sufficient to write $e^{{\rm i}n\theta}$. Contributions from its complex
conjugate $e^{-{\rm i}n\theta}$ would be redundant.

Substituting (29) and (30) for $\varphi_{2}$ into (16) yields

$\displaystyle\Big{(}\sum_{j=0}^{+\infty}G_{j}(\eta)\Big{)}(R+\eta)^{-|n|}e^{{\rm
i}n\theta}$ $\displaystyle=$ $\displaystyle\Big{[}(R+\eta)^{-|n|-1}|n|+{\rm
i}(\partial_{\theta}\eta)(R+\eta)^{-|n|-2}n\Big{]}e^{{\rm i}n\theta}\,,$
$\displaystyle\Big{(}\sum_{j=0}^{+\infty}G_{j}(\eta)\Big{)}\Big{(}1+\frac{\eta}{R}\Big{)}^{-|n|}e^{{\rm
i}n\theta}$ $\displaystyle=$ $\displaystyle
R^{-1}\Big{[}\Big{(}1+\frac{\eta}{R}\Big{)}^{-|n|-1}|n|-\Big{(}D\frac{\eta}{R}\Big{)}\Big{(}1+\frac{\eta}{R}\Big{)}^{-|n|-2}n\Big{]}e^{{\rm
i}n\theta}\,,$

where $D=-{\rm i}\,\partial_{\theta}$ and the absolute value $|n|$ is employed
to enforce exponent positivity for $n\in\mathbb{Z}$ (all integers over the
full spectrum), as motivated above. In terms of the binomial expansion

$\Big{(}1+\frac{\eta}{R}\Big{)}^{-n}=\sum_{\ell=0}^{+\infty}\Big{(}\frac{\eta}{R}\Big{)}^{\ell}C_{\ell}^{-n}\,,$

with coefficient

$C_{\ell}^{-n}=\begin{pmatrix}-n\\\
\ell\end{pmatrix}=(-1)^{\ell}\begin{pmatrix}n+\ell-1\\\
\ell\end{pmatrix}=(-1)^{\ell}\,\frac{n(n+1)(n+2)\cdots(n+\ell-1)}{\ell!}\,,$

this equation becomes

$\Big{(}\sum_{j=0}^{+\infty}G_{j}(\eta)\Big{)}\Big{(}\sum_{\ell=0}^{+\infty}\Big{(}\frac{\eta}{R}\Big{)}^{\ell}C_{\ell}^{-|n|}\Big{)}e^{{\rm
i}n\theta}=R^{-1}\Big{[}\sum_{\ell=0}^{+\infty}\Big{(}\frac{\eta}{R}\Big{)}^{\ell}C_{\ell}^{-|n|-1}|n|-\Big{(}D\frac{\eta}{R}\Big{)}\sum_{\ell=0}^{+\infty}\Big{(}\frac{\eta}{R}\Big{)}^{\ell}C_{\ell}^{-|n|-2}n\Big{]}e^{{\rm
i}n\theta}\,.$

By inspection, the first term for $j=0$ in (30) is given by

$G_{0}e^{{\rm i}n\theta}=\frac{|n|}{R}e^{{\rm i}n\theta}\,,$

which can be viewed as the Fourier symbol of the pseudo-differential operator

$G_{0}=\frac{|D|}{R}\,,$ (31)

acting on any function that has a Fourier series decomposition in $\theta$.
Recall that, by definition, the Fourier symbol associated with $D$ is $n$. For
$j>0$, collecting terms of the same degree in $\eta$ leads to

$\displaystyle G_{j}(\eta)e^{{\rm i}n\theta}$ $\displaystyle=$
$\displaystyle-\sum_{\ell=0}^{j-1}G_{\ell}(\eta)\Big{(}\frac{\eta}{R}\Big{)}^{j-\ell}C_{j-\ell}^{-|n|}e^{{\rm
i}n\theta}$
$\displaystyle+R^{-1}\Big{[}\Big{(}\frac{\eta}{R}\Big{)}^{j}C_{j}^{-|n|-1}|n|-\Big{(}\frac{\eta}{R}\Big{)}^{j-1}\Big{(}D\frac{\eta}{R}\Big{)}C_{j-1}^{-|n|-2}n\Big{]}e^{{\rm
i}n\theta}\,,$

which can also be expressed symbolically as

$\displaystyle G_{j}(\eta)$ $\displaystyle=$
$\displaystyle-\sum_{\ell=0}^{j-1}G_{\ell}(\eta)\Big{(}\frac{\eta}{R}\Big{)}^{j-\ell}C_{j-\ell}^{-|D|}$
$\displaystyle+R^{-1}\Big{[}\Big{(}\frac{\eta}{R}\Big{)}^{j}C_{j}^{-|D|-1}|D|-\Big{(}\frac{\eta}{R}\Big{)}^{j-1}\Big{(}D\frac{\eta}{R}\Big{)}C_{j-1}^{-|D|-2}D\Big{]}\,.$

By convention, unless parentheses are specified, any operator acts on all
functions to its right in the same term. Pseudo-differential operators like
$D$, $|D|$ or $C_{j}^{-|D|}$ are also called Fourier multipliers due to their
multiplicative action in the Fourier space (dual to the physical
$\theta$-space). Then recognizing that

$D\Big{(}\frac{\eta}{R}\Big{)}^{j}D=j\Big{(}\frac{\eta}{R}\Big{)}^{j-1}\Big{(}D\frac{\eta}{R}\Big{)}D+\Big{(}\frac{\eta}{R}\Big{)}^{j}D^{2}\,,\quad
D^{2}=|D|^{2}\,,$

together with the property

$C_{j}^{-|D|-1}=-C_{j-1}^{-|D|-2}\frac{|D|+1}{j}=-\frac{1}{j}C_{j-1}^{-|D|-2}|D|-\frac{1}{j}C_{j-1}^{-|D|-2}\,,$

for the binomial coefficient, we arrive at

$G_{j}(\eta)=-\frac{1}{j}\Big{[}D\Big{(}\frac{\eta}{R}\Big{)}^{j}\frac{D}{R}+\Big{(}\frac{\eta}{R}\Big{)}^{j}G_{0}\Big{]}C_{j-1}^{-|D|-2}-\sum_{\ell=0}^{j-1}G_{\ell}(\eta)\Big{(}\frac{\eta}{R}\Big{)}^{j-\ell}C_{j-\ell}^{-|D|}\,.$

Finally, we invoke the self-adjointness of the DNO to obtain

$G_{j}(\eta)=-\frac{1}{j}\Big{[}D\Big{(}\frac{\eta}{R}\Big{)}^{j}\frac{D}{R}+\Big{(}\frac{\eta}{R}\Big{)}^{j}G_{0}\Big{]}C_{j-1}^{-|D|-2}-\sum_{\ell=0}^{j-1}C_{j-\ell}^{-|D|}\Big{(}\frac{\eta}{R}\Big{)}^{j-\ell}G_{\ell}(\eta)\,,$
(32)

after reversing the sequence of application for the various operations in the
summation above [15; 45; 46]. The reason for doing so will be made more clear
when discussing the numerical scheme in a subsequent section.

We now repeat this calculation by using (29) with $\varphi_{1}$. For this
purpose, we slightly tweak the harmonic function as follows
$\varphi_{1}=\ln(R_{m}/r)$ where the constant radius $R_{m}$ is chosen
sufficiently large such that $R_{m}\gg R$. Because $\varphi_{1}$ only depends
on $r$, it may be associated with the zeroth Fourier mode $n=0$ in $\theta$.
Accordingly, Eq. (16) reduces to

$\displaystyle\Big{(}\sum_{j=0}^{+\infty}G_{j}(\eta)\Big{)}\mathbb{P}_{0}\ln\left(\frac{R_{m}}{R+\eta}\right)e^{{\rm
i}n\theta}$ $\displaystyle=$
$\displaystyle\mathbb{P}_{0}\frac{1}{R+\eta}e^{{\rm i}n\theta}\,,$
$\displaystyle\Big{(}\sum_{j=0}^{+\infty}G_{j}(\eta)\Big{)}\mathbb{P}_{0}\left[\ln\left(\frac{R_{m}}{R}\right)-\ln\left(1+\frac{\eta}{R}\right)\right]e^{{\rm
i}n\theta}$ $\displaystyle=$ $\displaystyle
R^{-1}\mathbb{P}_{0}\left(1+\frac{\eta}{R}\right)^{-1}e^{{\rm i}n\theta}\,,$

where $\mathbb{P}_{0}$ denotes the projection onto the zeroth mode $n=0$, i.e.
the DNO only affects this specific mode here. Then Taylor expanding about
$\eta=0$ gives

$\Big{(}\sum_{j=0}^{+\infty}G_{j}(\eta)\Big{)}\mathbb{P}_{0}\left[\ln\left(\frac{R_{m}}{R}\right)-\sum_{\ell=1}^{+\infty}\frac{(-1)^{\ell+1}}{\ell}\Big{(}\frac{\eta}{R}\Big{)}^{\ell}\right]e^{{\rm
i}n\theta}=R^{-1}\mathbb{P}_{0}\sum_{j=0}^{+\infty}(-1)^{j}\Big{(}\frac{\eta}{R}\Big{)}^{j}e^{{\rm
i}n\theta}\,.$

For $j=0$, we readily infer

$\displaystyle G_{0}\mathbb{P}_{0}\ln\left(\frac{R_{m}}{R}\right)e^{{\rm
i}n\theta}$ $\displaystyle=$ $\displaystyle R^{-1}\mathbb{P}_{0}e^{{\rm
i}n\theta}\,,$ $\displaystyle G_{0}\mathbb{P}_{0}e^{{\rm i}n\theta}$
$\displaystyle=$ $\displaystyle\frac{1}{R\ln(R_{m}/R)}\mathbb{P}_{0}e^{{\rm
i}n\theta}\,,$

which can be written symbolically as

$G_{0}\mathbb{P}_{0}=\frac{1}{R\ln(R_{m}/R)}\mathbb{P}_{0}\,.$ (33)

For $j>0$, by identifying terms of the same degree in $\eta$, we get

$G_{j}(\eta)\mathbb{P}_{0}e^{{\rm
i}n\theta}=\frac{1}{R\ln(R_{m}/R)}\mathbb{P}_{0}(-1)^{j}\Big{(}\frac{\eta}{R}\Big{)}^{j}e^{{\rm
i}n\theta}+\frac{1}{\ln(R_{m}/R)}\sum_{\ell=0}^{j-1}G_{\ell}(\eta)\mathbb{P}_{0}\frac{(-1)^{j-\ell+1}}{j-\ell}\Big{(}\frac{\eta}{R}\Big{)}^{j-\ell}e^{{\rm
i}n\theta}\,,$

and thus we can define its action symbolically as

$G_{j}(\eta)\mathbb{P}_{0}=\frac{1}{R\ln(R_{m}/R)}\mathbb{P}_{0}(-1)^{j}\Big{(}\frac{\eta}{R}\Big{)}^{j}+\frac{1}{\ln(R_{m}/R)}\sum_{\ell=0}^{j-1}\mathbb{P}_{0}\frac{(-1)^{j-\ell+1}}{j-\ell}\Big{(}\frac{\eta}{R}\Big{)}^{j-\ell}G_{\ell}(\eta)\,,$
(34)

after invoking again the self-adjointness property (to reverse the sequence of
operations in the inner summation). It now makes sense why $R_{m}\gg R$ was
introduced in such a way, otherwise we would end up with a possible
logarithmic singularity when $R_{m}=R$ (or equivalently when $R=1$ in some
dimensional or dimensionless units for the choice $\varphi_{1}=\ln r$). The
range $R_{m}\gg R$ is preferable over $0<R_{m}\ll R$ because it is inside the
fluid domain ${\cal D}$ and thus complies with the mathematical formulation
(1)–(5) for the fluid flow. Moreover, in that upper range, $\ln(R_{m}/R)$ is
positive and its variation is relatively mild.

Equations (31)–(32) or (33)–(34) provide recursion formulas to evaluate the
DNO in its series form (30) given $\eta$ and $\xi$. Depending on the physical
situation under consideration, each set may be used separately, i.e. either
formulas (31)–(32) or (33)–(34). More generally, they may be combined by
superposition as follows. The zeroth-order operator ($j=0$) takes the form

$G_{0}=G_{0}^{(1)}+G_{0}^{(2)}=\frac{|D|}{R}+\frac{1}{R\ln(R_{m}/R)}\mathbb{P}_{0}\,,$
(35)

from (31) and (33), while the higher-order contributions ($j>0$) read

$G_{j}(\eta)=G_{j}^{(1)}(\eta)+G_{j}^{(2)}(\eta)\,,$ (36)

with $G_{j}^{(1)}$ and $G_{j}^{(2)}$ determined by (32) and (34) respectively.
This superposition enables the coupling between radial pulsations and shape
distortions. Note in particular how $G_{0}^{(1)}$ and $G_{0}^{(2)}$ complement
each other in (35) over the full spectrum $n\in\mathbb{Z}$, considering that
$G_{0}^{(1)}=|D|/R$ has a trivial effect on the zeroth mode $n=0$, while
$G_{0}^{(2)}=R^{-1}\mathbb{P}_{0}/\ln(R_{m}/R)$ has a nontrivial effect there.
Compared to the water wave problem in the perturbed rectangular geometry where
the DNO has been commonly involved [25], the additional component
$G_{j}^{(2)}$ represents a distinctive new feature of the present formulation
in polar coordinates.

##  4 Approximate regimes

Several analytical approximations can be made in order to simplify the full
nonlinear problem.

###  4.1 Linearized problem

We first investigate the linearized problem for small-amplitude deformations
about $\eta=0$, relative to the equilibrium state $r=R$ (with $R$ fixed). As
can be deduced from (3), by setting all derivatives to zero, this equilibrium
corresponds to

$\frac{\sigma}{R}+\Delta p=0\,,$ (37)

in light of the expansion

$\kappa=\frac{1}{R}-\frac{1}{R^{2}}(\partial_{\theta}^{2}\eta+\eta)+O(\eta^{2})\,.$

Equation (37), also known as the Young–Laplace equation, states that the
pressure jump $\Delta p$ across the bubble boundary is balanced by the
constant curvature due to surface tension. In this case, Eqs. (1)–(5) via
their lower-dimensional form (21)–(22) simplify to

$\partial_{t}\eta=-G_{0}\xi\,,\quad\partial_{t}\xi=-\frac{\sigma}{\rho
R^{2}}(\partial_{\theta}^{2}\eta+\eta)\,,$ (38)

with only contributions from up to first order in $(\eta,\xi)$, which can be
combined as

$\partial_{t}^{2}\eta-\frac{\sigma}{\rho
R^{2}}G_{0}(\partial_{\theta}^{2}\eta+\eta)=0\,,$ (39)

and similarly for $\xi$. Having rotating wave solutions of the form
$\eta,\xi\sim e^{{\rm i}(n\theta-\omega_{0}t)}$ in mind here, it is sufficient
to choose $G_{0}=G_{0}^{(1)}$ from (31). Equation (39) then implies the linear
dispersion relation

$\omega_{0}^{2}=\frac{\sigma}{\rho R^{3}}|n|(n^{2}-1)\,,$ (40)

between the angular frequency $\omega_{0}$ and wavenumber $n$. For all
$n\in\mathbb{Z}$, we see that $\omega_{0}^{2}\geq 0$ and thus $\omega_{0}$ is
real. Assuming $n$, $\omega_{0}>0$, the angular phase speed can be derived
from (40) as

$c_{0}=\frac{\omega_{0}}{n}=\sqrt{\frac{\sigma(n^{2}-1)}{\rho R^{3}n}}\,.$
(41)

Dyachenko [13] obtained the same result using a different mathematical
formulation of this problem. The modes $n=\pm 1$ would correspond to a rigid
translation of the bubble as in the three-dimensional case [52], but this
mechanism is ignored here.

###  4.2 Second-order Stokes approximation for rotating waves

Extending this linear analysis to the weakly nonlinear regime, we focus our
attention on wave solutions of the form

$\eta(\theta,t)=\eta(\Theta)\,,\quad\xi(\theta,t)=\xi(\Theta)\,,\quad\Theta=\theta-c\,t\,,$
(42)

which rotate counter-clockwise at constant angular speed $c$ under static
pressure (37). Equations (21)–(22) then reduce to the nonlinear system of
ordinary differential equations (ODEs)

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
c\,\eta^{\prime}-G(\eta)\xi\,,$ (43) $\displaystyle 0$ $\displaystyle=$
$\displaystyle c\,\xi^{\prime}-\frac{1}{2(1+s^{-2}\eta^{\prime
2})}\Big{[}s^{-2}\xi^{\prime
2}+2s^{-2}\eta^{\prime}\xi^{\prime}G(\eta)\xi-\big{(}G(\eta)\xi\big{)}^{2}\Big{]}$
(44) $\displaystyle+\frac{\sigma}{\rho}\Big{(}\kappa-\frac{1}{R}\Big{)}\,,$

by the chain rule, where the primes denote differentiation with respect to
$\Theta$. Similarly, Fourier multipliers with respect to $\theta$ in $G(\eta)$
are replaced by their counterparts with respect to $\Theta$. The change of
variables (42) is equivalent to reformulating this problem in a reference
frame rotating at constant speed $c$.

Alternatively, from the Hamiltonian viewpoint, such rotating solutions can be
interpreted as fixed points of the renormalized Hamiltonian
$\widehat{H}=H-c\,I-\sigma Q/(\rho R)$ where $Q$ (mean surface level) and $I$
(angular momentum) are also invariants of motion as defined by (7) and (9). In
other words, these solutions satisfy

$\begin{pmatrix}0&-1\\\
1&0\end{pmatrix}\begin{pmatrix}\partial_{\eta}\widehat{H}\\\
\partial_{\xi}\widehat{H}\end{pmatrix}=\begin{pmatrix}0\\\ 0\end{pmatrix}\,,$

which coincide with (43)–(44). This equivalence can be checked directly by
invoking (25), (27) and (28), modulo the substitution of $\Theta$ for
$\theta$.

Along the lines of Stokes theory for the water wave problem [17], we seek
$(\eta,\xi)$ perturbatively via an asymptotic series

$\eta=\varepsilon\eta_{0}+\varepsilon^{2}\eta_{1}+\cdots\,,\quad\xi=\varepsilon\xi_{0}+\varepsilon^{2}\xi_{1}+\cdots\,,$
(45)

together with $c=c_{0}+\varepsilon c_{1}+\cdots$, where the perturbation
parameter $\varepsilon\ll 1$ is a dimensionless measure of the wave amplitude
(e.g. the wave steepness). Inserting (45) in (43)–(44) and equating
coefficients in factor of the same power in $\varepsilon$, we find

$c_{0}\eta_{0}^{\prime}-G_{0}\xi_{0}=0\,,\quad
c_{0}\xi_{0}^{\prime}-\frac{\sigma}{\rho
R^{2}}(\eta_{0}^{\prime\prime}+\eta_{0})=0\,,$

for $(\eta_{0},\xi_{0})$ at first order $O(\varepsilon)$, which is nothing but
the homogeneous linear system (38) in terms of $\Theta$ with $c_{0}$ given by
(41). We may take a general solution of the form

$\eta_{0}=a_{0}\cos(n\Theta)\,,\quad\xi_{0}=-a_{0}c_{0}R\sin(n\Theta)\,,\quad
n>1\,,$ (46)

which sets the amplitude and phase of this Stokes wave. As stated in the
previous section, only the contributions (31)–(32) are implemented here for
the DNO. At second order $O(\varepsilon^{2})$, we obtain the inhomogeneous
linear system

$\displaystyle c_{0}\eta_{1}^{\prime}-G_{0}\xi_{1}$ $\displaystyle=$
$\displaystyle-c_{1}\eta_{0}^{\prime}+G_{1}(\eta_{0})\xi_{0}\,,$ (47)
$\displaystyle c_{0}\xi_{1}^{\prime}-\frac{\sigma}{\rho
R^{2}}(\eta_{1}^{\prime\prime}+\eta_{1})$ $\displaystyle=$ $\displaystyle-
c_{1}\xi_{0}^{\prime}+\frac{1}{2}\Big{[}\frac{1}{R^{2}}\xi_{0}^{\prime
2}-(G_{0}\xi_{0})^{2}\Big{]}$ $\displaystyle-\frac{\sigma}{\rho
R^{3}}\Big{(}\eta_{0}^{2}+2\eta_{0}\eta_{0}^{\prime\prime}+\frac{1}{2}\eta_{0}^{\prime
2}\Big{)}\,,$

for $(\eta_{1},\xi_{1})$ with inhomogeneous terms depending on
$(\eta_{0},\xi_{0})$. This system can be merged into the single equation

$\displaystyle\left(1-\frac{\sigma G_{0}}{c_{0}^{2}\rho
R^{2}}\right)\eta_{1}^{\prime\prime}-\frac{\sigma G_{0}}{c_{0}^{2}\rho
R^{2}}\eta_{1}$ $\displaystyle=$
$\displaystyle-\frac{c_{1}}{c_{0}}\eta_{0}^{\prime\prime}+\frac{G_{1}(\eta_{0}^{\prime})}{c_{0}}\xi_{0}+\frac{G_{1}(\eta_{0})}{c_{0}}\xi_{0}^{\prime}$
$\displaystyle-\frac{G_{0}}{c_{0}^{2}}\left[c_{1}\xi_{0}^{\prime}-\frac{1}{2R^{2}}\xi_{0}^{\prime
2}+\frac{1}{2}(G_{0}\xi_{0})^{2}\right.$
$\displaystyle\left.+\frac{\sigma}{\rho
R^{3}}\Big{(}\eta_{0}^{2}+2\eta_{0}\eta_{0}^{\prime\prime}+\frac{1}{2}\eta_{0}^{\prime
2}\Big{)}\right]\,,$

for $\eta_{1}$ alone after eliminating $\xi_{1}$. Note that

$G_{0}=\frac{|D|}{R}\,,\quad
G_{1}(\eta)=|D|\frac{\eta}{R}G_{0}-D\frac{\eta}{R^{2}}D-\frac{\eta}{R}G_{0}\,,$

according to (31)–(32). Substitution of (46) for $(\eta_{0},\xi_{0})$ leads to

$\displaystyle\left(1-\frac{\sigma G_{0}}{c_{0}^{2}\rho
R^{2}}\right)\eta_{1}^{\prime\prime}-\frac{\sigma G_{0}}{c_{0}^{2}\rho
R^{2}}\eta_{1}$ $\displaystyle=$
$\displaystyle\frac{2c_{1}a_{0}n^{2}}{c_{0}}\cos(n\Theta)+\frac{a_{0}^{2}n^{2}}{R}\cos(2n\Theta)+\frac{a_{0}^{2}n^{3}}{R}\cos(2n\Theta)$
(48) $\displaystyle+\frac{5a_{0}^{2}n^{3}\sigma}{2c_{0}^{2}\rho
R^{4}}\cos(2n\Theta)-\frac{a_{0}^{2}n\sigma}{c_{0}^{2}\rho
R^{4}}\cos(2n\Theta)\,.$

In doing so, we have applied the identities

$D\begin{pmatrix}\cos(n\Theta)\\\ \sin(n\Theta)\end{pmatrix}={\rm
i}\,n\begin{pmatrix}\sin(n\Theta)\\\
-\cos(n\Theta)\end{pmatrix}\,,\quad|D|\begin{pmatrix}\cos(n\Theta)\\\
\sin(n\Theta)\end{pmatrix}=n\begin{pmatrix}\cos(n\Theta)\\\
\sin(n\Theta)\end{pmatrix}\,,$

for $n>0$. Because the second-order corrections are supposed to be bound to
the first-order components, we choose the particular form

$\eta_{1}=a_{1}\cos(2n\Theta)\,,$

as suggested by the right-hand side of (48). Then equating the coefficients in
factor of $\cos(n\Theta)$ and $\cos(2n\Theta)$ implies $c_{1}=0$ for
$a_{0}\neq 0$ (nontrivial solutions) and

$a_{1}=\frac{a_{0}^{2}(2n^{3}+7n^{2}-2n-4)}{4R(2n^{2}+1)}\,,$

respectively. Once $\eta_{1}$ is known, $\xi_{1}$ can be determined from (47),
namely

$G_{0}\xi_{1}=c_{0}\eta_{1}^{\prime}-G_{1}(\eta_{0})\xi_{0}=-c_{0}n\left(2a_{1}+\frac{a_{0}^{2}}{2R}\right)\sin(2n\Theta)\,,$
(49)

which yields

$\xi_{1}=-c_{0}\left(a_{1}R+\frac{a_{0}^{2}}{4}\right)\sin(2n\Theta)\,,$

as dictated by the right-hand side of (49). In principle, this procedure can
be pursued up to an arbitrary order at the expense that the equations become
increasingly more complicated (involving higher-order terms $G_{j}$) as the
level of approximation rises. Collecting all these expressions in (45) given
$a_{0}$ and $n$ provides an explicit weakly nonlinear estimate, accurate up to
second order in Stokes theory, for steadily rotating waves at angular speed
$c$. It may serve as a validation test for two-dimensional numerical solvers
including the present one, when applied to wave speeds near $c_{0}$. Numerical
illustrations will be shown in a subsequent section.

###  4.3 Rayleigh–Plesset model

The Rayleigh–Plesset (RP) equation and its variants have been an important
model to understand the dynamics of cavitation bubbles [16; 49]. For a single
bubble in an infinite body of incompressible fluid and under the assumption of
spherical symmetry, a second-order ODE can be derived for the time evolution
of its radius. Because this model has been mostly considered for three-
dimensional spherical applications in the literature, we find it suitable to
present its detailed derivation in the two-dimensional circular configuration
so as to make this paper sufficiently self-contained. We also restrict
ourselves to the inviscid limit, although the RP equation typically includes a
viscous term.

Based on the potential-flow formulation (1)–(5) in the special case with
circular symmetry, we take a harmonic function of the form

$\varphi=B\ln\left(\frac{R_{m}}{r}\right)\,,\quad\partial_{r}\varphi=-\frac{B}{r}\,,$
(50)

as in Sec. 3, and define the free surface by ${\cal
S}=\\{0\leq\theta<2\pi,r=s(t)\\}$. Here we prefer to employ the variable $s$
(full surface deformation) rather than $\eta$ (surface perturbation relative
to $r=R$) in order to comply with the typical formulation of the RP equation.
The kinematic boundary condition (2) then simplifies to
$ds/dt=\partial_{r}\varphi$ on ${\cal S}$, which implies that $B=-s\,ds/dt$.
The dynamic boundary condition (3) reduces to

$\partial_{t}\varphi=-\frac{1}{2}(\partial_{r}\varphi)^{2}+\frac{\sigma}{\rho
s}+\frac{1}{\rho}\Delta p\,,$ (51)

where

$\partial_{r}\varphi=\frac{ds}{dt}\,,\quad\partial_{t}\varphi=-\left[\left(\frac{ds}{dt}\right)^{2}+s\frac{d^{2}s}{dt^{2}}\right]\ln\left(\frac{R_{m}}{s}\right)\,,\quad\mbox{on}\quad{\cal
S}\,.$

To complete this boundary condition, we need to specify the pressure jump
$\Delta p=p_{\infty}-p_{B}$. The bubble pressure is assumed to obey the static
law

$p_{B}=\left(p_{\infty}+\frac{\sigma}{R}\right)\left(\frac{V_{0}}{V}\right)^{\gamma}\,,\quad\gamma\geq
0\,,$

for a polytropic process so that it satisfies (37) at equilibrium $r=R$, where
the circular bubble volume is given by $V=\pi s^{2}$ at any time $t$ according
to (10) and $V_{0}=\pi R^{2}$ denotes its equilibrium value [40]. With all
these substitutions, Eq. (51) becomes

$s\frac{d^{2}s}{dt^{2}}\ln\left(\frac{R_{m}}{s}\right)+\left(\frac{ds}{dt}\right)^{2}\left[\ln\left(\frac{R_{m}}{s}\right)-\frac{1}{2}\right]+\frac{\sigma}{\rho
s}+\frac{p_{\infty}}{\rho}-\frac{1}{\rho}\left(p_{\infty}+\frac{\sigma}{R}\right)\left(\frac{R}{s}\right)^{2\gamma}=0\,.$
(52)

This two-dimensional version of the RP equation is somewhat different from its
typical three-dimensional counterpart, due to the different forms of harmonic
function between these two geometries.

In three dimensions and in the absence of viscosity, the RP equation is known
to possess a canonical Hamiltonian structure assuming a uniform far-field
pressure $p_{\infty}$ [16]. Such a mathematical formulation is suitable for a
phase-plane stability analysis of solutions about fixed points when varying
the parameter $p_{\infty}$ [41]. We now reveal a similar Hamiltonian structure
for the two-dimensional model (52) by following Feng and Leal [16] with a key
adjustment of their change of variables. If we set

$\mathrm{q}=\ln\left(\frac{R_{m}}{s}\right)\,,\quad\mathrm{p}=-s^{3}\frac{ds}{dt}\ln\left(\frac{R_{m}}{s}\right)\,,$
(53)

as motivated by (50), then Eq. (52) is equivalent to the system of two first-
order ODEs

$\displaystyle\frac{d\mathrm{q}}{dt}$ $\displaystyle=$
$\displaystyle\frac{\mathrm{p}\,e^{4\mathrm{q}}}{\mathrm{q}\,R_{m}^{4}}\,,$
(54) $\displaystyle\frac{d\mathrm{p}}{dt}$ $\displaystyle=$
$\displaystyle-\frac{2\mathrm{p}^{2}\,e^{4\mathrm{q}}}{\mathrm{q}\,R_{m}^{4}}+\frac{\mathrm{p}^{2}\,e^{4\mathrm{q}}}{2\mathrm{q}^{2}\,R_{m}^{4}}+\frac{\sigma
R_{m}}{\rho}e^{-\mathrm{q}}$ (55)
$\displaystyle+\frac{p_{\infty}R_{m}^{2}}{\rho}e^{-2\mathrm{q}}-\frac{R_{m}^{2}}{\rho}\left(p_{\infty}+\frac{\sigma}{R}\right)\left(\frac{R}{R_{m}}\right)^{2\gamma}e^{-2(1-\gamma)\mathrm{q}}\,.$

The first equation is simply an identity from the definitions (53) of
$\mathrm{q}$ and $\mathrm{p}$. The second equation is an alternate form of
(52) in terms of these auxiliary variables. Moreover, this system can be
expressed as

$\frac{d\mathrm{q}}{dt}=\partial_{\mathrm{p}}{\cal
H}\,,\quad\frac{d\mathrm{p}}{dt}=-\partial_{\mathrm{q}}{\cal H}\,,$

which are canonical Hamiltonian equations for the two conjugate coordinates
$\mathrm{q}$ and $\mathrm{p}$, associated with the Hamiltonian

${\cal H}=\frac{\mathrm{p}^{2}\,e^{4\mathrm{q}}}{2q\,R_{m}^{4}}+\frac{\sigma
R_{m}}{\rho}e^{-\mathrm{q}}+\frac{p_{\infty}R_{m}^{2}}{2\rho}e^{-2\mathrm{q}}-\frac{R_{m}^{2}}{2\rho(1-\gamma)}\left(p_{\infty}+\frac{\sigma}{R}\right)\left(\frac{R}{R_{m}}\right)^{2\gamma}e^{-2(1-\gamma)\mathrm{q}}\,.$
(56)

All these relations can be verified by direct calculation given (53) and (56).
This Hamiltonian structure for (52) is of interest in its own right because it
is not directly related to the Hamiltonian formulation (11) of the full
system. It is not a straightforward extension of that in the three-dimensional
case due to the presence of the logarithm function $\ln(R_{m}/s)$ in this two-
dimensional problem. As a result, the new Hamiltonian (56) and associated
evolution equations (54)–(55) exhibit both rational and exponential
dependences on $q$.

Note that the possible singularity in (56) at $\gamma=1$ for an isothermal
process is consistent with the fact that such conditions typically occur when
a system is in contact with an outside thermal reservoir, allowing e.g. for
heat exchange in order to maintain a constant temperature. Accordingly, energy
conservation is not expected in this system. On the other hand, for a general
isentropic process ($\gamma\neq 1$) where the conditions are adiabatic and
reversible, energy conservation as represented by (56) is realizable.

The linearized problem for the RP equation (52) about equilibrium $r=R$ also
requires an examination. Substituting $s(t)=R+\eta(t)$ into (52) and retaining
terms of up to first order in $\eta$ gives

$\frac{d^{2}\eta}{dt^{2}}+\Omega_{0}^{2}\,\eta=0\,,$ (57)

for a harmonic oscillator at frequency

$\Omega_{0}^{2}=\frac{2\gamma\left(p_{\infty}+\frac{\sigma}{R}\right)-\frac{\sigma}{R}}{\rho
R^{2}\ln(R_{m}/R)}\,,$ (58)

if the far-field pressure $p_{\infty}$ is assumed to be uniform. This
expression for the linear fundamental frequency of radial oscillations shares
similarities with the three-dimensional version [47]. Not surprisingly, Eqs.
(57)–(58) coincide exactly with (39) for $\gamma=0$ and under circular
symmetry (i.e. $\theta$-invariance) when $G_{0}=G_{0}^{(2)}$ according to
(33). The fact that $\Omega_{0}^{2}<0$ in this situation ($\gamma=0$,
$R_{m}\gg R$) is indicative of an instability triggering bubble collapse due
to surface tension in the absence of internal pressure.

##  5 Numerical methods

We present a numerical scheme to compute the DNO and discretize (21)–(22) in
space and time. We emphasize that, while the DNO is given by a Taylor series
expansion, the full system (21)–(22) is directly solved for the time evolution
problem (as opposed to invoking a reduced asymptotic model).

###  5.1 Space discretization

Considering the periodic boundary conditions in $\theta$, we use a pseudo-
spectral method to discretize the DNO and equations of motion (21)–(22) in
space [5]. This is a natural choice for the computation of $G(\eta)$ because
each term in its Taylor series (30) is evaluated via recursion formulas
(35)–(36) involving concatenations of Fourier multipliers with powers of
$\eta/R$. More specifically, both functions $\eta$ and $\xi$ are expressed as
truncated Fourier series

$\begin{pmatrix}\eta_{j}\\\
\xi_{j}\end{pmatrix}=\sum_{n=-N/2}^{N/2-1}\begin{pmatrix}\widehat{\eta}_{n}\\\
\widehat{\xi}_{n}\end{pmatrix}e^{{\rm
i}n\theta_{j}}\,,\quad\theta_{j}=\frac{2\pi}{N}j\,,\quad j=0,\dots,N-1\,.$

Applications of spatial derivatives or Fourier multipliers are performed in
the Fourier space, while nonlinear products are calculated in the physical
space on a regular grid of $N$ collocation points. For example, if we wish to
apply the zeroth-order operator $G_{0}^{(1)}$ to a function $\xi$ in the
physical space, we implement it as follows

$G_{0}^{(1)}\xi={\cal F}^{-1}\left(\frac{|n|}{R}{\cal F}(\xi)\right)={\cal
F}^{-1}\left(\frac{|n|}{R}\widehat{\xi}_{n}\right)\,.$

Similarly, the projection $\mathbb{P}_{0}$ in $G_{j}^{(2)}$ can be computed as

$\mathbb{P}_{0}\xi={\cal F}^{-1}\big{(}\widehat{\xi}_{0}\big{)}\,,$

where ${\cal F}$ (resp. ${\cal F}^{-1}$) denotes the direct (resp. inverse)
Fourier transform. All operations from the physical to Fourier space and vice
versa are carried out via the fast Fourier transform (FFT).

The Taylor series of the DNO is also truncated to a finite number of terms

$G(\eta)\simeq G^{(M)}(\eta)=\sum_{j=0}^{M}G_{j}(\eta)\,,$ (59)

for which the choice of truncation order $M$ will be discussed in more detail
in a subsequent section when showing convergence tests. Due to the analyticity
property, a small number of terms would be sufficient (typically $M<10\ll N$)
to achieve highly accurate results. For this computation, the adjoint formulas
(32) and (34) are more efficient than the original ones (those before
reversing the sequence of inner operations) because they allow us to save and
re-use the $G_{\ell}$’s on $\xi$ as vectors, without having to re-compute
these operators at each order $j$ when applied to concatenations of Fourier
multipliers and powers of $\eta/R$. In this adjoint form, the computational
cost for evaluating (59) can be estimated to be $O(M^{2}N\ln N)$ via the FFT.
This point highlights another advantage of the present approach where the
Taylor series representation of the DNO together with the recursive
calculation of its constitutive terms avoids setting up and solving a large or
dense matrix system (related to the Laplace problem (1)–(5)), which contrasts
with the classical strategy in other numerical solvers such as boundary
integral methods or volumetric finite-difference/element methods [26].

###  5.2 Time integration

Time integration of (21)–(22) is performed in the Fourier space so that linear
terms can be solved exactly by the integrating factor technique. To do so, we
first split these equations into

$\partial_{t}\mbox{\boldmath$v$}={\cal L}\mbox{\boldmath$v$}+{\cal
N}(\mbox{\boldmath$v$})\,,$ (60)

for $\mbox{\boldmath$v$}=(\eta,\xi)^{\top}$, where the linear part is given by

${\cal L}\mbox{\boldmath$v$}=\begin{pmatrix}0&-G_{0}^{(1)}\\\
-\frac{\sigma}{\rho
R^{2}}(\partial_{\theta}^{2}+1)&0\end{pmatrix}\begin{pmatrix}\eta\\\
\xi\end{pmatrix}\,,$ (61)

according to (38), and the nonlinear part takes the form ${\cal
N}(\mbox{\boldmath$v$})=({\cal N}_{1},{\cal N}_{2})^{\top}$ with

$\displaystyle{\cal N}_{1}$ $\displaystyle=$
$\displaystyle-\big{(}G(\eta)-G_{0}^{(1)}\big{)}\xi\,,$ $\displaystyle{\cal
N}_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2\big{(}1+s^{-2}(\partial_{\theta}\eta)^{2}\big{)}}\Big{[}s^{-2}(\partial_{\theta}\xi)^{2}+2s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\xi)G(\eta)\xi-\big{(}G(\eta)\xi\big{)}^{2}\Big{]}$
$\displaystyle+\frac{\sigma}{\rho}\Big{[}\kappa+\frac{1}{R^{2}}(\partial_{\theta}^{2}\eta+\eta)\Big{]}+\frac{1}{\rho}\Delta
p\,.$

The pressure difference $\Delta p$ is also included in ${\cal N}_{2}$ to
tackle out-of-equilibrium regimes. The subtraction of $-G_{0}^{(1)}\xi$ in
${\cal N}_{1}$ and of $-\sigma(\partial_{\theta}^{2}\eta+\eta)/(\rho R^{2})$
in ${\cal N}_{2}$ is meant to compensate for their presence in the linear part
${\cal L}\mbox{\boldmath$v$}$. Curvature terms are known to cause stiffness in
the numerical simulation of capillarity-dominated flows [34], hence treating
their linear contributions exactly will help mitigate this issue.

More specifically, we take the Fourier transform of (60) and make the change
of variables

$\widehat{\mbox{\boldmath$v$}}_{n}(t)=\Phi_{n}(t)\widehat{\mbox{\boldmath$w$}}_{n}(t)\,,$
(62)

in terms of the integrating factor

$\Phi_{n}(t)=\begin{pmatrix}\cos(\omega_{0}t)&-\frac{|n|}{\omega_{0}R}\sin(\omega_{0}t)\\\
\frac{\sigma(n^{2}-1)}{\omega_{0}R^{2}}\sin(\omega_{0}t)&\cos(\omega_{0}t)\end{pmatrix}\,,$
(63)

where $\omega_{0}$ is given by (40). This integrating factor is the
fundamental matrix of the linear system

$\partial_{t}\widehat{\mbox{\boldmath$v$}}_{n}=\widehat{\cal
L}_{n}\widehat{\mbox{\boldmath$v$}}_{n}=\begin{pmatrix}0&-\frac{|n|}{R}\\\
\frac{\sigma}{\rho
R^{2}}(n^{2}-1)&0\end{pmatrix}\begin{pmatrix}\widehat{\eta}_{n}\\\
\widehat{\xi}_{n}\end{pmatrix}\,,$

associated with the linear part (61). It thus encodes the exact linear
solution of (60). Because $\omega_{0}$ vanishes at $n=\\{0,\pm 1\\}$, there
may be some indetermination at these modes for either one of the off-diagonal
entries in $\Phi_{n}(t)$, as suggested by (63). Via l’Hôpital’s rule, we find

$\Phi_{0}(t)=\begin{pmatrix}1&0\\\ -\frac{\sigma t}{R^{2}}&1\end{pmatrix}\,,$

for $n=0$, and

$\Phi_{\pm 1}(t)=\begin{pmatrix}1&-\frac{t}{R}\\\ 0&1\end{pmatrix}\,,$

for $n=\pm 1$.

This leads to a strictly nonlinear system

$\partial_{t}\widehat{\mbox{\boldmath$w$}}_{n}=\Phi_{n}(t)^{-1}\widehat{\cal
N}_{n}\big{[}\Phi_{n}(t)\widehat{\mbox{\boldmath$w$}}_{n}\big{]}\,,$

for the new variable $\widehat{\mbox{\boldmath$w$}}_{n}$, which is solved
numerically in time using the fourth-order Runge–Kutta method with constant
step $\Delta t$. The resulting scheme can be inverted back to
$\widehat{\mbox{\boldmath$v$}}_{n}$ by virtue of (62), yielding

$\widehat{\mbox{\boldmath$v$}}_{n}^{k+1}=\Phi_{n}(\Delta
t)\widehat{\mbox{\boldmath$v$}}_{n}^{k}+\frac{\Delta t}{6}\Phi_{n}(\Delta
t)(f_{1}+2f_{2}+2f_{3}+f_{4})\,,$

for the numerical solution at time $t_{k+1}=t_{k}+\Delta t$, where

$\displaystyle f_{1}$ $\displaystyle=$ $\displaystyle\widehat{\cal
N}_{n}(\widehat{\mbox{\boldmath$v$}}_{n}^{k})\,,$ $\displaystyle f_{2}$
$\displaystyle=$ $\displaystyle\Phi_{n}\left(-\frac{\Delta
t}{2}\right)\widehat{\cal N}_{n}\left[\Phi_{n}\left(\frac{\Delta
t}{2}\right)\left(\widehat{\mbox{\boldmath$v$}}_{n}^{k}+\frac{\Delta
t}{2}f_{1}\right)\right]\,,$ $\displaystyle f_{3}$ $\displaystyle=$
$\displaystyle\Phi_{n}\left(-\frac{\Delta t}{2}\right)\widehat{\cal
N}_{n}\left[\Phi_{n}\left(\frac{\Delta
t}{2}\right)\left(\widehat{\mbox{\boldmath$v$}}_{n}^{k}+\frac{\Delta
t}{2}f_{2}\right)\right]\,,$ $\displaystyle f_{4}$ $\displaystyle=$
$\displaystyle\Phi_{n}(-\Delta t)\,\widehat{\cal N}_{n}\Big{[}\Phi_{n}(\Delta
t)\big{(}\widehat{\mbox{\boldmath$v$}}_{n}^{k}+\Delta
t\,f_{3}\big{)}\Big{]}\,.$

Here we have exploited the fact that $\Phi_{n}(t)$ is a semigroup that
satisfies the properties

$\Phi_{n}(t)^{-1}=\Phi_{n}(-t)\,,\quad\Phi_{n}(t+\tau)=\Phi_{n}(t)\Phi_{n}(\tau)\,.$

These identities can be verified by direct calculation. Again, in this time-
integration process, the FFT makes it possible to go back and forth between
$v$ and $\widehat{\mbox{\boldmath$v$}}_{n}$ in an efficient manner.

###  5.3 De-aliasing and filtering

For pseudo-spectral methods applied to nonlinear problems, numerical errors
may stem from the aliasing phenomenon [5]. In the present algorithm, aliasing
may occur when evaluating the equations of motion and the DNO with the FFT.
The $j$th-order term $G_{j}(\eta)\xi$ involves nonlinearities of degree $j+1$
as indicated by (32) or (34), therefore aliasing may be severe for large $j$.
To deal with this issue, the zero-padding technique is a simple and effective
option. Typically, for a quadratic nonlinearity and given resolution $N$, this
amounts to doubling the size of the discretized spectra of $(\eta,\xi)$ and
setting the Fourier coefficients for the extra modes to zero, consistent with
the $2/3$ rule for de-aliasing in this configuration [5]. Because
nonlinearities are of polynomial or rational type in (21)–(22) and (30), we
accommodate each nonlinear term by breaking it up into successive products of
two functions and by applying the aforementioned de-aliasing technique at each
multiplicative step. We have successfully employed such a procedure in other
physical settings [28; 30; 55] where nonlinear surface waves were simulated.
It comes along with increased memory storage but this turns out not to be a
major concern in the present two-dimensional case via the FFT.

For shape distortions of appreciable steepness, we have found it necessary to
filter the numerical solution in order to stabilize the computation so that it
can run over a sufficiently long period of time. Otherwise, spurious high-
wavenumber oscillations tend to develop, eventually leading to the computation
breakdown. Possible causes for this phenomenon include numerical ill-
conditioning of the DNO in its series form (59) (which will be further
assessed in a subsequent section), or the ill-posed character of the governing
equations. A detailed mathematical analysis on these points is outside the
scope of this study. As a remedy, we apply a hyperviscosity-type filter of the
form

${\exp}\left(-36\left|\frac{n}{N/2}\right|^{36}\right)\,,$

to both Fourier coefficients $\widehat{\eta}_{n}$ and $\widehat{\xi}_{n}$ at
each time step, where $N/2$ is the largest wavenumber resolved by the spatial
discretization. Such a filter has been commonly adopted in numerical
simulations of nonlinear fluid flows by spectral methods [5; 33; 34; 27; 55].
Its smooth but steep behavior near $n=N/2$ ensures that only Fourier
coefficients at high wavenumbers are affected. Accordingly, for a sufficiently
high resolution, this filtering strategy should help suppress spurious
instabilities while preserving the overall solution. Combined with zero-
padding, it contributes further to control of aliasing errors.

##  6 Numerical results

We present numerical results to illustrate the performance of this numerical
model. These include convergence tests on the DNO computation and direct
simulations of time-dependent nonlinear solutions. For the latter, we explore
purely circular motions under transient excitation as well as shape
deformations that rotate steadily or oscillate periodically on the bubble
surface.

All variables are non-dimensionalized with respect to the characteristic
length $R$, mass $\rho R^{3}$ and time $\sqrt{\rho R^{3}/\sigma}$ so that
$R=1$, $\rho=1$ and $\sigma=1$ in the resulting dimensionless equations. For
convenience, we will use the same notations but it is now understood that
values of any variable or parameter are dimensionless.

###  6.1 Convergence tests on the DNO

Because we approximate the DNO in terms of a Taylor series as justified by its
analyticity property, it is suitable to assess its numerical convergence which
is expected to be exponential with respect to $M$ for sufficiently smooth
deformations $\eta$. This question is especially relevant considering that
such a representation has been shown to be ill-conditioned numerically in
related problems [30; 46; 55]. Indeed, the series expansion (59) relies on
cancellation of terms to guarantee convergence but in practice an exact
cancellation never happens due to round-off errors. As indicated by (32), the
higher the order $j$ of $G_{j}$, the stronger the exponentiation power of
$|D|$ in such Fourier multipliers as $C_{j}^{-|D|}$. Therefore, numerical
errors can be amplified dramatically through the recursive process as the
order $j$ increases in (59). Needless to say that this error amplification
will contaminate the numerical solution during the time integration of
(21)–(22), and may precipitate the accuracy deterioration or promote spurious
instabilities.

Here time is frozen and Eqs. (60)–(61) are not solved yet. Instead, we inspect
the numerical convergence of (59) by testing this approximation against an
exact expression that can be derived based on the harmonic function
$\varphi_{2}$ from (29) together with a specific profile of $\eta$. More
precisely, using the real-valued form

$\varphi_{2}=r^{-n}\sin(n\theta)\,,\quad n>0\,,$ (64)

and substituting it into the definition (16), we obtain

$\displaystyle G^{(E)}(\eta)\xi$ $\displaystyle=$
$\displaystyle-\partial_{r}\varphi_{2}+s^{-2}(\partial_{\theta}\eta)(\partial_{\theta}\varphi_{2})\big{|}_{r=s}\,,$
(65) $\displaystyle=$ $\displaystyle
n(R+\eta)^{-n-2}\Big{[}(R+\eta)\sin(n\theta)+(\partial_{\theta}\eta)\cos(n\theta)\Big{]}\,.$

We perform tests for two different surface profiles

$\eta=a_{1}\cos(n\theta)\,,$ (66)

and

$\eta=a_{2}\left[\theta^{4}(2\pi-\theta)^{4}-\frac{128\pi^{8}}{315}\right]\,,$
(67)

in combination with (59) (numerical approximation) or (65) (exact expression)
for the DNO. The choice (66) corresponds to a smooth boundary ($C^{\infty}$
profile) while (67) represents a rougher boundary (with finite smoothness)
[46]. Both profiles are prescribed in such a way that their mean value is
zero. Note also that, as opposed to (66) which is more symmetric, the profile
(67) tends to be more eccentric to the left, with a more ovate shape, as
$a_{2}$ is increased. The maximum distortion is $\|\eta\|_{\infty}=a_{1}$ (at
any crest or trough) for (66), while it is given by
$\|\eta\|_{\infty}=a_{2}187\pi^{8}/315$ for (67) and is achieved at
$\theta=\pi$ (left side of the bubble). Recalling that the DNO (13) is linear
in $\xi$ but depends nonlinearly on $\eta$, any constant coefficient in (64)
would be superfluous while the amplitudes $a_{1}$, $a_{2}$ in (66), (67) are
expected to be relevant.

In all our computations, we set the parameter $R_{m}=1000R$ so that $R_{m}\gg
R$ as mentioned earlier. We have checked that results are insensitive to the
specific value of $R_{m}$ in this range.

Figure 1(a) plots the relative $L^{2}$ error

${\rm
Error}=\frac{\|G^{(M)}(\eta)\xi-G^{(E)}(\eta)\xi\|_{2}}{\|G^{(E)}(\eta)\xi\|_{2}}\,,$
(68)

between (59) and (65) as a function of truncation order $M$ for varying
amplitude $a_{1}$ in the case of a smooth profile (66) with wavenumber $n=2$
and resolution $N=256$ (i.e. angular grid size $\Delta\theta=2\pi/N=0.024$). A
first observation concerns the convergence with respect to $a_{1}$ which is
demonstrated by the lower error curves for smaller $a_{1}$, spanning multiple
orders of magnitude. Given $a_{1}$, exponential convergence seems to take
place over only the first few values of $M$ (as indicated by the near-linear
slope in these semilog plots). Past $M\simeq 2$, the errors quickly stagnate,
though they remain small in general. For large amplitudes, say $a_{1}=0.9$
(which is comparable to the radius $R=1$ of the unperturbed circular bubble),
the errors are large and reach $100$%. After some stagnation, we see a
dramatic error growth past some critical value of $M$ (past $M\simeq 5$ for
$a_{1}=0.9$). This phenomenon is characteristic of the numerical ill-
conditioning for the DNO computation via its series form (59), as reported by
previous studies on free-surface flows in different geometric configurations
[24; 30; 46; 55]. To help the reader visualize the present geometry, Fig. 1(b)
portrays various possible bubble shapes when varying $a_{1}$. For $a_{1}=0.9$
(with $n=2$), the bubble is so strongly pinched along the central vertical
axis that two lobes are generated. The surface profile $\eta$ is clearly an
important factor influencing the accuracy of the DNO computation, as pointed
out earlier, and this is further illustrated in our next tests.

The fact that convergence promptly stagnates with respect to $M$ is likely
another outcome of this numerical ill-conditioning. A similar phenomenon has
been observed in the axisymmetric cylindrical case for surface deformations of
a ferrofluid jet [30], but it differs somewhat from results in the rectangular
geometry with Cartesian coordinates [24; 46; 55] where faster (i.e.
exponential) convergence is achieved over several values of $M$ before a
sudden deterioration occurs. This difference may be explained here by the
action of binomial-type Fourier multipliers like $C_{j}^{-|D|}$ (aside from
$G_{0}$), which are more complicated than the mere power-type $|D|^{j}$ in the
Cartesian coordinate setting. As a consequence, the present recursion formula
(32) for each $G_{j}$ contains more nonlinear terms (without counting
contributions from (33)), and their number increases with $j$. This promotes
stronger amplification of numerical errors, even at low values of $M$.

Figure 2 shows $L^{2}$ errors (68) versus $M$ for varying $n$ (with
$a_{1}=0.1$, $N=256$) and varying $N$ (with $a_{1}=0.1$, $n=2$). Here a
moderately large amplitude $a_{1}=0.1$ is selected as a compromise so we can
expect small errors while possibly capturing ill-conditioning effects. In
fact, quite large surface steepnesses $\varepsilon=(a_{1}n)/R$ are examined in
Fig. 2(a), e.g. $\varepsilon=\\{0.5,1.0,1.2\\}$ for $n=\\{5,10,12\\}$
respectively. Overall, we discern similar features to Fig. 1. The higher the
steepness of $\eta$ (i.e. the larger $a_{1}$ or $n$), the slower the
convergence and the sooner it deteriorates (i.e. the smaller the critical
value of $M$ at which a rapid loss of accuracy occurs). This issue is
especially pronounced when increasing $N$, which is compatible with a previous
statement that the higher the order $j$ of $G_{j}$ (combined with higher
resolution $N$), the more drastically Fourier multipliers such as
$C_{j}^{-|D|}$ can amplify round-off errors. This result suggests that
prescribing unnecessarily fine resolutions is not recommended.

From Fig. 3, the same observations can be made on error plots for the surface
profile (67). Because the bubble shape is oval-like in this case with $\eta$
being independent of $n$, we simply set $n=2$ in (65) and skip any convergence
test by varying $n$. As expected, for this rougher boundary, the numerical
issues appear to be more severe when increasing $a_{2}$ or $N$. For reference,
the maximum distortion associated with $a_{2}=\\{10^{-6},\dots,10^{-4}\\}$ is
$\|\eta\|_{\infty}=\\{5\times 10^{-3},\dots,5\times 10^{-1}\\}$ respectively,
as depicted in Fig. 3(a). Finally, we illustrate in Fig. 4 the effectiveness
of our de-aliasing procedure for both (66) and (67) by comparing the
corresponding $L^{2}$ errors with those obtained from aliased computations
(without zero-padding). Because aliasing typically arises when evaluating
nonlinear terms and affects the high-wavenumber tail of the discretized
spectrum, not surprisingly it accentuates ill-conditioning effects at large
values of $M$ and is more pronounced for higher surface steepnesses.
Representative cases are displayed in Fig. 4 for $(a_{1},n)=(0.1,10)$ and
$(a_{2},n)=(5\times 10^{-5},2)$ with $N=256$.

In summary, it may be inferred from these convergence tests that choosing a
value around $M=6$ (say, $4\leq M\leq 8$) would generally be a good compromise
between accuracy and efficiency for the DNO approximation. Similar values of
$M$ have been successfully implemented in other physical contexts [27; 28; 30;
55].

###  6.2 Comparison with the Rayleigh–Plesset equation

In this section, we further assess the performance of our boundary
perturbation approach by comparing time-dependent simulations of (60)–(61)
against predictions from the RP equation (52). For this purpose, the second-
order nonlinear ODE is transformed to a system of two first-order ODEs

$\displaystyle\frac{ds}{dt}$ $\displaystyle=$ $\displaystyle v\,,$ (69)
$\displaystyle\frac{dv}{dt}$ $\displaystyle=$
$\displaystyle\frac{1}{s\ln(R_{m}/s)}\left[v^{2}\left(\frac{1}{2}-\ln\left(\frac{R_{m}}{s}\right)\right)-\frac{\sigma}{\rho
s}-\frac{p_{\infty}}{\rho}+\frac{1}{\rho}\left(p_{\infty}+\frac{\sigma}{R}\right)\left(\frac{R}{s}\right)^{2\gamma}\right]\,,$
(70)

which is integrated numerically in time via the fourth-order Runge–Kutta
scheme as applied to (60)–(61). Equations (69)–(70) thus lead to a direct
numerical solver for this nonlinear problem, modulo the restriction to purely
circular motions. Accordingly, the present tests may be viewed as
complementary to those from the previous section in the sense that the focus
now is on (33)–(34) for $G^{(2)}$ as it is the relevant contribution to the
DNO. Note that Eqs. (69)–(70) are completely equivalent to their Hamiltonian
counterparts (54)–(55) but were preferred for our simulations because they
directly yield the variable $s$ as an output. For the sake of comparison, the
same time step $\Delta t$ is used when solving (60)–(61) and (69)–(70). The
same choice of non-dimensionalization is adopted for both systems of
equations.

In the following illustrative experiments, we set $\gamma=1.4$ under
isentropic conditions [11; 40] and assume a uniform far-field pressure, say
$p_{\infty}=1$ in dimensionless units, for simplicity. Again, the parameter
$R_{m}=1000R$ is prescribed sufficiently large to avoid any possible
logarithmic singularity in (70). Figure 5 shows the relative error

${\rm Error}=\left|\frac{R+\eta(0,t)-s(t)}{s(t)}\right|\,,$ (71)

between $\eta$ from (60)–(61) and $s$ from (69)–(70), as a function of time
$t$. Due to angular invariance in this case, data on $\eta$ from (60)–(61) can
be taken at any arbitrary angle $\theta$, say $\theta=0$, when evaluating
(71). For both systems of equations, the same initial conditions are specified
at $t=0$, namely

$\big{(}\eta(\theta,0),\xi(\theta,0)\big{)}=(a_{3},0)\,,\quad\big{(}s(0),v(0)\big{)}=(R+a_{3},0)\,,$
(72)

with $\Delta t=0.01$ to integrate in time. Errors (71) are plotted in Fig.
5(a) for varying $a_{3}$ (with $M=8$, $N=256$) and in Fig. 5(b) for varying
$M$ (with $a_{3}=0.3$, $N=256$). Convergence with respect to surface amplitude
(here the initial amplitude $a_{3}$) or truncation order $M$ is confirmed by
these results, i.e. the lower $a_{3}$ or the larger $M$, the smaller the error
at any time $t$. This contrasts with the quick stagnation of errors when
increasing $M$, as reported in Sec. 6.1 for $\theta$-dependent surface
profiles. Along the lines of our previous comment on this point, the better
outcome here may be attributed to $G^{(1)}$ being inactive and the absence of
such Fourier multipliers as $C_{j}^{-|D|}$ from (34), thus making $G^{(2)}$
better conditioned for purely circular applications. Even in the large-
deflection regime $a_{3}=1$ (being on the order of $R$), the errors still
remain acceptable (below $100$%). Because Eqs. (60)–(61) are tested against
the RP system (69)–(70) which is itself solved numerically (and thus is also
affected by truncation or round-off errors, though to a lesser extent), error
fluctuations over time are expected as displayed in Fig. 5. Beside the low
levels of these instantaneous errors, their overall stable behavior over time
attests further to the accuracy and effectiveness of our boundary perturbation
approach. The computations were not found to be particularly sensitive to $N$
in this regular $\theta$-independent geometry, so convergence tests with
varying $N$ are not shown here for convenience.

Graphs of the bubble radii predicted by these two models during the time
evolution are compared in Fig. 6 for $a_{3}=0.3$ and $1.0$ (with $M=8$,
$N=256$). The competition between $p_{\infty}$ and $p_{B}$ induces a cycle of
compression-dilatation around the equilibrium state $r=R$, with surface
deflections of maximum amplitude $a_{3}$ (as determined by the initial
amplitude) occurring at the dilatation peaks. This cycle of compression-
dilatation is asymmetric in the sense that the surface amplitude (relative to
$r=R$) at the compression dips is slightly smaller than $a_{3}$. Moreover, the
dynamics during compression seems to be faster than that during dilatation, as
suggested by the steeper dips and smoother peaks along the curves of Fig. 6. A
small discrepancy between these two estimated radii is discernible around
$t=14$ (near a compression dip) for $a_{3}=1.0$, which is consistent with the
larger error as revealed by Fig. 5(a) in this case. The cyclic period is
graphically inferred to be $\tau=8.0$ and $9.8$ for $a_{3}=0.3$ and $1.0$
respectively, which is comparable to the linear fundamental period
$\tau_{0}=2\pi/\Omega_{0}=7.7$ of radial oscillations, with $\Omega_{0}$ given
by (58) for $\gamma=1.4$ and $p_{\infty}=1$ based on the RP equation (52). The
slight deviation is attributable to nonlinear effects, considering that $\tau$
is found to increase with $a_{3}$. If the far-field pressure $p_{\infty}$ was
amplified, the bubble dynamics would retain similar features but the cyclic
period would be shorter, in accordance with the dependence of (58) on
$p_{\infty}$. We point out that filtering was not used in any of these
simulations, for either (60)–(61) or (69)–(70).

The canonical conjugate variables $(\mathrm{q},\mathrm{p})$ for the RP
equation (52) can be reconstructed from the simulated variables $(s,v)$
according to (53), so that the associated Hamiltonian (56) can be evaluated.
The conservation of this Hamiltonian ${\cal H}$ is confirmed by Fig. 7 which
plots the time evolution of the error

${\rm Error}=\left|\frac{{\cal H}-{\cal H}_{0}}{{\cal H}_{0}}\right|\,,$

relative to the initial value ${\cal H}_{0}$, for $a_{3}=0.3$ and $1.0$.
Recall that this Hamiltonian structure is specific to the RP equation (52) and
is not a direct consequence of the Hamiltonian property (11) exhibited by the
full system. In particular, ${\cal H}$ remains a conserved quantity under the
action of $\Delta p$ (with uniform $p_{\infty}$) as given by (56), while the
original Hamiltonian formulation (11) with (23) does not extend to such a
forced regime.

For computations in this purely circular setting, the RP model (69)–(70) is
clearly a better option than (60)–(61), owing to its simpler and more explicit
form. Nevertheless, this discussion helps validate (60)–(61) as a viable
alternative to the RP equation, highlighting the role of $G_{j}^{(2)}$ which
can be readily incorporated into the DNO (35)–(36) for the general problem
(1)–(5), to produce a nonlinear model capable of describing a variety of
bubble shapes.

###  6.3 Simulation of nonlinear rotating waves

Further examining the range of applicability of this boundary perturbation
approach, we now turn our attention to the simulation of nonlinear waves
rotating steadily on the bubble surface. Unlike the previous case, these shape
distortions exhibit angular dependence according to (42). Equations (43)–(44)
constitute a boundary value problem with periodic boundary conditions in
$\Theta\in[0,2\pi)$. Via a pseudo-spectral method as outlined in Sec. 5.1,
$2N$ discrete equations arise for $2N$ unknowns $(\eta_{j},\xi_{j})$. They are
solved iteratively after specifying $c$ along with the initial guess (45)
based on second-order Stokes theory. Through this initial guess, the parameter
$n$ selects the wavenumber of the nonlinear wave. The initial amplitude
$a_{0}$ may be tuned depending on the choice of $c$ to help find the fixed-
point solution. This task is accomplished by using the Matlab routine fsolve,
which is essentially a quasi-Newton method with a numerical approximation of
the Jacobian matrix. We have successfully employed this Matlab routine in
previous work [23; 29; 56] to obtain numerical predictions from nonlinear
algebraic or differential equations. Given a value of $n>1$, we look here for
a nonlinear wave rotating at speed $c>c_{0}$. By analogy with the water wave
problem [17], an asymptotic behavior of the form

$c\simeq c_{0}\left(1+\frac{1}{2}\varepsilon^{2}\right)\simeq
c_{0}\left[1+\frac{1}{2}\left(\frac{a_{0}n}{R}\right)^{2}\right]\,,$

may be anticipated according to Stokes theory up to third order (recall that
the second-order contribution $c_{1}=0$ as shown in Sec. 4.2), which offers
the following estimate

$a_{0}\simeq\frac{R}{n}\sqrt{2\left(\frac{c}{c_{0}}-1\right)}\,,$

for $a_{0}$ in terms of $c$ and $n$, to be prescribed in the initial guess
(45).

Figure 8 depicts the bubble shapes produced by such rotating waves for
wavenumber $n=3$ and angular speeds $c=\\{1.64,1.70,1.75,1.86\\}$. In this
case, the linear phase speed is $c_{0}=1.63$. These nonlinear solutions were
computed using $N=256$ and $M=4$ or $6$ for (59) with (31)–(32). Indeed, for
wave speeds higher than $c=1.70$, it was necessary to specify a lower
truncation order of the DNO ($M=4$ rather than $M=6$) in order to ensure
convergence of the fixed-point iterative scheme. This issue is likely related
to ill-conditioning of the DNO computation as discussed in Sec. 6.1. For high
wave speeds, we also needed to include intermediate steps (with intermediate
wave speeds) as part of the iterative process going from the initial $c_{0}$
to the target $c$, where the initial guess at each step is given by the
converged solution from the previous step. This strategy improves the
convergence as compared to a single fsolve search from $c_{0}$ to $c$. We can
discern in Fig. 8 that the higher $c$, the larger the wave amplitude and
steepness (with steeper crests vs. smoother troughs). More quantitatively, the
amplitude of steadily rotating waves can be evaluated as
$a=|\eta_{\max}-\eta_{\min}|/2$. For $c=\\{1.64,1.70,1.75,1.86\\}$, we find
$a=\\{0.02,0.07,0.09,0.12\\}$ respectively, hence wave steepness
$\varepsilon=(a\,n)/R=\\{0.07,0.22,0.28,0.37\\}$ which is a measure of the
nonlinearity strength. These features are clearly revealed by Fig. 10(a) where
the bubble’s surface deformations are plotted versus $\Theta$ in a rectangular
format. These graphs also indicate that the wave steepening is accompanied by
a decrease of mean surface level $Q$, hence a decrease of bubble volume $V$.

Similar observations can be made when inspecting computations of mode-4
($n=4$) rotating waves as portrayed in Figs. 9 and 10(b). For
$c=\\{1.95,2.00,2.10,2.29\\}>c_{0}=1.93$ in this setting, we find
$a=\\{0.02,0.04,0.07,0.09\\}$ and $\varepsilon=\\{0.09,0.19,0.29,0.39\\}$
respectively. While the bubble under $n=3$ wave excitation tends to exhibit a
triangular-like shape as $c$ increases, the limiting geometry for $n=4$ is
more square-like. In either case, we were not able to achieve convergence for
wave speeds higher than those considered in Fig. 10.

We further check the accuracy of our fixed-point iterative algorithm for
(43)–(44) by prescribing its converged solutions as initial conditions for
(60)–(61) and by solving these evolution equations via the time integration
method described in Sec. 5.2. Under the same conditions as in Sec. 4.2, their
implementation is restricted to (31)–(32) for $G(\eta)$ and to (37) for
$\Delta p$. The steady character of these rotating waves is confirmed by Fig.
11 where snapshots of the bubble deformation are shown at $t=0$, $20$ in the
two highly nonlinear regimes $c=1.86$ ($n=3$) and $c=2.29$ ($n=4$). For both
solutions, the surface profiles at $t=0$ and $t=20$ look indistinguishable,
modulo a phase shift due to the uniform wave propagation. While these
snapshots happen to be located close together in either plot, it should be
kept in mind that the wave has actually rotated multiple times across the
computational domain over the span of the simulation, because of the periodic
boundary conditions. The associated numerical parameters are $\Delta t=0.001$,
$M=4$ and $N=256$. For such high wave speeds, filtering was required to run
the time-dependent computations as explained earlier (see Sec. 5.3).

A more quantitative assessment of these rotating wave solutions is provided in
Fig. 12 which plots the relative errors

${\rm
Error}=\left|\frac{H-H_{0}}{H_{0}}\right|\,,\quad\left|\frac{Q-Q_{0}}{Q_{0}}\right|\,,\quad\left|\frac{V-V_{0}}{V_{0}}\right|\,,$

on energy $H$, mean surface level $Q$ and volume $V$ versus time $t$ for
$c=1.86$ ($n=3$) and $c=2.29$ ($n=4$) as presented by Fig. 11. The quantities
$H_{0}$, $Q_{0}$, $V_{0}$ denote the corresponding initial values at $t=0$.
The trapezoidal rule was employed to evaluate integrals in the definitions
(7), (10), (23) of these invariants. Overall, the minuscule values of these
errors (in particular near machine precision for $Q$ and $V$) as well as their
stable evolution over time confirms that $H$, $Q$, $V$ are very well conserved
numerically and that steadily rotating waves are simulated accurately. The
computation of $H$ being more prone to inaccuracy than $Q$ or $V$ is likely
due to its more complicated expression (23) which directly involves the DNO
approximation. Moreover, the slightly larger errors for $c=2.29$ ($n=4$,
$\varepsilon=0.39$) as compared to $c=1.86$ ($n=3$, $\varepsilon=0.37$) are
expected considering that our boundary perturbation approach is better suited
for milder bubble deformations.

###  6.4 Simulation of nonlinear standing waves

Another class of interesting nonlinear solutions is that of time-periodic
standing waves. Such solutions have been much investigated in the three-
dimensional axisymmetric configuration, especially by means of analytical
perturbation calculations [52]. Their direct numerical computation is
particularly challenging in the nonlinear regime, which may be tackled by
solving the associated constrained problem via an adjoint method, and thus is
envisioned for future work. Alternatively, a standing wave can be produced by
superimposing two identical constant-amplitude waves that rotate in opposite
directions.

Exploiting this idea here, we numerically solve (60)–(61) in time together
with (37) from initial conditions given by

$\eta(\theta,0)=\eta_{1}(\theta)+\eta_{2}(\theta)\,,\quad\xi(\theta,0)=\xi_{1}(\theta)-\xi_{2}(\theta)\,,$

at $t=0$, where both pairs $(\eta_{1},\xi_{1})$ and $(\eta_{2},\xi_{2})$
denote steadily rotating waves with the same speed $c>c_{0}$ and wavenumber
$n>1$, as discussed in the previous section, but with possibly a phase shift
between them. The minus sign in $\xi(\theta,0)$ above is to ensure that these
two initial waves rotate in opposite directions. Their interaction will create
a standing (non-rotating) wave pattern that exhibits the same wavenumber $n$
as the original rotating waves but oscillates in time at a specific frequency
$\omega$. Because of this superposition, the resulting standing wave reaches
an amplitude equal to the sum of the initial amplitudes, which makes it a more
nonlinear solution than each individual rotating wave.

Figure 13 illustrates this process with two initial rotating waves
corresponding to $c=1.26$ and $n=2$. They are initially positioned in phase to
interfere with each other in such a way that oblate deformations associated
with the mode-2 standing wave occur along the axis $\theta=0$ while prolate
deformations occur along the axis $\theta=\pi/2$. The bubble shapes at various
instants are displayed in Fig. 13 during the first cycle of oblate-prolate
oscillation. This cycle indeed repeats itself over time in our simulations.
Note the substantial surface deflection relative to the unit circle at the
times of maximum oblate ($t=0.00$, $2.52$) and prolate ($t=1.26$, $3.78$)
oscillations.

The favorable results on conservation of $H$, $Q$, $V$ as demonstrated by Fig.
14 attest to the computation accuracy in this standing wave case as well. The
idea that we synthetically superimpose two rotating waves at $t=0$ (instead of
prescribing an actual initial condition for the standing wave) together with
the fact that the resulting solution is highly nonlinear as pointed out
earlier may explain why we see higher error levels here than say, in Fig. 12
for a single rotating wave.

Finally, we make an attempt to compare our numerical results in this setting
with previous data from the literature. Following Tsamopoulos and Brown [52],
we plot in Fig. 15 estimates of the frequency shift
$(\omega-\omega_{0})/\omega_{0}$ relative to $\omega_{0}$, the linear angular
frequency (40) with $n=2$, versus the aspect ratio $L/W$ of the bubble’s major
axis $L$ to its minor axis $W$ at maximum prolate deformation. We calculate
$L/W$ from our simulations by tracking the bubble radius $s$ whereas, for the
frequency $\omega$, we first estimate the period $T$ of a prolate oscillation
and then evaluate $\omega=2\pi/T$. In this two-dimensional context, the oblate
shape dynamics turn out to exactly mirror the prolate shape dynamics,
exhibiting an identical cycle of oscillation with the same period $T$ and
aspect ratio $L/W$. Figure 15 shows a comparison with laboratory measurements
by Trinh and Wang [51] for almost neutrally buoyant drops of silicon oil and
carbon tetrachloride suspended in distilled water. Data sets for two different
drop radii $R=0.49$ cm and $0.62$ cm are reported. Numerical estimates by
Foote [18] for drops as well as asymptotic predictions by Tsamopoulos and
Brown [52] for bubbles are also included in Fig. 15. All these data were
extracted from Fig. 4 in [52]. Note that this assessment is only meant to be
qualitative because, on one hand, we compare two-dimensional simulations with
three-dimensional results for which quantitative differences are expected. In
particular, the numerical or asymptotic results of Foote [18] and Tsamopoulos
and Brown [52] assume a three-dimensional spherical geometry that is invariant
under azimuthal rotation. Secondly, our numerical estimates are based on
inviscid potential-flow theory while laboratory measurements are subject to
viscous effects. Furthermore, drops and bubbles may reveal different behaviors
depending on the situation, and it should also be kept in mind that our
computations were not set up under the same conditions as in the laboratory
experiments (regarding e.g. the wave excitation mechanism).

Our own data in Fig. 15 correspond to $c=\\{1.23,1.24,1.25,1.26,1.27\\}$
representing standing waves of amplitude $a=\\{0.08,0.15,0.19,0.22,0.25\\}$ at
maximum deformation. The respective maximum wave steepnesses are
$\varepsilon=\\{0.17,0.30,0.38,0.45,0.50\\}$ for $n=2$, which is indicative of
the highly nonlinear character of these computed solutions. To this aim, we
specified $M=4$, $N=256$ with a time step as small as $\Delta t=0.0002$ and
$0.00001$ for $c=1.26$ and $1.27$ respectively. Despite this fine resolution
in time, filtering was still required to stabilize the simulation of such
large-amplitude waves during their oscillatory evolution.

It is first remarked that our $L/W$ estimates fall within the range of values
observed by these previous authors. In our situation, the higher $c$, the
larger $a$ (or $\varepsilon$) and so is the aspect ratio $L/W$ of standing
waves. We were not able to obtain two-dimensional solutions for $L/W$ values
as high as those reported in three-dimensional configurations. However, a
striking difference from these previous studies is that we find a (positive)
frequency upshift $(\omega-\omega_{0})/\omega_{0}$ relative to $\omega_{0}$,
rather than a (negative) frequency downshift. Note that the negative ordinate
for our last data point (the one farthest to the right associated with
$c=1.27$) is likely a numerical artifact. Although our simulation managed to
be completed in this highly nonlinear case (by means of filtering together
with a small time step), producing an acceptable solution overall, it was
particularly subject to errors to the point that spurious noise was detectable
near the wave crests and troughs.

The qualitative agreement between the asymptotic predictions of Tsamopoulos
and Brown [52] and the other data sets suggests that the discrepancy with our
numerical results may be attributed to the two different geometries. An
argument in favor of this explanation stems from the different expression of
the linear angular frequency

$\omega_{0}^{2}=\frac{\sigma}{\rho R^{3}}(n+2)(n^{2}-1)\,,\quad n>1\,,$ (73)

for three-dimensional axisymmetric bubbles [52]. Clearly, $\omega_{0}$ from
(73) is larger than that from (40) for any $n>1$, so the frequency shift
$(\omega-\omega_{0})/\omega_{0}$ may be more prone to being negative in three
dimensions than in two dimensions. Another argument is related to the fact
that three-dimensional drops or bubbles undergoing mode-2 oscillations have
been observed to spend a longer part of each period in a prolate form than in
an oblate one, which may be due to the vertical direction ($\theta=\pi/2$)
being an axis of symmetry in the axisymmetric case [52]. By contrast, our two-
dimensional problem is isotropic, giving rise to identical oblate and prolate
shape dynamics as mentioned earlier. Accordingly, the typical period $T$ of a
nonlinear mode-2 standing wave in three dimensions may be longer than the
linear period, or equivalently the associated frequency $\omega$ may be
shorter than $\omega_{0}$, hence a negative frequency shift
$(\omega-\omega_{0})/\omega_{0}$. Finally, perhaps a simpler (more intuitive)
argument is that, because a bubble volume is less in two dimensions than in
three dimensions for a given radius $R$, it thus takes less effort to excite
the bubble surface at a certain mode $n$. As a consequence, the excitation
frequency $\omega$ should be larger in two dimensions than in three
dimensions, which promotes a frequency upshift in the former setting versus a
frequency downshift in the latter setting. Overall, this geometrical
explanation for the different dynamics depending on the spatial dimension is
consistent with e.g. recent laboratory measurements by Duplat [12] who noticed
that, during the late stages of a circular bubble’s collapse, its radius
vanishes like $s(t)\sim(t_{0}-t)^{1/2}$ as $t\to t_{0}$, which contrasts with
the faster collapse like $s(t)\sim(t_{0}-t)^{2/5}$ for a spherical bubble.

##  7 Conclusions

We consider the two-dimensional problem of free or forced deformations of a
gas bubble immersed in a liquid of infinite extent. This configuration has
drawn much less attention than the three-dimensional axisymmetric geometry,
especially in applications to bubble acoustics. Based on nonlinear theory for
potential flow in the presence of a moving boundary under surface tension
without gravity, and adopting the Hamiltonian formulation proposed by Benjamin
[1], we introduce the DNO to clarify the reduction in terms of surface
variables alone and to verify the canonical symplectic structure. Benjamin [1]
did not specifically treat this two-dimensional case in polar coordinates.

We propose a Taylor series representation of the DNO about a quiescent
circular state of the bubble. A recursion formula is devised to evaluate this
Taylor series up to an arbitrary order of nonlinearity, with each term having
two distinct components to describe shape distortions and radial pulsations as
well as their coupling. Each of these components is given by a sum of
concatenations of Fourier multipliers with powers of the surface displacement.

In this theoretical framework, various analytical approximations of the
problem are examined. We obtain a Stokes wave solution that is accurate up to
second order in wave steepness for steadily rotating shape oscillations
without volume change. Assuming circular symmetry, we also derive an inviscid
RP model for the time evolution of the bubble radius under the excitation of a
far-field pressure. Its independent Hamiltonian structure is established in
the case of uniform pressure forcing.

Exploiting the lower-dimensional form of the full governing equations together
with the Taylor expansion of the DNO, we develop a direct numerical solver
where each term in the DNO series is computed efficiently and accurately by a
pseudo-spectral method with the FFT. Via this Fourier decomposition, an
arbitrary number of shape modes can be specified to contribute to nonlinear
wave interactions, and the linear terms can be solved exactly by the
integrating factor technique. This approach is thus well suited for time-
dependent simulations in a forced or freely evolving regime.

We show extensive numerical tests on the convergence of the DNO as a function
of the truncation order $M$ for varying bubble shapes, surface steepnesses and
spatial resolutions. The significance of aliasing errors when computing
nonlinear products is revealed and a zero-padding technique is implemented to
remove them. While the convergence errors remain small overall, they are not
found to decrease exponentially with increasing $M$, in contrast to what might
be expected from the analyticity property of the DNO. Rather, the convergence
quickly stagnates past $M\simeq 2$ and even deteriorates past a critical value
of $M$ for severe shape distortions or fine spatial resolutions. This behavior
is consistent with the ill-conditioning of the DNO series, which has been
documented in other contexts [30; 46; 55].

Despite these numerical issues associated with the DNO, computations of bubble
deformations of moderately large amplitude or steepness can be performed with
satisfactory accuracy. We apply this algorithm to simulating cycles of
compression-dilatation for a purely circular bubble driven by a uniform
pressure field. We test these results against predictions by the RP model,
thus providing mutual validation for these two proposed approaches. Nonlinear
shape oscillations of a non-circular bubble in the unforced regime without
volume change are also investigated. We compute steadily rotating waves and
time-periodic standing waves for a few first modes of the bubble surface.
Focusing on mode-2 standing waves, we examine the shift in oscillation
frequency (relative to the linear value) as a function of the bubble’s maximum
aspect ratio, and we compare our numerical estimates to laboratory
measurements and other theoretical predictions. Our two-dimensional results
show a frequency upshift that accentuates with increasing aspect ratio, while
previous studies have reported a frequency downshift. This notable discrepancy
is attributed to the difference in geometry and highlights peculiarities of
the two-dimensional problem which may not be widely known, as opposed to the
three-dimensional case.

Possible future work includes a detailed phase plane analysis of solutions to
the new RP equation, exploiting its Hamiltonian structure and examining the
role of the external pressure field. In addition to mathematical results, the
present study explores the possibility of using a boundary perturbation method
(as represented by the Taylor expansion of the DNO) to simulate the nonlinear
surface dynamics of a gas bubble. The focus here is on the algorithmic
development and preliminary assessment of this numerical solver, in cases
where the processes of radial pulsations and shape distortions do not
interact. On the other hand, the nonlinear coupling between these different
modes under resonant or near-resonant conditions (e.g. under tuned harmonic
forcing) constitutes an interesting problem in its own right, leading to
complex phenomena [22; 43]. Its in-depth investigation via this numerical
scheme is envisioned for a subsequent paper.

Finally, we emphasize that the proposed mathematical formulation and
computational strategy are not restricted to the two-dimensional setting.
Provided a FFT version in spherical harmonics is available [44], we plan to
extend the present results to the more general problem on nonlinear
deformations of three-dimensional bubbles. Related work in the context of
three-dimensional water waves can be found in [10; 55].

## Acknowledgements

P. Guyenne is partially supported by the NSF under grant No. DMS-2307712.

## References

  * [1] T.B. Benjamin, Hamiltonian theory for motions of bubbles in an infinite liquid, J. Fluid Mech. 181 (1987) 349–379.
  * [2] T.B. Benjamin, A.T. Ellis, Self-propulsion of asymmetrically vibrating bubbles, J. Fluid Mech. 212 (1990) 65–80.
  * [3] M.P. Brenner, S. Hilgenfeldt, D. Lohse, Single-bubble sonoluminescence, Rev. Mod. Phys. 74 (2002) 425–484.
  * [4] A. Bu$\beta$mann, F. Riahi, B. Gökce, S. Adami, S. Barcikowski, N.A. Adams, Investigation of cavitation bubble dynamics near a solid wall by high-resolution numerical simulation, Phys. Fluids 35 (2023) 016115.
  * [5] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag Berlin, Heidelberg, 1988.
  * [6] R. Coifman, Y. Meyer, Nonlinear harmonic analysis and analytic dependence, Proc. Symp. Pure Math. 43 (1985) 71–78.
  * [7] W. Craig, P. Guyenne, C. Sulem, The water wave problem and Hamiltonian transformation theory, in Waves in Flows, eds. T. Bodnár, G. Galdi & S̆. Nec̆asová, Adv. Math. Fluid Mech., pp. 113–196, Birkhäuser, 2021.
  * [8] W. Craig, C. Sulem, Numerical simulation of gravity waves, J. Comput. Phys. 108 (1993) 73–83.
  * [9] D. Crowdy, Circulation-induced shape deformations of drops and bubbles: Exact two-dimensional models, Phys. Fluids 11 (1999) 2836–2845.
  * [10] R. de la Llave, P. Panayotaros, Gravity waves on the surface of the sphere, J. Nonlinear Sci. 6 (1996) 147–167.
  * [11] A.A. Doinikov, Translational motion of a bubble undergoing shape oscillations, J. Fluid Mech. 501 (2004) 1–24.
  * [12] J. Duplat, Dynamics of expansion and collapse of explosive two-dimensional bubbles, J. Fluid Mech. 859 (2019) 677–690.
  * [13] S.A. Dyachenko, Traveling capillary waves on the boundary of a fluid disc, Stud. Appl. Math. 148 (2021) 125–140.
  * [14] A. Esmaeeli, G. Tryggvason, Direct numerical simulations of bubbly flows. Part 2. Moderate Reynolds number arrays, J. Fluid Mech. 385 (1999) 325–358.
  * [15] Q. Fang, D.P. Nicholls, J. Shen, A stable, high-order method for three-dimensional, bounded-obstacle, acoustic scattering, J. Comput. Phys. 224 (2007) 1145–1169.
  * [16] Z.C. Feng, L.G. Leal, Nonlinear bubble dynamics, Annu. Rev. Fluid Mech. 29 (1997) 201–243.
  * [17] J.D. Fenton, A fifth-order Stokes theory for steady waves, J. Waterw. Port Coast. Ocean Eng. 111 (1985) 216–234.
  * [18] G.B. Foote, A numerical method for studying simple drop behavior: simple oscillation, J. Comput. Phys. 11 (1973) 507–530.
  * [19] D. Fuster, S. Popinet, An all-Mach method for the simulation of bubble dynamics problems in the presence of surface tension, J. Comput. Phys. 374 (2018) 752–768.
  * [20] F. Garoosi, T. Merabtene, T.-F. Mahdi, Numerical simulation of merging of two rising bubbles with different densities and diameters using an enhanced Volume-Of-Fluid (VOF) model, Ocean Eng. 247 (2022) 110711.
  * [21] A. Guckenberger, S. Gekle, A boundary integral method with volume-changing objects for ultrasound-triggered margination of microbubbles, J. Fluid Mech. 836 (2018) 952–997.
  * [22] M. Guédra, C. Inserra, Bubble shape oscillations of finite amplitude, J. Fluid Mech. 857 (2018) 681–703.
  * [23] P. Guyenne, Large-amplitude internal solitary waves in a two-fluid model, C. R. Mecanique 334 (2006) 341–346.
  * [24] P. Guyenne, A high-order spectral method for nonlinear water waves in the presence of a linear shear current, Comput. Fluids 154 (2017) 224–235.
  * [25] P. Guyenne, HOS simulations of nonlinear water waves in complex media, in Nonlinear Water Waves–An Interdisciplinary Interface, eds. D. Henry et al., Tutorials, Schools and Workshops in the Mathematical Sciences, pp. 53–69, Birkhäuser, 2019.
  * [26] P. Guyenne, S.T. Grilli, Numerical study of three-dimensional overturning waves in shallow water, J. Fluid Mech. 547 (2006) 361–388.
  * [27] P. Guyenne, D.P. Nicholls, A high-order spectral method for nonlinear water waves over moving bottom topography, SIAM J. Sci. Comput. 30 (2007) 81–101.
  * [28] P. Guyenne, E.I. Părău, Computations of fully nonlinear hydroelastic solitary waves on deep water, J. Fluid Mech. 713 (2012) 307–329.
  * [29] P. Guyenne, E.I. Părău, Finite-depth effects on solitary waves in a floating ice sheet, J. Fluids Struct. 49 (2014) 242–262.
  * [30] P. Guyenne, E.I. Părău, An operator expansion method for computing nonlinear surface waves on a ferrofluid jet, J. Comput. Phys. 321 (2016) 414–434.
  * [31] A.A. Harkin, T.J. Kaper, A. Nadim, Energy transfer between the shape and volume modes of a nonspherical gas bubble, Phys. Fluids 25 (2013) 062101.
  * [32] X. He, B. Wang, J. Meng, S. Zhang, S. Wang, How to prevent bubbles in microfluidic channels, Langmuir 37 (2021) 2187–2194.
  * [33] T.Y. Hou, R. Li, Computing nearly singular solutions using pseudo-spectral methods, J. Comput. Phys. 226 (2007) 379–397.
  * [34] T.Y. Hou, J.S. Lowengrub, M.J. Shelley, Removing the stiffness from interfacial flows with surface tension, J. Comput. Phys. 114 (1994) 312–338.
  * [35] J. Hua, J. Lou, Numerical simulation of bubble rising in viscous liquid, J. Comput. Phys. 222 (2007) 769–795.
  * [36] O.A. Hurricane et al., Fuel gain exceeding unity in an inertially confined fusion implosion, Nature 506 (2014) 343–348.
  * [37] K. Khoshmanesh, A. Almansouri, H. Albloushi, P. Yi, R. Soffe, K. Kalantar-Zadeh, A multi-functional bubble-based microfluidic system, Sci. Rep. 5 (2015) 9942.
  * [38] S. Li, Y. Saade, D. van der Meer, D. Lohse, Comparison of Boundary Integral and Volume-of-Fluid methods for compressible bubble dynamics, Int. J. Multiph. Flow 145 (2021) 103834.
  * [39] J.R. Lindner, Microbubbles in medical imaging: current applications and future directions, Nat. Rev. Drug Discov. 3 (2004) 527–533.
  * [40] M.S. Longuet-Higgins, Monopole emission of sound by asymmetric bubble oscillations. Part 1. Normal modes, J. Fluid Mech. 201 (1989) 525–541.
  * [41] J.T.S. Ma, P.K.C. Wang, Effect of initial air content on the dynamics of bubbles in liquids, IBM J. Res. Dev. 6 (1962) 472–474.
  * [42] N.K. McDougald, L.G. Leal, Numerical study of the oscillations of a non-spherical bubble in an inviscid, incompressible liquid. Part I: free oscillations from non-equilibrium initial conditions, Int. J. Multiph. Flow 25 (1999) 887–919.
  * [43] N.K. McDougald, L.G. Leal, Numerical study of the oscillations of a non-spherical bubble in an inviscid, incompressible liquid. Part II: the response to an impulsive decrease in pressure, Int. J. Multiph. Flow 25 (1999) 921–941.
  * [44] M.J. Mohlenkamp, A fast transform for spherical harmonics, J. Fourier Anal. Appl. 5 (1999) 159–184.
  * [45] D.P. Nicholls, N. Nigam, Exact non-reflecting boundary conditions on general domains, J. Comput. Phys. 194 (2004) 278–303.
  * [46] D.P. Nicholls, F. Reitich, Stability of high-order perturbative methods for the computation of Dirichlet–Neumann operators, J. Comput. Phys. 170 (2001) 276–298.
  * [47] M.S. Plesset, A. Prosperetti, Bubble dynamics and cavitation, Annu. Rev. Fluid Mech. 9 (1977) 145–185.
  * [48] C. Pozrikidis, Three-dimensional oscillations of rising bubbles, Eng. Anal. Bound. Elem. 28 (2004) 315–323.
  * [49] A. Prosperetti, Vapor bubbles, Annu. Rev. Fluid Mech. 49 (2017) 221–248.
  * [50] R.P. Tong, A new approach to modelling an unsteady free surface in boundary integral methods with application to bubble-structure interactions, Math. Comput. Simul. 44 (1997) 415–426.
  * [51] E. Trinh, T.G. Wang, Large amplitude drop shape oscillations, in Proc. 2nd Int. Colloq. on Drops and Bubbles, JPL Publ., pp. 82–98, Pasadena, 1982.
  * [52] J.A. Tsamopoulos, R.A. Brown, Nonlinear oscillations of inviscid drops and bubbles, J. Fluid Mech. 127 (1983) 519–537.
  * [53] Q.X. Wang, J.R. Blake, Non-spherical bubble dynamics in a compressible liquid. Part 1. Travelling acoustic wave, J. Fluid Mech. 659 (2010) 191–224.
  * [54] R. Wegmann, D. Crowdy, Shapes of two-dimensional bubbles deformed by circulation, Nonlinearity 13 (2000) 2131–2141.
  * [55] L. Xu, P. Guyenne, Numerical simulation of three-dimensional nonlinear water waves, J. Comput. Phys. 228 (2009) 8446–8466.
  * [56] B. Xu, P. Guyenne, Assessment of a porous viscoelastic model for wave attenuation in ice-covered seas, Appl. Ocean Res. 122 (2022) 103122.
  * [57] V.E. Zakharov, Stability of periodic waves of finite amplitude on the surface of a deep fluid, J. Appl. Mech. Tech. Phys. 9 (1968) 190–194.
  * [58] A.M. Zhang, Y.L. Liu, Improved three-dimensional bubble dynamics model based on boundary element method, J. Comput. Phys. 294 (2015) 208–223.

(a)  (b)

Figure 1: (a) Relative $L^{2}$ error on the DNO vs. truncation order $M$ for
$\eta$ given by (66). (b) Surface profile $s$ in polar coordinates
$(r,\theta)$. Angles are indicated in degrees. In each panel, graphs are shown
for varying amplitude $a_{1}$ with $n=2$ and $N=256$.

(a)  (b)

Figure 2: Relative $L^{2}$ error on the DNO vs. truncation order $M$ for
$\eta$ given by (66). Graphs are shown for (a) varying wavenumber $n$ (with
$a_{1}=0.1$, $N=256$) and (b) varying resolution $N$ (with $a_{1}=0.1$,
$n=2$).

(a)  (b)

Figure 3: Relative $L^{2}$ error on the DNO vs. truncation order $M$ for
$\eta$ given by (67). Graphs are shown for (a) varying amplitude $a_{2}$ (with
$n=2$, $N=256$) and (b) varying resolution $N$ (with $a_{2}=5\times 10^{-5}$,
$n=2$).

(a)  (b)

Figure 4: Relative $L^{2}$ error on the DNO vs. truncation order $M$ for
$\eta$ given by (a) profile (66) (with $a_{1}=0.1$, $n=10$) and (b) profile
(67) (with $a_{2}=5\times 10^{-5}$, $n=2$). In each panel, the aliased and de-
aliased results are compared. In all these cases, the resolution is $N=256$.

(a)  (b)

Figure 5: Relative error on bubble radius $s$ vs. time $t$ for far-field
pressure $p_{\infty}=1$ and initial conditions given by (72). Graphs are shown
for (a) varying initial amplitude $a_{3}$ (with $M=8$, $N=256$) and (b)
varying truncation order $M$ (with $a_{3}=0.3$, $N=256$).

(a) $a_{3}=0.3$ (b) $a_{3}=1$

Figure 6: Bubble radius $s$ vs. time $t$ for far-field pressure $p_{\infty}=1$
and initial conditions given by (72). Graphs are shown for (a) $a_{3}=0.3$ and
(b) $a_{3}=1.0$ (with $M=8$, $N=256$). For each $a_{3}$, predictions from
(60)–(61) (blue circles) and (69)–(70) (red line) are compared. Figure 7:
Relative error on ${\cal H}$ vs. time $t$ for far-field pressure
$p_{\infty}=1$ and initial conditions given by (72). Predictions from the RP
model (69)–(70) are shown for $a_{3}=0.3$ (blue curve) and $a_{3}=1$ (red
curve) with $\Delta t=0.01$.

(a) $c=1.64$ (b) $c=1.70$ (c) $c=1.75$ (d) $c=1.86$

Figure 8: Bubble shape $s$ in polar coordinates $(r,\Theta)$ for steadily
rotating waves with $n=3$ and $N=256$. The different panels correspond to
angular speeds (a) $c=1.64$, (b) $c=1.70$, (c) $c=1.75$, (d) $c=1.86$. As a
reference, the red circle represents the unperturbed bubble of radius $R=1$.
Angles are indicated in degrees.

(a) $c=1.95$ (b) $c=2.00$ (c) $c=2.10$ (d) $c=2.29$

Figure 9: Bubble shape $s$ in polar coordinates $(r,\Theta)$ for steadily
rotating waves with $n=4$ and $N=256$. The different panels correspond to
angular speeds (a) $c=1.95$, (b) $c=2.00$, (c) $c=2.10$, (d) $c=2.29$. As a
reference, the red circle represents the unperturbed bubble of radius $R=1$.
Angles are indicated in degrees.

(a) $n=3$ (b) $n=4$

Figure 10: Bubble deformation $s$ vs. angle $\Theta$ for steadily rotating
waves with (a) $n=3$ and (b) $n=4$. In each panel, graphs are shown for
varying speed $c$ with $N=256$.

(a) $n=3$ (b) $n=4$

Figure 11: Bubble deformation $s$ vs. angle $\theta$ at $t=0$ (red line) and
$t=20$ (blue line) during the propagation of steadily rotating waves for (a)
$c=1.86$, $n=3$ and (b) $c=2.29$, $n=4$ (with $M=4$, $N=256$).

(a) $n=3$ (b) $n=4$

Figure 12: Relative errors on $H$, $Q$, $V$ vs. time $t$ during the
propagation of steadily rotating waves for (a) $c=1.86$, $n=3$ and (b)
$c=2.29$, $n=4$ (with $M=4$, $N=256$).

(a) $t=0.00$ (b) $t=0.34$ (c) $t=0.50$ (d) $t=0.62$ (e) $t=0.76$ (f) $t=0.94$
(g) $t=1.26$ (h) $t=2.52$ (i) $t=3.78$

Figure 13: Bubble shape $s$ in polar coordinates $(r,\theta)$ for a time-
periodic standing wave due to the superposition of two steadily rotating waves
with $c=1.26$ and $n=2$ ($M=4$, $N=256$). The different panels correspond to
times (a) $t=0.00$, (b) $t=0.34$, (c) $t=0.50$, (d) $t=0.62$, (e) $t=0.76$,
(f) $t=0.94$, (g) $t=1.26$, (h) $t=2.52$, (i) $t=3.78$. As a reference, the
red circle represents the unperturbed bubble of radius $R=1$. Angles are
indicated in degrees. (a)  Figure 14: Relative errors on $H$, $Q$, $V$ vs.
time $t$ during the evolution of a time-periodic standing wave due to the
superposition of two steadily rotating waves with $c=1.26$ and $n=2$ ($M=4$,
$N=256$). (a)  Figure 15: Frequency shift $(\omega-\omega_{0})/\omega_{0}$ vs.
aspect ratio $L/W$ at maximum deformation for standing waves corresponding to
$n=2$ and $c=\\{1.23,1.24,1.25,1.26,1.27\\}$ ($M=4$, $N=256$). Our numerical
results are compared to experimental data for drops by Trinh and Wang [51],
numerical estimates for drops by Foote [18] and asymptotic predictions for
bubbles by Tsamopoulos and Brown [52].