# A molecular dynamics simulation of the abrupt changes in the thermodynamic
properties of water after formation of nano-bubbles / nano-cavities induced by
passage of charged particles

Ramin Abolfath1,5†, Niayesh Afshordi2,3,4, Sohrab Rahvar5, Adri van Duin6,
Martin Rädler7, Reza Taleei8, Katia Parodi7, Julie Lascaud7, Radhe Mohan1
1Department of Radiation Physics and Oncology, University of Texas MD Anderson
Cancer Center, Houston, TX, 75031, USA
2Department of Physics and Astronomy, University of Waterloo, Waterloo, ON,
N2L 3G1, Canada
3Waterloo Centre for Astrophysics, University of Waterloo, Waterloo, ON, N2L
3G1, Canada
4Perimeter Institute for Theoretical Physics, 31 Caroline St. N., Waterloo,
ON, N2L 2Y5, Canada
5Physics Department, Sharif University of Technology, P.O.Box 11365-9161,
Azadi Avenue, Tehran, Iran
6Department of Mechanical and Nuclear Engineering, Pennsylvania State
University, Pennsylvania 16802, USA
7Department of Medical Physics, Ludwig-Maximilians-Universität München,
Munich, Germany
8Department of Radiation Oncology, Jefferson University, Philadelphia PA, USA

###### Abstract

Purpose: To present a multi-scale formalism that accounts for the formation of
nano-scale bubbles/cavities owing to a burst of water molecules after the
passage of high energy charged particles that leads to the formation of “hot”
non-ionizing excitations or thermal spikes (TS).

Methods: We constructed amorphous track structures to account for the
formation of TS in a step-by-step Monte Carlo (MC) simulation of ionizing
radiation in liquid water. Subsequently, we simulate sudden expansion and
collective motion of water molecules employing a molecular dynamics (MD)
simulation that allows computation of ${\cal O}(10^{6})$ particle trajectories
and breaking/forming of chemical bonds on the fly using a reactive force
field, ReaxFF. We calculate the fluctuations of thermodynamic variables before
and after TS formation to model the macroscopic abrupt changes of the system,
possibly the occurrence of a first-order phase transition, and go beyond the
accessible simulation times by engaging fluid dynamic equations with
appropriate underlying symmetries and boundary conditions.

Results: We demonstrate the coexistence of a rapidly growing condensed state
of water and a hot spot that forms a stable state of diluted water at high
temperatures and pressures, possibly at a supercritical phase. Depending on
the temperature of TS, the thin shell of a highly dense state of water grows
by three to five times the speed of sound in water, forming a thin layer of
shock wave (SW) buffer, wrapping around the nano-scale cylindrical symmetric
cavity. The cavity contains a low density of water molecules; hence, it forms
a bubble-like state of water embedded in a room-temperature background of
water. The stability of the cavity as a result of the incompressibility of
water at ambient conditions and the surface tension allows the transition of
supersonic SW to a subsonic contact discontinuity and dissipation to thermo-
acoustic sound waves. Thus, TS gradually decays to acoustic waves, a channel
of deexcitation that competes with the spontaneous emission of photons, and a
direct mechanism for water luminescence. We further study the mergers of
nanobubbles that lead to fountain-like or jet-flow structures at the collision
interface. We introduce a time delay in the nucleation of nano-bubbles, a
novel mechanism, responsible for the growth and stability of much larger or
even micro-bubbles, possibly relevant to FLASH ultra-high dose rate (UHDR).

Conclusions: The formation of fast-growing nano-bubbles is an indication of a
much higher diffusion constant currently implemented in the reaction-diffusion
of the Monte Carlo models in the radiation physics community. The current
study is potentially significant in modeling overlap integrals and inter-track
recombination processes among the radiation tracks at UHDRs and is useful for
the FLASH radiotherapy data analysis. Our analysis predicts that thermal
radiation from the transient supercritical state of water localized in nano-
cavities wrapping around the track of charged particles can be manifested in
the (indirect) water luminescence spectrum.

###### pacs:

##  I Introduction

The formation of microscopic bubbles [1, 2, 3] with the size of sub-microns in
water has attracted interest due to a potential paradigm shift in the
treatment of cancer, such as hypothermia [4], and sonodynamic therapy [5].
Understanding a potential connection with other recently developed treatment
techniques, such as spatially fractionated radiotherapy [6], and FLASH
radiotherapy [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], may impact
the clinical applications of these modalities. However, the topic is more
directly relevant to various imaging techniques designed for the detection of
charged particle tracks by thermo-acoustic [21, 20, 22, 23], optical
(scintillation) [24] and nuclear magnetic resonance (NMR) beam topography or,
by magnetic resonance imaging (MRI) [25, 26], including a long list of other
scientific and engineering applications from nuclear physics [27] to cold
fusion and sonoluminescence effect [28], astrophysics, high-energy physics and
space radiation [29].

The passage of ionizing radiation results in the formation of ionizing and
non-ionizing energy deposition [30] to molecules within a nanometer radial
distance from the center of the track. Non-ionizing excitations deposited into
such a cloud of molecules release their internal energy either in a slow
nanosecond process by spontaneous emission of photons or in a fast femtosecond
process through the host nuclear vibrations. The latter further transforms the
release of excitation energies into heat exchange with the surrounding
molecules. Statistically, the spatial dimension of these cylindrical symmetric
clouds of molecules averaged over a large ensemble of tracks, is a rapidly
decreasing function of distance from the center of the track, $b$. See Eqs.
(2) and (LABEL:eq24x), analogous to a quadratic falling ($1/b^{2}$) of the
dose profile shown in [33]. Under such an ensemble average, the spatial
details of individual tracks that exhibit the inchoate distribution of energy
transfer [31] are washed out, converging into an amorphous track structure.

Although for several decades this energy transfer has been known to lead to
the formation of thermal spikes (TS) [34] and microscopic gas bubbles in
liquids [1], the mechanism and its possible role in advanced radiotherapy
applications have not been fully understood. In particular, TS’s wrap around
the cold ions in the medium, boosting the initial diffusion of species. This
acts as a driver for the outward burst of ions, and consequently, it can
potentially impact the dynamics of the fast chemical reaction following
ionizing radiation [30]. In radio/particle therapy, fast chemical and early
biological processes can be further influenced by the shock waves emerging
from TS [44], as well as bubble expansion and collapse [42]. These phenomena
may be of greater importance for treatment outcomes in novel radiotherapy
schemes such as FLASH radiotherapy, in which the dose is administered at an
ultrahigh dose rate (UHDR $\geq$ 40 Gy/s) in less than a second to reduce
toxicity in healthy tissues [7, 8, 9, 10, 13, 15]. In these conditions of
ultra-high fluence, the average distance between particle tracks is such that
they are close enough to interact. This has already been shown to influence
ultrafast chemical processes and is assumed to be one of the main reasons for
the decrease in toxicity in FLASH-RT [18]. However, current simulation
frameworks fail to reproduce experimental observations, highlighting the need
for more accurate modeling engines that account for physical phenomena
neglected in current state-of-the-art models.

Monte Carlo (MC) samplings of chemical species generated by radiation tracks
are the state-of-the-art computational model for the analysis and computer
simulations of the physico-chemical processes [35, 36]. The latest version of
these toolkits has added a series of chemical reaction rates based on
empirical data at low dose rates of homogeneous $\gamma$-rays. They, however,
do not completely model the collective dynamics of reactive species and change
of medium, the processes relevant to interplay between change in the
heterogeneous chemical changes of the tracks and the dynamics of TS’s at nano-
scales.

As noted above, current MC models cannot account for the FLASH-UHDR
experimental observations. Thus, they need substantial refinements to become
reliable computational tools for planning FLASH-UHDR radiotherapy. Our
technical aim is to go beyond this simple framework by simulating the
evolution and interaction of TS’s, produced in the wake of particle tracks,
using large-scale molecular dynamics (MD) at UHDR to demonstrate the
collective flow of molecules and introduce a mechanism that constitutes
elevation of reaction and diffusion constants of ionized species.

In this work, we present our theoretical predictions of the formation of nano-
bubbles, pico-seconds after generation of TS’s. Our motivation in this work is
to present the MD-based computational tool for the characterization of the
nano-bubbles and investigate the interplay between the TS’s and the inter-
track correlations such as recombination processes of chemical species that
are potentially substantial for the radio-biological effects and translational
applications of FLASH radiotherapy and their predictive outcomes.

Our calculation shows a rapid evolution from shock to contact discontinuity
after the creation of a thermal spike, where the pressure gradients inside and
outside the bubble cancel each other out. They act at two different locations
(inside and outside of the bubble interface) with opposite gradients with zero
curvature at the interface, which makes the bubble dynamically stable. Any
fluctuation inside the bubble that causes an increase in the negative slope of
the pressure (making it more negative) breaks the symmetry/balance in the
pressure profile and drags the molecules toward the center of the bubble (like
creating a vacuum in the center), which eventually enforces the bubble to
collapse suddenly.

##  II Method

Our methods to investigate the effects of TS’s and their thermal-diffusion
processes in a medium consist of two parts: (1) molecular dynamics (MD)
simulations using LAMMPS-ReaxFF with geometries and symmetries obtained from
amorphous track structure, and (2) solving a set of non-linear fluid dynamics
partial differential equations with appropriate boundary conditions, a coarse-
grained continuous model.

###  II.1 Amorphous track structure calculation

A model based on the amorphous track structure calculated by MC or other
analytical formulations [32, 33] allows the construction of initial conditions
in MD that start from the following assumptions: (1) sub-femtosecond passage
of a charged particle (e.g., proton) through nm-thick MD slab, (2)
instantaneous propagation of electromagnetic (EM)/photon fields to reach the
nm-size edge of a cylindrical symmetric TS, (3) sparse and relatively slow
delta-rays crossing out the MD box thus we neglect their effects in forming TS
inside of our simulation box. Violating these assumptions does not limit us to
constructing the initial condition in MD but reduces complex initial
conditions to simplistic ones. For ease of calculation and proof of principle,
we chose to present our results based on these simple initial conditions.

We neglect the contribution of the delta-rays with a range beyond the nm
lateral dimension of the TS confined in an MD computational box. They do not
contribute significantly to temperature increases and molecular thermal
excitations surrounding the charged particles. Note that 60% - 70% of the
energy deposition by the charged particles is spent on the creation of ions
and delta-rays. However, delta-rays originating from the position of the high
energy (low LET) charged particles travel far beyond the nm-scale size of our
computational box. In contrast, the delta-rays generated from low-energy (high
LET) charged particles carry low kinetic energies (KE) but with high numbers
and compactness. They deposit their energies within nm range from the position
they were created. Therefore, in a first assumption, the charged-particle
fields instantaneously create much denser molecular excitations than the slow-
moving and highly sparse energy deposition by delta-rays. The considered
uniform distribution does not lower the generalities of our propositions,
which could later integrate a more complex distribution of the thermal
energies among atoms.

In nanometer scales, relevant to our MD simulation geometries, the
instantaneous energy deposition within the subfemtosecond time scale, and the
time frame for sequential steps in MD is justifiable. To illustrate, let us
consider a fast proton passing through a 1 nm thick water slab in 0.1 fs, the
spatial and temporal scales in our MD simulations. To fulfill this condition,
the charged particle is required to pass through with minimum kinetic energy
calculated by $KE=(\gamma_{L}-1)Mc^{2}$ where $\gamma_{L}=(1-\beta^{2})^{1/2}$
and $\beta=v/c$. For a proton with a rest mass of $M=938.272$, MeV this
requires $KE\geq$ 0.52 MeV. For simplicity in calculating this threshold in
proton KE, we disregarded the deceleration of the proton along 1 nm distance
and the retardation of the electric field for the energy transform to that
volume. Taking into account these factors, we may end up with a slightly
higher KE threshold. Note that this energy is much lower than the nominal
energy of protons used in radiotherapy, i.e. 70 to 230 MeV. Considering the
range of protons as a function of KE, only protons at the very end of the
track violate our condition for generating cylindrical symmetric TS. In
amorphous track structure approximation protons below these minimum kinetic
energies tend to generate conical TS’s with circular symmetries, i.e., a
graded temperature profile. However, the rise in temperature for such low-
energy particles is not significant. We finally note that the energy transfer
in the lateral dimension takes place through the propagation of EM fields
below the time scales 0.001 to 0.01 fs, assuming a uniform spread of energy
transfer in a cylinder with a 1 nm diameter. Thus, the formation of molecular
excitations due to electron-proton impulsive Coulomb interaction is much
faster than the MD time step (0.1 fs or 100 as). Therefore, it would be a good
approximation to consider molecular excitations as an initial condition for
the simulation of the TS propagation.

