# The Kirkwood-Bethe hypothesis for bubble dynamics, cavitation and underwater
explosions

Fabian Denner<EMAIL_ADDRESS>Department of Mechanical Engineering,
Polytechnique Montréal, Montréal, H3T 1J4, Québec, Canada

###### Abstract

Pressure-driven bubble dynamics is a major topic of current research in fluid
dynamics, driven by innovative medical therapies, sonochemistry, material
treatments, and geophysical exploration. First proposed in 1942, the Kirkwood-
Bethe hypothesis provides a simple means to close the equations that govern
pressure-driven bubble dynamics as well as the resulting flow field and
acoustic emissions in spherical symmetry. The models derived from the
Kirkwood-Bethe hypothesis can be solved using standard numerical integration
methods at a fraction of the computational cost required for fully resolved
simulations. Here, the theoretical foundation of the Kirkwood-Bethe hypothesis
and contemporary models derived from it are gathered and reviewed, as well as
generalized to account for spherically symmetric, cylindrically symmetric, and
planar one-dimensional domains. In addition, the underpinning assumptions are
clarified and new results that scrutinize the predictive capabilities of the
Kirkwood-Bethe hypothesis with respect to the complex acoustic impedance
experienced by curved acoustic waves and the formation of shock waves are
presented. Although the Kirkwood-Bethe hypothesis is built upon simplifying
assumptions and lacks some basic acoustic properties, models derived from it
are able to provide accurate predictions under the specific conditions
associated with pressure-driven bubble dynamics, cavitation and underwater
explosions.

Copyright (2024) Author(s). This article is distributed under a Creative
Commons Attribution (CC BY) License.

##  I Introduction

Pressure-driven bubble dynamics are encountered and play an important role in
countless engineering applications and natural phenomena (Plesset and
Prosperetti, 1977; Brennen, 1995; Lauterborn and Kurz, 2010). The formation
and collapse of vapor bubbles (cf. hydrodynamic cavitation) are well known to
cause material erosion of propellers and hydro turbines (Blake and Gibson,
1987; Reuter, Deiter, and Ohl, 2022; Dular and Ohl, 2023), and the related
acoustic emissions are actively studied and monitored in situ (Johansen  _et
al._ , 2017; Gaisser, Kirschner, and Riedelbauch, 2023; Haskell  _et al._ ,
2023). However, the past decades have seen a paradigm shift in how pressure-
driven bubble dynamics and, especially, cavitation are studied, away from
treating cavitation solely as a hazard, towards leveraging its mechanical and
thermal effects. This has spawned innovative new applications, such as
contrast-enhanced medical imaging (Wan, Feng, and ter Haar, 2015),
intensifying chemical reactions (Meroni  _et al._ , 2022), breaking down
kidney and urinary stones (Maeda  _et al._ , 2018; Sapozhnikov, Maxwell, and
Bailey, 2021), treating tumors (Yeats  _et al._ , 2023), using cavitation as a
sonic scalpel for precision surgery (Vogel  _et al._ , 1986; Lee  _et al._ ,
2017), targeted drug delivery by ultrasound-driven micro- and nanobubbles
(Shakya  _et al._ , 2024), 3D printing with sound (Habibi  _et al._ , 2022),
actuating smart materials (Athanassiadis  _et al._ , 2022), and enhancing
bioadhesion (Ma  _et al._ , 2022). The energy focused by a strong bubble
collapse can be utilized to pierce biological cells (Ohl  _et al._ , 2006;
Helfield  _et al._ , 2016), to destruct bacteria for wastewater treatment
(Šarc  _et al._ , 2018), to clean surfaces and membranes (Chahine  _et al._ ,
2016; Reuter  _et al._ , 2017), to synthesize nanoparticles and nanostructured
materials (Barcikowski  _et al._ , 2019), and to manipulate hardened
materials, such as metals (Soyama and Korsunsky, 2022) and glass (Gutiérrez-
Hernández  _et al._ , 2023). Ferns use cavitation-driven catapults to eject
their spores (Llorens  _et al._ , 2016) and pistol shrimps use the shock waves
generated by collapsing bubbles to stun their prey (Koukouvinis, Bruecker, and
Gavaises, 2017). The pressure and shock waves of underwater explosions are
studied for defense purposes as a major component of marine warfare (Cole,
1948; Wu  _et al._ , 2020; Yu, Song, and Choi, 2023) and the pressure
transients produced by the generation and collapse of underwater bubbles are
used for seismic surveys in marine and geophysical exploration (de Graaf,
Penesis, and Brandner, 2014; MacGillivray, 2019; Athanassiadis, 2019). The
advent of modern monitoring and control techniques for bubble dynamics, many
of which now leverage machine learning algorithms, largely rely on acoustic
signals to identify and classify bubble activity (O’Reilly and Hynynen, 2012;
Novell  _et al._ , 2020; Doan, Huynh-The, and Kim, 2022; Rom, 2023; Gaisser,
Kirschner, and Riedelbauch, 2023; Klapcsik  _et al._ , 2023; Mallik, Jaiman,
and Jelovica, 2024). This, in turn, requires high-fidelity models of the
complex bubble dynamics and nonlinear acoustics to classify the related
acoustic signatures and to train machine learning algorithms.

The Kirkwood-Bethe hypothesis (Kirkwood and Bethe, 1942) is an important tool
for fundamental and theoretical research on pressure-driven bubble dynamics,
as well as the related fluid dynamics and acoustics. It stipulates that in
spherical symmetry the quantity

$g=r\left(h-h_{\infty}+\frac{u^{2}}{2}\right),$ (1)

emitted at the gas-liquid interface of a gas cavity or bubble, is constant
along outgoing characteristics and propagates with speed $c+u$, where $r$ is
the radial coordinate, $c$ is the liquid speed of sound, $u$ is the flow
velocity, and $h-h_{\infty}$ is the difference in specific enthalpy between
radial location $r$ and the ambient state denoted with subscript $\infty$. The
power of the Kirkwood-Bethe hypothesis lies in the fact that it enables, in
conjunction with a suitable equation of state, to close the set of equations
governing the pressure-driven behavior of gas bubbles and the resulting flow
of the surrounding liquid, such that mathematical models for the velocity and
pressure generated by a cavitation bubble or an underwater explosion can be
derived. Computer programs developed in the 1960s using models for pressure-
driven bubble dynamics and acoustics based on the Kirkwood-Bethe hypothesis
were also among the first open-source software tools in fluid dynamics (Ivany,
1965; Lilliston, 1966). Despite the widespread availability of simulation
tools that can predict two- and three-dimensional pressure-driven bubble
dynamics in complex scenarios with increasing accuracy and robustness (Denner,
Xiao, and van Wachem, 2018; Fuster and Popinet, 2018; Schmidmayer  _et al._ ,
2020; Bryngelson  _et al._ , 2021; Saade, Lohse, and Fuster, 2023), the
Kirkwood-Bethe hypothesis and the models derived from it, henceforth
collectively referred to as KBH models, still play an important role for the
prediction of pressure-driven bubble phenomena. Predicting the dynamic
behavior of collapsing bubbles and the pressure transients they produce
remains challenging and KBH models have emerged as robust benchmarks for
testing and validating these numerical methods (Koch  _et al._ , 2016; Denner,
Evrard, and van Wachem, 2020; Gonzalez-Avila, Denner, and Ohl, 2021; Saade,
Lohse, and Fuster, 2023; Wang  _et al._ , 2024a). In addition, acoustic
emissions can be modulated by different mechanisms (Crighton  _et al._ , 1992;
Christov and Christov, 2017; Denner, 2024) and, consequently, interpreting
acoustic signals at a distance remains a complex task. To this end, the
Kirkwood-Bethe hypothesis has seen a renewed interest over the past years as
the basis for gaining an improved understanding of the interaction between
acoustic emissions and fluid dynamics (Lai  _et al._ , 2022; Liang  _et al._ ,
2022; Denner and Schenke, 2023a; Wen  _et al._ , 2023; Denner, 2024),
primarily in relation to pressure-driven bubble dynamics.

Although the Kirkwood-Bethe hypothesis and models are well established and
widely used today, much of their theoretical foundation is scattered over many
institutional reports (Gilmore, 1952; Flynn, 1957; Hickling and Plesset, 1963;
Lilliston, 1966; Beyer, 1974; Flynn, 1978), declassified military reports
(Rice and Ginell, 1943, 1944; Whitham, 1953; Snay, 1966; Naugol’nykh and Roy,
1971), textbooks (Cole, 1948; Cooper, 1958; Rozenberg, 1971; Lauterborn,
1980), PhD theses (Schneider, 1949; Ivany, 1965), and papers (Hickling and
Plesset, 1964; Ivany and Hammitt, 1965; Mellen, 1956; Akulichev  _et al._ ,
1968; Kedrinskii, 1972; Ebeling, 1978; Vokurka, 1986) published in the years
and decades following the report of Kirkwood and Bethe (1942). Many of these
formidable works are difficult to obtain and are, thus, often not credited
correctly or even overlooked in the modern literature. Further adding to the
difficulties associated with gaining a comprehensive picture of the Kirkwood-
Bethe hypothesis and models, the original report of Kirkwood and Bethe
(Kirkwood and Bethe, 1942) does not appear to be indexed in scientific
databases, meaning that a list of subsequent publications referencing this
report is not readily available, and no review of the scientific literature
related to the Kirkwood-Bethe hypothesis has been published.

The aim of this review is to gather, for the first time, the theoretical
foundation of the Kirkwood-Bethe hypothesis, as well as models derived from it
that are still relevant today, in a single document. Moreover, the description
of the flow field arising from the Kirkwood-Bethe hypothesis is further
generalized and the analysis of the limitations and assumptions underpinning
the Kirkwood-Bethe hypothesis is extended.

As a first step, Section II provides a short historical overview of major
developments and applications related to the Kirkwood-Bethe hypothesis.
Subsequently, the Kirkwood-Bethe hypothesis is derived from the conservation
laws governing mass, momentum and energy in Section III, paying particular
attention to the underpinning assumptions. The application of different
equations of state describing the properties of the liquid is presented in
Section IV, the equation of motion of the gas-liquid interface of a pressure-
driven gas bubble or cavity is presented in Section V, and models to describe
the flow generated by oscillating or collapsing bubbles and cavities are
derived in Section VI. The model hierarchy with respect to the Mach number
following from the Kirkwood-Bethe hypothesis is briefly reviewed in Section
VII and numerical methods to solve KBH models for the bubble dynamics, flow
field and shock waves are discussed in Section VIII. In Section IX, the
assumptions and capabilities of the Kirkwood-Bethe hypothesis are further
scrutinized and new results are presented. Conclusions are drawn and open
questions are discussed in Section X.

##  II Retrospective

The Kirkwood-Bethe hypothesis was proposed as part of a technical report
(Kirkwood and Bethe, 1942) commissioned by the research office of the United
States federal government on predicting the pressure wave produced by an
underwater explosion, completed in 1942 by John G. Kirkwood and Hans A. Bethe.
This report focused on obtaining an approximation for the pressure amplitude
at some distance from an underwater explosion and, more precisely, the decay
constant describing an approximately exponential decay of the peak pressure.
To this end, Kirkwood and Bethe focused on the propagation of the shock front,
the conditions in the liquid behind the shock front, the rapidly expanding gas
bubble of burnt detonation products, and the rarefaction wave forming inside
the gas bubble. The work of Kirkwood and Bethe was first brought to the
attention of a wider audience through the textbook of Cole (1948).

###  II.1 Theoretical foundations

Shortly after its publication, Rice and Ginell (1943, 1944) extended the work
of Kirkwood and Bethe (1942) to an explosive charge in the shape of an
infinite cylinder by considering cylindrical instead of spherical symmetry.
Although they found that the Kirkwood-Bethe hypothesis is less accurate in
cylindrical symmetry than it is in spherical symmetry, it provided a first
glimpse into the dynamics of underwater explosions that are not spherically
symmetric. The work of Rice and Ginell (1943, 1944) later became the
foundation of the studies of Naugol’nykh and Roy (1971) on the hydrodynamics
of electrical discharges in liquids and of Kedrinskii (1972) on cylindrically
symmetric underwater explosions.

In 1952, Gilmore (1952) studied the Kirkwood-Bethe hypothesis further,
focusing, however, on cavitation bubble dynamics rather than on underwater
explosions. Aside from clarifying some assumptions and exploring the
theoretical limits of the Kirkwood-Bethe assumption, Gilmore presented a
second-order differential equation of the radial motion of the gas-liquid
interface of a spherical gas bubble in a compressible fluid as a result of
ambient pressure differences, now known as the Gilmore model, which has since
become a staple of fundamental research on cavitation. As it contains
important terms that are second order with respect to the Mach number of the
gas-liquid interface motion, the Gilmore model presents an extension to the
models of Herring (1941) and Trilling (1952) that are bound to first-order
acoustics. In addition, Gilmore (1952) also derived models for the flow
velocity and pressure induced by a bubble or cavity under different
assumptions. Among these models is an ordinary differential equation (ODE) for
the flow velocity with respect to the radial coordinate along outgoing
characteristics based on the Kirkwood-Bethe hypothesis, which allows to
reconstruct the flow and acoustics resulting from a bubble collapse or
oscillation. Generalizing the Gilmore model for spherical bubbles, Kedrinskii
and Kuzavov (1980) proposed a model of the radial dynamics of bubbles in
planar and cylindrical symmetry.

Hickling and Plesset (1963) published an institutional report on a study of
the collapse and rebound of an adiabatically compressed gas cavity in 1963,
using the Kirkwood-Bethe hypothesis in conjunction with the method of
characteristics. A shortened version of their report was published as a paper
(Hickling and Plesset, 1964) the following year. To describe the flow of the
liquid surrounding the gas cavity, Hickling and Plesset presented ODEs for
both the flow velocity and pressure with respect to time along outgoing
characteristics based on the Kirkwood-Bethe hypothesis, although they did not
present a detailed derivation of these ODEs. Under the assumptions associated
with the Kirkwood-Bethe hypothesis, the ODE for the flow velocity with respect
to time ($\text{d}u/\text{d}t$) proposed by Hickling and Plesset (1963) is
equivalent to the ODE with respect to the radial coordinate
($\text{d}u/\text{d}r$) of Gilmore (1952).

Ivany (1965) revisited the theory on cavitation bubble dynamics and acoustic
emissions in a PhD thesis on the collapse of gas cavities, completed in 1965.
In this thesis, which was in parts published in a subsequent paper (Ivany and
Hammitt, 1965), Ivany presented a detailed derivation of an ODE for the flow
velocity in the liquid starting from the momentum and continuity equations by
using the Kirkwood-Bethe hypothesis and, in fact, arrived at the same ODE as
Hickling and Plesset (1963) two years prior. However, while the model
presented by Hickling and Plesset required to also solve an ODE for the
pressure, Ivany employed an explicit expression for the pressure along
outgoing characteristics obtained directly from the Kirkwood-Bethe hypothesis
using the Tait (1888) equation of state (EoS), as previously suggested by
Gilmore (1952). In the appendix of the thesis, Ivany (1965) carefully explains
how the bubble dynamics and the flow field are solved and, remarkably,
provides listings of the developed computer codes, written in the Michigan
Algorithm Decoder programming language, making it perhaps the first open-
source software for cavitation research. Shortly after, in 1966, Lilliston
(1966) provided a computer code written in Fortran IV for bubble dynamics and
the resulting flow field based on the work of Gilmore (1952) in the appendix
of a report on the collapse of a gas-filled cavity.

Akulichev  _et al._ (1968) proposed an explicit approximation for the pressure
emitted by oscillating bubbles using Riemann invariants and the Tait EoS on
the basis of the Kirkwood-Bethe hypothesis, which does not require explicit
knowledge of the generated flow velocity. Moreover, while previous studies
(Gilmore, 1952; Hickling and Plesset, 1963; Ivany, 1965) had deliberately
avoided the multivalued solutions arising at shock fronts, Akulichev  _et al._
(1968) used the rule of equal areas (Landau and Lifshitz, 1959) to treat such
multivalued solutions. As a result, Akulichev and co-workers were able to
demonstrate that when the emitted pressure wave forms a shock front, the
pressure may be severely overpredicted if the additional dissipation arising
at the shock front is not taken into account.

All the research reported in the publications discussed above employ the
(modified) Tait equation of state (Tait, 1888) to describe the liquid
properties, starting with the original work of Kirkwood and Bethe (1942). The
main shortcoming of the Tait EoS is that it drastically overpredicts the
temperature due to isentropic compression (Le Métayer and Saurel, 2016;
Radulescu, 2020; Denner, 2021). However, the Kirkwood-Bethe hypothesis is not
limited to a specific EoS for the liquid. Flynn (1957, 1978) used a tabulated
EoS based on the work of Rice and Walsh (1957) to define the properties of the
liquid. Unfortunately, however, a copy of one of Flynn’s reports (Flynn, 1957)
on this matter could not be located (the reference here relies on the account
of Ivany (1965)), while the other report (Flynn, 1978) does not present the
methodology in detail. Recently, Denner (2021) extended the Gilmore model
using the Noble-Abel stiffened-gas (NASG) EoS and, subsequently, Denner and
Schenke (2023a) used the NASG EoS successfully to model the flow and acoustics
induced by cavitation bubbles in conjunction with the models of Gilmore
(1952), Hickling and Plesset (1963), and Ivany (1965).

These milestones form the foundation of KBH models today. Concurrent and
subsequent research efforts have primarily focused on applying these models to
various engineering problems, as well as to elucidate the physics of pressure-
driven bubble dynamics and acoustics.

###  II.2 Applications