Fig. 1 illustrates a schematic geometry of the amorphous track structure
considered in this work. The circles in Fig. 1 represent the radius of isoline
energy deposition of such amorphous tracks in our MD simulations.

###  II.2 Molecular dynamics

We employ a million particles MD simulation to calculate the trajectories of
water molecules, using LAMMPS. We chose ReaxFF as a model for inter-particle
interactions to be able to capture the molecular dissociation including
breaking and forming chemical bonds during the time evolution of particles.

ReaxFF [37] is a general bond-order-dependent potential that provides accurate
descriptions of bond breaking and bond formation. It is a bond-order-dependent
empirical force field method, which includes a polarizable charge function,
enabling application to a wide range of materials and accurate reproduction of
reaction energies and barriers. Recent simulations on a number of hydrocarbon-
oxygen systems [38, 39] and DNA [40] showed that ReaxFF reliably provides
energies, transition states, reaction pathways, and reactivity trends in
agreement with QM calculations and experiments.

For all the simulation results presented thereafter, the MD time step is
constant for a given run and chosen to be within 0.1-0.25 fs depending on the
desired total simulation time and temporal resolution. Time steps greater than
these values lead to numerical instability and early termination of the
simulation. Therefore, the physical simulation time can be easily calculated
by multiplying the total steps and the time step.

Information on the type of chemical reactions and the time evolution of damage
is collected by running the MD up to 50 ps, where the rearrangement of the
atomic coordinates has been deduced from a dynamical trajectory calculated by
ReaxFF. These simulations were performed using periodic and/or free boundary
conditions in canonical moles, pressure, and temperature (NPT), as well as a
constant energy and volume (NVE) ensemble with a Nose-Hoover thermostat for
temperature control.

In MD, NPT refers to a system of particles, equivalent to a canonical
ensemble, with a constraint on the values of the parameters, N, P, and T.
Practically a thermostat with virtual/auxiliary degrees of freedom is added to
apply a fluctuation-dissipation dynamics to the entire system, such as
coupling with a thermal bath following a series of Langevin equations to
enforce convergence of the thermodynamic parameters to the user-defined
values. Because of the constant P, the volume of the system, V, is a variable
in time. Similarly, for NVE the constraints apply to constant parameters V and
total energy E. Here, P is a free parameter and a time-dependent variable. The
system of particles as such is equivalent to a microcanonical ensemble.

We initially allow the entire system to run for several thousand steps to
reach thermal equilibrium at room temperature. Once the fluctuations in
temperature and pressure are damped, the passage of a charged particle and
induction of all ionizations and excitations scored within the MD time step
(e.g., 0.1 fs) that separates the state of the system before and after the
passage of fast-charged particles were performed. In amorphous track structure
construction, the ionizations and production of secondary electrons and delta
rays are neglected as they are sparse in their numbers and dimensions,
reaching outside of the core of TS domains. In contrast, the effects of non-
ionizing excitations are lumped and integrated into an abrupt change in the
velocities of a subset of water molecules forming TS’s. These new velocities
can be considered an initial condition for the rest of the simulation. We then
allow the system to evolve spontaneously and capture the trajectories of the
millions of particles.

Using the create style in LAMMPS, we set an initial condition to add TS’s by
changing the velocities of a group of atoms in a cylinder. The water molecules
outside the cylinder have a spectrum of kinetic energies that fluctuate
thermally around room temperature and ambient pressure with a macroscopic mass
density of room-temperature water (997 kg / m3). This method provides an
ensemble of velocities from a uniform distribution, using a random number
generator with the specified seed at the specified temperature. It selects
velocities randomly within a range of minimum to maximum values and
subsequently scales the velocities to match kinetic energies to the assigned
temperature. In our setup, the updated velocities are the sum of new and old
velocities (calculated under equilibrium conditions at room temperature). The
total velocities are subjected to two constraints to force the system of water
molecules in the TS to contain zero total linear and angular momentum.

###  II.3 MD computational box

The radius of cylindrical TS’s can be extracted from an amorphous track
structure calculation using MC toolkits such as Geant4 and Geant4-DNA [35, 36]
or using analytical formulation [33]. Hereby, we consider a cylinder with a
mean radius $\approx 1$ nm based on previous works [41, 42, 30].

Considering a single track was penetrating in an MD computational box. The
lateral and longitudinal size of the box can be determined by the radiation
deposited dose in Gy (J/kg) and the charged particle LET in eV/nm
(keV/$\mu$m), respectively.

For an illustration of the problem, let us consider a single track of 12C
close to its Bragg peak with LET of $8.45\times 10^{2}$ eV/nm creating a
cylindrical track with a radius of 1 nm and a length of 4 nm. The mass of
water in this cylinder, $m=\rho_{m}V=13\times 10^{-21}$ g, where $\rho_{m}=$1
g/cm3 and $V=\pi(1nm)^{2}\times 4(nm)\approx 13nm^{3}$ is the volume of the
cylinder. Considering that the energy deposited by 12C uniformly distributed
along the axis of the cylinder, we find $\Delta E={\rm LET}\times
L=(8.45\times 10^{2}{\rm eV/nm})\times 4{\rm nm}=34\times 10^{2}eV=5408\times
10^{-19}J$. The deposited dose to this volume is then $D=5408\times
10^{-19}J/(13\times 10^{-24}kg)=41.6\times 10^{6}$ Gy.

This is, of course, a very large dose, that is, 6 to 8 orders of magnitude
larger than the radiation dose in medical applications. However, by extending
the volume to an entire MD computational box, and including the mass of water
molecules in the denominator of $D$, the numerical value of the dose relative
to the calculated dose deposited in the cylinder is substantially lower. The
same is true for the increase in temperature, $\Delta T$, using the water-
specific heat, $C_{V}$, where the increase in temperature throughout the MD
computational box is substantially lower than $\Delta T$ calculated only
within the cylinder.

In this type of calculation, one needs to remove high-energy delta-rays and
any secondary particles that deposit their energies outside of the MD box.
However, we consider low-energy secondary electrons. The ones deposit energy
in nm distance from the primary particle. For simplicity in our numerical
illustration, because the majority of the secondary electrons in low energies
are generated close to the Bragg peak of a primary charged particle, we simply
do not take into account an escape of energy from the MD box due to delta-
rays.

As the deposited dose scales down linearly by increasing the lateral dimension
of the computational box (considering the entire mass of water in the box for
the calculation of the mass, in the denominator of the dose), we tend to
consider an MD box with a size of $10^{3}\times 10^{3}{\rm nm}^{2}(1\times
1\mu m^{2})$, or even $10^{4}\times 10^{4}{\rm nm}^{2}(10\times 10\mu m^{2})$.
The latter corresponds to a uniform lateral distribution of the single tracks
per single-cell nucleus.

In MD, using ReaxFF, a computational box with sizes of $10^{3}\times
10^{3}{\rm nm}^{2}$ and $10^{4}\times 10^{4}{\rm nm}^{2}$ is considered very
large. To put this in perspective, consider a box with lateral and
longitudinal dimensions of $8\times 8{\rm nm}^{2}$ and 4 nm. We can fit this
box with 24,399 hydrogen and oxygen atoms to construct 8,133 H2O molecules
with an average mass density of water under ambient conditions (1 g/cm3). To
scale up this number to a system of water molecules in a box of $8,000\times
8,000~{}{\rm nm}^{2}$ we need to nearly scale up 24,399 atoms to 24,399
million atoms, which is a huge number of particles. This is beyond current
computational capabilities.

To be able to perform MD simulations, we, therefore, need to limit ourselves
to smaller dimensions and truncate the number of particles up to a few million
fitted to up to $100\times 100~{}{\rm nm}^{2}$ computational box, which is
possibly the largest dimensions our current brute-force computational platform
can handle.

By truncating the number of atoms, we may encounter some computational
artifacts such as the rate of temperature decay of the entire system or the
rate of expansion of the TS wavefront (using periodic boundary conditions)
that we will discuss in detail below.

Note that for a given dose deposited by lighter charged particles than 12C,
such as the proton and electron, the minimum size of the computational box is
significantly smaller than the one we estimated above. This value and the
temperature of the core of TS scale up like $Z^{2}$ where $Z$ is the charged
particle atomic number.

##  III Results

Our results consist of the derivation of analytical equations, numerical
calculations, and MD simulations as follows.

###  III.1 Thermal-spike formation

We divide energy deposition into ionization and molecular excitation. The
molecular excitations that constitute TS are localized in space, within nm
distance from the core of the track where the charged particle passes through.
The charged particle induces an impulse to valence electrons that transfers to
their host nuclei due to electron-nucleus interaction that in turn converts to
the excess in the nuclear recoil kinetic energy, a mechanism that is faster by
at least five orders of magnitude than spontaneous emission of optical photons
[62]. Thus, it is plausible to assume that the TS’s are instantaneously
created following the passage of charged particles. The equivalent increase in
local temperature as a function of distance $b$ by a charged particle moving
with velocity $v$ can be calculated after integration over all normal modes of
these energy depositions. It can be calculated directly by the average energy
excitation on a molecule

$\displaystyle\Delta T(b,v)=\frac{1}{k_{B}}\langle
U_{\gamma}\rangle_{\tau}(b,v).$ (1)

or with a macroscopic mean value calculated by

$\displaystyle\left\langle\Delta T\right\rangle(b,v)=\frac{1}{c_{V}}\frac{2\pi
L}{m^{\prime}(b)}\int_{0}^{b}b^{\prime}db^{\prime}\langle
U_{\gamma}\rangle_{\tau}(b^{\prime},v).$ (2)

Here $c_{V}$ is the water-specific heat capacity at constant volume and
$m^{\prime}(b)$ is the mass of water embedded in a cylinder with radius $b$,
and longitudinal length, $L$. In the calculation of the stopping power, $b$ is
the impact factor. $\langle\cdots\rangle_{\tau}$ denotes the time averaging
over a time-dependent distribution function, embedded in the time integration
over the atomic excitation power-loss. $\langle U_{\gamma}\rangle_{\tau}$
describes the energy transferred to medium as a function of distance $b$ by a
charged particle moving with velocity $v$. In MD, Eq.(1) can be used for the
initial condition and distribution of KE among the atoms within TS. As a
function of time, this initial excessive KE decays and stores randomly into PE
given by the FF among all atoms including the ones surrounding TS, following
the equipartition theorem in interacting many particle systems. Whereas the
usage of water heat capacity in Eq.(2) assumes such relaxation and storage of
the initial thermal energy equally into KE and PE of the particles in the TS.

In a simple harmonic oscillator (SHO) model of an atom, the rate of change in
the kinetic and potential energies of a valence electron, the energy stored by
the charged particle can be dissipated by the rate $\gamma$ to the recoil
nucleus KE. This is the dark channel of energy dissipation that is five orders
of magnitude faster than the bright version of that channel (optical photon
emission), according to our QM calculation. Thus, the former excitation
converts to heat in the environment at a rate five orders of magnitude higher
than the latter [62].

It is a straightforward calculation to derive analytically the following
equation by taking into account the dissipation term perturbatively (small
$\gamma$ limit), using the solutions of SHO coupled with EM fields of a
relativistic moving charged particle [see, for example, Ref. 61]

$\displaystyle\langle U_{\gamma}\rangle_{\tau}(b,v)$ $\displaystyle=$
$\displaystyle\frac{m}{\sqrt{2\pi}}\frac{2}{\pi}\left(\frac{Zze^{2}}{mv^{2}}\right)^{2}$
$\displaystyle\int_{-\infty}^{\infty}d\omega\frac{\omega^{4}}{(\omega^{2}_{0}-\omega^{2})^{2}+\gamma^{2}\omega^{2}}\frac{\gamma}{\gamma^{4}_{L}(\omega)}$
$\displaystyle\left[K^{2}_{0}\left(\frac{\omega
b}{\gamma_{L}(\omega)v}\right)+\gamma^{2}_{L}(\omega)K^{2}_{1}\left(\frac{\omega
b}{\gamma_{L}(\omega)v}\right)\right],$

where $m$ and $z$ are the atom’s rest mass and the atomic number,
respectively. $Z$ and $v$ are the charge and velocity of the charged particle.
$\gamma_{L}(\omega)=1/\sqrt{1-v^{2}/(c/n(\omega))^{2}}$ and $K_{0}(x)$ and
$K_{1}(x)$ are the modified Bessel functions of the first kind. Here
$n(\omega)$ is the refractive index of the medium; hence $c/n(\omega)$
represents the speed of light in a dispersive medium. This equation is given
in the CGS unit and is derived simply by integration over the SHO power loss,
$\langle
U_{\gamma}\rangle_{\tau}=\int_{-\infty}^{\infty}dtP_{\gamma}=m\gamma\int_{-\infty}^{\infty}dt\left(\frac{d\vec{r}}{dt}\right)^{2}$,
using a uniform time-dependent distribution function in integration. Here,
$\vec{r}$ is the position of the atom, which gains excitation energy from the
charged particle.

It is worth noting that the points are drawn from Eqs. (2) and (LABEL:eq24x),
as $\Delta T$ scales by $Z^{2}$, heavy ions such as 12C, 56Fe and even 197Au,
near their Bragg peaks, cause an increase in temperature in a nanometer TS
from 104 to 105, and even 106 K.

Figure 1:  Sketch of a cylindrical symmetric TS surface which is a sharp
discontinuity interface separated inside and outside of a fluid denoted by
labels 2 and 1 respectively. (A) and (B) correspond to a view of an observer
in a frame of references at TS (living on a flat tangent space of the
cylinder) and the Lab. $v_{1}$ and $v_{2}$ give the rate of expansion of the
TS surface in (A) and (B) respectively. The TS rest frame is from an
observer’s view at rest on a co-moving tangent space of a cylinder.  Figure 2:
MD simulation of half-million water molecules in the presence of a cylindrical
symmetric single TS with 1 nm diameter and 0.1 fs MD time step. The initial
temperature profile that scales up proportionally to $Z^{2}$, corresponds to a
particle with LET$=8\times 10^{3}$ eV/nm, e.g., 56Fe, near the Bragg peak. The
shock front propagates initially with a speed five times the speed of sound in
water. It leaves a hole behind the shock and forms a thick layer of water with
a density larger than the average density of water at ambient conditions. The
shock front expands to approximately 10 nm in 2 ps, with a speed equivalent of
$v_{s}=5,000$ m/s.  Figure 3: Average temperature of the entire MD box,
equivalent of mean kinetic energy per atom, as a function of MD steps for two
types of charged particle, 12C (red) and proton, 1H (green). In this
simulation, each step is 0.25 fs time step, thus this simulation was
terminated at $t=5$ ps. For both 12C and 1H we used the same computational box
with PBC and 8133 water molecules and a size of $8\times 8\times 4$ nm. The TS
starts where the sharp peak is seen. The drop from the temperature peak can be
explained by the equal partitioning of the initial energy of atoms into their
kinetic and potential energies.

###  III.2 Molecular dynamics

Fig. 2 shows a typical state of water at a time frame after the sudden
expansion due to the initial rise in the temperature of water molecules within
the track of a particle, formed as a TS. It shows the time evolution of
approximately a million water molecules, a representation of a simulation for
a system of over a million particles. The interaction among the particles has
been given by the parameters of ReaxFF with an MD time step of 0.1 fs using
LAMMPS.

Induction of molecular thermal excitations and vibrations subsequent to the
passage of charged particles in a medium results in a rapid increase of the
kinetic energy and the thermal fluctuations within the nano-scale volumes,
wrapping around the charged particle trajectory. In the literature of
radiation physics, this phenomenon is known as TS formation.

This, in turn, results in an abrupt variation of the local thermodynamic
variables in the material embedded in a cylindrical shape, formed as an
amorphous track structure. Under high linear energy transfer, the nano-scale
temperature and pressure gradients drive the molecules and chemical species in
the track out of thermodynamic equilibrium relative to the rest of the
background molecules, leading to a sudden expansion, burst, and generation of
the chemicals in the TS. This includes the formation of the new chemicals
under extreme thermodynamic conditions just because of intensive temperature
and pressure in atomic scales, in addition to ionization events and ejection
of the secondary charged particles and delta-rays.

Macroscopically, the fluid dynamic equations with appropriate boundary
conditions govern the time evolution of the thermodynamic variables of the
system (atoms and molecules embedded in the track) coupled with the
environment (rest of fluid). For the boundary conditions, we note that the
profile of energy deposition calculated by either Monte Carlo and/or amorphous
track structure, predicts a gradual decrease from the center of the charged
particle track.

In the present model, we however choose slightly different geometry and
boundary conditions. Without loss of generality, and for the sake of
simplicity, we replace the gradual distribution of the thermodynamic variables
inside the track with a surface of discontinuity due to abrupt change in
pressure ($P$), number density ($\rho$), specific volume ($V=1/\rho$), and
temperature ($T$). We consider a uniform distribution of temperature inside
and outside of the track, denoted by $T_{2}$ and $T_{1}$, respectively.
Initially, at $t=0$ we consider $T_{2}>T_{1}$ where $T_{1}=300K$ is the room
temperature of the fluid. We calculate $T_{2}$ by balancing the internal
energy of the molecules inside the track with radius $R$ per unit length,
parallel to the direction of the charged particle track, and the charged
particle LET, hence $\Delta T={\rm LET}/(A\rho_{m}c_{V})$. Here $A$ is the
cross-section area of the track. $\rho_{m}$ and $c_{V}$ are the mass density
and the specific heat capacity at constant volume of unperturbed material
under radiation. We consider the ambient condition as the external condition
of the material before radiation exposure.

Alternatively one may consider the balance of LET with the kinetic energy of
the atoms in the cylinder using $\Delta KE=\frac{3}{2}Nk_{B}\Delta T$ where
$\Delta KE=LET\times L$ ($L$ is the height of cylinder in MD). Note that using
this formulation doubles $\Delta T$ as the energy stored in the form of PE (as
presumably included in $C_{V}$) is initially neglected. For illustration, let
us consider the red line in Fig. 3. It shows a rise in temperature after the
passage of a single track of Carbon, 12C with LET $=$ 845 eV/nm. For this
calculation, we used $C_{V}$ to estimate the initial rise in temperature.
Considering the alternative approach, the same temperature profile may be
applied to 12C with half of LET, i.e., LET = 422.5 eV/nm.

Going back to Fig. 2, we used $\Delta T$ (regardless of which formulation was
used for the calculation of LET corresponding to that $\Delta T$) to set up
the initial conditions of the TS in the MD calculation in a cylinder with a 1
nm diameter. The simulation shows the shock front propagates initially with a
speed higher than the speed of sound in water (1450 m/s). As seen, it leaves a
hole behind the shock and forms a thick layer of water with a density larger
than the average density of water at ambient conditions. In this calculation,
the shock front expands to approximately 10 nm in 2 ps, with a speed
equivalent of $v_{s}=5,000$ m/s.

The simulation was performed using periodic boundary conditions (PBC) in all
three directions. In PBC, there is an artifact of an interaction of the TS
with its own four lateral images via compression of water molecules between
the TS and the four sides of the computational box. Thus PBC underestimates
the rate of expansion of the TS’s compared with a similar computation using
free boundary conditions in lateral dimensions of the box and PBC in the
vertical direction (parallel to the axis of the cylinder), as pointed out in
Ref. [41]. The speed of propagation of TS calculated in our simulations has
been compared and verified with similar analytical and simulation models [41,
42, 43, 44, 45], using the same parameters as given in Fig. 4 in Ref. [41].

Regarding the comparison of the present work with early study from de Vera et
al. [41, 42], we remark that de Vera et al. have reported usage of both the
standard CHARMM-FF and a modification to reactive CHARMM-FF in their
publications for MD simulations of thermal spikes or shock waves in water in
the presence of DNA. The standard CHARMM-FF is not appropriate for simulating
chemical reactions. The CHARMM-FF potential (and AMBER-FF, GROMACS-FF that
have been used vastly for the simulation of protein folding) is a series of
simple harmonic oscillators without a transition point needed for breaking and
forming bonds. Therefore the FF as such can not simulate the formation of
radicals and other chemical species, in addition, it lacks simulating DNA
damage (either single-strand or double-strand breaks). That is the main
rationale for using ReaxFF (trained by a large series of quantum mechanical
calculations) in this work and comparing with reactive CHARMM-FF reported by
de Vera et al. [42] for the TS effects on the hydrodynamics of water.

###  III.3 Interface dynamics

Following the above conventions, we denote $P_{1},\rho_{1},T_{1}$, the
thermodynamic variables at a point outside of the track interface,
infinitesimally close to the surface of discontinuity as shown in Fig. 1.
Similarly $P_{2},\rho_{2},T_{2}$ denote the same set of thermodynamic
variables, at a similar point but on the opposite side of the interface,
inside of TS with an infinitesimal distance from the surface of discontinuity.
Note $V_{1}=1/\rho_{1}$ and $V_{2}=1/\rho_{2}$. We match these two sets of
variables using conventional fluid dynamic equations [46].

In the frame of reference of the TS interface (and moving laboratory frame of
reference with velocity in opposite direction) the conservation of mass,
energy, and momentum flux is given by

$\displaystyle{\rm mass:}~{}~{}~{}~{}~{}~{}\rho_{1}v_{1}=\rho_{2}v_{2}\equiv
J,$ (4) $\displaystyle{\rm
energy:}~{}~{}~{}~{}~{}~{}H_{1}+\frac{1}{2}v^{2}_{1}=H_{2}+\frac{1}{2}v^{2}_{2},$
(5) $\displaystyle{\rm
momentum:}~{}~{}~{}~{}~{}~{}P_{1}+\rho_{1}v^{2}_{1}=P_{2}+\rho_{2}v^{2}_{2}.$
(6)

Here $v_{1}$ and $v_{2}$ are the components of the velocities normal to the
thermal-spike cylinder, and $J$ is the conserved fluid current number density.
$H=\tilde{\gamma}/(\tilde{\gamma}-1)P/\rho$ is the enthalpy, independent of
the equation of the state of the fluid. Here $\tilde{\gamma}=C_{P}/C_{V}$. We
emphasize that this expression for $H$ is not based on the ideal gas equation
of state.

The system of coupled equations for the conservation laws can be solved
together by straightforward algebraic manipulations in terms of the Mach
number, $M=v_{1}/c_{s}$ and $\tilde{\gamma}$. As mentioned above, the side
denoted by 1, is in front of the TS interface, on the unperturbed side of the
fluid. The side denoted by 2, behind the TS interface and inside the TS
volume, may undergo a phase transition and form a hot and dilute gas with
different chemical compositions.

Using the continuity equations, and calling
$H=\tilde{\gamma}/(\tilde{\gamma}-1)(P/\rho)$ with
$\tilde{\gamma}=C_{P}/C_{V}$, assuming $\tilde{\gamma}_{1}=\tilde{\gamma}_{2}$
one can relate the thermodynamic variables right behind the TS surface,
$(\rho_{2},P_{2})$, to the same thermodynamics variables right in front of the
TS surface, $(\rho_{1},P_{1})$, two points adjacent to each other but in
opposite sides of the TS surface.

$\displaystyle\frac{\rho_{2}}{\rho_{1}}=\frac{v_{1}}{v_{2}}=\frac{(\tilde{\gamma}+1)\tilde{M}^{2}}{(\tilde{\gamma}+1)+(\tilde{\gamma}-1)(\tilde{M}^{2}-1)}\approx\frac{\tilde{\gamma}+1}{\tilde{\gamma}-1},$
(7)
$\displaystyle\frac{P_{2}}{P_{1}}=\frac{(\tilde{\gamma}+1)+2\tilde{\gamma}(\tilde{M}^{2}-1)}{\tilde{\gamma}+1}\approx\frac{2\tilde{\gamma}}{\tilde{\gamma}+1}\tilde{M}^{2}.$
(8)