The application in the context of which the Kirkwood-Bethe hypothesis was
originally conceived are underwater explosions, whereby the expansion of burnt
gaseous detonation products and the emitted shock wave in the water are of
primary interest. Early work on this subject, starting with the study of
Kirkwood and Bethe (1942), focused on predicting the evolution and decay of
the peak pressure of the emitted shock wave in spherical (Kirkwood and Bethe,
1942; Arons, 1954; McGrath, 1966; Poché, 1972; Rogers, 1977; Best, 1991; Shan
_et al._ , 2021; Zhang  _et al._ , 2021a) and cylindrical symmetry (Rice and
Ginell, 1943, 1944; Kedrinskii, 1972; Likhterov, 2000). To this end, the main
aim was to find an accurate definition for the decay constant $\theta$ of the
approximately exponential decay, $e^{-t/\theta}$, of the peak pressure. In the
21st century, the focus has shifted away from approximating the peak pressure
in favor of using KBH models to compute and study the full range of dynamics
arising in underwater explosions, including the associated bubble dynamics
(Kwak  _et al._ , 2012; Wang  _et al._ , 2021; Zhang  _et al._ , 2021b; Li
_et al._ , 2023), pressure and shock waves (Kwak  _et al._ , 2012; Sheng  _et
al._ , 2023), and fluid-structure interaction (Mavaleix-Marchessoux, 2020).
These models have repeatedly been found to be of very high predictive quality
(Li  _et al._ , 2023; Zhang  _et al._ , 2021b; Wang  _et al._ , 2021; Sheng
_et al._ , 2023). The research into underwater explosions is closely related
to subsequent applications of the Kirkwood-Bethe hypothesis and, especially,
the Gilmore model to predict the bubble dynamics of seismic air guns for
geophysical exploration and the shock testing of defense vessels (Ziolkowski,
1970, 1998; de Graaf, Penesis, and Brandner, 2014; Bing-Shou, Hao, and Nan,
2021). For instance, Landrø and Sollie (1992) proposed a methodology based on
the Kirkwood-Bethe hypothesis to determine the pressure transients produced by
an array of air guns, the results of which Laws, Landrø, and Amundsen (1998)
found to be in very good agreement with measurements in different experimental
scenarios. Similar methods to predict the pressure field produced by air guns
were proposed by Ziolkowski (1970, 1998). KBH models were also leveraged to
study the influence of rough seas (Sertlek and Blacquiere, 2019) and shallow
water (Wehner, Svensson, and Landrø, 2020) on the acoustic signatures of
seismic air guns.

Other applications in which a strong local energy deposition produces a
gaseous bubble that emits pressure waves, and subsequently collapses and
rebounces, are cavitation bubbles induced by an electrical discharge or a
focused laser beam. Especially laser-induced cavitation has established itself
as a widely-used tool to study cavitation bubble dynamics (Akhatov  _et al._ ,
2001; Vogel  _et al._ , 2008; Supponen  _et al._ , 2016; Gonzalez-Avila,
Denner, and Ohl, 2021; Bokman  _et al._ , 2023; Reuter  _et al._ , 2023; Zhao
and Coutier-Delgosha, 2023; Mur  _et al._ , 2023; Wang  _et al._ , 2023; Mnich
_et al._ , 2024; Fu  _et al._ , 2024) under controllable, precise, and
reproducible conditions, both experimentally and computationally. Mellen
(1956) was an early adopter of the Kirkwood-Bethe hypothesis in the 1950s, in
the form of the Gilmore model, to provide a theoretical basis for an
experimental study of the evolution of a spherical cavity produced by an
electrical discharge, reporting an encouraging agreement between experimental
measurements and computational results. With respect to laser-induced
cavitation, KBH models have been used, aside from studying the resulting
bubble dynamics, to understand and quantify the formation and attenuation of
shock waves (Vogel, Busch, and Parlitz, 1996; Byun and Kwak, 2004; Oh  _et
al._ , 2018; Geng  _et al._ , 2021; Lai  _et al._ , 2022; Liang  _et al._ ,
2022; Agrež  _et al._ , 2023; Wen  _et al._ , 2023; Yang  _et al._ , 2023) and
the energy partitioning (Liang  _et al._ , 2022; Wen  _et al._ , 2023). The
results produced by these models have been shown to be in excellent agreement
with experimental measurements (Oh  _et al._ , 2018; Geng  _et al._ , 2021;
Vassholz  _et al._ , 2021; Lai  _et al._ , 2022; Wen  _et al._ , 2023).

A major driver of research into pressure-driven bubble dynamics over the past
decades have been diagnostic and therapeutic medical applications. In 1984,
Vacher, Gimenez, and Goutte (1984) used the Gilmore model to show that
ultrasound-driven bubbles emit significant energy at the second harmonic of
the driving frequency, now a widely used indicator to detect and classify
cavitation in medical ultrasound treatments (Versluis  _et al._ , 2020). Since
then, KBH models have been used to study tissue ablation (Chavrier  _et al._ ,
2000; Berlinda Law and Zhou, 2019; Wang, Lei, and Zhou, 2019), histotripsy
(Pahk  _et al._ , 2015, 2018, 2019) and lithotripsy (Sokolov, 2002)
treatments, laser surgery (Byun and Kwak, 2004; Yang  _et al._ , 2023), and
ultrasound-driven encapsulated microbubbles, including their dynamic behavior
(Ayme and Carstensen, 1989; Aymé-Bellegarda, 1990; Kreider  _et al._ , 2011;
Berlinda Law and Zhou, 2019; Gümmer, Schenke, and Denner, 2021; Li, Zou, and
Qin, 2022), bubble-bubble interactions (Qin  _et al._ , 2023a) and the
associated acoustic emissions (Kreider  _et al._ , 2011; Gümmer, Schenke, and
Denner, 2021). To this end, Gümmer, Schenke, and Denner (2021) showed a
favorable comparison of the frequency spectrum predicted by the Gilmore model
for lipid-coated microbubbles with the experimental measurements of Song,
Moldovan, and Prentice (2019). To account for the viscoelastic rheology of
tissue, Zilonova, Solovchuk, and Sheu (2018, 2019) extended the Gilmore model
to include the linear elastic solid model, also known as the Zener model,
resulting in a markedly different bubble behavior (Zilonova, Solovchuk, and
Sheu, 2018, 2019; Li, Zou, and Qin, 2022).

The Kirkwood-Bethe hypothesis has also been utilized to study other
applications of pressure-driven bubble dynamics, including food processing
(Bredihin  _et al._ , 2021), the spall formation in aluminium (Glam  _et al._
, 2014), the treatment of materials (Sonde  _et al._ , 2018; Zhang  _et al._ ,
2024), material synthesis (Barcikowski  _et al._ , 2019; Huang  _et al._ ,
2021, 2023), the preferred nucleation of crystals during solidification
(Rakita and Han, 2017), surface cleaning (Minsier and Proost, 2008), as well
as the dynamics of a collapsing cylindrical cavity driven by an annular piston
(Kedrinskii, 2022). Studies on sonoluminescence, whereby ultrasound-driven
bubbles emit a light pulse upon collapse (Brenner, Hilgenfeldt, and Lohse,
2002), have also directly benefited from KBH models, for instance to estimate
the conditions insight the bubble during collapse (Gimenez, 1982; Nazari-
Mahroo  _et al._ , 2018; Hoeppe  _et al._ , 2024) and to correlate acoustic
measurements to bubble dynamics (Lee  _et al._ , 1997; Holzfuss, Rüggeberg,
and Billo, 1998; Holzfuss, 2010). Aside from specific applications, KBH models
have been used extensively to study and elucidate fundamental pressure-driven
bubble dynamics and acoustics. Some examples of studies using the Gilmore
model to study pressure-driven bubble dynamics include studies on the natural
bubble frequency (Shima, 1970), self-similar bubble behavior (Hunter, 1960;
Akulichev, 1971), period doubling of ultrasound-driven bubbles (Akulichev,
1971), the excitation of bubbles with multiple frequencies (Moholkar, Rekveld,
and Warmoeskerken, 2000), as well as the influence of the liquid
compressibility (Fuster, Dopazo, and Hauke, 2011), heat transfer (Mahdi,
Shams, and Ebrahimi, 2010), the vapor content (Preso  _et al._ , 2024), and a
narrow confinement (Zhang  _et al._ , 2023a) on the bubble dynamics. Moreover,
the Kirkwood-Bethe hypothesis has been used to gain quantitative insight into
the acoustic emissions associated with pressure-driven bubble dynamics, such
as the acoustic noise spectrum generated by pressure-driven bubbles close to a
free surface (Likhterov, 1998; Berman and Likhterov, 2001), the acoustic
translation of bubbles (Kim  _et al._ , 2023), and the pressure and shock
waves emitted during a strong bubble collapse (Holzfuss, 2010; Liang  _et al._
, 2022; Denner and Schenke, 2023a; Wen  _et al._ , 2023; Lai  _et al._ ,
2022).

##  III The Kirkwood-Bethe hypothesis

The Kirkwood-Bethe hypothesis can be derived using basic fluid dynamics, and
the arising definitions and expressions then allow us to close the system of
governing equations. In the following, we start in Section III.1 by bringing
the conservation laws governing a fluid flow into a convenient form.
Subsequently, the Kirkwood-Bethe hypothesis is derived in Section III.2 and
its major limitations are discussed in Section III.3.

###  III.1 Governing conservation laws

The conservation of mass, momentum and energy governing a fluid flow are given
as

$\displaystyle\frac{\text{D}\rho}{\text{D}t}+\rho\nabla\cdot\mathbf{u}$
$\displaystyle=$ $\displaystyle 0$ (2)
$\displaystyle\frac{\text{D}\mathbf{u}}{\text{D}t}+\frac{1}{\rho}\nabla p$
$\displaystyle=$ $\displaystyle\frac{1}{\rho}\,\nabla\cdot\boldsymbol{\tau}$
(3)
$\displaystyle\frac{\text{D}h}{\text{D}t}-\frac{1}{\rho}\frac{\text{D}p}{\text{D}t}$
$\displaystyle=$
$\displaystyle\frac{1}{\rho}\left(\boldsymbol{\tau}:\nabla\mathbf{u}-\nabla\cdot\dot{q}\right),$
(4)

where $\rho$ is the fluid density, $\mathbf{u}$ is the flow velocity, $p$
denotes the pressure, $h$ is the specific enthalpy, $\boldsymbol{\tau}$
represents viscous stresses, $\dot{q}$ is the heat flux due to thermal
conduction, and

$\frac{\text{D}}{\text{D}t}=\frac{\partial}{\partial
t}+\left(\mathbf{u}\cdot\nabla\right)$ (5)

denotes the material derivative operator. The thermodynamic relation
$\text{d}h=T\,\text{d}s+\text{d}p/\rho$ allows us to reformulate the energy
equation in terms of the specific entropy $s$, such that

$\frac{\text{D}s}{\text{D}t}=\frac{1}{\rho
T}\left(\boldsymbol{\tau}:\nabla\mathbf{u}-\nabla\cdot\dot{q}\right),$ (6)

where $T$ is the temperature. This system of equations is closed with a
suitable equation of state (EoS) that defines a relationship between the
thermodynamic quantities. Since the subsequent derivation is agnostic to the
applied EoS, we proceed without assuming a specific EoS and revisit commonly
used EoS in Section IV.

Assuming that the considered flow is ideal, without viscous dissipation
($\boldsymbol{\tau}=0$) and heat conduction ($\dot{q}=0$), we obtain the Euler
equations:

$\displaystyle\frac{\text{D}\rho}{\text{D}t}+\rho\nabla\cdot\mathbf{u}$
$\displaystyle=$ $\displaystyle 0$ (7)
$\displaystyle\frac{\text{D}\mathbf{u}}{\text{D}t}+\frac{1}{\rho}\nabla p$
$\displaystyle=$ $\displaystyle 0$ (8)
$\displaystyle\frac{\text{D}s}{\text{D}t}$ $\displaystyle=$ $\displaystyle 0.$
(9)

This set of governing equations, with which Kirkwood and Bethe (1942)
commenced their work, assumes that changes to the state of the fluid of a
sufficiently smooth flow are adiabatic and reversible. Consequently, the flow
is isentropic, as imposed by Eq. (9). Especially for water, this turns out to
be a good assumption even for pressure waves with a large amplitude, because
the resulting change in entropy is small. The previously suggested upper limit
for the validity of this approximation is
$p=\mathcal{O}(\text{GPa})$.(Kirkwood and Bethe, 1942; Hickling and Plesset,
1964) Assuming that the flow is isentropic does, however, not imply that all
fluid elements are on the same isentrope, because the entropy may change at a
shock front, which is further discussed in Section VIII.3.

The differential relation of the specific enthalpy for an isentropic fluid is

$\mathrm{d}h=\frac{\mathrm{d}p}{\rho}=c^{2}\,\frac{\mathrm{d}\rho}{\rho},$
(10)

with the speed of sound $c$ defined as

$c^{2}=\left(\frac{\partial p}{\partial\rho}\right)_{s},$ (11)

where subscript $s$ implies a constant entropy. In addition, to account for a
time-varying pressure in the far field, we introduce the ambient enthalpy
$h_{\infty}$ and redefine, in a slight abuse of notation, the specific
enthalpy difference as $h\rightarrow h-h_{\infty}$. The ambient enthalpy
$h_{\infty}$ is spatially invariant ($\nabla h_{\infty}=0$), such that

$\frac{\text{D}h_{\infty}}{\text{D}t}=\frac{\partial h_{\infty}}{\partial t}.$
(12)

Inserting Eq. (10) in the form $\text{d}\rho=\rho\,\text{d}h/c^{2}$ alongside
the ambient enthalpy $h_{\infty}$ into the continuity equation, Eq. (7),
yields

$\frac{1}{c^{2}}\,\frac{\text{D}h}{\text{D}t}+\nabla\cdot\mathbf{u}=\frac{1}{c^{2}}\,\frac{\partial
h_{\infty}}{\partial t}.$ (13)

For convenience, we also separate the kinetic energy and the rotation of the
flow by reformulating the convective velocity derivative as

$\left(\mathbf{u}\cdot\nabla\right)\mathbf{u}=\nabla\left(\frac{\mathbf{u}^{2}}{2}\right)-\mathbf{u}\times(\nabla\times\mathbf{u}),$
(14)

with which the momentum equation, Eq. (8), becomes

$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+\nabla\left(\frac{\mathbf{u}^{2}}{2}\right)-\mathbf{u}\times(\nabla\times\mathbf{u})+\nabla
h$ $\displaystyle=$ $\displaystyle 0.$ (15)

Note that, since the density and pressure have been replaced by the specific
enthalpy and the speed of sound on the basis of the flow being isentropic, the
entropy equation, Eq. (9), is now obsolete.

###  III.2 Hypothesis

The Kirkwood-Bethe hypothesis is typically derived and applied under the
assumption of spherical symmetry. However, it can be derived more generally
for one-dimensional problems, e.g. exploiting cylindrical symmetry (Rice and
Ginell, 1943), to explore the influence of the shape of a bubble on its
dynamics and the induced flow field. Here, we consider an arbitrary (quasi)
one-dimensional domain, illustrated in Figure 1, characterized by the
dimensionality coefficient $\alpha$ that defines the effective dimensions of
the system: $\alpha=0$ describing a planar case, $\alpha=1$ describing
cylindrical symmetry, and $\alpha=2$ describing the conventionally assumed
spherical symmetry.

Figure 1: Schematic illustration of the considered (quasi) one-dimensional
domains: one-dimensional planar ($\alpha=0$), cylindrical symmetry
($\alpha=1$), and spherical symmetry ($\alpha=2$). The one-dimensional spatial
coordinate is denoted with $r$ and the width or radius of the gas cavity or
bubble is denoted with $R$.

Considering a (quasi) one-dimensional flow by exploiting the symmetry of the
system, we know that the flow is irrotational, with
$\nabla\times\mathbf{u}=0$. The continuity equation, Eq. (13), and the
momentum equation, Eq. (15), along the one-dimensional spatial coordinate $r$
are, therefore, given as

$\displaystyle\frac{1}{c^{2}}\left(\frac{\partial h}{\partial
t}+u\,\frac{\partial h}{\partial
r}\right)+\frac{1}{r^{\alpha}}\,\frac{\partial}{\partial
r}\left(r^{\alpha}u\right)$ $\displaystyle=$
$\displaystyle\frac{1}{c^{2}}\,\frac{\partial h_{\infty}}{\partial t}$ (16)
$\displaystyle\frac{\partial u}{\partial t}+\frac{\partial}{\partial
r}\left(\frac{u^{2}}{2}\right)+\frac{\partial h}{\partial r}$ $\displaystyle=$
$\displaystyle 0.$ (17)

Since the flow is irrotational, the velocity may be expressed in terms of the
potential $\psi$,

$u=-\frac{\partial\psi}{\partial r},$ (18)

with which the continuity and momentum equations can be written as

$\displaystyle\frac{1}{c^{2}}\left(\frac{\partial h}{\partial
t}-\frac{\partial h_{\infty}}{\partial t}+u\,\frac{\partial h}{\partial
r}\right)-\frac{1}{r^{\alpha}}\frac{\partial}{\partial
r}\left(r^{\alpha}\frac{\partial\psi}{\partial r}\right)\\!\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!0$ (19)
$\displaystyle-\frac{\partial}{\partial t}\frac{\partial\psi}{\partial
r}+\frac{\partial}{\partial r}\left(\frac{u^{2}}{2}\right)+\frac{\partial
h}{\partial r}\\!\\!\\!$ $\displaystyle=$ $\displaystyle\\!\\!0.$ (20)

Integrating Eq. (20) from $r$ to $\infty$, assuming the potential $\psi$
vanishes at infinity, yields the transient Bernoulli equation,

$-\frac{\partial\psi}{\partial t}+\frac{u^{2}}{2}+h-h_{\infty}=0,$ (21)