Here $\tilde{M}=v_{1}/(\lambda_{1}c_{s})=v_{1}/\tilde{c}_{s}$ where
$\lambda_{1}=\sqrt{P_{1}\beta_{T,1}}$ is the speed of sound correction factor
due to incompressibility of the fluids, $\beta_{T,1}=V^{-1}(\partial
V/\partial P)_{T}$ is the fluid compressibility at point 1, and
$\tilde{c}_{s}=\lambda_{1}c_{s}=\sqrt{\frac{\tilde{\gamma}_{1}P_{1}}{\rho_{1}}}$
that is the equation for the speed of sound in the ideal gases, but because of
different values of the thermodynamic parameters using the incompressible
fluids at side 1 in this equation, $\tilde{c}_{s}$ doesn’t match with the
empirical values of speed of sound in ideal gases. Thus

$\displaystyle\tilde{M}=\frac{1}{\lambda_{1}}M$ (9)

does not represent Mach number unless the fluid inside 1 would be an ideal
gas.

To relate the temperatures $T_{2}$ to $T_{1}$, we must know the equation of
state of the fluid. Here, to proceed with the derivation of the analytical
solutions, let us assume the perfect gas equation of the state, $PV=RT$
($V=1/\rho$). This simplification of approximating water to an ideal gas
(which is not correct) allows us to imply
$H=\tilde{\gamma}/(\tilde{\gamma}-1)(P/\rho)=\tilde{\gamma}/(\tilde{\gamma}-1)k_{B}T$.

$\displaystyle\frac{T_{2}}{T_{1}}$ $\displaystyle=$
$\displaystyle\frac{[(\tilde{\gamma}+1)+2\tilde{\gamma}(\tilde{M}^{2}-1)][(\tilde{\gamma}+1)+(\tilde{\gamma}-1)(\tilde{M}^{2}-1)]}{(\tilde{\gamma}+1)^{2}\tilde{M}^{2}}$
(10)
$\displaystyle\approx\frac{2\tilde{\gamma}(\tilde{\gamma}-1)}{(\tilde{\gamma}+1)^{2}}\tilde{M}^{2}.$

Here $\tilde{M}=M$ as we have used the ideal gas equation of state to derive
Eq.(10). We recall that $R=N_{A}k_{B}$ therefore $\rho^{-1}=V/N_{A}$ is the
gas-specific volume (here $V$ is the total volume of the gas). Finally, by
substituting the density with specific volume, $V$, where $\rho=1/V$ in Eq.
(7), hence, $\rho_{2}/\rho_{1}=V_{1}/V_{2}$, we can check the consistency of
Eqs. (7-10) after multiplying Eqs. (7) to (8) which leads to Eq. (10), hence

$\displaystyle\frac{T_{2}}{T_{1}}=\frac{P_{2}}{P_{1}}\frac{V_{2}}{V_{1}}.$
(11)

From the RHS of Eqs. (7), (8), and (10), a supersonic motion of particles
localized on the TS interface is predicted in the limit of high LET, where
$\tilde{M}>>1$. In this limit, the pressure inside the thermal-spike volume
grows quadratically by $\tilde{M}$. In addition, $P_{2}>P_{1}$,
$\rho_{2}>\rho_{1}$, $T_{2}>T_{1}$, hence $v_{2}<v_{1}$ in TS frame of
reference. See Fig. 1.

In this limit, the rate of expansion of a cylindrical TS assuming a perfect
gas equation of the state of the fluid can be calculated analytically. To this
end, from Eq. (8) we find
$1/T_{1}\propto\tilde{M}^{2}=v_{1}^{2}/\tilde{c}_{s}^{2}$ (here
$\tilde{c}_{s}=c_{s}$ thus $\tilde{M}=M$). As $v_{1}=\dot{R}$, where $R$ is
the radius of the cylinder that encompasses the TS surface, we obtain
$1/T_{1}\propto\dot{R}^{2}$. For ideal gases, $T_{1}\propto V_{1}$ where
$V_{1}$ is the volume embedded the cylinder, $V_{1}\propto R_{1}^{2}$, thus
$1/T_{1}\propto 1/V_{1}\propto 1/R^{2}$. Combining these we obtain
$\dot{R}^{2}\propto 1/R^{2}$ and $RdR\propto dt$ which lead to $R^{2}\propto
t$ or $R\propto\sqrt{t}$. These are solutions of strong explosion in
cylindrical symmetry, calculated using scaling arguments and dimensional
analysis [43]. The same calculation, as described here, predicts $R\propto
t^{2/5}$ for a spherical surface as $V_{1}\propto R^{3}$ originally calculated
by Sedov and von Neumann [52, 53, 54].

As a function of time, the dissipation mechanisms of supersonic motion prevail
and temperature-pressure gradients damp out rapidly until $M$ drops to unity
asymptotically where $\rho_{2}=\rho_{1}$, $v_{2}=v_{1}$, $P_{2}=P_{1}$, and
$T_{2}=T_{1}$. From that time, the interface that has evolved into a thin
volume propagates like an amplitude of perturbed wave in pressure and
temperature with the speed of sound without transporting the molecules
collectively. We will evaluate the speed of sound of such a system by
performing MD in the next section.

For an imperfect gas or a liquid with $1/(P\beta_{T})>>1$, the condition
$\tilde{M}=1$ where $v_{1}=\tilde{c}_{s}$ in Eqs. (7) and (8) leads to the
following thermodynamic conditions $\rho_{2}=\rho_{1}$, $v_{2}=v_{1}$, and
$P_{2}=P_{1}$, respectively. For a perfect gas, these conditions can be
reached at a Mach number equal to 1, that is, as soon as $v_{1}$ reaches the
speed of sound in the gas. But for incompressible fluids, such as water, the
TS dissipates to a sound wave at the wavefront speed and continues to expand
even with a speed below the speed of sound, as in $c_{s}>\tilde{c}_{s}$.

###  III.4 Thermodynamic variables

The dynamical formulation of TS expansion in a continuous fluid, presented in
the preceding section assumes the surface of TS’s, a sharp moving
discontinuity. However, at the microscopic level, the TS interface consists of
particles with a finite size thickness. Away from the interface, inside TS, we
expect the formation of a diluted and/or gas phase cavity with a structure
sketched in Fig. 2 where the density drops below $\rho_{1}$ with temperature
and pressure above $T_{1}$ and $P_{1}$.

In general, the time evolution of the equation of the state of the fluid is
unknown. Many authors simplify the model by considering the fluid as a perfect
gas. To avoid such simplifications and to make the model as realistic as
possible, we perform an MD simulation on a large number of molecules to be
able to handle uncertainties in the equation of the state of the water in the
presence of TSs, presumably induced by the passage of charged particles. A
representative example of the simulation is shown in Fig. 2. To capture the
dynamical expansion of a typical TS in real-time, we have uploaded samples of
our simulations and videos, available on YouTube: [55, 56, 58].

In MD, as a function of time and space, we calculate the local thermodynamic
variables such as mass density ($\rho$), temperature ($T$), pressure ($P$),
specific heat capacities, $c_{P}=C_{P}/m$, $c_{V}=C_{V}/m$, isentropic index
$\gamma_{S}$ where the index $S$ refers to entropy, isothermal compressibility
$\beta_{T}$ and its inverse, the bulk modules, and finally the speed of sound
$c_{s}$ in the fluid. We denote $\lambda^{2}=P\beta_{T}$, a correction factor
in the transformation between ideal gas and real fluids, needed for TS
boundary conditions in the equations for conservation of mass, momentum, and
energy, as described above.

We perform this calculation in two ways: (1) ensemble averaging over the
entire system and taking all atoms into account and (2) partitioning the
computational box into small volume elements, and voxels, and using them as
statistical bins. These volumes are equivalent to grand canonical ensembles
where the exchange of energy and particles occurs in their boundaries.

Numerically it is straightforward to calculate the average mass ($m$) and
kinetic energy ($K.E.$) in each volumetric element, at each time slice, to
obtain $\rho(\vec{r},t)=\langle m(\vec{r},t)\rangle/V$, and $\langle
K.E.(\vec{r},t)\rangle=3k_{B}T(\vec{r},t)/2$.

In MD the local pressure ($P$) and its deviation from the ideal gas (system of
non-interacting particles) is calculated using the Virial equation [47, 48]

$\displaystyle P(\vec{r},t)$ $\displaystyle=$
$\displaystyle\rho(\vec{r},t)k_{B}T(\vec{r},t)+\frac{1}{3V}\left\langle\sum_{i<j}^{N}\vec{f}_{ij}(\vec{r},t)\cdot\vec{r}_{ij}(\vec{r},t)\right\rangle.$

Here $\vec{f}_{ij}$ is a pair-wise force acting between particles $i$th and
$j$th with separation $\vec{r}_{ij}=\vec{r}_{i}-\vec{r}_{j}$. The state of an
ideal gas with negligible $\vec{f}_{ij}$ is given by the first term in Eq.
(LABEL:eq00). Note that $\vec{r}$ refers to the location of volumetric
elements.

It would be useful to present a review of the relations between thermodynamic
variables, relevant to our simulations. The enthalpy and the internal energy
of the system are given by $H=E+PV$ and U = E = K. E. + P. E., thus by
definition, $C_{V}=\left(\frac{\partial U}{\partial
T}\right)_{V}=\left(\frac{\partial U}{\partial
T}\right)_{P}+P\left(\frac{\partial V}{\partial T}\right)_{p}$, and
$C_{P}=\left(\frac{\partial H}{\partial T}\right)_{P}=\left(\frac{\partial
U}{\partial T}\right)_{P}+P\left(\frac{\partial V}{\partial T}\right)_{p}$. At
constant $P$, the internal energy is only a function of $T$ and $V$, thus
$U=U(T,V)$ then $dU=\left(\frac{\partial U}{\partial
T}\right)_{V}dT+\left(\frac{\partial U}{\partial V}\right)_{T}dV$. Applying a
chain rule of partial differentiation we obtain $\left(\frac{\partial
U}{\partial T}\right)_{P}=\left(\frac{\partial U}{\partial
T}\right)_{V}+\left(\frac{\partial U}{\partial
V}\right)_{T}\left(\frac{\partial V}{\partial T}\right)_{P}$. Hence
$C_{P}=\left(\frac{\partial U}{\partial T}\right)_{P}+P\left(\frac{\partial
V}{\partial T}\right)_{p}=\left(\frac{\partial U}{\partial
T}\right)_{V}+\left(\frac{\partial U}{\partial
V}\right)_{T}\left(\frac{\partial V}{\partial
T}\right)_{P}+P\left(\frac{\partial V}{\partial
T}\right)_{P}=C_{V}+\left[P+\left(\frac{\partial U}{\partial
V}\right)_{T}\right]\left(\frac{\partial V}{\partial T}\right)_{P}$. In non-
interacting ideal gas, $U=KE$ and $PE=0$. Thus the mean value of $U$ is only a
linear function of temperature, $T$, and it is independent of volume $V$,
hence $\left(\frac{\partial U}{\partial V}\right)_{T}=0$. Applying the
equation of the state for ideal fluids, $PV=nRT$ gives
$C_{P}-C_{V}=P\left(\frac{\partial V}{\partial T}\right)_{P}=nR$ with
$n=N/N_{A}$. $N$ and $N_{A}$ are the number of particles and Avagadro’s number
respectively. Defining per molar heat capacities, $C_{P}/n$ and $C_{V}/n$,
gives $C_{P}-C_{V}=R$.

For a system of interacting particles, one can show (using thermodynamic’s
free energies and Maxwell relations)
$C_{P}=C_{V}+VT\frac{\alpha^{2}_{P}}{\beta_{T}}$, where
$\alpha_{P}=\frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_{P}$ and
$\beta_{T}=-\frac{1}{V}\left(\frac{\partial V}{\partial P}\right)_{T}$ are the
isobaric expansion coefficient and isothermal compressibility respectively.
For incompressible fluids (including water), $\alpha$ and $\beta$ are
negligible thus $C_{P}/C_{V}\approx 1$ as we found in this study.

The isentropic (adiabatic) index, was calculated using
$\gamma_{S}=\left(\partial H/\partial
U\right)_{S}=(C_{P}/C_{V})/(P\beta_{T})$. For ideal gas $P\beta_{T}=1$ hence
$\gamma_{S}=\tilde{\gamma}=C_{P}/C_{V}$ where $\tilde{\gamma}$ used to
describe the TS wavefront in the preceding section. For non-ideal fluids,
$\lambda^{2}=P\beta_{T}$, calculated in MD, is the correction factor that must
be applied to $\gamma_{S}$.

To calculate the thermodynamic variables listed above, we use well-established
relations that correlate the fluctuations of the variables calculated in a
finite MD system and their thermodynamic counterparts. In NVE microcanonical
and/or NPT ensembles. For more details, see page 65 [49], eq.(2.70) in [50],
and [51] (page. 333 on a general concept of Fluctuations in thermodynamics)
where

$\displaystyle C_{V}=\frac{\langle E^{2}\rangle-\langle
E\rangle^{2}}{k_{B}T^{2}}$ (13)

where $E$ is the atomistic total energy, the sum of the particles’ kinetic and
potential energies, $E=KE+PE$. Similarly to calculate $C_{P}$ and $\beta_{T}$
(Eqs. 2.85 and 2.86 in [50]) in NPT ensemble we use

$\displaystyle C_{P}$ $\displaystyle=$ $\displaystyle\frac{\langle
H^{2}\rangle-\langle H\rangle^{2}}{k_{B}T^{2}}$ (14) $\displaystyle=$
$\displaystyle\frac{\langle(E+PV)^{2}\rangle-\langle
E+PV\rangle^{2}}{k_{B}T^{2}}$ $\displaystyle=$ $\displaystyle\frac{\langle
E^{2}\rangle-\langle E\rangle^{2}}{k_{B}T^{2}}+P^{2}\frac{\langle
V^{2}\rangle-\langle V\rangle^{2}}{k_{B}T^{2}}$ $\displaystyle+$
$\displaystyle 2P\frac{\langle EV\rangle-\langle E\rangle\langle
V\rangle}{k_{B}T^{2}}$

and

$\displaystyle\beta_{T}=\frac{1}{k_{B}T}\frac{\langle V^{2}\rangle-\langle
V\rangle^{2}}{V}$ (15)

Note that the the total volume in the NTP ensemble fluctuates whereas in NVE
is constant. In MD, we calculate Voronoi as a metric for the inter-atomic
occupation space and the specific volume per atom, $V_{\rm vor}$, a variable
as a function of time in both NPT and NVE ensembles. It satisfies a
geometrical sum rule, $V=\sum_{i=1}^{N}V_{{\rm vor},i}$, where $N$ is the
number of atoms in the computational box. Hence, $\delta V^{2}=\langle
V^{2}\rangle-\langle V\rangle^{2}$ where $\langle
V^{2}\rangle=\sum_{i=1}^{N}V_{\rm vor,i}^{2}/N$ and $\langle
V\rangle=\sum_{i=1}^{N}V_{\rm vor,i}/N$. The results of this calculation for
two ensembles of NPT and NVE before and after the introduction of a single TS
are shown in Fig. 4. It is interesting to note that the transition before and
after the addition of TS is continuous in the atomic-specific volumes. We also
observe that fluctuations in $V$ are substantially smaller in NVE compared to
the NPT ensemble, as expected.

As shown in Fig. 5, our MD calculation of the specific heat capacities of the
bulk water with ReaxFF at ambient conditions, and before creation of the TS’s,
converges to $C_{V}=18~{}{\rm Cal/mol/K}$ and $C_{P}=20~{}{\rm Cal/mol/K}$.
The initial fluctuations at the beginning of the simulations are the steps in
MD before reaching the system to the thermal equilibrium at room temperature.

After the creation of TS’s, a diluted state of water, mixed with chemical
species that are the products of dissociation of water molecules due to high-
temperature fluctuations, are formed, and an abrupt drop in heat capacities is
seen as an indication of a first-order phase transition.

These values are in agreement with the reported experimental data from NIST
[59]. Thus we obtain $\tilde{\gamma}=C_{P}/C_{V}=1.11$, hence the condition
for the strong explosion [46] can be obtained if the ratio of internal to
external pressure in TS’s is satisfied
$P_{2}/P_{1}>\rho_{2}/\rho_{1}=(\tilde{\gamma}+1)/(\tilde{\gamma}-1)=19$. See,
e.g., Eqs. (7) and (8) in the limit of $\tilde{M}>>1$). Our calculated value
for $\tilde{\gamma}=1.11$ is close to the value given by de Vera et al. [41],
$\tilde{\gamma}=1.22$. However, the predicted strong explosion criteria using
this value reduces the lower bound of the pressure ratio to 10.

To compare the numerical values with the prediction of adiabatic expansion we
use $c_{s}=1450$ m/s in water and the experimental value of the specific heat
that predicts $\gamma_{S}\approx 2.9$ which is three times larger than
reported $C_{P}/C_{V}\approx 1.11$. With these considerations, one may expect
deviations from Sedov’s solutions [52, 43] obtained analytically through ideal
gas assumptions.

Figure 4: STD of volumes per atom computed by Voronoi-method in LAMMPS.
Figure 5: Heat capacities calculated for the entire computational box at
constant pressure and volume using NPT and NVE ensembles. A sharp drop in
$C_{P}$ and $C_{V}$ occurs at the onset of TS formation. This is an indication
of a first-order phase transition and the emergence of a dilute gas of water
molecules mixed with other chemical species inside the nano-cavity. The
predicted values of heat capacities at room temperature (300 K) show a good
match with published NIST data.

Fig. 6 shows the result of our MD calculation for the isothermal
incompressibility in the NPT ensemble as a function of pressure in the water
and MD steps (time). Before TS generation, $P\beta_{T}=(P/V)(\partial
V/\partial P)_{T}$ approaches a small value due to the incompressibility of
water. The small fluctuations are seen because of the local pressure
controlling the Nosé–Hoover thermostat and barostat to equilibrate the system
to the requested $T$ and/or $P$. After TS’s though we initially observe strong
fluctuations in $P\beta_{T}$ which quickly saturates to an asymptotic value of
$P\beta_{T}\approx 0.4$. This is because of a diluted gas of water inside TS,
mixed with a condensed state of water on the shell and outside of TS. We note
that for an ideal gas $P\beta_{T}=1$.

Finally, we combine all relevant thermodynamic variables to calculate the
speed of sound as a local discontinuity due to the formation of TS. As seen in
Fig. 7, after initial fluctuations, the speed of water converges slowly to a
correct value of $\approx 1450$ m/s. Once the TS was added the speed of sound
of mixed phases drops to that value and converges to $\approx 1200$ m/s.

Figure 6: Top: Calculation of $P\beta_{T}=(P/V)(\partial V/\partial P)_{T}$ as
a function of pressure for normal water is shown. As
$\gamma_{S}=(C_{P}/C_{V})/(P\beta_{T})$, strong deviation from ideal gas due
to incompressibility of liquid water is expected at low pressures. Bottom:
$P\beta_{T}$ vs. MD steps. It seems the system of TS fluctuates around
$P\beta_{T}=0.4$, lower than the corresponding value for ideal fluids
($P\beta_{T}=1$ for ideal gas).  Figure 7: Calculated speed of sound using
$c_{s}=\sqrt{\tilde{\gamma}/\rho\beta_{T}}$ as a function of MD steps for two
NTP and NVE ensembles using the entire atoms in the MD simulation. In normal
water at the simulated pressure, the speed of sound is calculated close to
$1600m/s$ in the NPT ensemble. With an increase of MD steps in normal water
(before the formation of TS) the numerical value slowly converges to $1450m/s$
consistent with the reported values. We truncated the simulation at $1600m/s$
to speed up the simulations. The value of $c_{s}$ after TS formation is close
to SW velocity calculated from the rate of nano-cavity radius.

###  III.5 Chemical changes

Depending upon the initial temperature of TS, it may cause changes in the
local chemical composition of the track and its surrounding environment.

Dissociation of $\ce{{}^{.}OH}-\ce{{}^{.}H}$ bond of a water molecule requires
118.8 kcal/mol equivalent of $\Delta T=$ 59,776 K. Whereas electron ionization
of water molecule requires $\approx$ 13 eV, equivalent of 150,858 K in the
core of TS’s.

Thus isolated TS’s with high enough kinetic energy, presumably generated by
high LET charged particles (heavier elements than C, such as Fe) may act alone
as a source of ionizing radiation. Due to strong thermal fluctuations, bounded
electrons and/or atoms to water molecules vibrate severely, overcome the
potential energy barriers, and go over a transition state that turns the
constituted electrons and/or atoms into unbounded from their host and form a
local plasma. A dissociation mechanism, similar to the photoelectric and/or
Compton effect in the presence of $\gamma$-rays.

To illustrate this, during our MD simulations, we have picked isolated atoms
and singled out events in which we have observed the separation of hydrogen
atoms, hence the formation of ^.OH-radicals, incorporated with the transient
formation and subsequently dissociation of $H_{2}$ molecules as a function of
time, as shown in Figs. 8-10.

Figure 8: Time-evolution of H-OH bond distance, corresponding to water
molecules in SW ring/layer of nano-cavity. The initial compression of covalent
bonds due to the propagation of highly pressured shock-front stores a large
potential energy. Release of stored energy into kinetic energy of H
dissociates water molecule and forms an atomistic H and ^.OH-radical. This
simulation demonstrates that the nano-cavity is not only mechanically active
and its shock-front may destroy the biomolecules but it is also chemically
active and it is full of ROS.  Figure 9: Time-evolution of two randomly
selected hydrogen atoms separated from their host water molecules in the nano-
cavity ring. During MD simulation, they come close to each other, form a
hydrogen atom, and become separated due to high-energy collisions by other
molecules.  Figure 10: Time-evolution of a randomly selected water molecule in
the ring of nano-cavity. Due to highly energetic collisions hydrogen
dissociates from the host and leaves an ^.OH-radical behind.  Figure 11:
Kinetic energy (KE) profile of water molecules scored in $0.5\times 0.5$ nm
bins along the dashed line as depicted in the inset, corresponding to $t=0.5$
ps and $t=2$ ps of MD simulation in Fig. 2. The cusps in the KE are
reminiscent of the compressed ring in the density of water as shown in the
inset and in Fig. 2. The gaps in KE, behind and in front of the shoulders
(inside and outside of the cylinder) can be interpreted as discontinuities in
the local thermodynamic parameters. The decay in the amplitude of the gap as a
function of time which is an indication of SW crossing to a perturbative
thermo-acoustic wave is visible. The displacement of the cusps as a function
of time provides a rough estimate of the SW front speed, $v_{s}=5,000$ m/s.
Figure 12: Mass density profile of water molecules scored as in Fig. 11. The
initial mean density of water molecules after thermal equilibration (up to
0.15 ps), fluctuates around the bulk mass density (green dashed line) even
slightly after initiation of the TS. It shows a time delay in response of
molecular density to the rise in temperature and pressure of the TS. A drop in
density in the core of SW is a result of TS expansion in the form of SW.
Figure 13:  Sketch of an ideal cavity inside a fluid formed by pushing out the
particles outward. The hole represents the cavity and the shaded ring
represents a nano-scale shock-wave shell. The cavity grows with a rate
proportional to the TS expansion velocity, $u_{s}$, that determines the shell
thickness and the SW speed, $v_{s}$. The velocities are given in the Lab frame
of reference.

###  III.6 Nano cavities/bubbles

In this section, we seek an extension of the analytical equations for the
growth of the cavity separated by a wall thickness from the rest of the water
molecules, consistent with the results obtained from MD simulations.

For the illustration of the problem, we have calculated a profile of KE and
the mass density as a function of time as depicted in Figs. 11 and 12, using
our MD simulation. The cusps, interpolating the valleys inside of the cavity
to the outside are the positions of the TS interface with a finite width.