based on which Kirkwood and Bethe (1942) defined the specific kinetic enthalpy
as

$\Omega=\frac{\partial\psi}{\partial t}=h-h_{\infty}+\frac{u^{2}}{2}.$ (22)

Rearranging Eq. (22) for the enthalpy difference $h-h_{\infty}$ and inserting
the resulting expression into the continuity equation, Eq. (19), yields

$\frac{1}{c^{2}}\left(\frac{\partial^{2}\psi}{\partial
t^{2}}+u\frac{\partial}{\partial r}\frac{\partial\psi}{\partial
t}-\frac{1}{2}\frac{\partial u^{2}}{\partial t}-\frac{1}{2}u\frac{\partial
u^{2}}{\partial r}\right)\\\ -\frac{1}{r^{\alpha}}\frac{\partial}{\partial
r}\left(r^{\alpha}\frac{\partial\psi}{\partial r}\right)=0.$ (23)

Noting that the last term on the left-hand side is the Laplacian of the
velocity potential,

$\nabla^{2}\psi=\frac{1}{r^{\alpha}}\frac{\partial}{\partial
r}\left(r^{\alpha}\frac{\partial\psi}{\partial r}\right),$ (24)

Eq. (23) can be rewritten as a wave equation with a nonlinear correction of
order $\mathcal{O}(u^{2})$, (Kirkwood and Bethe, 1942; Keller and Kolodner,
1956)

$\frac{1}{c^{2}}\,\frac{\partial^{2}\psi}{\partial
t^{2}}-\nabla^{2}\psi=\frac{1}{c^{2}}\left(\frac{\partial u^{2}}{\partial
t}+\frac{1}{2}u\frac{\partial u^{2}}{\partial r}\right),$ (25)

suggesting that the velocity potential $\psi$ propagates like a wave. As
discussed in detail by Prosperetti (1984), the finite propagation speed is
unimportant in the near field, where the flow behaves similar to an
incompressible fluid ($c\rightarrow\infty$), with $\nabla^{2}\psi=0$, while
the perturbations are small in the far field, such that terms of order
$\mathcal{O}(u^{2})$ become negligible (Akulichev  _et al._ , 1968), and Eq.
(25) reduces to the classic form of the wave equation,

$\frac{1}{c^{2}}\,\frac{\partial^{2}\psi}{\partial t^{2}}-\nabla^{2}\psi=0.$
(26)

Comparing the magnitude of the right-hand side of Eq. (25) to the magnitude of
the second-order time derivative of the velocity potential on the left-hand
side may be used to quantify the local error of using a wave equation to
propagate acoustic information in this case (Mavaleix-Marchessoux, 2020).

The work of Kirkwood and Bethe (1942) rests on a reformulated velocity
potential in spherical symmetry, defined as $\phi=r\psi$. With this new
potential, the wave equation, Eq. (25), becomes

$\frac{1}{c^{2}}\,\frac{\partial^{2}\phi}{\partial
t^{2}}-\frac{\partial^{2}\phi}{\partial
r^{2}}=\frac{r}{c^{2}}\left(\frac{\partial u^{2}}{\partial
t}+\frac{1}{2}u\frac{\partial u^{2}}{\partial r}\right)$ (27)

and the transient Bernoulli equation, Eq. (21), yields

$\frac{\partial\phi}{\partial t}=r\,\frac{\partial\psi}{\partial
t}=r\left(h-h_{\infty}+\frac{u^{2}}{2}\right),$ (28)

which is, in fact, the quantity defined in Eq. (1). Adopting this principle
for the general one-dimensional domain considered here, with $\alpha$ defining
the symmetry of the problem, the reformulated velocity potential is defined as
$\phi=r^{\alpha/2}\psi$, originally proposed by Rice and Ginell (1943). The
wave equation then reads as

$\frac{1}{c^{2}}\,\frac{\partial^{2}\phi}{\partial
t^{2}}-\frac{\partial^{2}\phi}{\partial
r^{2}}=\frac{r^{\alpha/2}}{c^{2}}\left(\frac{\partial u^{2}}{\partial
t}+\frac{1}{2}u\frac{\partial u^{2}}{\partial r}\right)$ (29)

and the transient Bernoulli equation yields

$\frac{\partial\phi}{\partial t}=r^{\alpha/2}\,\frac{\partial\psi}{\partial
t}=r^{\alpha/2}\left(h-h_{\infty}+\frac{u^{2}}{2}\right).$ (30)

The validity of this definition of the velocity potential is further discussed
together with the limitations presented by the dimensionality coefficient
$\alpha$ in Section III.3.

In principle, the method of characteristics may be applied to solve the wave
equation governing the velocity potential and any of its derivatives, even
though an analytical solution is not readily obtainable due to the
contributions of $\mathcal{O}(u^{2})$.(Akulichev  _et al._ , 1968) The
Kirkwood-Bethe hypothesis stipulates that the quantity
$g={\partial\phi}/{\partial t}$, as given by Eq. (30), is constant along an
outgoing characteristic, satisfying (Kirkwood and Bethe, 1942)

$\left.\frac{\mathrm{d}g}{\mathrm{d}t}\right|_{\mathrm{c}}=0,$ (31)

where the derivative operator along an outgoing characteristic is defined as

$\left.\frac{\mathrm{d}}{\mathrm{d}t}\right|_{\mathrm{c}}=\frac{\partial}{\partial
t}+\left.\frac{\mathrm{d}r}{\mathrm{d}t}\right|_{\mathrm{c}}\frac{\partial}{\partial
r}.$ (32)

Observing that small disturbances to the velocity potential of the fluid
propagate like a wave according to Eq. (25), i.e. with speed $c$ relative to
the flow, Kirkwood and Bethe (1942) proposed to approximate the propagation
speed along an outgoing characteristic by

$\left.\frac{\mathrm{d}r}{\mathrm{d}t}\right|_{\mathrm{c}}=c+u.$ (33)

###  III.3 Limitations

The assumption of a constant value of $g=r^{\alpha/2}(h-h_{\infty}+u^{2}/2)$
along outgoing characteristics can be justified by the conservation of
specific energy when dissipative effects are neglected, as it follows from the
transient Bernoulli equation. However, the assumptions that the acoustic
potential is $\phi=r^{\alpha/2}\psi$ and propagates with speed $c+u$ may
entail errors.

A plane finite-amplitude acoustic wave moving in one direction (cf. simple
wave) propagates with speed $c+u$, (Rayleigh, 1910) and $c\pm u$ are the
eigenvalues of the Euler equations for plane waves,(LeVeque, 2008) meaning
that plane waves travel with the local speed of sound $c$ relative to the flow
(Landau and Lifshitz, 1959). Moreover, Eq. (33) reduces to the asymptotic
speed $c$ in the far field, since $u\rightarrow 0$ for $r\rightarrow\infty$,
and is consistent with the Chapman-Jouguet condition (Chapman, 1899; Jouguet,
1905) for detonation waves. However, the propagation speed $c+u$ is not exact.
This can be seen by considering the Riemann function $\varsigma$ for an
isentropic flow (Riemann, 1860),

$\varsigma=\int_{p_{\infty}}^{p}\frac{\text{d}p}{\rho
c}=\int_{h_{\infty}}^{h}\frac{\text{d}h}{c}.$ (34)

Introducing $\varsigma$ into the continuity and momentum equations, Eqs. (16)
and (17), yields

$\displaystyle\frac{\partial}{\partial
t}(\varsigma+u)+(c+u)\frac{\partial}{\partial r}(\varsigma+u)$
$\displaystyle=$ $\displaystyle-\frac{\alpha uc}{r}$ (35)
$\displaystyle\frac{\partial}{\partial
t}(\varsigma-u)+(c-u)\frac{\partial}{\partial r}(\varsigma-u)$
$\displaystyle=$ $\displaystyle-\frac{\alpha uc}{r},$ (36)

with the quantities $\varsigma+u$ and $\varsigma-u$ being so-called Riemann
invariants (Riemann, 1860). For a plane wave ($\alpha=0$), the Riemann
invariant $\varsigma+u$ propagates with constant value along the right-going
characteristic with speed $c+u$, whereas $\varsigma-u$ propagates with
constant value along the left-going characteristic with speed $c-u$. When the
right-hand side, which accounts for the geometric divergence of, for instance,
cylindrical ($\alpha=1$) or spherical ($\alpha=2$) domains, is not zero, the
Riemann invariants $\varsigma+u$ and $\varsigma-u$ do not propagate with $c+u$
and $c-u$, respectively. Hence, quantities that are invariant along outgoing
characteristics of a curved wave do not propagate with speed $c+u$, rendering
Eq. (33) an approximation that becomes more accurate as a curved wave moves
outwards or when the velocity is small with respect to the speed of sound
(Gilmore, 1952). Assuming the (dominant) fluid flow is the result of the
bubble motion and the associated acoustic emissions only, the velocity is
$u\propto r^{-\alpha/2}$ in the far field and, hence, the right-hand side of
Eq. (35) decays rapidly with $r^{-n}$, where $n\geq 1+\alpha/2$.

Considering a plane acoustic wave, the acoustic particle velocity
$u_{1}=p_{1}/(\rho c)$ is in phase with the acoustic pressure $p_{1}$, whereby
$p_{1}$ and $u_{1}$ follow from the expansions $p=p_{0}+p_{1}$ and
$u=u_{0}+u_{1}$, and where subscript $0$ denotes the reference or ambient
properties of the flow field. However, $u_{1}$ and $p_{1}$ of a curved wave
are not in phase (Kinsler, 2000). As a consequence, the specific acoustic
impedance $z=p_{1}/u_{1}$ experienced by a curved wave is complex, given for a
spherical harmonic wave as (Kinsler, 2000; Randall, 2005)

$z=\rho c\,\frac{(kr)^{2}+ikr}{1+(kr)^{2}},$ (37)

which is dependent on the wavenumber $k$ and radial coordinate $r$, and where
$i$ denotes the imaginary unit. For large $r$, the imaginary part (specific
acoustic reactance) vanishes and the real part (specific acoustic resistance)
approaches the specific acoustic impedance of a plane wave, $z=\rho c$,
whereas the imaginary part dominates for small $r$. KBH models do not account
for this complex acoustic impedance and the associated error is investigated
in Section IX.2.

The dimensionality coefficient $\alpha$ defines the dimension and symmetry of
the considered problem. Inserting $\phi=r^{\alpha/2}\psi$ in the corresponding
wave equation, Eq. (29), reads as

$\frac{1}{c^{2}}\,\frac{\partial^{2}(r^{\alpha/2}\psi)}{\partial
t^{2}}-\frac{\partial^{2}(r^{\alpha/2}\psi)}{\partial
r^{2}}=\frac{r^{\alpha/2}}{c^{2}}\left(\frac{\partial u^{2}}{\partial
t}+\frac{1}{2}u\frac{\partial u^{2}}{\partial r}\right),$ (38)

which, after expanding the left-hand side and dividing by $r^{\alpha/2}$,
becomes

$\frac{1}{c^{2}}\,\frac{\partial^{2}\psi}{\partial
t^{2}}-\underbrace{\left(\frac{\partial^{2}\psi}{\partial
r^{2}}+\frac{\alpha}{r}\frac{\partial\psi}{\partial
r}\right)}_{\nabla^{2}\psi,\text{\leavevmode\nobreak\ Eq.\leavevmode\nobreak\
\eqref{eq:laplacian}}}-\frac{\alpha(\alpha-2)}{4r^{2}}\psi\\\
=\frac{1}{c^{2}}\left(\frac{\partial u^{2}}{\partial
t}+\frac{1}{2}u\frac{\partial u^{2}}{\partial r}\right).$ (39)

The last term on the left-hand side vanishes for plane ($\alpha=0$) and
spherical ($\alpha=2$) waves, such that Eq. (39) is identical to Eq. (25).
This means that the potential $\psi$ and the shape of plane and spherical
waves of small amplitude, where $u^{2}\ll c^{2}$, remain unchanged as they
propagate away from the emitter. However, the corresponding wave equation for
cylindrical waves ($\alpha=1$),

$\frac{1}{c^{2}}\,\frac{\partial^{2}\psi}{\partial
t^{2}}-\nabla^{2}\psi+\frac{\psi}{4r^{2}}=\frac{1}{c^{2}}\left(\frac{\partial
u^{2}}{\partial t}+\frac{1}{2}u\frac{\partial u^{2}}{\partial r}\right),$ (40)

contains a spurious term proportional to $1/r^{2}$ and provides an accurate
solution only in the asymptotic limit $r\rightarrow\infty$. (Rice and Ginell,
1943; Kedrinskii, 2005) In fact, in the interval $0\leq\alpha\leq 2$, the
spurious term in Eq. (39) is largest for $\alpha=1$, and Rice and Ginell
(1943, 1944) reported more accurate results for cylindrical underwater
explosions with $\alpha=0.8$.

In summary, the assumptions underpinning the Kirkwood-Bethe hypothesis are,
strictly speaking, valid only when either the wavelength is relatively small,
such that the plane wave assumption is valid (Blackstock, 1964), or when the
Mach number of the induced flow is small. The various modelling errors of the
Kirkwood-Bethe hypothesis vanish in the asymptotic limits
$r\rightarrow\infty$, representing the far field, or $u\rightarrow 0$,
representing the linear acoustic regime. More generally, the Kirkwood-Bethe
hypothesis can be regarded as an extrapolation of the linear acoustic theory
(Prosperetti, 1984) and, as briefly revisited in Section VII, is part of a
consistent hierarchy of models (Denner and Schenke, 2023a; Zhang  _et al._ ,
2023b). Nevertheless, as further discussed in Section IX, the Kirkwood-Bethe
hypothesis provides accurate predictions under the specific conditions
associated with pressure-driven bubble dynamics, cavitation and underwater
explosions.

##  IV Liquid equation of state

The Kirkwood-Bethe hypothesis does not assume a specific equation of state
(EoS) for the liquid, even though it is almost always presented and used in
conjunction with the modified Tait EoS (Tait, 1888), which Kirkwood and Bethe
(1942) used in their original work. Similarly, the influential works of Cole
(1948) on underwater explosions and Gilmore (1952) on cavitation also used the
modified Tait EoS. The reasons for the popularity of this EoS are rather
straightforward: it has a simple analytical form and it provides an accurate
(isentropic) relationship between the pressure $p$, the specific enthalpy $h$,
the density $\rho$, and the speed of sound $c$ of water, even when compared
against modern standards, such as the IAPWS R6-95(2018) standard (Wagner and
Pruß, 2002).

Figure 2: Density $\rho$, speed of sound $c$, and temperature $T$ of water
obtained with the modified Tait EoS, using the model constants of Cole (1948)
($n=7.15$ and $B=3047\times 10^{5}\,\text{Pa}$), compared against the IAPWS
R6-95(2018) standard. The reference values are $p_{0}=10^{5}\,\text{Pa}$,
$\rho_{0}=997\,\text{kg/m}^{3}$, and $T_{0}=298.3\,\text{K}$.

While the original Tait EoS is formulated along isotherms (Tait, 1888),

$-\left(\frac{\partial v}{\partial
p}\right)_{T}=\frac{1}{n\left[p+B(T)\right]},$ (41)

where $v=1/\rho$ is the specific volume, $n$ is the polytropic exponent, and
$B$ is a pressure constant, the modified Tait EoS used widely in conjunction
with the Kirkwood-Bethe hypothesis is formulated along isentropes (Cole, 1948)

$-\frac{1}{v}\left(\frac{\partial v}{\partial
p}\right)_{s}=\frac{1}{n\left[p+B(s)\right]},$ (42)

which is, in essence, the isentropic (or, more generally, polytropic) form of
the stiffened-gas EoS (Harlow and Amsden, 1971; Le Métayer, Massoni, and
Saurel, 2004). The modified Tait EoS defines the fluid properties based on the
pressure $p$ as

$\displaystyle h$ $\displaystyle=$
$\displaystyle\frac{n}{n-1}\frac{p+B}{\rho}$ (43) $\displaystyle\rho$
$\displaystyle=$
$\displaystyle\rho_{0}\left(\frac{p+B}{p_{0}+B}\right)^{\frac{1}{n}}$ (44)
$\displaystyle c$ $\displaystyle=$ $\displaystyle\sqrt{n\,\frac{p+B}{\rho}},$
(45)

where subscript $0$ defines a reference state. For water, commonly used values
are $n\approx 7$ and $B\approx 3\times 10^{8}\,\text{Pa}$, or
$B\approx\rho_{0}c_{0}^{2}/n$.(Kirkwood and Bethe, 1942) The reference
density, pressure and speed of sound of water are typically assumed to be
$\rho_{0}\approx 1000\,\text{kg/m}^{3}$, $p_{0}\approx 10^{5}\,\text{Pa}$, and
$c_{0}\approx 1500\,\text{m/s}$, respectively. Figure 2 shows the density,
speed of sound and temperature predicted by the Tait EoS in comparison to the
IAPWS R6-95(2018) standard.

The polytropic exponent for an isentropic flow represents the ratio of
specific heats, which is approximately unity for water and other liquids.
Thus, the large polytropic exponent for water, $n\approx 7$, points to the
main shortcoming of the Tait EoS: it severely overpredicts changes in
temperature as a result of adiabatic compression or expansion (Radulescu,
2020; Denner, 2021),

$T=T_{0}\left(\frac{p+B}{p_{0}+B}\right)^{(n-1)/n},$ (46)

as seen in Figure 2. In fact, during a violent bubble collapse, the Tait EoS
may predict the temperature of the liquid in the vicinity of the gas-liquid
interface to be higher than the temperature inside the gas bubble (Denner,
2021), which is unphysical on account of the considerably larger heat capacity
of liquids compared to gases.

Table 1: Model constants of the NASG EoS for water proposed in recent years. |  Le Métayer and Saurel (2016) |  Chandran and Salih (2019) |  Denner and Schenke (2023a)
---|---|---|---
$n$ | 1.19 | 1.19 | 1.11
$B\ [\text{Pa}]$ | $7028\times 10^{5}$ | $6218\times 10^{5}$ | $6480\times 10^{5}$
$b\ [\text{m}^{3}/\text{kg}]$ | $6.61\times 10^{-4}$ | $6.72\times 10^{-4}$ | $6.80\times 10^{-4}$
$\rho_{0}\ [\text{kg/m}^{3}]$ | $957.7$ | $997.0$ | $997.0$
$p_{0}\ [\text{Pa}]$ | $1.0453\times 10^{5}$ | $10^{5}$ | $10^{5}$

To rectify this shortcoming of the Tait EoS and, similarly, the stiffened-gas
EoS, Le Métayer and Saurel (2016) proposed the Noble-Abel stiffened-gas (NASG)
EoS, which had already been postulated in similar (isothermal) form by Tammann
(1912). The NASG EoS extends the Tait and stiffened-gas EoS by introducing the
co-volume $b$ that represents the finite volume occupied by the fluid
molecules, rendering the NASG EoS a simplified version of the van-der-Waals
EoS (Kontogeorgis, Privat, and Jaubert, 2019) with constant coefficients $B$
and $b$. For an isentropic flow the specific enthalpy, density and speed of
sound are defined as (Denner, 2021)

$\displaystyle h$ $\displaystyle=$
$\displaystyle\frac{n}{n-1}\frac{p+B}{\rho}-\frac{n\,b}{n-1}\,(p+B)+b\,p$ (47)
$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{K\,(p+B)^{\frac{1}{n}}}{1+b\,K\,(p+B)^{\frac{1}{n}}}$ (48)
$\displaystyle c$ $\displaystyle=$
$\displaystyle\sqrt{n\,\frac{p+B}{\rho-b\rho^{2}}},$ (49)

where $K=\rho_{0}/[(p_{0}+B)^{{1/n}}\ (1-b\,\rho_{0})]$ describes a constant
reference state. Just like the Tait EoS, to which it reduces for $b=0$, the
NASG EoS is unconditionally convex (Le Métayer and Saurel, 2016), making it an
equally robust alternative to the Tait EoS. The temperature as a result of
adiabatic compression and expansion is also for the NASG EoS given by Eq.
(46). Different model constants for the NASG EoS of water have been proposed
in recent years, as listed in Table 1 and shown in Figure 3. In particular the
model constants for isentropic water proposed by Denner and Schenke (2023a)
exhibit a very good agreement with the IAPWS R6-95(2018) standard. The NASG
EoS was first used in conjunction with the Kirkwood-Bethe hypothesis in the
Gilmore-NASG model for bubble dynamics by Denner (2021) and, using the theory
derived in Section VI for spherical bubbles, to model the full flow field
introduced by cavitation bubbles by Denner and Schenke (2023a).

Figure 3: Density $\rho$, speed of sound $c$, and temperature $T$ of water
obtained with the NASG EoS, using the model constants listed in Table 1,
compared against the IAPWS R6-95(2018) standard. The reference temperature is
$T_{0}=298.3\,\text{K}$.

Other equations of state, such as the van-der-Waals EoS or any tabulated EoS,
may also be applied in conjunction with the Kirkwood-Bethe hypothesis.
However, the simplicity and robustness of the Tait EoS and, more recently, the
NASG EoS have prevailed so far.

##  V Bubble dynamics

The propagation of the invariant $g$ along outgoing characteristics, as
stipulated by the Kirkwood-Bethe hypothesis, can be readily employed to derive
the equation of motion of the gas-liquid interface of a gas bubble or cavity.

Expanding Eq. (31) under consideration of the time-dependent ambient enthalpy
defined by Eq. (12) yields

$r\frac{\text{D}h}{\text{D}t}-r\frac{\text{D}h_{\infty}}{\text{D}t}+ru\frac{\text{D}u}{\text{D}t}+rc\left(\frac{\partial
h}{\partial r}+u\frac{\partial u}{\partial r}\right)\\\
+\frac{\alpha}{2}\left(h-h_{\infty}+\frac{u^{2}}{2}\right)\left.\frac{\mathrm{d}r}{\mathrm{d}t}\right|_{\mathrm{c}}=0.$
(50)

Expressions for the spatial derivatives of the flow velocity and specific
enthalpy follow by rearranging the continuity and momentum equations, Eqs.
(16) and (17), as

$\displaystyle\frac{\partial u}{\partial r}$ $\displaystyle=$
$\displaystyle-\frac{1}{c^{2}}\left(\frac{\text{D}h}{\text{D}t}-\frac{\text{D}h_{\infty}}{\text{D}t}\right)-\frac{\alpha
u}{r}$ (51) $\displaystyle\frac{\partial h}{\partial r}$ $\displaystyle=$
$\displaystyle-\frac{\text{D}u}{\text{D}t},$ (52)

respectively. Inserting Eqs. (33), (51), and (52) into Eq. (50), we obtain

$r\left(\frac{\text{D}h}{\text{D}t}-\frac{\text{D}h_{\infty}}{\text{D}t}\right)+r\frac{\text{D}u}{\text{D}t}(u-c)-\frac{ru}{c}\left(\frac{\text{D}h}{\text{D}t}-\frac{\text{D}h_{\infty}}{\text{D}t}\right)\\\
-\alpha{cu^{2}}+\frac{\alpha}{2}\left(h-h_{\infty}+\frac{u^{2}}{2}\right)(c+u)=0.$
(53)

Evaluating this equation at the gas-liquid interface ($r=R$) yields a second-
order ODE describing the change in radius $R$ of a gas bubble or cavity, given
as

$\left(1-\frac{\dot{R}}{c_{\mathrm{L}}}\right)R\ddot{R}+\frac{3}{4}\alpha\left(1-\frac{\dot{R}}{3c_{\mathrm{L}}}\right)\dot{R}^{2}\\\
=\frac{\alpha}{2}\left(1+\frac{\dot{R}}{c_{\mathrm{L}}}\right)H+\left(1-\frac{\dot{R}}{c_{\mathrm{L}}}\right)\frac{R\dot{H}}{c_{\mathrm{L}}},$
(54)

with $H=h_{\mathrm{L}}-h_{\infty}$ and
$\dot{H}=\dot{h}_{\mathrm{L}}-\dot{h}_{\infty}$, where
$\dot{\lozenge}={\text{D}\lozenge}/{\text{D}t}$ denotes the material
derivative at the gas-liquid interface and the subscript $\mathrm{L}$ denotes
quantities of the liquid at the gas-liquid interface. In this general form,
Eq. (54) was first presented by Kedrinskii and Kuzavov (1980), while in
spherical symmetry ($\alpha=2$), Eq. (54) is the well-known Gilmore equation
(Gilmore, 1952)

$\left(1-\frac{\dot{R}}{c_{\mathrm{L}}}\right)R\ddot{R}+\frac{3}{2}\left(1-\frac{\dot{R}}{3c_{\mathrm{L}}}\right)\dot{R}^{2}\\\
=\left(1+\frac{\dot{R}}{c_{\mathrm{L}}}\right)H+\left(1-\frac{\dot{R}}{c_{\mathrm{L}}}\right)\frac{R\dot{H}}{c_{\mathrm{L}}}.$
(55)

Note that the assumption of a cylindrical or spherical bubble, as well as the
assumption of a spatially invariant ambient enthalpy $h_{\infty}$, are only
valid if $\lambda_{\infty}\gg R$, (Neppiras, 1980; Tilmann, 1980) where
$\lambda_{\infty}$ is the wavelength of $h_{\infty}(t)$.

The specific enthalpy $h_{\mathrm{L}}$ incorporates the kinematic boundary
conditions at the gas-liquid interface. Considering a Newtonian liquid and a
clean gas-liquid interface without surface-active substances, the liquid
pressure $p_{\mathrm{L}}$ at the gas-liquid interface is defined as (Ilinskii
_et al._ , 2012)

$p_{\mathrm{L}}=p_{\text{G}}-\frac{\alpha\sigma}{R}-2\alpha\,\frac{\mu\dot{R}}{R},$
(56)

where $\sigma$ is the surface tension coefficient and $\mu$ is the dynamic
viscosity of the liquid. The Gilmore model has also been used in conjunction
with viscoelastic materials surrounding the gas bubble (Zilonova, Solovchuk,
and Sheu, 2018; Denner and Schenke, 2023b; Qin  _et al._ , 2023a; Li, Zou, and
Qin, 2022) and considering a coated gas-liquid interface (Gümmer, Schenke, and
Denner, 2021; Denner and Schenke, 2023a; Li, Zou, and Qin, 2022; Qin  _et al._
, 2023a).

The gas pressure inside the bubble, $p_{\text{G}}$, is defined by an
appropriate EoS as a function of the bubble size. Assuming the gas to be
ideal, the pressure in the bubble to be spatially uniform and the compression
to be adiabatic, the gas pressure of a cavitation bubble under consideration
of the dimensionality parameter $\alpha$ is given as

$p_{\text{G}}=p_{\text{G},0}\left(\frac{R_{0}^{\alpha+1}}{R^{\alpha+1}}\right)^{\gamma},$
(57)

where $p_{\text{G},0}$ is the gas pressure at initial radius $R_{0}$ and
$\gamma$ is the ratio of specific heats (or, more generally, the polytropic
exponent) of the gas. The finite volume occupied by the gas molecules, which
becomes important during a strong bubble collapse, can be taken into account
by the van-der-Waals hardcore radius $r_{\text{hc}}$,(Brenner, Hilgenfeldt,
and Lohse, 2002) such that

$p_{\text{G}}=p_{\text{G},0}\left(\frac{R_{0}^{\alpha+1}-r_{\text{hc}}^{\alpha+1}}{R^{\alpha+1}-r_{\text{hc}}^{\alpha+1}}\right)^{\gamma}.$
(58)

Alternatively, the same can be achieved with the Noble-Abel (NA) EoS
(Johnston, 2005; Toro, 2009), a simplification of the NASG EoS discussed in
Section IV, with which the gas pressure and density as a result of adiabatic
compression are defined as (Denner, 2021)

$\displaystyle p_{\text{G}}$ $\displaystyle=$ $\displaystyle
p_{\text{G},0}\left[\frac{\rho_{\text{G}}(1-b\,\rho_{\text{G,0}})}{\rho_{\text{G},0}(1-b\,\rho_{\text{G}})}\right]^{\gamma}$
(59) $\displaystyle\rho_{\text{G}}$ $\displaystyle=$
$\displaystyle\rho_{\text{G,0}}\left(\frac{R_{0}}{R}\right)^{\alpha+1},$ (60)

respectively, where $\rho_{\text{G,0}}$ is the gas density at pressure
$p_{\text{G},0}$ and $b$ is the co-volume of the gas. More elaborate gas
models that dispense with the polytropic relation or use more complex gas
equations of state (Prosperetti, Crum, and Commander, 1988; Qin  _et al._ ,
2023b; Nazari-Mahroo  _et al._ , 2020), account for heat transfer
(Prosperetti, 1991; Stricker, Prosperetti, and Lohse, 2011; Zhou and
Prosperetti, 2020; Zhou, 2021; Samiei, Shams, and Ebrahimi, 2010), mass
transfer (de Graaf, Penesis, and Brandner, 2014; Byun and Kwak, 2004; Samiei,
Shams, and Ebrahimi, 2010; Qin  _et al._ , 2023b; Nazari-Mahroo  _et al._ ,
2020), and external damping effects (de Graaf, Penesis, and Brandner, 2014)
can also be readily applied in conjunction with the Gilmore model. For
underwater explosions, the Jones-Wilkins-Lee (JWL) EoS (Lee, Hornig, and Kury,
1968; Baudin and Serradeill, 2010) is a popular choice in conjunction with the
Chapman-Jouguet isentrope (Chapman, 1899; Jouguet, 1905) to describe the
properties of the gaseous detonation products constituting the bubble and to
provide initial conditions for the velocity of the gas-liquid interface
(Riley, 2010; Kwak  _et al._ , 2012; Zhang  _et al._ , 2021a, b; Wang  _et
al._ , 2021; Gao  _et al._ , 2022; Jia  _et al._ , 2022; Sheng  _et al._ ,
2023).

The Gilmore equation for spherical bubbles, Eq. (55), has been shown in many
studies to yield remarkably accurate results, which Prosperetti and Lezzi
(Prosperetti and Lezzi, 1986; Lezzi and Prosperetti, 1987) attributed to the
fact that the specific enthalpy of the liquid is retained in the equation of
motion. Even though the Gilmore equation contains some high-order terms, it is
not consistently second-order accurate with respect to the Mach number of the
gas-liquid interface. Benjamin (1958) pointed to missing second-order terms in
the Kirkwood-Bethe hypothesis, which was met with some skepticism by Plesset
(see the discussion on p. 232-233 in the book edited by Cooper (1958)),
although similar second-order errors in the Gilmore model (Gilmore, 1952) were
later identified by Tilmann (1980) and Prosperetti and Lezzi (1986). Lezzi and
Prosperetti (1987) proposed a more elaborate equation of motion for pressure-
driven dynamics of spherical bubbles that is second-order accurate with
respect to the Mach number, which was, however, derived under the explicit
assumption that the velocity of the induced flow is small. Recent studies
(Zhang  _et al._ , 2021a; Wang  _et al._ , 2021; Sheng  _et al._ , 2023) have
used the model of Lezzi and Prosperetti (1987) in conjunction with the
Kirkwood-Bethe hypothesis successfully to predict the pressure transients and
shock waves generated by underwater explosions.

##  VI Flow field around a bubble

The flow field generated by an oscillating or collapsing bubble, or an
underwater explosion, can be predicted with expressions derived from the
Kirkwood-Bethe hypothesis in conjunction with the appropriate conservation
laws. To this end, the original work of Kirkwood and Bethe (1942), as well as
subsequent studies by Rice and Ginell (1943, 1944), Akulichev  _et al._
(1968), and Kedrinskii (1972) focused on deriving models for the evolution of
the peak pressure of a shock wave emitted by an underwater explosion or a
pressure-driven bubble, leveraging the Kirkwood-Bethe hypothesis together with
the Rankine-Hugoniot conditions. Since the Rankine-Hugoniot conditions follow
from the governing conservation laws at shock discontinuities, these peak-
pressure models can be seen a special case of the more general descriptions of
the flow field devised by Gilmore (1952), Hickling and Plesset (1963), and
Ivany (1965).

The enthalpy of the fluid is readily obtained by rearranging Eq. (30) as

$h=h_{\infty}+\frac{g}{r^{\alpha/2}}-\frac{u^{2}}{2},$ (61)

where the invariant $g$ is defined by Eq. (30), using the conditions under
which it originates at the gas-liquid interface ($R=r$), as

$g=R^{\alpha/2}\,\Omega_{\text{L}}=R^{\alpha/2}\left(h_{\text{L}}-h_{\infty}+\frac{\dot{R}^{2}}{2}\right).$
(62)

Inserting Eq. (62) into Eq. (61), the specific enthalpy at radial location $r$
and time $t$ is

$h(r,t)=h_{\infty}(t)+\left[\frac{R(\tau)}{r(t)}\right]^{\alpha/2}\Omega_{\text{L}}(\tau)-\frac{u(r,t)^{2}}{2},$
(63)

where $\tau=t-r/(c+u)$ is the retarded time that corresponds to the time at
which the invariant $g$ or, likewise, the specific kinetic enthalpy $\Omega$
are emitted at the gas-liquid interface. This definition of the specific
enthalpy satisfies the boundary condition $h(R,\tau)=h_{\text{L}}(\tau)$ at
the gas-liquid interface ($r=R$). Considering, in addition, that $h_{\infty}$
is an explicitly imposed ambient enthalpy and the radial position $r$ can be
obtained by integrating Eq. (33), the flow velocity $u$ is the only remaining
unknown.

Eqs. (33) and (63), in conjunction with either the explicit expression, Eq.
(69), or the differential equation, Eq. (75), for the flow velocity derived
below, describe the flow field induced by a bubble or cavity with varying
size, assuming the induced flow is isentropic.

###  VI.1 Explicit expression for the flow velocity

By directly inserting the velocity potential $\psi=\phi/r^{\alpha/2}$ into the
definition of the velocity given by Eq. (18), similar to the derivation of the
quasi-acoustic description of the flow field presented by Gilmore (1952) and
Trilling (1952), we obtain

$u=\frac{\alpha\,\phi}{2\,r^{1+\alpha/2}}-\frac{1}{r^{\alpha/2}}\frac{\partial\phi}{\partial
r}.$ (64)

Making the assumption that the potential $\phi$ propagates with constant value
along outgoing characteristics,

$\left.\frac{\text{d}\phi}{\text{d}t}\right|_{\text{c}}=\frac{\partial\phi}{\partial
t}+\left.\frac{\text{d}r}{\text{d}t}\right|_{\text{c}}\frac{\partial\phi}{\partial
r}=0,$ (65)

allows us to write

$\frac{\partial\phi}{\partial
r}=-\left.\frac{\text{d}r}{\text{d}t}\right|_{\text{c}}^{-1}\frac{\partial\phi}{\partial
t}.$ (66)

Inserting this expression into Eq. (64) with the propagation velocity given by
Eq. (33) and $g=\partial\phi/\partial t$, the flow velocity is defined as

$u=\frac{\alpha\,\phi}{2\,r^{1+\alpha/2}}+\frac{g}{r^{\alpha/2}(c+u)}.$ (67)