An ideal cavity inside a fluid can be approximated with an empty space, like a
void. However, in the following section, we discuss a very rich internal
structure of the cavity which exhibits oscillatory patterns in temperature,
pressure, and density profiles. As shown in Fig. 13, at the interface between
unperturbed fluid and cavity, internal pressure due to high temperature in the
cavity pushes out the fluid particles outward. Thus the cavity can be treated
like a rigid cylinder with a radius growing outward with speed, $u_{s}$. A
rapid expansion of the cavity compresses the fluid in the shell, even if the
unperturbed fluid is incompressible, like liquid water.

The compressed fluid forms a nano-meter size ring with a layer thickness which
is determined by the in-flux of fluid particles and the speed of expansion,
$u_{s}$. The outer surface of the cylinder expands with a speed, higher than
$u_{s}$, and can even reach a supersonic velocity at high temperatures, hence
$v_{s}>c_{s}$.

To analyze this ideal situation, we neglect the velocity of the fluid in
region 1 (see Fig. 13). We further assume the compressed fluid within the TS
shell moves with the cavity velocity, thus $v_{1}=0$ and $v_{2}=u_{s}$.

In the rest frame of TS wavefront, where $v_{s}=0$, the radial velocities
across the TS surface are transformed as $v_{1}=-v_{s}$ and
$v_{2}=u_{s}-v_{s}=u_{s}+v_{1}$. This corresponds to an observer sitting at
the outer surface of TS, as sketched in Fig. 13. In terms of the conserved
mass-flux density, $J$, the velocities can be rearranged as $v_{1}=JV_{1}$,
$v_{2}=JV_{2}$ where $V_{1}=1/\rho_{1}$, $V_{2}=1/\rho_{2}$ are the
corresponding specific volumes. The conservation of momentum across the outer
cavity interface in this TS frame of reference can be written

$\displaystyle P_{1}+J^{2}V_{1}=P_{2}+J^{2}V_{2},$ (16)

which gives

$\displaystyle J=\pm\sqrt{\frac{P_{2}-P_{1}}{V_{1}-V_{2}}}.$ (17)

Considering the minus sign in Eq. (17) for $J$ which describes the flow of
particles from 1 to 2 across the cavity interface as it grows, yields

$\displaystyle-u_{s}$ $\displaystyle=$ $\displaystyle
v_{1}-v_{2}=J(V_{1}-V_{2})=-\sqrt{(P_{2}-P_{1})(V_{1}-V_{2})}$ (18)
$\displaystyle=-\sqrt{\left(\frac{\rho_{2}-\rho_{1}}{\rho_{1}\rho_{2}}\right)(P_{2}-P_{1})}.$

Knowing that in the rest frame $v_{1}>v_{2}$, we expect the layer thickness of
TS shell to grow up as a function of time. This is due to the positive influx
rate of the particles, accumulated behind the cavity interface in region 2,
entering from the interface front in region 1.

The conservation of mass, momentum, and energy flux leads to

$\displaystyle\frac{\rho_{1}}{\rho_{2}}=\frac{v_{2}}{v_{1}}=\frac{u_{s}+v_{1}}{v_{1}}=\frac{P_{1}(\tilde{\gamma}+1)+P_{2}(\tilde{\gamma}-1)}{P_{1}(\tilde{\gamma}-1)+P_{2}(\tilde{\gamma}+1)},$
(19)

in which it can be algebraically manipulated to give

$\displaystyle\frac{P_{2}}{P_{1}}$ $\displaystyle=$ $\displaystyle
1+\frac{\tilde{\gamma}(\tilde{\gamma}+1)}{4}\left(\frac{u_{s}}{\tilde{c}_{s}}\right)^{2}$
(20) $\displaystyle+$
$\displaystyle\tilde{\gamma}\left(\frac{u_{s}}{\tilde{c}_{s}}\right)\sqrt{1+\frac{(\tilde{\gamma}+1)^{2}}{16}\left(\frac{u_{s}}{\tilde{c}_{s}}\right)^{2}}.$

With this ratio and with the help of Eq.(19) it is straightforward to solve
for $v_{s}$ in the Lab frame

$\displaystyle
v_{s}=\frac{\tilde{\gamma}+1}{4}u_{s}+\tilde{c}_{s}\sqrt{1+\frac{(\tilde{\gamma}+1)^{2}}{16}\left(\frac{u_{s}}{\tilde{c}_{s}}\right)^{2}}.$
(21)

Once the system reaches the equilibrium, the difference between the internal
and external density and pressure vanishes, i.e., $\rho_{2}=\rho_{1}$ and
$P_{2}=P_{1}$. Applying these conditions in Eqs. (19-21), we find $u_{s}=0$,
hence $v_{s}=\tilde{c}_{s}$. The SW front reduces to a regular sound wave that
travels with speed $\tilde{c}_{s}$.

A similar equation to Eq. (21) has been proposed empirically by Sing et al.
[60], such that $v_{s}=a+bu_{s}$. Here $a=1.364\times 10^{3}$ m/s and
$b=2.128$ where $a$ and $b$ are the fitting parameters obtained after
interpolation to the experimental data. At the limit of $u_{s}=0$, Sing et al.
[60] empirical results predict $v_{s}=1.364\times 10^{3}$ m/s, which is
slightly below speed of sound in water. Interestingly, in this limit our
analytical Eq. (21) predicts $v_{s}=a=\tilde{c}_{s}<c_{s}$ where
$\tilde{c}_{s}=\sqrt{\tilde{\gamma}P/\rho}$, and
$c_{s}=\sqrt{\tilde{\gamma}/\rho\beta_{T}}$. Recall that $c_{s}$ is the
analytical expression for the speed of sound in water, calculated by our MD
calculation (in the absence of TS), and was converged to 1450 m/s. Moreover,
as pointed out before, our MD simulation predicts $\tilde{\gamma}=1.1$.
Comparing Sing et al.’s results with ours, $b=2.128$ is almost four times
higher than our analytical expression, $(\tilde{\gamma}+1)/4\approx 0.5$ for
$u_{s}=0$, whereas in the limit of $u_{s}>>\tilde{c}_{s}$, where Eq.(21)
approaches to $v_{s}=u_{s}(\tilde{\gamma}+1)/2\approx u_{s}$, and our
prediction gives, $b=1$, which is closer to the empirical value of the Sing et
al. results.

Figure 14: Time-evolution of temperature, mass density, and pressure profiles
in a series of concentric cylindrical bins. The initial TS was generated based
on 12C parameters. Symbols corresponding to numbers are the sequence of the
time. the black line corresponding to the symbol, 0, shows the state of water
at the time before TS entrance. The oscillations in time in the first bin are
visible. Across the interface, on the sixth bin, a sharp variation in
temperature and mass density is seen in contrast to a smooth variation in
pressure, an indication of a cross-over from shock to contact discontinuity.
The lateral thickness of bins is $\Delta r=0.2$ nm. The longitudinal dimension
of the bins spans over the entire 4 nm height of the computational box. The
time difference between sequential symbols corresponds is $\Delta t=0.025$ ps
thus the last curve (shown by $+$), was sampled at $t=0.25$ ps after the
generation of the TS. The unit of pressure is the Atmosphere (1 Atm = 101.325
kPa).

###  III.7 Transition from shock to contact discontinuity

As can be seen in Fig. 11, the KE of the TS dissipates significantly within
the first few ps after the creation of a TS. However, the cusp in KE that
separates inside and outside of the TS persists. Similarly, the density
difference between in and out of TS persists to exhibit a cusp, as seen in
Fig. 12. However, one may expect the initial sharp difference in the pressure,
inside and outside of the surface of discontinuity, to deform into a smooth
curve. We interpret this behavior as a persistent jump in the density and
temperature at the surface of discontinuity but a smooth/continuous transition
in the pressure difference, a time evolution and a cross-over from shock to
contact discontinuity [46].

Beyond this transition time, the thermal spike is expected to expand following
standard heat transfer from a source of the circular symmetric high-
temperature expanding cloud. The mass transfer follows a similar diffusion
equation [57] but with a time-varying diffusion constant that decays to the
equilibrium value at room temperature.

To examine this hypothesis, we construct a slightly different type of binning
using a series of concentric cylinders to sample pressure, temperature, and
mass density as a function of time, as shown in Fig. 14. Here 20 concentric
cylindrical volumes from the center of the axis of the cylinder to a 4 nm
radius, are used. These volumes, which are stationary at the Lab frame of
reference, span over the lateral dimension within the following range of
$0\leq r\leq 4$ nm. Thus $\Delta r=0.2$ nm is the lateral thickness of the
cylindrical shells. The height of each bin spans over the entire height of the
computational box, in this case, 4 nm. As the thickness of each shell is fixed
to 0.2 nm, the bin volume and the number of particles increase by the bin
radius.

In Fig. 14, the symbols correspond to the time dependence scoring of the
temperature, mass density, and pressure in each 100-time step, hence $\Delta
t=100\times 0.25$ fs $=0.025$ ps. The black line with the symbol, 0,
corresponds to the state of water, 0.025 ps before the initiation of the TS.
The last curve presented in Fig. 14, corresponding to the symbol ($+$), was
sampled at $t=0.25$ ps after the generation of the TS.

The bin volume can be calculated by $L\times$ lateral area which is given by
$\pi(r_{2}^{2}-r_{1}^{2})=\pi(r_{2}-r_{1})(r_{2}+r_{1})=\pi(4{\rm
nm})(r_{2}+r_{1})\Delta r=251.2\AA^{2}\times(r_{2}+r_{1})$. Here $r_{1}$ and
$r_{2}$ are the inner and outer radius of the cylindrical bins. For small
$\Delta r=r_{2}-r_{1}$, we can approximate $(r_{2}+r_{1})$ to $2r$. As the bin
volume grows linearly with the cylinder radius, $r$, the number of particles
grows linearly, given the constant density of water (1g/cm3). Hence, as it is
visible in Fig. 14, the statistical fluctuations of the calculated variables
lower down with an increase in $r$.

Evidently, within ps’s time frame after the creation of the TS, during the
formation of the nano-cavity, our analysis of temperature and pressure
profiles reveals a novel internal periodic structure and an oscillatory
pattern between low and high temperatures and pressures of the gas inside the
cavity. An internal periodic “beating” structure of the gas motion, resembles
the standing normal modes in a cylindrical acoustic pipe with expanding
boundaries. After induction of a large pressure inside the TS, the sudden
expansion of the gas forms a local vacuum in the middle of the cavity. Because
the flow of the fluid goes along with a negative pressure gradient, the
positive gradient in pressure, in the middle of the cavity pushes the escaped
particles back to the center of the cavity. In this cycle, a rise in
temperature and pressure occurs once again.

Furthermore, because the direction of the gradient in pressure changes signs
at the interface of the cavity, the stability of the cavity can be guaranteed
due to the symmetry in pressure gradients with respect to the cavity
interface, as seen in Fig. 14. As a function of time, the curvature at the
pressure peak, on the cavity interface, lowers and the slope approaches zero.
It is however possible that internal fluctuations disturb the balance between
positive and negative pressure gradients. If the internal positive pressure
gradient prevails over the external negative pressure, the cavity collapses
suddenly, and a phenomenon similar to the sonoluminescence effect takes place.
These fluctuations can be driven by classical and/or quantum zero-point
oscillations of the electromagnetic field around the plasmon modes of the gas
polarization in the cavity that breaks the symmetry of pressure gradients
spontaneously, i.e., the Casimir effect [64]. Our hypothesis on the sudden
collapse of TS nano-cavity may be considered as an alternative mechanism that
explains the production of the Cherenkov light spectrum, with the emitted
light intensity $\propto 1/\lambda$ (where $\lambda$ is the wavelength of the
emitted light) during luminescence of water than prompt-$\gamma$ and emission
of light at lower energy than the Cherenkov light threshold hypothesized
recently in Ref. [65].

The transient state of water in the cavity at high temperatures and pressures
but low densities is a manifestation of super-critical phase of water [66].
Thus we may also interpret the observed water luminescence as thermal
radiation (e.g., black-body radiation, with gray-body corrections at longer
wavelength in the optically thick regime) from the transient supercritical
state of water localized in nanocavities wrapping around the track of charged
particles. We suggest the measurement of the spectrum, fitted to thermal
radiation, where the intensity of emitted light scales like $1/\lambda^{4}$
for large $\lambda$’s may determine the average temperature of the spikes as a
function of the charged particle penetration length, and LET.