This explicit expression of the flow velocity can already be found in the
original report of Kirkwood and Bethe (1942) for spherical symmetry
($\alpha=2$) as well as in the report of Rice and Ginell (1944) for
cylindrical symmetry ($\alpha=1$). Evaluating Eq. (67) at the gas-liquid
interface ($R=r$), we find the invariant potential $\phi$ to be defined as

$\displaystyle\phi=\frac{2}{\alpha}\left(R^{1+\alpha/2}\dot{R}-R\,\frac{g}{c_{\text{L}}+\dot{R}}\right).$
(68)

Inserting Eqs. (62) and (68) into Eq. (67), the velocity at radial location
$r$ and time $t$ is

$u(r,t)=\left[\frac{R(\tau)}{r(t)}\right]^{1+\alpha/2}\left[\dot{R}(\tau)-\frac{\Omega_{\text{L}}(\tau)}{c_{\text{L}}(\tau)+\dot{R}(\tau)}\right]\\\
+\left[\frac{R(\tau)}{r(t)}\right]^{\alpha/2}\frac{\Omega_{\text{L}}(\tau)}{c(r,t)+u(r,t)},$
(69)

which satisfies the boundary condition $u(R,\tau)=\dot{R}(\tau)$ at the gas-
liquid interface ($r=R$).

The explicit expression of the velocity given by Eq. (69) separates the flow
velocity into three distinct contributions (Denner and Schenke, 2023a): the
first term is the incompressible displacement of the fluid by the moving gas-
liquid interface, the second term represents the compression/expansion of the
fluid resulting from its displacement, and the third term is the particle
velocity associated with the emitted acoustic wave. Hence, the first two terms
represent hydrodynamic contributions and the third term represents the
acoustic contribution to the flow field. However, the assumption that the
potential $\phi$ propagates as an outgoing wave with constant amplitude is
valid only if the nonlinear terms on the right-hand side of Eq. (29) are
negligible. Hence, Eq. (67) is accurate only if the flow velocity is
sufficiently small, $u^{2}\ll c^{2}$.

###  VI.2 Differential equation for the flow velocity

Alternative to the explicit expression for the flow velocity presented in the
previous section, it is also possible to derive an ODE describing the flow
velocity, as proposed by Gilmore (1952) and Hickling and Plesset (1963),
without the assumption that $\phi$ retains a constant value as it propagates
outward. To this end, the continuity equation, Eq. (16), is added to the
momentum equation, Eq. (17), to obtain

$\frac{1}{c}\left(\frac{\partial h}{\partial t}+u\frac{\partial h}{\partial
r}\right)+c\frac{\partial u}{\partial r}+\frac{\alpha uc}{r}+\frac{\partial
u}{\partial t}\\\ +u\frac{\partial u}{\partial r}+\frac{\partial h}{\partial
r}=\frac{1}{c}\,\frac{\partial h_{\infty}}{\partial t}.$ (70)

This equation may be rearranged as

$\frac{1}{c}\left[\frac{\partial h}{\partial t}+(c+u)\frac{\partial
h}{\partial r}\right]+\frac{\partial u}{\partial t}+(c+u)\frac{\partial
u}{\partial r}\\\ =-\frac{\alpha uc}{r}+\frac{1}{c}\,\frac{\partial
h_{\infty}}{\partial t},$ (71)

which is, considering that the ambient specific enthalpy $h_{\infty}$ is
assumed to be spatially invariant, equivalent to the equation for the outgoing
Riemann invariant, see Eq. (35). Thus, the derivation of a differential
equation for the flow velocity starts from the equation governing a curved
Riemann wave in an isentropic flow. This further shows that the propagation
speed $c+u$, albeit only approximately accurate, is inferred by the governing
conservation laws.

The enthalpy derivative along outgoing characteristics follows from Eq. (71)
as

$\displaystyle\left.\frac{\mathrm{d}h}{\mathrm{d}t}\right|_{\text{c}}=-c\left.\frac{\mathrm{d}u}{\mathrm{d}t}\right|_{\text{c}}-\frac{\alpha
uc^{2}}{r}+\frac{\partial h_{\infty}}{\partial t}.$ (72)

Taking the derivative of Eq. (30) along an outgoing characteristic reads as

$\left.\frac{\mathrm{d}g}{\mathrm{d}t}\right|_{\mathrm{c}}=\frac{\alpha
g}{2r}\left.\frac{\mathrm{d}r}{\mathrm{d}t}\right|_{\mathrm{c}}+r^{\alpha/2}\left.\frac{\mathrm{d}h}{\mathrm{d}t}\right|_{\mathrm{c}}\\\
-r^{\alpha/2}\left.\frac{\mathrm{d}h_{\infty}}{\mathrm{d}t}\right|_{\mathrm{c}}+r^{\alpha/2}u\left.\frac{\mathrm{d}u}{\mathrm{d}t}\right|_{\mathrm{c}}=0,$
(73)

where $(\mathrm{d}g/\mathrm{d}t)_{\mathrm{c}}=0$ by virtue of the Kirkwood-
Bethe hypothesis, and
$(\mathrm{d}h_{\infty}/\mathrm{d}t)_{\mathrm{c}}=\partial h_{\infty}/\partial
t$, since $h_{\infty}$ is spatially invariant by definition. Inserting Eq.
(73) into Eq. (72),

$\frac{\partial h_{\infty}}{\partial t}-\frac{\alpha
g}{2r^{(1+\alpha/2)}}\left.\frac{\mathrm{d}r}{\mathrm{d}t}\right|_{\mathrm{c}}-u\left.\frac{\mathrm{d}u}{\mathrm{d}t}\right|_{\mathrm{c}}\\\
=-c\left.\frac{\mathrm{d}u}{\mathrm{d}t}\right|_{c}-\frac{\alpha
uc^{2}}{r}+\frac{\partial h_{\infty}}{\partial t},$ (74)

the derivative of the flow velocity $u$ along an outgoing characteristic is
given as

$\left.\frac{\mathrm{d}u}{\mathrm{d}t}\right|_{c}=\frac{\alpha}{r(c-u)}\bigg{[}\frac{g}{2\,r^{\alpha/2}}\left.\frac{\mathrm{d}r}{\mathrm{d}t}\right|_{\mathrm{c}}-uc^{2}\bigg{]},$
(75)

or, by inserting Eq. (33),

$\left.\frac{\mathrm{d}u}{\mathrm{d}t}\right|_{c}=\frac{\alpha}{r(c-u)}\bigg{[}\frac{g}{2\,r^{\alpha/2}}\,(c+u)-uc^{2}\bigg{]}.$
(76)

In spherical symmetry ($\alpha=2$), Eq. (76) corresponds to the ODE first
presented by Hickling and Plesset (1963),

$\left.\frac{\mathrm{d}u}{\mathrm{d}t}\right|_{c}=\frac{1}{r(c-u)}\left[\frac{g}{r}\,(c+u)-2uc^{2}\right].$
(77)

The derivation of this ODE for the flow velocity in the way presented here
follows the derivation presented for the special case of spherical symmetry
($\alpha=2$) by Ivany (1965).

An ODE for the flow velocity with respect to the radial coordinate may be
obtained by dividing Eq. (76) by Eq. (33),

$\left.\frac{\mathrm{d}u}{\mathrm{d}r}\right|_{c}=\left.\frac{\mathrm{d}u}{\mathrm{d}t}\right|_{c}\left.\frac{\mathrm{d}r}{\mathrm{d}t}\right|_{c}^{-1}=\frac{\alpha}{r(c-u)}\bigg{[}\frac{g}{2\,r^{\alpha/2}}-\frac{uc^{2}}{c+u}\bigg{]},$
(78)

an ODE first proposed by Gilmore (1952) for spherical symmetry ($\alpha=2$)
using the modified Tait EoS.

###  VI.3 Liquid pressure

With an EoS for the liquid at hand, the pressure of the liquid can be computed
based on the specific enthalpy $h$. Using the modified Tait EoS, the pressure
of the liquid follows by inserting Eq. (44) into Eq. (43) to obtain

$p=\left[\frac{(n-1)\rho_{0}}{n\left(p_{0}+B\right)^{1/n}}\,h\right]^{\frac{1}{1-1/n}}-B$
(79)

or, using Eqs. (43) and (61), the pressure can be computed based on the
invariant $g$,

$p=\rho\left[\frac{p_{\infty}+B}{\rho_{\infty}}+\frac{n-1}{n}\left(\frac{g}{r^{\alpha/2}}-\frac{u^{2}}{2}\right)\right]-B,$
(80)

as introduced for spherical symmetry ($\alpha=2$) by Ivany (1965) and, later,
Ebeling (1978). The pressure amplitude $\Delta p=p-p_{\infty}$ is dominated by
the term $g/r^{\alpha/2}$, especially in the far field where $u\rightarrow 0$,
thus recovering the well-known amplitude decay with $1/r$ in spherical
symmetry ($\alpha=2$) and $1/\sqrt{r}$ in cylindrical symmetry
($\alpha=1$).(Landau and Lifshitz, 1959) With the NASG EoS, the liquid
pressure follows from Eqs. (47) and (61) in a similar manner as

$p=\frac{(n-1)\rho h-(1-b\rho)nB}{n-b\rho}.$ (81)

Since $p$ and $\rho$ depend explicitly on each other, Eqs. (80) and (81) may
require an iterative evaluation (Denner and Schenke, 2023a), in order to reach
a desired accuracy.

Alternatively, inserting Eq. (76) into Eq. (73), and rearranging for the
enthalpy derivative yields

$\left.\frac{\text{d}h}{\text{d}t}\right|_{\text{c}}=\frac{\alpha}{r(c-u)}\left[u^{2}c^{2}-\frac{g}{2r^{\alpha/2}}(c+u)c\right]+\frac{\partial
h_{\infty}}{\partial t}.$ (82)

By applying Eq. (10), the ODE for the pressure along outgoing characteristics
follows as

$\left.\frac{\text{d}p}{\text{d}t}\right|_{\text{c}}=\frac{\rho\alpha}{r(c-u)}\left[u^{2}c^{2}-\frac{g}{2r^{\alpha/2}}(c+u)c\right]+\frac{\partial
p_{\infty}}{\partial t},$ (83)

which may be solved together with the ODEs for the radial position, Eq. (33),
and the velocity, Eq. (76). Assuming spherical symmetry in conjunction with
the modified Tait EoS, and neglecting ambient pressure changes, Eq. (83)
simplifies to

$\left.\frac{\text{d}p}{\text{d}t}\right|_{\text{c}}=\frac{\rho_{0}}{r(c-u)}\left(\frac{p+B}{p_{0}+B}\right)^{1/n}\left[2u^{2}c^{2}-\frac{g}{r}(c+u)c\right],$
(84)

which has been used frequently instead of Eq. (80).

##  VII Model hierarchy

KBH models for pressure-driven bubble dynamics and acoustics are part of a
consistent hierarchy of models (Naugol’nykh and Roy, 1971; Prosperetti, 1984;
Denner and Schenke, 2023a; Zhang  _et al._ , 2023b) that can be established by
successive simplifications of the Kirkwood-Bethe hypothesis.

Assuming $u^{2}\ll c_{0}^{2}$ and constant fluid properties ($\rho=\rho_{0}$,
$c=c_{0}$), such that the EoS of the liquid becomes obsolete, the potential
$\phi=r^{\alpha/2}\psi$ is governed by the wave equation

$\frac{1}{c_{0}^{2}}\frac{\partial^{2}\phi}{\partial t^{2}}-\nabla^{2}\phi=0,$
(85)

the propagation speed along outgoing characteristics is constant,

$\left.\frac{\text{d}r}{\text{d}t}\right|_{\text{c}}=c_{0},$ (86)

and the specific enthalpy is defined by the first-order approximation
$h=p/\rho_{0}$. The pressure, consequently, follows from Eq. (61) as

$p=p_{\infty}+\rho_{0}\left(\frac{g}{r^{\alpha/2}}-\frac{u^{2}}{2}\right),$
(87)

which neglects terms proportional to $c_{0}^{-2}$ and higher. Considering Eq.
(86), the explicit expression for the flow velocity, Eq. (67), simplifies to

$u=\frac{\alpha\,\phi}{2\,r^{1+\alpha/2}}+\frac{g}{r^{\alpha/2}c_{0}}.$ (88)

Based on Eqs. (87) and (88), the invariants $g$ and $\phi$ are defined at the
gas-liquid interface ($r=R$), such that the boundary conditions are satisfied,
as

$\displaystyle g$ $\displaystyle=$ $\displaystyle
R^{\alpha/2}\left(\frac{p_{\text{L}}-p_{\infty}}{\rho_{0}}+\frac{\dot{R}^{2}}{2}\right)$
(89) $\displaystyle\phi$ $\displaystyle=$
$\displaystyle\frac{2}{\alpha}\left(R^{1+\alpha/2}\dot{R}-R\,\frac{g}{c_{0}}\right).$
(90)

This presents, in fact, the quasi-acoustic model of Gilmore (1952) and
Trilling (1952), generalized with respect to the dimensionality coefficient
$\alpha$. Although this quasi-acoustic model incorporates a finite wave speed,
it cannot describe the nonlinear distortion of acoustic waves and the
formation of shock fronts, since all parts of a wave propagate with the same
constant speed $c_{0}$. Applying the same simplifications to the ODE for the
radial bubble dynamics, Eq. (54), and neglecting the term proportional to
$c_{0}^{-2}$, yields a generalized form of the Keller-Miksis equation,(Keller
and Miksis, 1980)

$\left(1-\frac{\dot{R}}{c_{0}}\right)R\ddot{R}+\frac{3}{4}\alpha\left(1-\frac{\dot{R}}{3\,c_{0}}\right)\dot{R}^{2}\\\
=\frac{\alpha}{2}\left(1+\frac{\dot{R}}{c_{0}}\right)\frac{p_{\mathrm{L}}-p_{\infty}}{\rho_{0}}+R\,\frac{\dot{p}_{\mathrm{L}}-\dot{p}_{\infty}}{\rho_{0}\,c_{0}}.$
(91)

The quasi-acoustic model and the Keller-Miksis equation are first-order
accurate with respect to the Mach number and, thus, limited to small Mach
numbers, $(\dot{R}/c_{0})^{2}\ll 1$.

Going one step further by assuming the fluid is incompressible, with
$c_{0}\rightarrow\infty$, the velocity potential is governed by
$\nabla^{2}\phi=0$ and the propagation speed along outgoing characteristics is
infinite. The pressure in the liquid is given by Eq. (87) and the flow
velocity follows from Eq. (88) as

$u=\left(\frac{R}{r}\right)^{1+\alpha/2}\dot{R}.$ (92)

With $c_{0}\rightarrow\infty$, the ODE for the radial bubble dynamics, Eq.
(91), further reduces to a generalized form of the Rayleigh-Plesset equation
(Plesset and Prosperetti, 1977; Ilinskii  _et al._ , 2012; Wang  _et al._ ,
2024b),

$R\ddot{R}+\frac{3}{4}\alpha\dot{R}^{2}=\frac{\alpha(p_{\mathrm{L}}-p_{\infty})}{2\rho_{0}}.$
(93)

##  VIII Solution methods

The equations governing the bubble dynamics and the equations describing the
flow field can be solved independent of each other. From the viewpoint of the
bubble dynamics, this is possible because the Kirkwood-Bethe hypothesis
connects the kinematic conditions at the gas-liquid interface to the spatially
invariant ambient condition at an assumed infinite distance. The flow field,
in turn, is defined by discrete pieces of information that evolve along a
single spatial dimension with speed $c+u$ and are connected to their initial
conditions originating at the gas-liquid interface by the invariant $g$. The
solution methods presented in this section are implemented in the open-source
software library APECSS (Denner and Schenke, 2023b).

###  VIII.1 Bubble dynamics

In order to solve the equation of motion of the gas-liquid interface, Eq.
(54), we require expressions for the specific enthalpy difference $H$ and its
derivative $\dot{H}$. The specific enthalpy difference can be computed using a
suitable EoS in conjunction with the liquid pressure at the gas-liquid
interface $p_{\text{L}}$ and the ambient pressure $p_{\infty}$.

Considering a Newtonian liquid, $p_{\text{L}}$ is given by Eq. (56). Since
$\text{d}h=\text{d}p/\rho$ for the considered isentropic flow, the derivative
of the specific enthalpy difference is defined as

$\dot{H}=\frac{\dot{p}_{\text{L}}}{\rho_{\text{L}}}-\frac{\dot{p}_{\infty}}{\rho_{\infty}}=\frac{1}{\rho_{\text{L}}}\left(\dot{p}_{\text{G}}+\frac{\alpha\sigma}{R^{2}}\dot{R}\right.\\\
\left.+2\alpha\mu\frac{\dot{R}^{2}}{R^{2}}-2\alpha\mu\frac{\ddot{R}}{R}\right)-\frac{\dot{p}_{\infty}}{\rho_{\infty}}.$
(94)

Inserting Eq. (94) into Eq. (54) and gathering the $\ddot{R}$-terms on the
left-hand side yields

$\left(c_{\text{L}}R+\frac{2\alpha\mu}{\rho_{\text{L}}}\right)\ddot{R}=\frac{\alpha}{2}\left(\frac{c_{\text{L}}+{\dot{R}}}{c_{\text{L}}-{\dot{R}}}\,H-\frac{3c_{\text{L}}-{\dot{R}}}{c_{\text{L}}-{\dot{R}}}\,\frac{\dot{R}^{2}}{2}\right)+\\\
\frac{R}{\rho_{\text{L}}}\left(\dot{p}_{\text{G}}+\frac{\alpha\sigma}{R^{2}}\dot{R}+2\alpha\mu\frac{\dot{R}^{2}}{R^{2}}\right)-\frac{R\dot{p}_{\infty}}{\rho_{\infty}}.$
(95)