###  III.8 Nano-bubble stability and inter-bubble interaction

To further examine the stability of the TSs discussed in the preceding
section, we consider two TSs created in the vicinity of each other. Upon their
expansion, they may collide and form ”molecular crowding” at their interface,
one of the hypothetical track configurations at UHDRs as shown in Figs. 15 and
16.

Fig. 15(a) shows a snapshot of two thermal spikes created simultaneously 10 nm
apart from each other at 2 ps after their creations, whereas in (b) there is a
0.5 ps time delay to their creation time. A time delay beyond a threshold,
which is a function of the TS’s inter-separation distance and radius, allows
penetration of the second TS into the first one. In contrast, two TS’s with
time delay below the threshold, collide like two rigid objects with a high-
density wall between them, like two stable bubbles touching each other at
their interface. The surface tension of the interface keeps the bubbles
impenetrable.

Fig. 17 shows the sequence of the selected events and how the time delay
allows their penetration. The circles in (a) are the initial location and
dimension of the TS’s. The radius of each of them is 1 nm and the separation
of 5 nm. In (b) the first thermo-acoustic waveforms with a radius of 5 nm when
the second TS enters the wall and expands (c) and finally they combine and
form a large cavity with a diameter of 20 nm.

The flow of water molecules through the interface to outside TS’s interfaces
looks like a fountain-like and/or a jet stream with a side-wise direction to
the outside of cavities. The direction of the flows as shown by the arrows in
Fig. 16 (d) forms a saddle point in the center. Such severe flow of water
molecules makes them chemically active after they break apart through their
covalent bonds losing electrons and leaving ions behind.

The creation of these two thermo-acoustic waves has been uploaded in the form
of a video [58] to help the readers visualize their time evolution and the
inter-track couplings. It shows the overall 5 ps running time.

Figure 15: Time-evolution of a system of two TS’s created simultaneously (a)
and with a time delay of 500 fs (b).  Figure 16: Time-evolution of a system of
two TS’s created simultaneously (a) and with a time delay of 500 fs (b).
Figure 17: Time-evolution of a system of two TS’s created with a time delay of
500 fs.

###  III.9 Summary and remarks

TS’s induced by charged particles generate a nano-scale cylindrical symmetric
gradient in the temperature profile along the center of the track structure.
Assuming a one-to-one correspondence among the thermodynamic variables,
specified by the equation of the state of the medium across the interface with
finite width, $\ell$, a temperature gradient,
$\vec{\nabla}T\approx(T_{2}-T_{1})/\ell$, produces a pressure gradient
$\vec{\nabla}P\approx(P_{2}-P_{1})/\ell$ normal to the cylindrical interface
with volume, $\pi R^{2}L$, with the positive direction toward outside of this
volume.

The location of the interface relative to the center of the TS, in the Lab
frame of reference, is a dynamical variable and changes as a function of time
with a characteristic rate given by $v_{1}$. This is the speed of unperturbed
fluid relative to TS surface of discontinuity, from the perspective of an
observer stationary relative to the track interface. In this system of
coordinates where the Lab frame is moving in the opposite direction, $v_{1}$
represents the speed of the interface. The location of the interface in the
Lab frame, i.e., the radius of the cylinder can be calculated by
$R(t)=\int_{0}^{t}dt^{\prime}v_{1}(t^{\prime})$.

Another important variable is the speed of sound in the unperturbed side of
fluid which is given by the following equation
$c_{s}=\sqrt{(C_{P}/C_{V})/(\rho\beta_{T})}=\sqrt{(C_{P}/C_{V})(P/\rho)/(P\beta_{T})}$.
Note that $P\beta_{T}=1$ for an ideal gas where
$c_{s}=\sqrt{\tilde{\gamma}P/\rho}$ and $\tilde{\gamma}=C_{P}/C_{V}$ with
numerical value corresponding to $c_{s}\approx 300m/s$ in air at the standard
condition. For incompressible fluids $P\beta_{T}<<1$. The inclusion of
$P\beta_{T}$ in the equation for the speed of sound corrects and matches
$c_{s}$ with their empirical values. For water at room temperature, our MD
calculation predicts the correct value of $c_{s}\approx 1450$ m/s.

For large enough gradients, $v_{1}$ may exceed the speed of sound. We denote
$M$, the Mach’s number of the unperturbed fluids. The supersonic flow of TS
interfaces expands like shock waves (SW) of explosives with $M=v_{1}/c_{s}>1$.

Our analysis has revealed that variabilities in the molecular mass and charge
distribution and in general chemical, physical, and mechanical composition of
the micro-environment of cells may dramatically alter the diffusion of ^.OH-
radicals within nano-scales.

##  IV Discussion and conclusion

In this work, we have demonstrated the coexistence of a rapidly growing
condensed state of water with a fast-growing hot spot, forming a stable state
of water mixed with free radicals and chemical species generated due to bond
breaking at high temperatures. The thin shell of a highly dense state of water
grows by three to five folds of the speed of sound in water, forming a SW
buffer, and wrapping around the nano-scale cylindrical symmetric cavity. The
cavity consists of low density of water molecules, hence it forms a bubble-
like state of water embedded in a room-temperature background of water.
Eventually, SW amplitude dissipates to a sound wave. We have also studied the
interaction of SWs and mapped the trajectory of water molecules into a jet
stream along a vector field and identified a stable saddle point normal to the
direction of SW propagation. We have demonstrated a random time-elapse
mechanism among the TS’s accounted for the nucleation, growth, and stability
of the large bubbles.

The formation of supersonic growth of nano-bubbles / nano-cavities has been
shown as an indication of a much higher diffusion constant currently predicted
by the reaction-diffusion of the MC models in radiation physics. The current
study is potentially significant in the accurate calculation of overlap
integrals among the radiation tracks at ultra-high dose rates and the data
analysis in FLASH radiotherapy that requires detailed knowledge of the spatial
and temporal distribution of the particles in microscopic-scale pulsed beams.
The more accurate simulation framework proposed hereby should help elucidate
the intriguing mechanisms involved in FLASH radiotherapy. It may be used to
define the dose fractionation requirements (i.e., instantaneous dose rate,
irradiation time or beam scanning parameters) to maintain a FLASH effect while
minimizing potential adverse effects [23]. Therefore, it may be a critical
step toward safe clinical deployment of FLASH radiotherapy.

Recent efforts on developing of ionizing radiation acoustic imaging (iRAI)
technique for real-time dosimetric measurements for FLASH radiotherapy [20] or
more broadly radiotherapy with pulsed ionizing sources may also benefit from
our simulation framework. Novel theoretical approaches may open up an in-depth
understanding of the first principle mechanisms consistent with the picture
explored in this work, i.e. starting from the formation of TSs and ending with
the acoustic pulses used to reconstruct the dose in iRAI. It may foster the
development of a new model-based reconstruction technique accounting for TS
and nano-bubble generation to retrieve information not only on the dose
distribution, as currently implemented, but also on the linear energy transfer
(LET). Hence, it may unlock the in vivo assessment of LET known to be more
biologically relevant [63] but which cannot be imaged with any methods
currently available.

The new knowledge may also help improving various imaging techniques designed
for the detection of charged particle tracks beyond radiotherapy, such as the
optical emission of photons a mechanism known as scintillation in organic,
inorganic, and in particular in liquid scintillation materials [24]. In
addition to bubble formation in water in nuclear plants, NMR and MRI imaging
beam topography [25, 26], and a many other applications in fundamental science
and industry [27, 28, 29].

Embedded in these results is a connection to the latest experiment on water
luminescence, reported in [65]. The authors have estimated the photon
intensity by the Cherenkov light production from prompt gamma photons using MC
simulations. However, in this paper, we suggest an alternative mechanism. We
propose that the indirect photon emission in the observed water luminescence
may result either from the thermal radiation of the hot gas bubbles or the
collapse of nano-bubbles because of strong fluctuations that drive the
internal negative pressure gradient dominating the opposite pressure gradient
outside of the bubble interface. One may speculate classical/quantum
fluctuations that force bubbles (induced by high-energy charged particles) to
collapse, to produce light, a mechanism similar to the Casimir effect [64] in
the context of the sonoluminescence phenomenon [28]. More characterization of
the spectrum is needed to identify which channel of luminescence is
responsible for the observation of the optical lights in the experiments [28].

In our model, for every charged particle passing through water, a nano-bubble
is formed. They statistically collapse which results in photon emission from
the location of the primary charged particles, hence their intensity is
expected to be greater at the high LET regions, near the Bragg peak, and
within the spots of low energy secondaries. This mechanism can explain the
tail of photon production from secondary protons beyond the range of the
primary beam of the proton, distal to Bragg peak as illustrated in Ref. [65].

Our simulations presented in this study have revealed a transient state of
water in the cavity at high temperatures and pressures but low densities are
formed within a radius of nano-meter from the charged particles. The
calculated mass density profile of water inside these nano-cavities, shown in
Figs. 12 and 14, matches with the thermodynamically stable state of
supercritical water [66]. We therefore interpret the water luminescence
observed and reported in Ref. [65] as a manifestation of the thermal radiation
from the transient supercritical state of water localized in nano-cavities
wrapping around the track of charged particles.

###### Acknowledgements.

The authors would like to thank useful scientific communications with Arash
Darafsheh, Alberto Fraile, Jorge Kohanoff, Eugene Surdutovich, and Andrey
Solov’yov. The authors gratefully acknowledge the computational and data
resources provided by the Leibniz Supercomputing Centre (www.lrz.de) and
support from George Dedes. NA is supported by the Natural Sciences and
Engineering Research Council of Canada, the University of Waterloo, and the
Perimeter Institute for Theoretical Physics. Research at Perimeter Institute
is supported in part by the Government of Canada through the Department of
Innovation, Science and Economic Development Canada and by the Province of
Ontario through the Ministry of Colleges and Universities. MR, KP, and JL
acknowledge the support from the Centre for Advanced Laser Applications and
LMUexcellent, funded by the Federal Ministry of Education and Research (BMBF)
and the Free State of Bavaria under the Excellence Strategy of the Federal
Government and the Länder.

Authors contributions: RA: wrote the main manuscript, prepared figures, and
performed mathematical derivations and computational steps. NA, SR, MR, RT,
JL: contributed to developing scientific ideas, computational analysis of the
problem, and writing the manuscript. AvD: provided ReaxFF, data analysis, and
wrote the manuscript. KP, JL, and RM: wrote the main manuscript, proposed the
scientific problem, , co-supervised the project, and provided critical
feedback.

Competing financial interest: The authors declare no competing financial
interests.

Corresponding Author:
†<EMAIL_ADDRESS>
## References

  * [1] A. Norman, P. Spiegler, Radiation Nucleation of Bubbles in Water, Nucl. Sci. Eng. 16, 213-217 (1963).
  * [2] M. Vedadi, A. Choubey, K. Nomura, R. K. Kalia, A. Nakano, P. Vashishta, and A. C. T. van Duin, Structure and Dynamics of Shock-Induced Nanobubble Collapse in Water, Phys. Rev. Lett. 105, 014503 (2010).
  * [3] S. H. Min, M. L. Berkowitz, Bubbles in water under stretch-induced cavitation, J. Chem. Phys. 150, 054501 (2019).
  * [4] M. Mallory, E. Gogineni, G. C. Jones, L. Greer, C. B. Simone II, Therapeutic hyperthermia: The old, the new, and the upcoming, Critic. Rev. Onc. Hemato., 97, 56-64 (2016).
  * [5] G.-Y. Wan, Y. Liu, B.-W. Chen, Y.-Y. Liu, Y.-S. Wang, N. Zhang, Recent advances of sonodynamic therapy in cancer treatment Cancer Biol Med 2016. doi: 10.20892/j.issn.2095-3941.2016.0068
  * [6] W. Yan, M. K. Khan, X. Wu, C. B. Simone, J. Fan, E. Gressen, X. Zhang, C. L. Limoli, H. Bahig, S. Tubin, and W. F. Mourad, Spatially fractionated radiation therapy: History, present and the future, Clin. Transl. Radiat. Oncol. 20, 30–38 (2020).
  * [7] Favaudon V, Caplier L, Monceau V, et al. Ultrahigh dose-rate FLASH irradiation increases the differential response between normal and tumor tissue in mice. Sci Transl Med. 2014;6:1-9.
  * [8] Montay-Gruel P, Bouchet A, Jaccard M, et al. X-rays can trigger the FLASH effect: Ultra-high dose-rate synchrotron light source prevents normal brain injury after whole brain irradiation in mice. Radiother Oncol. 2018;129(3):582-588.
  * [9] Vozenin MC, De Fornel P, Petersson K, et al. The advantage of FLASH radiotherapy confirmed in mini-pig and catcancer patients. Clin Cancer Res 2019 Jan 1;25(1):35-42.
  * [10] Montay-Gruel P, Acharya MM, Petersson K, et al. Long-term neurocognitive benefits of FLASH radiotherapy driven by reduced reactive oxygen species. Proc Natl Acad Sci USA. 2019;116(22):10943-10951.
  * [11] Buonanno M, Grilj V, Brenner DJ. Biological effects in normal cells exposed to FLASH dose rate protons. Radiother Oncol. 2019;139:51-55.
  * [12] Vozenin MC, Baumann M, Coppes RP, Bourhis J. FLASH radiotherapy international workshop. Radiother Oncol. 2019;139:1-3.
  * [13] Darafsheh A, Hao Y, Zwart T, Wagner M, Catanzano D, Williamson JF, Knutson N, Sun B, Mutic S, Zhao T. Feasibility of proton FLASH irradiation using a synchrocyclotron for preclinical studies. Med Phys. 2020; doi: 10.1002/mp.14253.
  * [14] Spitz DR, Buettner GR, Petronek MS, et al. An integrated physico-chemical approach for explaining the differential impact of FLASH versus conventional dose rate irradiation on cancer and normal tissue responses. Radiother Oncol. 2019; 139:23-27.
  * [15] Koch CJ. Re: Differential impact of FLASH versus conventional dose rate irradiation. Radiother Oncol. 2019;139:62-63.
  * [16] R. Abolfath, D. Grosshans, R. Mohan, Oxygen depletion in FLASH ultra-high-dose-rate radiotherapy: A molecular dynamics simulation, Med. Phys. 47, 6551-6561 (2020).
  * [17] J. Jansen, J. Knoll, E. Beyreuther, J. Pawelke, R. Skuza, R. Hanley, S. Brons, F. Pagliari, J. Seco, Does FLASH deplete oxygen? Experimental evaluation for photons, protons, and carbon ions, Med. Phys. 48, 3982 (2021).
  * [18] R. Abolfath, A. Baikalov, S. Bartzsch, E. Schuler, R. Mohan, A stochastic reaction-diffusion modeling investigation of FLASH ultra-high dose rate response in different tissues, Special issue on Multidisciplinary Approaches to the FLASH radiotherapy, Front. Phys., 11, 1060910 (2023).
  * [19] A. Baikalov, R. Abolfath, E. Schuler, R. Mohan, Jan Wilkens, S. Bartzsch, A Intertrack interaction at ultra-high dose rates and its role in the FLASH effect, Front. Phys. 11, 1215422 (2023).
  * [20] I. Oraiqat, W. Zhang, D. Litzenberg, K. Lam, N. B. Sunbul, J. Moran, K. Cuneo, P. Carson, X. Wang, I. E. Naqa, An ionizing radiation acoustic imaging (iRAI) technique for real-time dosimetric measurements for FLASH radiotherapy, Med. Phys. 47, 5090-5101 (2020).
  * [21] A.I. Kalinichenko, V.T. Lazurik, and I.I. Zalyubovsky, The Physics and Technology of Particle and Photon Beams Volume 9: Introduction to Radiation Acoustics. The Netherlands: Harwood Academic Publishers, 2001.
  * [22] S. Hickling, L. Xiang, K. C. Jones, K. Parodi, W. Assmann, S. Avery, M. Hobson, I. El Naqa, Ionizing radiation-induced acoustics, Med. Phys. 45, e707 (2018).
  * [23] J. Lascaud and K. Parodi, On the potential biological impact of radiation-induced acoustic emissions during ultra-high dose rate electron radiotherapy: a preliminary study, Phys. Med. Biol. 68, 05LT01 (2023).
  * [24] A. S. Beddar, and L Beaulieu, Scintillation Dosimetry. Boca Raton, FL: CRC Press (2016).
  * [25] H. Odeen, D. L. Parker, Magnetic resonance thermometry and its biological applications – Physical principles and practical considerations, Prog. Nucl. Magn. Reson. Spectrosc. 110 34–61 (2019).
  * [26] S. Gantza, L. Karscha, J. Pawelkea , J. Petera, S. Schellhammer, J. Smeets, Erik van der Kraaij, A. Hoffmann, Direct visualization of proton beam irradiation effects in liquids by MRI, PNAS 120, 23 e2301160120 (2023).
  * [27] G. E. Winter, C. M. Cooling, M. D. Eaton, Linear energy transfer of fission fragments of 235U and nucleation of gas bubbles in aqueous solutions of uranyl nitrate Annals of Nuclear Energy 142, 107379 (2020).
  * [28] M. P. Brenner, S. Hilgenfeldt, D. Lohse, Single-bubble sonoluminescence, Rev. Mod. Phys., 74, 2, (2002).
  * [29] M. Durante, F. A. Cucinotta, Physical basis of radiation protection in space travel, Rev. Mod. Phys. 83, 1245 (2011).
  * [30] R. Abolfath, A. Baikalov, S. Bartzsch, N. Afshordi, and R. Mohan, “The effect of non-ionizing excitations on the diffusion of ion species and inter-track correlations in flash ultra-high dose rate radiotherapy,” Phys. Med. Biol. 67, 105 (2022).
  * [31] A. M. Kellerer, Fundamentals of Microdosimetry the Dosimetry of Ionizing Radiation. Vol. 1, Kase KR, et al. London: Academic, 77–161 (1985).
  * [32] J. Fain, M. Monnin, M. Montret, Spatial Energy Distribution Around Heavy-Ion Path, Radiat. Res. 57, 379-389 (1974).
  * [33] H. Wang, O. N Vassiliev, Radial dose distributions from protons of therapeutic energies calculated with Geant4-DNA, Phys. Med. Biol. 59 (2014) 3657–3668.
  * [34] J. A. LaVerne, Track Effects of Heavy Ions in Liquid Water, Radiat. Res. 153, 487–496 (2000).
  * [35] S. Agostinelli et. al., Nucl. Instrum. Meth. A 506, 250 (2003).
  * [36] Incerti S. et. al., The GEANT4-DNA project, Int. J. Modelling Simul. Sci. Comput. 1, 157–78 (2010).
  * [37] van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A., ReaxFF: A Reactive Force Field for Hydrocarbons, J. Phys. Chem. A 2001, 105, 9396.
  * [38] M Yusupov, E C Neyts, U Khalilov, R Snoeckx, A C T van Duin, and A Bogaerts, Atomic-scale simulations of reactive oxygen plasma species interacting with bacterial cell walls, New Journal of Physics 14 (2012) 093043.
  * [39] Verlackt, C.C., Neyts, E.C., Jacob, T., Fantauzzi, D., Golkaram, M., Shin, Y.K., van Duin, A.C.T. and Bogaerts, A. (2015) Atomic scale pathways in DNA oxidation by hydroxyl radicals: a reactive molecular dynamics study for plasma oncology. New Journal of Physics 17, 103005-103001.
  * [40] R. M. Abolfath, A. C. T. van Duin, T. Brabec, it Reactive molecular dynamics study on the first steps of DNA damage by free hydroxyl radicals, J. Phys. Chem. A 115, 11045 (2011). See the real-time simulations and movies at: http://qmsimulator.wordpress.com/
  * [41] P. de Vera, E. Surdutovich, A. V. Solov’yov, The role of shock waves on the biodamage induced by ion beam radiation, Cancer Nano 10, 5, (2019). https://doi.org/10.1186/s12645-019-0050-3
  * [42] P. de Vera, E. Surdutovich, N. J. Mason, F. J. Currell, A. V. Solov’yov, Simulation of the ion-induced shock waves effects on the transport of chemically reactive species in ion tracks Eur. Phys. J. D 72, 147 (2018).
  * [43] E. Surdutovich, A. V. Solov’yov, Shock wave initiated by an ion passing through liquid water, Phys Rev E. 82, 051915 (2010).
  * [44] E. Surdutovich, A.V. Solov’yov, Multiscale approach to the physics of radiation damage with ions, Eur. Phys. J. D 68, 353 (2014).
  * [45] A. Fraile, M. Smyth, J. Kohanoff, A. V. Solov’yov, First principles simulation of damage to solvated nucleotides due to shock waves, J. Chem. Phys. 150, 015101 (2019); https://doi.org/10.1063/1.5028451
  * [46] Landau L, Lifshitz E. Fluid dynamics, vol. 6. 2nd ed. Oxford: Reed-Elsevier; 1987.
  * [47] K. Huang, Statistical Mechanics, 2Ed (Wiley, 1987).
  * [48] T. W. Lion, R. J. Allen, Computing the local pressure in molecular dynamics simulations, J. Phys.: Condens. Matter 24, 284133 (2012).
  * [49] D. Chandler, Introduction to modern statistical, Mechanics. Oxford University Press, Oxford, UK (1987).
  * [50] M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press (1989).
  * [51] L.D. Landau, and E.M. Lifshitz, Statistical Physics. Part 1, third ed., Pergamon Press, Oxford, UK (1980).
  * [52] Sedov LI. Prikl Mat Mek (SSSR). 10, 241 (1946).
  * [53] von Neumann J. In: Fuchs K, Hirschfelder JO, Magee JL, Peierls R, von Neumann J, editors. The point source solution. Blast wave: Los Alamos; 1947.
  * [54] Zel’dovich Y, Raiser Y. Physics of shock waves and high-temperature hydrodynamic phenomena, vol I. New York: Oxford; 1996.
  * [55] R. Abolfath, Lipid Bilayer destruction by ionizing radiation, TS’s and shock waves, available online: https://www.youtube.com/watch?v=0GYgT8sIpUk
  * [56] R. Abolfath, Nano-bubbles / cavities induced by charged particles simulated by MD-ReaxFF with $\approx$ 150,000 water molecules. The running time up to 2 ps shows the cavity expands up to 10 nm in diameter. The video is available online: https://www.youtube.com/watch?v=Sh$\\_$6s34WzNE
  * [57] S. Ljunggren, J. C. Eriksson, The lifetime of a colloid-sized gas bubble in water and the cause of the hydrophobic attraction Colloids and Surfaces A: Physicochemical and Engineering Aspects 129-130, 151-155 (1997).
  * [58] R. Abolfath, Molecular crowding, and inter-track interaction The video is available online: https://www.youtube.com/watch?v=lTGNbH2npME
  * [59] https://webbook.nist.gov/cgi/cbook.cgi?ID=C7732185
&Units=CAL&Type=JANAFL&Table=on

  * [60] V. P. Singh, A. K. Madan, H. R. Suneja, D. Chand, Propagation of spherical shock waves in water, Prec. Indian Acad. Sci. (Engg. Sei.) 3, 169-175 (1980).
  * [61] J. D. Jackson, Classical Electrodynamics, Wiley Third Edition (1999).
  * [62] R. Abolfath et al., M. Rädler et al. (unpublished).
  * [63] Abolfath, R.; Helo, Y.; Bronk, L.; Carabe, A.; Grosshans, D.; Mohan, R. Renormalization of radiobiological response functions by energy loss fluctuations and complexities in chromosome aberration induction: Deactivation theory for proton therapy from cells to tumor control. _Eur. Phys. J. D_ 2019, 73, 64.
  * [64] S. Lamoreaux, The Casimir force: Background, experiments, and applications. Reports on Progress in Physics. 68, 201–236 (2005).
  * [65] S. Yamamoto, Discovery of the luminescence of water during irradiation of radiation at a lower energy than the Cherenkov light threshold, Radiol. Phys. Technol. 14, 16–24 (2021).
  * [66] H. Zheng, T. Yu, C. Qu, W. Li, Y. Wang, Basic Characteristics and Application Progress of Supercritical Water, Earth and Environmental Science 555, 012036 (2020).