We can solve this second-order ODE by separating it into two first-order ODEs
for $\dot{R}$ and $\dot{U}=\ddot{R}$, such that at each time instance $t_{j}$
we integrate the equations

$\displaystyle\dot{R}_{j+1}=U_{j}$ (96)

and

$\dot{U}_{j+1}=\frac{\alpha}{2\mathcal{A}_{j}}\left(\frac{c_{\text{L},j}+{\dot{R}_{j}}}{c_{\text{L},j}-{\dot{R}_{j}}}H_{j}-\frac{3c_{\text{L},j}-{\dot{R}_{j}}}{c_{\text{L},j}-{\dot{R}_{j}}}\,\frac{\dot{R}_{j}^{2}}{2}\right)\\\
+\frac{R_{j}}{\mathcal{A}_{j}}\left(\frac{\dot{p}_{\text{G},j}}{\rho_{\text{L},j}}+\frac{\alpha\sigma\dot{R}_{j}}{R_{j}^{2}\rho_{\text{L},j}}+\frac{2\alpha\mu\dot{R}_{j}^{2}}{R^{2}_{j}\rho_{\text{L},j}}-\frac{\dot{p}_{\infty,j}}{\rho_{\infty,j}}\right),$
(97)

where

$\mathcal{A}_{j}=c_{\text{L},j}R_{j}+\frac{2\alpha\mu}{\rho_{\text{L},j}}.$
(98)

The embedded fifth-order Runge-Kutta scheme of Dormand and Prince (1980),
which adapts the numerical time step $\Delta t_{j}=t_{j}-t_{j-1}$ based on the
estimated solution error, is a common choice to integrate this system of ODEs,
as it is, for instance, the basis of the ode45 function in Matlab, the default
solver of the solve_ivp function of the SciPy package, and used in the solver
of the open-source cavitation software library APECSS (Denner and Schenke,
2023b).

###  VIII.2 Flow field

To solve the flow velocity and specific enthalpy (or pressure), the
information originating at the gas-liquid interface is tracked along the
outgoing characteristics, whereby each parcel of information propagates on its
own characteristic with velocity $c+u$. To this end, the bubble radius
$R_{j}$, the speed of the gas-liquid interface $\dot{R}_{j}$, and the specific
enthalpy of the liquid at the gas-liquid interface $h_{\text{L},j}$ obtained
by solving Eqs. (96) and (97) serve as the initial conditions.

The invariant $g$ of information parcel $i$ is defined upon emission at the
gas-liquid interface, at emission time $\tau_{i}$, as

$g_{i}=R(\tau_{i})^{\alpha/2}\left[h_{\text{L}}(\tau_{i})-h_{\infty}(\tau_{i})+\frac{\dot{R}(\tau_{i})^{2}}{2}\right].$
(99)

The radial position $r_{i,j}$ and flow velocity $u_{i,j}$ of each information
parcel $i$ at time $t_{j}$ result from Eqs. (33) and (75), by integrating

$\displaystyle\left.\frac{\mathrm{d}r}{\mathrm{d}t}\right|_{\text{c}(i,j)}=c_{i,j}+u_{i,j}$
(100)
$\left.\frac{\mathrm{d}u}{\mathrm{d}t}\right|_{\text{c}(i,j)}=\frac{\alpha}{r_{i,j}(c_{i,j}-u_{i,j})}\Bigg{[}\frac{g_{i}}{2\,r_{i,j}^{\alpha/2}}\left.\frac{\mathrm{d}r}{\mathrm{d}t}\right|_{\mathrm{c}(i,j)}-u_{i,j}c_{i,j}^{2}\Bigg{]}$
(101)

using a suitable integration scheme for ODEs together with the initial
conditions $R(\tau_{i})$ and $\dot{R}(\tau_{i})$. Once the new radial position
$r_{i,j}$ and flow velocity $u_{i,j}$ of information parcel $i$ are
determined, its specific enthalpy is readily given by Eq. (61) as

$h_{i,j}=h_{\infty,j}+\frac{g_{i}}{r_{i,j}^{\alpha/2}}-\frac{u_{i,j}^{2}}{2},$
(102)

from which the local pressure and speed of sound can be computed using the
employed EoS.

As an alternative to integrating the velocity using Eq. (101), although not
used for the results presented in Section IX, the flow velocity may also be
determined explicitly using Eq. (69), such that

$u_{i,j}=\frac{\alpha\,\phi_{i}}{2\,r_{i,j}^{1+\alpha/2}}+\frac{g_{i}}{r_{i,j}^{\alpha/2}(c_{i,j}+u_{i,j})},$
(103)

with

$\displaystyle\phi_{i}=\frac{2}{\alpha}\left[R(\tau_{i})^{1+\alpha/2}\dot{R}(\tau_{i})-R(\tau_{i})\,\frac{g_{i}}{c_{\text{L}}(\tau_{i})+\dot{R}(\tau_{i})}\right]$
(104)

if sufficiently small flow velocities are expected.

Ivany (1965) proposed to solve the flow velocity and pressure along outgoing
characteristics originating at the gas-liquid interface at selected time
instances, after the solution of the bubble dynamics had been fully obtained.
Denner and Schenke (2023a) presented a Lagrangian wave tracking algorithm,
whereby the flow field is solved together with the bubble dynamics and the
information originating from the gas-liquid interface is tracked along
outgoing characteristics. To this end, Eqs. (100) and (101) are integrated
together as a coupled system of ODEs by a conventional fourth-order Runge-
Kutta scheme, using the same time step $\Delta t_{j}$ as for the integration
of Eqs. (96) and (97) (Denner and Schenke, 2023b).

###  VIII.3 Shock fronts

Because both the particle velocity and the speed of sound are monotone-
increasing functions of the pressure (Fay, 1931; Hamilton and Blackstock,
1988), the parts of an acoustic wave with a higher pressure propagate faster
than the parts of the same wave with a lower pressure. This pressure
dependency of the propagation speed drives a progressive steepening of the
emitted waves and, eventually, may lead to the formation of a shock front for
waves with sufficiently large initial amplitude. These cumulative nonlinear
effects may become dominant for $Mkr\gtrsim 1$, where $k$ is the wavenumber of
the emitted wave, irrespective of the Mach number $M$; (Naugol’nykh, 1971) in
the absence of dissipation, any plane wave forms a shock wave in finite time.
At the shock front, a frequency-dependent attenuation prevents fluid particles
from overtaking the shock front (Fay, 1931) and, as postulated by Rudnick
(1952), the stabilization of the shock front is independent of the origin of
the acting dissipation process, such as viscous stresses or thermal
conduction. This can be observed in numerical simulations, where a stable
shock front can form even in inviscid and non-conducting flows, if sufficient
numerical dissipation is present (Laney, 1998; LeVeque, 2008; Schenke  _et
al._ , 2022).

Figure 4: Example of the different treatments of a shock front forming after
the collapse of a cavitation bubble produced by a femtosecond laser, as
studied by Liang  _et al._ (2022). (a) Selected outgoing wave characteristics
are shown in a time-space plot, where the characteristic associated with the
minimum radius is seen to intersect previously emitted characteristics. The
thick black line shows the bubble radius. (b) The resulting multivalued
density profile $50\,\text{ns}$ after the bubble attained its minimum radius.
(c) Close up of the density profile, to which two different shock treatments
are applied. The clipping of the profile as proposed by Denner and Schenke
(2023a) is applied at run time, whereas the rule of equal areas (Landau and
Lifshitz, 1959) is applied as a post-processing step.

When solving the flow field by propagating the relevant information along
outgoing characteristics using a KBH model, the formation of a shock front is
signified by the intersection of characteristics, as seen in Figure 4(a). This
intersection of characteristics yields a multivalued solution, shown in Figure
4(b), that is not physical. Akulichev  _et al._ (1968) applied the rule of
equal areas (Landau and Lifshitz, 1959) to post process the multivalued
solutions resulting from intersecting characteristics. Taking a multivalued
density profile as shown in Figure 4(b), the corresponding weak solution is
obtained graphically by eliminating the multivalued solution in such a way,
that the area under the curve remains unchanged, as illustrated in Figure
4(c). Thus, the integral of the multivalued solution and the integral of the
discontinuous weak solution of the density profile are the same, such that
mass is conserved (Whitham, 1999). This implies that the wave on either side
of the shock front remains a simple wave traveling in one direction, without
reflections introduced at the discontinuity (Naugol’nykh, 1971). This method
has been used, for instance, in studies on laser-induced cavitation (Vogel,
Busch, and Parlitz, 1996; Lai  _et al._ , 2022; Liang  _et al._ , 2022; Wen
_et al._ , 2023) and sonoluminescence (Holzfuss, 2010) to model the formation
and attenuation of shock waves using the Kirkwood-Bethe hypothesis.

Figure 5: Sketch of a discrete wave profile (left) with a multivalued solution
and (right) after treating this multivalued solution by clipping the wave
profile, as proposed by Denner and Schenke (2023a).

Utilizing the Lagrangian wave tracking algorithm mentioned in the previous
section, Denner and Schenke (2023a) proposed to simply discard any information
that overtakes the information propagating ahead of it, as illustrated in
Figure 5, to account for the formation of shock fronts and the associated
dissipation. Naugol’nykh (1971) described a similar procedure as “clipping of
the vertices of its profile” when referring to the decay of plane sawtooth
waves. This rather simple procedure aims to exploit the argument of Rudnick
(1952) that the stabilization of the shock front is independent of the origin
of dissipation. However, even though the results obtained with this method are
in reasonably good agreement with the result obtained by applying the rule of
equal areas for shock waves produced by collapsing bubbles (Denner and
Schenke, 2023a), discarding information by clipping the density profile does
not redistribute mass. Hence, this simple shock treatment yields an error in
mass conservation and overpredicts the attenuation of the shock front, as
observed in Figure 4(c).

##  IX Further investigation

KBH models have been tested and validated extensively (Oh  _et al._ , 2018;
Geng  _et al._ , 2021; Vassholz  _et al._ , 2021; Lai  _et al._ , 2022; Wen
_et al._ , 2023; Li  _et al._ , 2023; Zhang  _et al._ , 2021b; Wang  _et al._
, 2021; Sheng  _et al._ , 2023; Denner and Schenke, 2023a), utilizing high-
fidelity experiments and fully resolved simulations, and have been found to
produce accurate and reliable results in a variety of scenarios. However, some
fundamental questions regarding the validity and applicability of the
Kirkwood-Bethe hypothesis remain. To further scrutinize and improve KBH
models, this section focuses on the prediction of the propagation and
attenuation of shock fronts at simulation run time (Section IX.1) and the
complex impedance of curved acoustic waves (Section IX.2). The presented
results were produced with version v1.6 of the open-source software library
APECSS (Denner and Schenke, 2023b), available at
https://doi.org/10.5281/zenodo.10981878.

###  IX.1 Predicting shock fronts at simulation run time

Figure 6: Sketch of a discrete wave profile (left) with a multivalued solution
and (right) after treating this multivalued solution by the proposed averaging
procedure, Eqs. (105)-(107). Figure 7: Pressure amplitude $\Delta p(r,t)$,
normalized by the excitation pressure amplitude $\Delta p_{\text{a}}$,
generated by an oscillating planar emitter as a function of the dimensionless
distance $(r-R_{0})/r_{\text{sh}}$, where $R_{0}$ is the initial size of the
emitter and $r_{\text{sh}}$ is the shock formation distance of a plane
harmonic wave in an inviscid fluid, Eq. (109). The analytical solution for the
amplitude of the decaying shock wave of Fay (1931), Eq. (108), is shown for
reference.

The robust and accurate treatment of multivalued solutions associated with
shock fronts at simulation run time presents challenges. Firstly, the full
wave profile is required for an accurate treatment of a multivalued solution.
This can be addressed when solving the flow along the outgoing characteristics
by reconstructing the pressure or density profile in space. Secondly, however,
the multivalued solution should also be appropriately resolved in space.
Assuming that we follow the discrete pieces of information that originate at
the gas-liquid interface along the outgoing characteristic, e.g. as described
in Section VIII.2, the spatial resolution of the wave profile is directly
dependent on the time step applied to solve the bubble dynamics.

While the application of the rule of equal areas provides an accurate means to
rectify multivalued solutions, it is usually applied as a post-processing
step. The clipping method of Denner and Schenke (2023a) can treat multivalued
solutions at run time, but does not conserve mass and tends to overestimate
the attenuation at shock fronts. To remedy this shortcoming, a run-time shock
treatment based on a simple averaging procedure of the multivalued solution is
proposed, as illustrated in Figure 6. With the overtaking information parcel
$i+1$ tagged to be deleted, the location $r_{i}$ and invariant $g_{i}$ of
information parcel $i$ are redefined as

$\displaystyle r_{i,j}$ $\displaystyle=$
$\displaystyle\frac{r_{i+1,j}+r_{i,j}}{2}$ (105) $\displaystyle g_{i}$
$\displaystyle=$ $\displaystyle\frac{g_{i+1}+g_{i}}{2}.$ (106)

The average of the invariant $g$ follows from the conservation of specific
kinetic enthalpy and presumes that the applied time step is sufficiently small
so that the invariant $g$ varies linearly in the close vicinity of the shock
front. This leaves the enthalpy $h_{i}$ and the velocity $u_{i}$ undefined.
Considering that the velocity has two solutions that satisfy a given $g$, the
robust choice is to average the velocity as

$u_{i,j}=\frac{u_{i+1,j}+u_{i,j}}{2}.$ (107)

The specific enthalpy $h_{i}$ is then readily defined by Eq. (102).
Furthermore, preliminary tests showed that the choice of quantities that are
averaged (i.e., $g_{i}$ and $u_{i}$, $h_{i}$ and $u_{i}$, or $\rho_{i}$ and
$u_{i}$) has a negligible impact on the results. After the quantities
associated with information parcel $i$ have been redefined, the overtaking
information parcel $i+1$ is discarded.

As a first test case, we consider the shock formation of a planar wave emitted
in an inviscid fluid by an oscillating planar emitter ($\alpha=0$). The
oscillating emitter generates an harmonic acoustic wave that gradually
steepens to eventually form a shock front, at which point it starts to develop
into a sawtooth wave that decays due to the attenuation acting at the shock
front. For such a sawtooth wave, Fay (1931) derived an analytical solution for
the pressure amplitude $\Delta p$, given in function of the one-dimensional
coordinate $r$ as

$\Delta p(r)=\pm\frac{\pi r_{\text{sh}}}{r+r_{\text{sh}}}\,\Delta
p_{\text{a}},$ (108)

where the shock forming distance of a plane harmonic wave is (Blackstock,
1966)

$r_{\text{sh}}=\frac{\rho_{0}c_{0}^{3}}{2\pi\beta f_{\text{a}}\Delta
p_{\text{a}}},$ (109)

where $f_{\text{a}}$ and $\Delta p_{\text{a}}$ are the initial frequency and
pressure amplitude of the emitted acoustic wave. The nonlinearity coefficient
$\beta$ for an ideal gas or a Tait fluid is defined as (Hamilton and
Blackstock, 1988; Lauterborn, Kurz, and Akhatov, 2007) $\beta=(\gamma+1)/2$.
Note that the solution of Fay, Eq. (108), is only applicable to “perfect”
sawtooth waves (Blackstock, 1966) for $r\gg r_{\text{sh}}$ and that the
excitation amplitude $\Delta p_{\text{a}}$ ought to be sufficiently small such
that the nonlinear dependency of the fluid properties on the pressure (see
Section IV) is negligible. Figure 7 shows the results obtained using Eqs.
(100)-(102) in conjunction with the clipping-based shock treatment of Denner
and Schenke (2023a) and the averaging-based shock treatment proposed above,
for an initially harmonic wave in water modeled using the modified Tait EoS
($\beta=4.075$). At $r=r_{\text{sh}}+R_{0}$, the wave profile forms vertical
fronts, where $\partial p/\partial r\rightarrow\infty$, as predicted by the
theory of Blackstock (1966). Subsequently, the wave evolves into a fully
developed sawtooth wave and starts to decay, with its amplitude asymptotically
approaching the solution of Fay (1931). The clipping procedure of Denner and
Schenke (2023a), whereby a multivalued solution is treated by clipping the
“overtaking” part of the wave profile, fails to predict this sawtooth wave
correctly; only the wave peaks reduce, whereas the wave troughs remain
unchanged. This shortcoming is evidently less of an issue for the pressure
pulses emitted by a strong bubble collapse for which the clipping procedure
has previously been tested (Denner and Schenke, 2023a). The proposed averaging
procedure is able to predict the decay of both the peaks and troughs of the
sawtooth wave accurately.

Figure 8: Radius $R(t)$, normalized by the initial radius $R_{0}$, and gas
pressure $p_{\text{G}}(t)$, normalized by the initial gas pressure
$p_{\text{G,0}}$, as a function of time $t$, normalized by the Rayleigh
collapse time $t_{\text{c}}$, of the Rayleigh collapse of a spherical gas
bubble in water, obtained with a Navier-Stokes solver and the Gilmore
equation. Figure 9: Temporal evolution of the pressure amplitude $\Delta
p(r,t)=p(r,t)-p_{\infty}$, normalized by the constant ambient pressure
$p_{\infty}$, and flow velocity $u(r,t)$, normalized by the characteristic
collapse velocity $u_{\text{c}}=R_{0}/t_{\text{c}}$, emitted by the Rayleigh
collapse of a spherical gas bubble in water at different radial locations
$r=\\{0.2,0.5,1\\}R_{0}$, obtained with the KBH model described in Section
VIII.2. Time $t$ is normalized by the Rayleigh collapse time $t_{\text{c}}$.
The multivalued solutions at the shock front are treated either by clipping
the wave profile or by applying the proposed averaging procedure. Results
obtained with a Navier-Stokes solver are shown for reference. Figure 10:
Temporal evolution of the radius $R(t)$, bubble wall Mach number
$\dot{R}(t)/c_{\text{L}}(t)$, and gas pressure $p_{\text{G}}(t)$ of the
sonoluminescence bubble of Holzfuss (2010), obtained with the Gilmore model.
Figure 11: Spatial evolution of the pressure amplitude $\Delta
p(r,t)=p(r,t)-p_{\infty}(t)$ and the flow velocity $u(r,t)$ of the acoustic
wave emitted by the primary collapse of the sonoluminescence bubble of
Holzfuss (2010), obtained with the KBH model described in Section VIII.2. The
multivalued solutions at the shock front are treated either by clipping the
wave profile or by applying the proposed averaging procedure. The expected
pressure and velocity decay in the far field, proportional to $r^{-1}$, is
indicated by the gray dashed lines.

To assess the capabilities of the run-time shock treatment in a more
representative scenario, the Rayleigh collapse of a spherical gas bubble
($\alpha=2$) in water is simulated. The bubble has an initial radius $R_{0}$
and its collapse is driven by an initial gas pressure of
$p_{\mathrm{G,0}}=p_{\infty}/1000$ compared to a constant ambient pressure of
$p_{\infty}=10^{5}\,\mathrm{Pa}$. The gas is modeled using the NA EoS, see
Eqs. (59) and (60), with $\rho_{\mathrm{G},0}=1.2\,\mathrm{kg/m}^{3}$ at
$p_{\infty}$, $\kappa=1.4$, and $b=1.5\times
10^{-3}\,\mathrm{m}^{3}/\mathrm{kg}$, whereas water is modeled using the
modified Tait EoS, with $n=7.15$, $B=3.046\,\times 10^{8}\,\text{Pa}$,
$\rho_{0}=997\,\text{kg/m}^{3}$ at $p_{\infty}$, and $\mu=0.001\,\mathrm{Pa\
s}$. The surface tension of the gas-liquid interface is neglected. For
validation purposes, the results obtained based on the Kirkwood-Bethe
hypothesis are compared to a fully resolved simulation using a state-of-the-
art Navier-Stokes solver (Denner, Xiao, and van Wachem, 2018; Denner, Evrard,
and van Wachem, 2020). For this simulation, the pressure in the liquid is
initialized as $p(r)=p_{\infty}+R_{0}(p_{\mathrm{G,0}}-p_{\infty})/r$, the
bubble is resolved with $1000$ mesh cells per initial bubble radius, and the
time step applied at the time when the bubble reaches its minimum radius is
$\Delta t\simeq 5\times 10^{-8}t_{\mathrm{c}}$, where $t_{\mathrm{c}}\simeq
0.915\,R_{0}\sqrt{\rho_{\infty}/p_{\infty}}$ is the Rayleigh collapse time
(Rayleigh, 1917). Figure 8 shows a close agreement of the bubble radius and
gas pressure obtained with the Navier-Stokes solver and with the Gilmore
equation, Eq. (55). The temporal evolution of the pressure and velocity in the
vicinity of the bubble is shown in Figure 9, at $r=\\{0.2,0.5,1\\}R_{0}$. The
emitted acoustic wave forms a shock front, with a vertical left slope of the
pressure and velocity profiles. Note that in the temporal-evolution plots
shown in Figure 9, the left-hand side of the profile reaches the measurement
location first. The results obtained with both methods to treat the shock
front, either clipping the wave profile or applying the proposed averaging
procedure, yield a good agreement with the Navier-Stokes solution for this
specific test case. However, clipping the wave profile overpredicts the
attenuation of the shock front in the far field, which manifests in a smaller
peak pressure and velocity at $r=R_{0}$.

Figure 12: Temporal evolution of the radius $R(t)$, bubble wall Mach number
$\dot{R}(t)/c_{\text{L}}(t)$, and gas pressure $p_{\text{G}}(t)$ of the laser-
induced bubble of Liang  _et al._ (2022), obtained with the Gilmore model.
Figure 13: Wave profiles of the pressure amplitude $\Delta
p(r,t)=p(r,t)-p_{\infty}(t)$ and the flow velocity $u(r,t)$ at selected time
instances of the acoustic wave emitted by the first collapse of the laser-
induced bubble of Liang  _et al._ (2022), obtained with the KBH model
described in Section VIII.2. The solid line shows the results obtained by
applying the proposed averaging procedure at the shock front, the dotted line
shows the results obtained by clipping the wave profile, and the colored dots
show the corresponding results of Liang  _et al._ (2022). Time $t_{0}$ is the
time at which the peak pressure is emitted.

The strong collapse and shock formation of an argon bubble that produces
sonoluminescence (Holzfuss, 2010) and of a laser-induced bubble (Liang  _et
al._ , 2022) have previously been studied with respect to the emission of
shock waves using KBH models. In both studies (Holzfuss, 2010; Liang  _et al._
, 2022), the rule of equal areas was applied to rectify the multivalued
solutions at the shock front. The sonoluminescence bubble considered by
Holzfuss (2010) has an initial radius of $R_{0}=5\,\mu\text{m}$ and is excited
with a frequency of $f_{\text{a}}=23.5\,\text{kHz}$ and a pressure amplitude
of $\Delta p_{\text{a}}=145\,\text{kPa}$, whereas the laser-induced bubble
studied by Liang  _et al._ (2022) is generated by a laser pulse with
$20\,\mu\text{J}$ energy and $265\,\text{fs}$ duration, and reaches a maximum
bubble radius of approximately $36\,\mu\text{m}$. In both cases, water is
modeled by the modified Tait EoS. The exact parameters of these cases can be
found in the respective publications, or the recent work of Denner and Schenke
(2023a). For orientation, the radius $R$, bubble wall Mach number
$\dot{R}/c_{\text{L}}$, and gas pressure $p_{\text{G}}$ of the
sonoluminescence and laser-induced bubbles are shown as a function of time in
Figures 10 and 12, respectively. Note that the discontinuous changes in gas
pressure of the laser-induced bubble are intentional and account for the
complex mass-transfer dynamics of this bubble, as proposed by Liang  _et al._
(2022). In Figure 11, the pressure and velocity amplitude of the acoustic wave
emitted by the primary collapse of the sonoluminescence bubble are shown,
where subtle differences can be observed dependent on the treatment of the
developing shock front. The proposed averaging procedure exhibits a tendency
to larger pressure and velocity values compared to simply clipping the wave
profile. Nevertheless, the results obtained with both methods to treat the
shock front at run time compare very well with the results of Holzfuss (2010).
The spatial profiles of the acoustic wave produced by the collapse of the
laser-induced bubble, see Figure 13, exhibit a comparable trend: the proposed
averaging procedure yields a faster shock front with larger pressure and
velocity values compared to both the results of Liang  _et al._ (2022) and the
results obtained by clipping the wave profile. The observed differences are
confined to the immediate vicinity of the shock front and are largest in the
region in which the shock dissipation is strongest.

###  IX.2 Complex impedance of curved acoustic waves

The acoustic impedance experienced by a curved acoustic wave is complex, see
Eq. (37), and depends on the radius of curvature of the acoustic wave relative
to its wavenumber. Only for relatively short waves ($kr\gg 1$) does the simple
definition $z=\rho c$ of the specific acoustic impedance hold in good
approximation.

As discussed in Section VI.1, the explicit velocity expression given by Eq.
(69) is the sum of the displacement of the fluid by the moving gas-liquid
interface and the acoustic particle velocity associated with the emitted
acoustic wave. According to Eq. (69), the acoustic particle velocity $u_{1}$
is, thus, given as

$u_{1}(r,t)=\left[\frac{R(\tau)}{r(t)}\right]^{\alpha/2}\frac{h_{\text{L}}(\tau)-h_{\infty}(\tau)}{c(r,t)+u(r,t)}.$
(110)

Applying the isentropic relation $\text{d}h=\text{d}p/\rho$ and assuming the
acoustic wave has a small pressure amplitude, such that $u\ll c$,
$\rho=\rho(r,t)=\rho_{\text{L}}(\tau)$ and $c=c(r,t)=c_{\text{L}}(\tau)$, the
velocity is

$u_{1}(r,t)=\left[\frac{R(\tau)}{r(t)}\right]^{\alpha/2}\frac{p_{\text{L}}(\tau)-p_{\infty}(\tau)}{\rho
c},$ (111)

or, using the notation used in the acoustic expansions, with $p=p_{0}+p_{1}$,

$u_{1}(r,t)=\left[\frac{R(\tau)}{r(t)}\right]^{\alpha/2}\frac{p_{1}(\tau)}{\rho
c}.$ (112)

This expression is evidently not dependent on the relative curvature of the
emitted wave and, therefore, suggests that velocity expressions derived from
the Kirkwood-Bethe hypothesis do not account for the complex acoustic
impedance of curved waves.

Figure 14: The pressure (left) and velocity (right) of acoustic waves emitted
by a pulsating spherical emitter with (a) $k_{\text{a}}R_{0}=0.01$, (b)
$k_{\text{a}}R_{0}=1$, and (c) $k_{\text{a}}R_{0}=100$, predicted by the
Kirkwood-Bethe hypothesis. The analytical solution for curved acoustic waves,
Eqs. (115) and (117), and the particle velocity for plane acoustic waves are
shown for reference.

To test this, we consider a spherical emitter situated in water that
oscillates and emits acoustic waves at angular frequency $\omega_{\text{a}}$,
chosen such that a specific dimensionless wavenumber $k_{\text{a}}R_{0}$ is
obtained, with $k_{\text{a}}=\omega_{\text{a}}/c$. The pressure at the emitter
wall and the time-dependent velocity of the emitter wall are defined as

$\displaystyle p_{\text{L}}(t)$ $\displaystyle=$ $\displaystyle p_{0}+\Delta
p_{\text{a}}\,\cos\left(\omega_{\text{a}}t-\frac{\pi}{2}\right)$ (113)
$\displaystyle\dot{R}(t)$ $\displaystyle=$ $\displaystyle\frac{\Delta
p_{\text{a}}}{\rho c}\,\cos\left(\omega_{\text{a}}t-\frac{\pi}{2}\right),$
(114)

respectively. The considered pressure amplitude $\Delta p_{\text{a}}$ is
sufficiently small such that the pressure-dependency of the liquid density and
speed of sound are negligible. At the emitter wall, pressure and velocity are,
therefore, in phase and there is no relative velocity between the emitter wall
and the fluid. The water surrounding the emitter is described using the
modified Tait EoS, with $n=7.15$, $B=3.046\times 10^{8}\,\text{Pa}$, and
$\rho_{0}=997\,\text{kg/m}^{3}$. Following Eq. (37), the analytical solutions
for the acoustic pressure and particle velocity are given as(Kinsler, 2000)

$p_{1}(r,t)=\frac{R}{r}\,\Delta
p_{\text{a}}\,\cos\left(\omega_{\text{a}}t-\frac{\pi}{2}-k(r-R)\right)$ (115)
$u_{1}(r,t)=\frac{R}{r}\,\frac{\Delta p_{\text{a}}}{\rho
c\cos(\theta)}\,\cos\left(\omega_{\text{a}}t-\frac{\pi}{2}-k(r-R)-\theta\right),$
(116)

respectively, where the complex impedance manifests in the phase shift
$\theta=\pi/2-\arctan(kr)$. Taking the oscillating emitter into account by
inserting the analytical particle velocity defined by Eq. (116) into the flow
velocity given by Eq. (67), the analytical solution for the flow velocity is

$u(r,t)=\frac{R^{2}}{r^{2}}\dot{R}(t)-\frac{R}{r}u_{1}(r,t)+u_{1}(r,t).$ (117)

Figure 14 shows the spatial profiles of the pressure and velocity generated by
the spherical emitter ($\alpha=2$) for $k_{\text{a}}R_{0}\in\\{0.01,1,100\\}$,
predicted by the Kirkwood-Bethe hypothesis, using Eqs. (100)-(102). For all
three cases, the pressure predicted by the Kirkwood-Bethe hypothesis is in
virtually perfect agreement with the analytical solution, Eq. (115). However,
the velocity profile predicted by the Kirkwood-Bethe hypothesis agrees with
the analytical solution for a curved wave, Eq. (117), only for the case with
$k_{\text{a}}R_{0}=10^{2}$, i.e. for a short wave with $kr\gg 1$. As the
dimensionless wavenumber decreases, i.e. longer waves, the accuracy of the
velocity predicted by the Kirkwood-Bethe hypothesis deteriorates near the
emitter. Despite this shortcoming, the Kirkwood-Bethe hypothesis predicts the
velocity accurately after the wave has traveled for more than approximately
one wavelength ${\lambda_{\text{a}}}$, which corresponds to
$k_{\text{a}}(r-R_{0})\gtrsim 2\pi$. We can, therefore, conclude that the
Kirkwood-Bethe hypothesis does not account for the complex acoustic impedance
of curved waves and, as a consequence, fails to predict the particle velocity
of long waves ($kr\lesssim 1$) close to the emitter correctly. In fact, the
velocity ODE arising from the Kirkwood-Bethe hypothesis, Eq. (101), predicts
$u_{1}\simeq p_{1}/(\rho c)$ for all three cases, as observed in Figure 14,
irrespective of the dimensionless wavenumber $kr$. Nevertheless, the
predictions provided by the Kirkwood-Bethe hypothesis are accurate for
$r\gtrsim R+\lambda_{\text{a}}$.

Recent experimental and numerical studies that reported accurate predictions
of KBH models stand in apparent contradiction to the shortcoming associated
with the complex nature of the acoustic impedance. For instance, the
experimental measurements of the pressure and velocity of shock waves emitted
by laser-induced cavitation bubbles of Lai  _et al._ (2022) and the fully
resolved numerical results of the Rayleigh collapse of a spherical gas bubble
shown in Figure 9 are in very good agreement with the corresponding
predictions provided on the basis of the Kirkwood-Bethe hypothesis. However,
in these cases, which are representative for pressure-driven bubble dynamics,
including underwater explosions as well as medical and sonochemical
applications, the length of the emitted pressure pulse is short, with
$kR=\mathcal{O}(1)$ or larger Lai  _et al._ (2022); Denner and Schenke
(2023a); Wen  _et al._ (2023); Liang  _et al._ (2022). In addition, the
velocity very close to the bubble is dominated by the motion of the gas-liquid
interface and even though the particle velocity is nominally large, it is
typically less than $10\%$ of the total propagation velocity $c+u$.

##  X Conclusions

Originally developed to estimate the peak pressure of an underwater explosion
at some distance, the Kirkwood-Bethe hypothesis has given rise to a collection
of models that describe pressure-driven bubble dynamics, cavitation and
underwater explosions, as well as their acoustic emissions, with remarkable
accuracy. By using standard numerical integration methods and modern
computational equipment, these models allow to predict complex processes of
fluid dynamics and nonlinear acoustics, such as the laser-induced cavitation
bubble discussed in Section IX.1, in a fraction of a second. Here, the
theoretical foundations of the Kirkwood-Bethe hypothesis and contemporary
models derived from it have been reviewed, as well as generalized from the
typically considered spherical symmetry, to account for spherically symmetric,
cylindrically symmetric, and planar one-dimensional domains. In addition, a
new method to treat multivalued solutions associated with the formation of
shock fronts at simulation run time has been proposed, which further improves
the predictive capabilities of the models based on the Kirkwood-Bethe
hypothesis. Even though the Kirkwood-Bethe hypothesis fails to account for the
curvature of acoustic waves, as they are, for instance, emitted by a bubble
collapse, it has been shown to provide accurate predictions under the specific
conditions associated with pressure-driven bubble dynamics, cavitation and
underwater explosions.

Despite the extensive and formidable work achieved by the research community
over more than eight decades, models based on the Kirkwood-Bethe hypothesis
still lack some basic capabilities that can further extend their applicability
and improve their accuracy. For instance, although the Kirkwood-Bethe
hypothesis can produce accurate pressure and velocity predictions under
realistic scenarios, it does not account for heat transfer. As a first step,
the recent combination of the Kirkwood-Bethe hypothesis with the Noble-Abel
stiffened-gas equation of state for bubble dynamics (Denner, 2021) and
acoustic emissions (Denner and Schenke, 2023a) provides a consistent and
accurate basis for predicting the temperature in the liquid surrounding a gas
bubble or cavity. While adiabatic temperature changes can be readily predicted
in this way, thermal conduction and viscous heating, such as the heat
introduced by the energy dissipation at a shock front (Liang  _et al._ ,
2022), are not accounted for. In this respect, medical ultrasound stands out
as a field that requires high-fidelity temperature predictions if
computational tools are to play a more prominent role in the design of new
treatments, since temperature changes of only $1\,\text{K}$ can be safety
relevant (ter Haar, 2011). More generally, bubbles are neither stationary nor
do they appear in isolation or in a liquid of quasi-infinite extent in most
engineering applications, which contradicts some of the main assumptions
underpinning the available models based on the Kirkwood-Bethe hypothesis.
Recent advances for methods developed based on the (quasi-)acoustic assumption
point to possible solutions that may be incorporated into models based on the
Kirkwood-Bethe hypothesis, to account for the forces on bubbles as they move
through a liquid (Zhang  _et al._ , 2023b), the acoustic interaction of
multiple bubbles (Fan, Li, and Fuster, 2021), and the influence of confinement
and nearby walls (Fu  _et al._ , 2023).

## Conflict of interest

The author has no conflicts that could have influenced this work.

## Data availability

All presented results were produced, and can be readily reproduced, with the
open-source software library APECSS, as referenced in Section IX.

## References

## References

  * Plesset and Prosperetti (1977) M. S. Plesset and A. Prosperetti, “Bubble Dynamics and Cavitation,” Annual Review of Fluid Mechanics 9, 145–185 (1977).
  * Brennen (1995) C. E. Brennen,  _Cavitation and Bubble Dynamics_ , Oxford Engineering Science Series No. 44 (Oxford University Press, New York, 1995).
  * Lauterborn and Kurz (2010) W. Lauterborn and T. Kurz, “Physics of bubble oscillations,” Reports on Progress in Physics 73, 106501 (2010).
  * Blake and Gibson (1987) J. R. Blake and D. C. Gibson, “Cavitation Bubbles Near Boundaries,” Annual Review of Fluid Mechanics 19, 99–123 (1987).
  * Reuter, Deiter, and Ohl (2022) F. Reuter, C. Deiter, and C.-D. Ohl, “Cavitation erosion by shockwave self-focusing of a single bubble,” Ultrasonics Sonochemistry 90, 106131 (2022).
  * Dular and Ohl (2023) M. Dular and C. D. Ohl, “Bulk material influence on the aggressiveness of cavitation – Questioning the microjet impact influence and suggesting a possible way to erosion mitigation,” Wear , 205061 (2023).
  * Johansen  _et al._ (2017) K. Johansen, J. H. Song, K. Johnston, and P. Prentice, “Deconvolution of acoustically detected bubble-collapse shock waves,” Ultrasonics 73, 144–153 (2017).
  * Gaisser, Kirschner, and Riedelbauch (2023) L. Gaisser, O. Kirschner, and S. Riedelbauch, “Cavitation detection in hydraulic machinery by analyzing acoustic emissions under strong domain shifts using neural networks,” Physics of Fluids 35, 027128 (2023).
  * Haskell  _et al._ (2023) S. C. Haskell, N. Lu, G. E. Stocker, Z. Xu, and J. R. Sukovich, “Monitoring cavitation dynamics evolution in tissue mimicking hydrogels for repeated exposures via acoustic cavitation emissions,” The Journal of the Acoustical Society of America 153, 237–247 (2023).
  * Wan, Feng, and ter Haar (2015) M. Wan, Y. Feng, and G. ter Haar, eds.,  _Cavitation in Biomedicine_ (Springer Netherlands, Dordrecht, 2015).
  * Meroni  _et al._ (2022) D. Meroni, R. Djellabi, M. Ashokkumar, C. L. Bianchi, and D. C. Boffito, “Sonoprocessing: From Concepts to Large-Scale Reactors,” Chemical Reviews 122, 3219–3258 (2022).
  * Maeda  _et al._ (2018) K. Maeda, A. D. Maxwell, T. Colonius, W. Kreider, and M. R. Bailey, “Energy shielding by cavitation bubble clouds in burst wave lithotripsy,” The Journal of the Acoustical Society of America 144, 2952–2961 (2018).
  * Sapozhnikov, Maxwell, and Bailey (2021) O. A. Sapozhnikov, A. D. Maxwell, and M. R. Bailey, “Maximizing mechanical stress in small urinary stones during burst wave lithotripsy,” The Journal of the Acoustical Society of America 150, 4203–4212 (2021).
  * Yeats  _et al._ (2023) E. Yeats, N. Lu, J. R. Sukovich, Z. Xu, and T. L. Hall, “Soft Tissue Aberration Correction for Histotripsy Using Acoustic Emissions From Cavitation Cloud Nucleation and Collapse,” Ultrasound in Medicine & Biology 49, 1182–1193 (2023).
  * Vogel  _et al._ (1986) A. Vogel, W. Hentschel, J. Holzfuss, and W. Lauterborn, “Cavitation Bubble Dynamics and Acoustic Transient Generation in Ocular Surgery with Pulsed Neodymium:YAG Lasers,” Ophthalmology 93, 1259–1269 (1986).
  * Lee  _et al._ (2017) T. Lee, W. Luo, Q. Li, H. Demirci, and L. J. Guo, “Laser-Induced Focused Ultrasound for Cavitation Treatment: Toward High-Precision Invisible Sonic Scalpel,” Small 13, 1701555 (2017).
  * Shakya  _et al._ (2024) G. Shakya, M. Cattaneo, G. Guerriero, A. Prasanna, S. Fiorini, and O. Supponen, “Ultrasound-responsive microbubbles and nanodroplets: A pathway to targeted drug delivery,” Advanced Drug Delivery Reviews 206, 115178 (2024).
  * Habibi  _et al._ (2022) M. Habibi, S. Foroughi, V. Karamzadeh, and M. Packirisamy, “Direct sound printing,” Nature Communications 13, 1800 (2022).
  * Athanassiadis  _et al._ (2022) A. G. Athanassiadis, Z. Ma, N. Moreno-Gomez, K. Melde, E. Choi, R. Goyal, and P. Fischer, “Ultrasound-Responsive Systems as Components for Smart Materials,” Chemical Reviews 122, 5165–5208 (2022).
  * Ma  _et al._ (2022) Z. Ma, C. Bourquard, Q. Gao, S. Jiang, T. De Iure-Grimmel, R. Huo, X. Li, Z. He, Z. Yang, G. Yang, Y. Wang, E. Lam, Z.-h. Gao, O. Supponen, and J. Li, “Controlled tough bioadhesion mediated by ultrasound,” Science 377, 751–755 (2022).
  * Ohl  _et al._ (2006) C.-D. Ohl, M. Arora, R. Ikink, N. de Jong, M. Versluis, M. Delius, and D. Lohse, “Sonoporation from Jetting Cavitation Bubbles,” Biophysical Journal 91, 4285–4295 (2006).
  * Helfield  _et al._ (2016) B. Helfield, X. Chen, S. C. Watkins, and F. S. Villanueva, “Biophysical insight into mechanisms of sonoporation,” Proceedings of the National Academy of Sciences 113, 9983–9988 (2016).
  * Šarc  _et al._ (2018) A. Šarc, J. Kosel, D. Stopar, M. Oder, and M. Dular, “Removal of bacteria Legionella pneumophila, Escherichia coli, and Bacillus subtilis by (super)cavitation,” Ultrasonics Sonochemistry 42, 228–236 (2018).
  * Chahine  _et al._ (2016) G. L. Chahine, A. Kapahi, J.-K. Choi, and C.-T. Hsiao, “Modeling of surface cleaning by cavitation bubble dynamics and collapse,” Ultrasonics Sonochemistry 29, 528–549 (2016).
  * Reuter  _et al._ (2017) F. Reuter, S. Lauterborn, R. Mettin, and W. Lauterborn, “Membrane cleaning with ultrasonically driven bubbles,” Ultrasonics Sonochemistry 37, 542–560 (2017).
  * Barcikowski  _et al._ (2019) S. Barcikowski, A. Plech, K. S. Suslick, and A. Vogel, “Materials synthesis in a bubble,” MRS Bulletin 44, 382–391 (2019).
  * Soyama and Korsunsky (2022) H. Soyama and A. M. Korsunsky, “A critical comparative review of cavitation peening and other surface peening methods,” Journal of Materials Processing Technology 305, 117586 (2022).
  * Gutiérrez-Hernández  _et al._ (2023) U. J. Gutiérrez-Hernández, H. Reese, F. Reuter, C.-D. Ohl, and P. A. Quinto-Su, “Nano-Cracks and Glass Carving from Non-Symmetrically Converging Shocks,” Advanced Physics Research , 2300030 (2023).
  * Llorens  _et al._ (2016) C. Llorens, M. Argentina, N. Rojas, J. Westbrook, J. Dumais, and X. Noblin, “The fern cavitation catapult: Mechanism and design principles,” Journal of The Royal Society Interface 13, 20150930 (2016).
  * Koukouvinis, Bruecker, and Gavaises (2017) P. Koukouvinis, C. Bruecker, and M. Gavaises, “Unveiling the physical mechanism behind pistol shrimp cavitation,” Scientific Reports 7, 13994 (2017).
  * Cole (1948) R. H. Cole,  _Underwater Explosions_ (Princeton University Press, Princeton, New Jersey, 1948).
  * Wu  _et al._ (2020) W. Wu, A.-M. Zhang, Y.-L. Liu, and M. Liu, “Interaction between shock wave and a movable sphere with cavitation effects in shallow water,” Physics of Fluids 32, 016103 (2020).
  * Yu, Song, and Choi (2023) W. Yu, S. Song, and J.-I. Choi, “Numerical simulations of underwater explosions using a compressible multi-fluid model,” Physics of Fluids 35, 106102 (2023).
  * de Graaf, Penesis, and Brandner (2014) K. de Graaf, I. Penesis, and P. Brandner, “Modelling of seismic airgun bubble dynamics and pressure field using the Gilmore equation with additional damping factors,” Ocean Engineering 76, 32–39 (2014).
  * MacGillivray (2019) A. O. MacGillivray, “An Airgun Array Source Model Accounting for High-Frequency Sound Emissions During Firing—Solutions to the IAMW Source Test Cases,” IEEE Journal of Oceanic Engineering 44, 582–588 (2019).
  * Athanassiadis (2019) A. G. Athanassiadis,  _Optical Breakdown Acoustics: Transduction and Sensing Underwater_ , Ph.D. thesis, Massachusetts Institute of Technology, Cambridge, MA, USA (2019).
  * O’Reilly and Hynynen (2012) M. A. O’Reilly and K. Hynynen, “Blood-Brain Barrier: Real-time Feedback-controlled Focused Ultrasound Disruption by Using an Acoustic Emissions–based Controller,” Radiology 263, 96–106 (2012).
  * Novell  _et al._ (2020) A. Novell, H. A. S. Kamimura, A. Cafarelli, M. Gerstenmayer, J. Flament, J. Valette, P. Agou, A. Conti, E. Selingue, R. Aron Badin, P. Hantraye, and B. Larrat, “A new safety index based on intrapulse monitoring of ultra-harmonic cavitation during ultrasound-induced blood-brain barrier opening procedures,” Scientific Reports 10, 10088 (2020).
  * Doan, Huynh-The, and Kim (2022) V.-S. Doan, T. Huynh-The, and D.-S. Kim, “Underwater Acoustic Target Classification Based on Dense Convolutional Neural Network,” IEEE Geoscience and Remote Sensing Letters 19, 1500905 (2022).
  * Rom (2023) M. Rom, “Physics-informed neural networks for the Reynolds equation with cavitation modeling,” Tribology International 179, 108141 (2023).
  * Klapcsik  _et al._ (2023) K. Klapcsik, B. Gyires-Tóth, J. M. Rosselló, and F. Hegedűs, “Position control of an acoustic cavitation bubble by reinforcement learning,” (2023), arxiv:2312.05674 [cs, eess] .
  * Mallik, Jaiman, and Jelovica (2024) W. Mallik, R. Jaiman, and J. Jelovica, “Deep neural network for learning wave scattering and interference of underwater acoustics,” Physics of Fluids 36, 017137 (2024).
  * Kirkwood and Bethe (1942) J. G. Kirkwood and H. A. Bethe, “The pressure wave produced by an underwater explosion I,” Tech. Rep. Report No. 588 (Office of Scientific Research and Development, 1942).
  * Ivany (1965) R. Ivany,  _Collapse of a Cavitation Bubble in Viscous Compressible Liquid : Numerical and Experimental Analysis_ , Ph.D. thesis, University of Michigan, Ann Arbour, MI, USA (1965).
  * Lilliston (1966) R. R. Lilliston, “Calculations on the collapse of a spherical gas-filled cavity in a compressible liquid,” Tech. Rep. 2223 (SP 6-0002) (David Taylor Model Basin, Washington, DC, USA, 1966).
  * Denner, Xiao, and van Wachem (2018) F. Denner, C.-N. Xiao, and B. van Wachem, “Pressure-based algorithm for compressible interfacial flows with acoustically-conservative interface discretisation,” Journal of Computational Physics 367, 192–234 (2018).
  * Fuster and Popinet (2018) D. Fuster and S. Popinet, “An all-Mach method for the simulation of bubble dynamics problems in the presence of surface tension,” Journal of Computational Physics 374, 752–768 (2018).
  * Schmidmayer  _et al._ (2020) K. Schmidmayer, F. Petitpas, S. Le Martelot, and É. Daniel, “ECOGEN: An open-source tool for multiphase, compressible, multiphysics flows,” Computer Physics Communications 251, 107093 (2020).
  * Bryngelson  _et al._ (2021) S. H. Bryngelson, K. Schmidmayer, V. Coralic, J. C. Meng, K. Maeda, and T. Colonius, “MFC: An open-source high-order multi-component, multi-phase, and multi-scale compressible flow solver,” Computer Physics Communications 266, 107396 (2021).
  * Saade, Lohse, and Fuster (2023) Y. Saade, D. Lohse, and D. Fuster, “A multigrid solver for the coupled pressure-temperature equations in an all-Mach solver with VoF,” Journal of Computational Physics 476, 111865 (2023).
  * Koch  _et al._ (2016) M. Koch, C. Lechner, F. Reuter, K. Köhler, R. Mettin, and W. Lauterborn, “Numerical modeling of laser generated cavitation bubbles with the finite volume and volume of fluid method, using OpenFOAM,” Computers & Fluids 126, 71–90 (2016).
  * Denner, Evrard, and van Wachem (2020) F. Denner, F. Evrard, and B. van Wachem, “Modeling Acoustic Cavitation Using a Pressure-Based Algorithm for Polytropic Fluids,” Fluids 5, 69 (2020).
  * Gonzalez-Avila, Denner, and Ohl (2021) S. R. Gonzalez-Avila, F. Denner, and C.-D. Ohl, “The acoustic pressure generated by the cavitation bubble expansion and collapse near a rigid wall,” Physics of Fluids 33, 032118 (2021).
  * Wang  _et al._ (2024a) S.-P. Wang, H. Geng, S. Zhang, and S.-W. Wang, “Pressure waves from air gun bubbles: A numerical analysis based on the finite volume method,” Physics of Fluids 36, 013345 (2024a).
  * Crighton  _et al._ (1992) D. G. Crighton, A. P. Dowling, J. E. Ffowcs Williams, M. Heckl, and F. G. Leppington,  _Modern Methods in Analytical Acoustics_ (Springer London, London, 1992).
  * Christov and Christov (2017) I. C. Christov and C. I. Christov, “On mechanical waves and Doppler shifts from moving boundaries,” Mathematical Methods in the Applied Sciences 40, 4481–4492 (2017).
  * Denner (2024) F. Denner, “Acoustic pressure modulation driven by spatially non-uniform flow,” Journal of the Acoustical Society of America 155, 984–989 (2024).
  * Lai  _et al._ (2022) G. Lai, S. Geng, H. Zheng, Z. Yao, Q. Zhong, and F. Wang, “Early Dynamics of a Laser-Induced Underwater Shock Wave,” Journal of Fluids Engineering 144, 011501 (2022).
  * Liang  _et al._ (2022) X.-X. Liang, N. Linz, S. Freidank, G. Paltauf, and A. Vogel, “Comprehensive analysis of spherical bubble oscillations and shock wave emission in laser-induced cavitation,” Journal of Fluid Mechanics 940, A5 (2022).
  * Denner and Schenke (2023a) F. Denner and S. Schenke, “Modeling acoustic emissions and shock formation of cavitation bubbles,” Physics of Fluids 35, 012114 (2023a).
  * Wen  _et al._ (2023) H. Wen, Z. Yao, Q. Zhong, Y. Tian, Y. Sun, and F. Wang, “Energy partitioning in laser-induced millimeter-sized spherical cavitation up to the fourth oscillation,” Ultrasonics Sonochemistry 95, 106391 (2023).
  * Gilmore (1952) F. R. Gilmore, “The growth or collapse of a spherical bubble in a viscous compressible liquid,” Tech. Rep. Report No. 26-4 (California Institute of Technology, Pasadena, California, USA, 1952).
  * Flynn (1957) H. G. Flynn, “Collapse of a transient cavity in a compressible liquid,” Technical Memorandum 38 (Acoustics Research Laboratory, Havard University, Cambridge, MA, USA, 1957).
  * Hickling and Plesset (1963) R. Hickling and M. S. Plesset, “The collapse of a spherical cavity in a compressible liquid,” Tech. Rep. 85-24 (California Institute of Technology, Pasadena, California, USA, 1963).
  * Beyer (1974) R. T. Beyer, “Nonlinear Acoustics,” Technical Report (Brown University, Providence, RI, USA, 1974).
  * Flynn (1978) H. G. Flynn, “Cavitation dynamics: IV. Collapse of transient cavities,” Tech. Rep. (University of Rochester, Rochester, NY, USA, 1978).
  * Rice and Ginell (1943) O. K. Rice and R. Ginell, “The pressure wave produced by an underwater explosion: Part VI: The case of cylindrical symmetry,” Tech. Rep. 2023 (Office of Scientific Research and Development, 1943).
  * Rice and Ginell (1944) O. K. Rice and R. Ginell, “The pressure wave produced by an underwater explosion, Part VIII. The case of cylindrical symmetry, II,” Tech. Rep. 3950 (Office of Scientific Research and Development, 1944).
  * Whitham (1953) G. B. Whitham, “Spherical waves and shocks,” Tech. Rep. IMM-NYU 190 (Courant Institute of Mathematical Sciences, New York University, New York, NY, USA, 1953).
  * Snay (1966) H. Snay, “Hydrodynamic concepts: Selected topics for underwater nuclear explosions,” Technical Report NOLTR 65-52 (US Naval Ordnance Laboratory, White Oak, MD, USA, 1966).
  * Naugol’nykh and Roy (1971) K. A. Naugol’nykh and N. A. Roy, “Electrical discharges in water: A hydrodynamic description,” Tech. Rep. (Academy of Sciences USSR, Acoustic Institute, Moscow, 1971).
  * Cooper (1958) R. Cooper, ed.,  _Second Symposium on Naval Hydrodynamics: Hydrodynamic Noise, Cavity Flow_ (National Academy of Sciences–National Research Council, 1958).