Modeling emission of acoustic energy during bubble expansion in PICO bubble chambers

TL;DR

This paper develops a physical model to understand and improve acoustic discrimination between alpha particles and nuclear recoils in bubble chambers used for dark matter detection, enhancing background rejection.

Contribution

It introduces a novel physical model of bubble growth dynamics that explains acoustic signals for the first time in PICO bubble chambers.

Findings

The model accurately predicts the acoustic energy distribution for alpha and recoil events.

Simulation results match experimental data, validating the model.

Improved discrimination enhances dark matter search sensitivity.

Abstract

The PICO experiment uses bubble chambers filled with superheated C$_3$F$_8$ for spin-dependent WIMP dark matter searches. One of the main sources of background in these detectors is alpha particles from decays of environmental $^{222}\mathrm{Rn}$, which nucleate bubbles that are visually indistinguishable from WIMP candidate events. Alpha-induced bubbles can be discriminated acoustically, because the signal from alpha events is consistently larger in magnitude than that from nuclear recoil/WIMP-like events. By studying the dynamics of bubbles nucleated by these two types of ionizing radiation from the first stages of their growth, we present a physical model for the acoustic discrimination for the first time. The distribution of acoustic energies that we generate for a simulated sample of bubble nucleations by alpha particles and nuclear recoils is compared directly to the experimental data.

Tables

Table 1. TABLE 1: Select thermodynamic properties of C 3 F 8 at 14   °C times 14 celsius 14\text{\,}\mathrm{\SIUnitSymbolCelsius} and 207   kPa times 207 kPa 207\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}  ( 30   psia times 30 psia 30\text{\,}\mathrm{p}\mathrm{s}\mathrm{i}\mathrm{a} ) extracted from REFPROP Lemmon et al. ( 2018 ) . When appropriate, we also list the respective values found for the C 3 F 8 -like LJ fluid and use a dagger † to denote these cases.

Property	Value	Unit
$h$: enthalpy of vaporization	$82.58$	$\mathrm{kJ}{\mathrm{kg}}^{-1}$
$c_l$: liquid heat capacity	$1130.5$	$J{\mathrm{kg}}^{-1}{K}^{-1}$
$a_l$: liquid thermal diffusivity	$0.030$	${\mathrm{mm}}^2{s}^{-1}$
${\rho }_v$: vapor density	$59.4$	$\mathrm{kg}{m}^{-3}$
${\rho }_l$: liquid density	$1379.0$	$\mathrm{kg}{m}^{-3}$
$T_{\mathrm{sat}}$: saturation temperature	253.57	K
$A$ (Eq. 2)	$17.7$	$m{s}^{-1}$
$B$ (Eq. 3)	$0.0036$	$m{s}^{-1/2}$
	$0.0086\pm {0.0002}^{†}$	$m{s}^{-1/2}$
${\nu }_l$: kinematic viscosity	$0.141$	${\mathrm{mm}}^2{s}^{-1}$
	$18.3\pm {0.4}^{†}$	${\mathrm{mm}}^2{s}^{-1}$
$\gamma$: surface tension	$4.9$	$\mathrm{mN}{m}^{-1}$
	$3.667\pm {0.007}^{†}$	$\mathrm{mN}{m}^{-1}$

Table 2. TABLE 2: Estimates of the time t s subscript 𝑡 s t_{\mathrm{s}} when α 𝛼 \alpha -induced bubbles become spherical and their radii R s subscript 𝑅 s R_{\mathrm{s}} at that moment. The values of t s subscript 𝑡 s t_{\mathrm{s}} and R s subscript 𝑅 s R_{\mathrm{s}} are evaluated with Eqs.   5 and  6 respectively for alphas from the Rn 222 superscript Rn 222 {}^{222}\mathrm{Rn} decay chain ( II.4 ). The reported l cyl subscript 𝑙 cyl l_{\mathrm{cyl}} values correspond to the mean track ranges of α 𝛼 \alpha particles in the superheated C 3 F 8 ( ρ l = 1379 ​ kg ​ m − 3 subscript 𝜌 l 1379 kg superscript m 3 \rho_{\mathrm{l}}=1379\,\mathrm{kg\,m^{-3}} ), projected with SRIM Ziegler et al. ( 2010 ) for each E α subscript 𝐸 𝛼 E_{\alpha} .

$E_{\alpha }\left[\mathrm{MeV}\right]$	$l_{\mathrm{cyl}}\left[\mu m\right]$	$R_s\left[\mu m\right]$	$t_s\left[\mu s\right]$
5.6	42.6	$16.2\pm 1.2$	$3.5\pm 0.3$
6.1	48.3	$18.4\pm 1.4$	$4.0\pm 0.3$
7.8	70.0	$26.6\pm 2.1$	$5.7\pm 0.4$

Table 3. TABLE 3: Liquid slab experiment results: MD measurements of surface tension and liquid/vapor densities at saturation ( T = 14   °C 𝑇 times 14 celsius T=$14\text{\,}\mathrm{\SIUnitSymbolCelsius}$ ) in LJ units as compared to REFPROP values in SI units. For each quantity, the last column shows σ 𝜎 \sigma that makes the measured LJ value match the REFPROP-predicted value. The uncertainties on REFPROP values are derived from the uncertainties on simulation-extracted temperatures for each phase ( T l = 13.95 ± 0.05 ∘ ​ C subscript 𝑇 l plus-or-minus 13.95 superscript 0.05 C T_{\mathrm{l}}=13.95\pm 0.05\,^{\circ}\mathrm{C} , T v = 14.4 ± 0.5 ∘ ​ C subscript 𝑇 v plus-or-minus 14.4 superscript 0.5 C T_{\mathrm{v}}=14.4\pm 0.5\,^{\circ}\mathrm{C} ).

	MD simulation	REFPROP	$\sigma \left[\mathrm{nm}\right]$
$\gamma$	$0.201\pm 0.001\frac{ϵ}{{\sigma }^2}$	$4.82\pm 0.05\frac{\mathrm{mN}}{m}$	$0.461\pm 0.003$
${\rho }_l$	$0.6715\pm 0.0001\frac{m}{{\sigma }^3}$	$1385.1\pm 0.3\frac{\mathrm{kg}}{m^3}$	$0.53293\pm 0.00005$
${\rho }_v$	$0.03983\pm 0.00009\frac{m}{{\sigma }^3}$	$62\pm 1\frac{\mathrm{kg}}{m^3}$	$0.586\pm 0.003$
$P_{\mathrm{eq}}$	$0.02240\pm 0.00003\frac{ϵ}{{\sigma }^3}$	$6.5\pm 0.1\mathrm{bar}$	$0.560\pm 0.003$

Equations

$$\begin{array}{rl}\frac{\mathrm{d}R}{\mathrm{d}t}=-\left[\frac{A^2\sqrt{t-t_{\mathrm{s}}}}{B}+\frac{2{\nu }_{\mathrm{l}}}{R}\right]& \\ & +\sqrt{A^2-\frac{2\gamma }{{\rho }_{\mathrm{l}}R}+{\left(\frac{2{\nu }_{\mathrm{l}}}{R}+\frac{A^2\sqrt{t-t_{\mathrm{s}}}}{B}\right)}^2}.\end{array}$$

$$A=\sqrt{\frac{2}{3}\frac{h{\rho }_{\mathrm{v}}\Delta T}{{\rho }_{\mathrm{l}}{T}_{\mathrm{sat}}}},$$

$$B=\sqrt{\frac{12}{\pi a_{\mathrm{l}}}}\text{Ja},$$

$$\tilde{e}=\sqrt{\left|1-\frac{R_{90}^2}{R_0^2}\right|},$$

$$t_{\mathrm{s}}\text{leavevmode}[\mathrm{ns}]\approx (0.082\pm 0.006)l_{\mathrm{cyl}}\text{leavevmode}[\mathrm{nm}]$$

$$\begin{array}{rl}{R}_{\mathrm{s}}\text{leavevmode}[\mathrm{nm}]\approx (0.38\pm 0.03)l_{\mathrm{cyl}}[\mathrm{nm}]& \\ & +(217\pm 27)\frac{E_{\mathrm{dep}}}{l_{\mathrm{cyl}}}\left[\frac{\mathrm{keV}}{\mathrm{nm}}\right].\end{array}$$

$${}^{222}\mathrm{Rn}\to {}^{218}\mathrm{Po}+\alpha +5.6\mathrm{MeV};$$

$${}^{218}\mathrm{Po}\to {}^{214}\mathrm{Po}+\alpha +6.1\mathrm{MeV};$$

$${}^{214}\mathrm{Po}\to {}^{210}\mathrm{Pb}+\alpha +7.8\mathrm{MeV}.$$

$$\mathcal{P}(t)=\frac{{\rho }_{\mathrm{l}}{\ddot{V}}^2}{4\pi c},$$

$$E_{\mathrm{ac}}($1300\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$)\equiv\int_{1\,\mathrm{kHz}}^{300\,\mathrm{kHz}}\left|\mathcal{A}(f)\right|^{2}\mathrm{d}f,$$

Full text

I Introduction

Bubble chambers played an important role in experimental advances of particle physics in the 1950s and 1960s Glaser and Rahm (1955). Initially, they were operated at very high superheats to allow reliable bubble nucleation by minimally ionizing particles. Over the last decade bubble chambers have made a resurgence in dark matter searches, where the chambers are operated at more moderate superheats that explicitly prevent nucleation from minimally ionizing particles. This way, the chambers are only sensitive to nuclear recoils and other highly ionizing particle interactions.

The PICO experiment searches for WIMP dark matter candidates using bubble chamber technology Amole et al. (2017, 2019a); Mitra (2018). To calibrate the detector, neutron radiation from a source of known activity, typically ${}^{252}\mathrm{Cf}$ or ${}^{241}\mathrm{Am}/^{9}\mathrm{Be}$, is used. Neutron source calibrations allow the observation of nuclear recoil events in the range of keV recoil energies, creating an acoustic signature that is identical to the expected recoil signal from dark matter interactions.

Alpha radiation from decays of environmental ${}^{222}\mathrm{Rn}$ is one of the dominant sources of background. Visually, bubbles nucleated by $\alpha$ particles are indistinguishable from nuclear recoil/WIMP-like events. It was discovered in the PICASSO experiment that the acoustic signature of a nucleation from an $\alpha$ event is more powerful than the acoustic signal of a nuclear recoil event Aubin et al. (2008). So far this has been explained with a model that lacked specific predictability. To reveal the origin of the observed differences in acoustic signal magnitudes, we connect the initial thermodynamic conditions of bubble nucleation by $\alpha$ particles and nuclear recoils with the specifics of bubble evolution. We approach this problem with a molecular dynamics simulation, which lets us vary the energy deposited by an incident ion and the shape of the initially vaporized region to reproduce $\alpha$-like and nuclear recoil-like scenarios. We further study the influence of these variables on the dynamics of a nucleated nonspherical bubble and analytically link the first expansion stages to the well-described spherical bubble growth Plesset and Zwick (1954); Rayleigh (1917). Finally, we reconstruct the acoustic pressure signal from the spherical bubble growth histories for both $\alpha$-induced and nuclear recoil-induced bubbles. The results of the developed model are compared to recent acoustic data from the PICO-60 detector Amole et al. (2017, 2019a), which used superheated C3F8 as the target fluid.

II Inputs to the acoustic emission model

II.1 Molecular dynamics (MD) simulations

To model the nucleation and subsequent expansion of bubbles in superheated C3F8, we make use of the Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) Plimpton (1995). Inside this package, we represent C3F8 as a system of atoms interacting according to the Lennard-Jones (LJ) potential with characteristic energy $\epsilon$ and characteristic distance $\sigma$ (see Appendix A for details). This method neglects interatomic interactions within the molecules themselves, following the approach of Denzel et al. (2016) to model C2ClF5 used in the SIMPLE experiment Felizardo et al. (2012). The Lennard-Jones potential description of any fluid is complete when the $\epsilon$ and $\sigma$ parameters are specified. Setting $\epsilon=$0.0318\text{,}\mathrm{e}\mathrm{V} makes the critical temperature of an LJ fluid match that of C3F8. We also aim to simulate the specific conditions of the PICO-60 run at temperature $T_{0}\simeq$14\text{\,}\mathrm{\SIUnitSymbolCelsius} and liquid pressure $P_{\mathrm{l}}\simeq$207\text{,}\mathrm{k}\mathrm{P}\mathrm{a}$\,($30\text{,}\mathrm{p}\mathrm{s}\mathrm{i}\mathrm{a}$)$. This corresponds to the C3F8 liquid density of $\rho_{\mathrm{l}}=$1379\text{,}\mathrm{k}\mathrm{g}$\,\mathrm{m^{-3}}$, which we can reproduce with $\sigma=$0.533\text{,}\mathrm{n}\mathrm{m}$$. At $14\text{\,}\mathrm{\SIUnitSymbolCelsius}$, such a choice results in a $16\%$ discrepancy between the REFPROP Lemmon et al. (2018) value of equilibrium pressure at saturation and the Lennard-Jones value measured from our MD simulations; $33\%$ between REFPROP-extracted and simulated vapor densities; and $25\%$ between the surface tension values. Appendix A summarizes the steps we took to arrive at these results.

To bring the LJ fluid approximating C3F8 to the specified superheated state, we follow the procedure described in Denzel et al. (2016) and outlined in Appendix B. As in Denzel et al. (2016), we proceed with simulating the heat spike by depositing energy $E_{\mathrm{dep}}$ within a narrow and long cylindrical region in the metastable liquid, mimicking the track of an ionizing particle. Practically, this is done by increasing the velocities of the atoms falling within the cylindrical region to the corresponding temperature. According to the Seitz model Seitz (1958), our chosen thermodynamic conditions require a minimum of $\mathcal{T}_{\mathrm{min}}=3.3\,\mathrm{keV}$ to be deposited within the critical bubble radius $R_{\mathrm{c}}$ for the nucleated bubble to continue growing stably. To imitate nuclear recoils, we perform several simulations with the cylinder length $l_{\mathrm{cyl}}$ equal to integer multiples of $R_{\mathrm{c}}$ ($40.58\sigma\approx$ $21.63\text{\,}\mathrm{n}\mathrm{m}$ at our conditions) and energy deposits close to $\mathcal{T}_{\mathrm{min}}$. The radius of the cylindrical energy deposition region is fixed at $r_{\mathrm{cyl}}=2\sigma\approx 1.066\,$\mathrm{n}\mathrm{m}. To simulate $\alpha$-like bubble nucleations, we deposit energies of order $10\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ along $\text{\,}\mathrm{\SIUnitSymbolMicro m}$-long cylindrical tracks, with $r_{\mathrm{cyl}}=3\sigma\approx 1.599\,$\mathrm{n}\mathrm{m}. We justify these geometry choices in Appendix B.

The simulated bubbles are allowed to grow under NPT (constant particle number, pressure, and temperature) constraints through several major stages. These include rapid initial bubble expansion, followed by slight decrease in the bubble volume over a short time interval, and subsequent transition to the inertial growth phase. In this phase, the effective bubble radius $R_{\mathrm{eff}}$ increases nearly linearly with time, being controlled by the forces of liquid inertia. When the bubble, which grew out of a long cylindrical vapor region, becomes spherical, we may predict its further behavior by solving the differential equation describing spherical bubble growth Plesset and Zwick (1954); Rayleigh (1917). This technique allows us to simulate only the first few tens of nanoseconds of bubble expansion in LAMMPS, feeding the MD-based $R(t)$ history into the familiar expressions for $R$ at later times. In Sec. II.4, we show that $\alpha$-induced bubbles transition from nonspherical to spherical in a similar fashion, meaning that the same approach is applicable to predicting their growth.

II.2 Post-simulation processing: Bubble surface tracking

From each molecular dynamics simulation, we save 2D distributions of liquid densities in the XY plane running through the middle of the simulation volume (see Appendix C for details). With energy deposited along the $x$-axis, the XY plane fully captures the 3D bubble growth, given that the $z$-axis is physically no different from the $y$-axis. These density distributions are then used as an input to the marching squares algorithm Lorensen et al. (1987); Walt et al. (2014), which allows us to extract the contours of constant density values. As per Appendix C, the longest contour of $\rho_{\mathrm{l}}\approx$830\text{,}\mathrm{k}\mathrm{g}$\,\mathrm{m^{-3}}$ accurately reproduces the liquid-vapor boundary at each $t$. The $(x,y)$ coordinates of this boundary let us evaluate the bubble radius as a function of polar angle $\phi\equiv\tan^{-1}(y/x)$. In the present work, we take advantage of this to track the eccentricity of bubbles and thereby define the time elapsed after the heat spike when each bubble can be considered nearly spherical. As discussed in Sec. II.3, specifying this time and the corresponding value of the bubble radius makes for an initial condition sufficient to predict the subsequent bubble evolution. In Fig. 1, we give an example of bubble boundary extraction results at several representative times after the heat spike.

II.3 Predicting evolution of a spherical bubble in the linear and thermal growth phases

Molecular dynamics simulations, even for liquid volumes as small as those considered so far, are quite computationally expensive: for instance, \sim$$1\text{\,}\mathrm{n}\mathrm{s} of a typical bubble evolution in an NPT ensemble takes about 720 CPU-hours to run in LAMMPS. At the same time, radiation of acoustic energy as the result of bubble growth spans much longer time scales. These are roughly defined by the time of transition from the inertial bubble growth phase to the thermal growth phase, which features much slower expansion and consequently smaller amplitudes of the emitted sound waves Landau and Lifshitz (2013). For a bubble growing in superheated C3F8, the characteristic inertial growth phase duration is of order $100\text{\,}\mathrm{\SIUnitSymbolMicro s}$ Mitra (2018). This makes simulating the whole acoustic emission practically impossible with molecular dynamics methods alone. Such a complication calls for a way to fully predict bubble evolution based on a short history of its growth. The method has to be further generalized to the case of an arbitrary energy deposition, letting us perform MD simulations for only a few different $E_{\mathrm{dep}}$ values and extrapolate the results to other energies.

To proceed, we need to make a connection between bubble nucleation conditions, growth of the resulting highly nonspherical bubbles, and well-described spherical bubble dynamics. By combining the Rayleigh-Plesset equation Rayleigh (1917) for the inertial growth phase and the Plesset-Zwick equation Plesset and Zwick (1954) for the thermal growth phase as in Mikic et al. (1970), we arrive at the bubble growth rate

[TABLE]

In Eq. 1, the constant $A$ is a characteristic speed of expansion in the inertial phase. It is given by

[TABLE]

where $h$ is the latent heat of vaporization, $\Delta T$ is the liquid superheat above the saturation temperature $T_{\mathrm{sat}}$ for a given liquid density $\rho_{\mathrm{l}}$, and $\rho_{\mathrm{v}}$ is the vapor density inside the bubble. The constant $B$ in (1) may be expressed as

[TABLE]

where $a_{\mathrm{l}}$ is the thermal diffusivity of the liquid and Ja is its Jakob number ($\textit{Ja}=\frac{\Delta Tc_{\mathrm{l}}\rho_{\mathrm{l}}}{h\rho_{\mathrm{v}}}$, with $c_{\mathrm{l}}$ being the liquid heat capacity.) The other quantities appearing in Eq. 1 are the liquid viscosity $\nu_{\mathrm{l}}$ and the surface tension $\gamma$, whose contributions are non-negligible for the small bubble sizes simulated in LAMMPS. In Table 1, we provide the REFPROP values of these properties for C3F8 at $14\text{\,}\mathrm{\SIUnitSymbolCelsius}$ and $207\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}$ ($30\text{\,}\mathrm{p}\mathrm{s}\mathrm{i}\mathrm{a}$), as well as the $A$ and $B$ constants evaluated according to Eqs. 2 and 3. For our C3F8-like LJ fluid, we may also find $A,B,\gamma$, and $\nu_{\mathrm{l}}$ by fitting that part of the simulation-extracted $R(t)$ data where the bubble is spherical. To determine when this happens, we examine the time evolution of eccentricity, defined as

[TABLE]

with $R_{0}$ being the bubble radius in the direction of energy deposition ($\phi=0^{\circ}$) and $R_{90}$ being that in the orthogonal direction ($\phi$ = 90*∘*). The bubble eccentricity as a function of time corresponding to $R(t)$ from Fig. 2a is shown in Fig. 2b. We choose to use the first minimum in bubble eccentricity as an approximation for the time $t_{\mathrm{s}}$ required for the bubble to become spherical. Each $t_{\mathrm{s}}$ corresponds to a certain bubble radius value $R_{\mathrm{s}}$. $R(t_{\mathrm{s}})=R_{\mathrm{s}}$ serves as the initial condition for Eq. 1, which we fit simultaneously to the four simulated data sets as described in Appendix D. With parameter $A$ fixed at its value predicted for C3F8 with Eq. 2, we find a satisfactory agreement between the best-fitting $B$ and $\gamma$ and their C3F8 values (see Table 1). The viscosity of C3F8 differs from the value we extract for our LJ fluid by more than two orders of magnitude. This is a natural result since the Lennard-Jones potential completely neglects the intramolecular interactions within C3F8 and cannot reproduce all of its thermodynamic properties at once. We choose to proceed with the simulation-extracted rather than REFPROP parameter values in order to develop a self-consistent description of bubble expansion in the C3F8-like Lennard-Jones fluid. We stress that the bubble dynamics and acoustic emission model is therefore constructed for such an LJ fluid and not C3F8 itself. However, we find that replacing the best-fitting LJ viscosity value with the REFPROP value for C3F8 would eventually result in $\mathcal{O}(1)$ scaling of the integrated acoustic power in our frequency range of interest ($1300\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$; see Sec. II.5) for both nuclear recoil- and $\alpha$-like bubble nucleations. Since both types of simulated particle interaction are affected in the same way, the relative behavior between the two cases remains well described.

With all parameters in Eq. 1 determined for our LJ fluid, we can predict the evolution of bubbles nucleated within it for all $t\geq t_{\mathrm{s}}$. In Fig. 3, we give an example of such a prediction for a bubble nucleated from $E_{\mathrm{dep}}=$5\text{,}\mathrm{k}\mathrm{e}\mathrm{V}$$ deposited along $l_{\mathrm{cyl}}=2R_{\mathrm{c}}$. However, it remains to be shown how exactly the initial condition $R(t_{\mathrm{s}})=R_{\mathrm{s}}$ varies with $E_{\mathrm{dep}}$ and $l_{\mathrm{cyl}}$. As per Appendix D, we find

[TABLE]

and

[TABLE]

Together, Eqs. 5 and 6 provide the initial condition for Eq. 1 given deposited energies and lengths of ion tracks, which are known for both nuclear recoils and $\alpha$ particles.

II.4 Predicting the growth of $\alpha$-induced bubbles

Three populations of alpha-decay events have been observed in the PICO bubble chambers:

[TABLE]

These events belong to the decay chain of environmental ${}^{222}\mathrm{Rn}$, which creates a background of $\alpha$ particles depositing their energy within the superheated C3F8. Once bubbles nucleated as the result of such an energy deposition grow to a visible size, they cannot be distinguished from neutron-induced bubbles based on visual appearance. Instead, the acoustic parameter (AP) has been introduced as a measure of total acoustic power radiated during bubble expansion Amole et al. (2017). In the 1-300 kHz frequency range, the experimental AP provides excellent discrimination against $\alpha$-induced bubbles, which have 5-20 times larger AP values than those measured for nuclear recoil events, depending on the $\alpha$ energy. The dynamics of $\alpha$-induced bubbles has not been fully investigated before, which led to the lack of understanding of the origin of the $\alpha$ signal amplitude. Here, we will examine the dynamics of the bubbles nucleated by the ${}^{222}\mathrm{Rn}$ decay chain $\alpha$ particles specifically.

In Fig. 4, we plot $\frac{\mathrm{d}E}{\mathrm{d}x}$ lost to ionization as a function of $\alpha$ penetration depth in C3F8 of $\rho_{\mathrm{l}}=1379\,\mathrm{kg\,m^{-3}}$. These data are extracted with the Stopping Range of Ions in Matter (SRIM) package Ziegler et al. (2010) for $\alpha$ particles with energies equal to the full energies released in the decays. The recoil energies of the heavy daughter nuclei are of order $100\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ and are deposited over tens of nanometers in close proximity to the significantly longer $\alpha$-induced vapor regions. For that reason, our model assumes that the total energy released in each decay is deposited continuously over the length of the respective $\alpha$ track.

To understand how different portions of the long $\alpha$-induced vapor tubes evolve, we simulate 10, 30, 50, 100, and $150\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ energy depositions over $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in five separate simulations. The details of the simulation box and energy deposition region geometry are given in Appendix B. The resulting tube of vapor is allowed to evolve in an NPT ensemble with periodic boundary conditions on all sides. This means that each micron-long vapor chunk is evolving exactly as if it was surrounded by the same vapor chunks on both sides, imitating long alpha-particle tracks quite well. We observe that among the tested configurations, only the vapor tube produced as the result of a 10\,\mathrm{\frac{keV}{\text{,}\mathrm{\SIUnitSymbolMicro m}}} energy deposition is separated into several protobubbles, all of which eventually collapse, as shown in Fig. 5. All other tubes keep expanding irreversibly as in Fig. 6. In Fig. 7, we show the evolution of the vapor tube volumes resulting from $E_{\mathrm{dep}}\geq 30\,\mathrm{keV}$. This observation implies that the bubbles induced by the three ${}^{222}\mathrm{Rn}$ chain $\alpha$ particles (Sec. II.4) will predominantly not split along $\alpha$ ionization tracks, as the per-$\text{\,}\mathrm{\SIUnitSymbolMicro m}$ energy depositions are smaller than $30\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ only over the last few microns to the right of the Bragg peaks (see Fig. 4).

This leads us to the conclusion that any bubble nucleated as the result of energy deposition by an $\alpha$ particle will evolve similarly to a bubble resulting from a nuclear recoil considered in Secs. II.2 and II.3. If we further approximate the Bragg curves from Fig. 4 by uniform depositions of $E_{\alpha}$ (5.6, 6.1, and $7.8\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$) over the respective penetration depths, the evolution of the produced bubbles will be described, as in the nuclear recoil case, by Eq. 1. The latter will apply at all times when an $\alpha$-induced bubble may be considered spherical, i.e. $t\geq t_{\mathrm{s}}$, where $t_{\mathrm{s}}$ is to be estimated using Eq. 5 for $l_{\mathrm{cyl}}=42.6,48.3,$ and 70.0\,\mathrm{\text{,}\mathrm{\SIUnitSymbolMicro m}} (given in the order of increasing $E_{\alpha}$). In Table 2, we list the estimated times following the instant of the $\alpha$-induced heat spike at which the resulting bubbles acquire spherical shapes, as well as the corresponding bubble radii.

II.5 Extraction of acoustic signal

Using Eq. 1, we modeled the evolution of bubbles nucleated as the result of energy depositions by both nuclear recoils (Sec. II.3) and $\alpha$ particles (Sec. II.4), starting from the time when these bubbles become spherical. Given $R(t)$ in each case, we follow Landau and Lifshitz (2013) to evaluate the total power radiated in sound waves in all directions:

[TABLE]

where $\ddot{V}$ is the second time derivative of the bubble volume $V(t)\equiv\frac{4}{3}\pi R^{3}(t)$ and $c\approx 333\,\mathrm{m\,s^{-1}}$ is the speed of sound in the liquid phase of superheated C3F8 at $T=$ $14\text{\,}\mathrm{\SIUnitSymbolCelsius}$. In Fig. 8, we plot the acoustic power evaluated from Eq. 8 as a function of time for fluorine recoils with energies $\mathcal{T}=5.5,50.5,$ and $500.5\,\mathrm{keV}$, with corresponding $l_{\mathrm{cyl}}$ extracted from SRIM Ziegler et al. (2010), as well as the three $\alpha$-induced bubbles from Table 2.

To find the frequency spectrum of the extracted acoustic signal, we compute discrete Fourier transform (DFT) of the time-dependent acoustic pressure ($p_{\mathrm{ac}}\propto\ddot{V}$). Fig. 9 gives the resulting frequency spectra for the cases considered in Fig. 8. In this study, we focus on the frequency range of $1300\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ to compare our distributions of acoustic energies $E_{\mathrm{ac}}$ to those observed in the PICO-60 calibration runs Mitra (2018); Amole et al. (2017). Specifically, we evaluate

[TABLE]

with $A(f)$ being the complex amplitude of the DFT at a frequency $f$ in the returned samples and $\left|\mathcal{A}(f)\right|^{2}$ giving the respective power spectrum.

From Fig. 8a, we see that at times $t\geq t_{\mathrm{s}}$, only a short interval covers $\mathcal{P}(t)$ monotonically increasing with $t$ for nuclear recoil events. Near $t\sim 1\,$\text{,}\mathrm{\SIUnitSymbolMicro s}$$, a maximum in $\mathcal{P}(t)$ is observed, after which the bubble growth slows down in the thermal phase and $\mathcal{P}(t)$ starts to drop. In Fig. 8b, we find a similar pattern in the acoustic power radiation due to expansion of $\alpha$-induced bubbles. A sudden drop in $\mathcal{P}(t)$ indicates an inflection point characteristic of the times close to $t_{\mathrm{s}}$. At all times prior to $t_{\mathrm{s}}$, the bubbles are not yet spherical and hence the exact $R(t)$ behavior cannot be predicted with the current model. However, using the $R(t)$ data obtained directly from LAMMPS simulations, we found that $\mathcal{P}(t)$ at $t\leq t_{\mathrm{s}}$ is nearly six orders of magnitude smaller than $\mathcal{P}(t)$ at $t\geq t_{\mathrm{s}}$ in the $1300\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ frequency range, which lets us neglect the contribution from the former.

III Generation of acoustic parameter distribution

In the PICO experiment, a typical data analysis reports the distribution of AP (acoustic parameter) values. The latter are derived directly from the magnitude of the acoustic signal recorded during the expansion of all bubbles observed as single events in the fiducial volume during the run. Neutron-induced nuclear recoils form a low-AP peak in the AP distribution. Expansion of $\alpha$-induced bubbles, which are background events in both calibration and dark matter search runs, is accompanied by radiation of larger acoustic energies. Accordingly, $\alpha$ events feature AP values higher than those characteristic of nuclear recoil (NR) events, and form distinguishable peaks to the right of the latter. The collected data are typically normalized so that center of the NR peak has an AP value of 1.

In Sec. II.5, we described a general technique for evaluating the acoustic energy $E_{\mathrm{ac}}$ emitted in $1300\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ frequency range for both nuclear recoil and $\alpha$-like energy depositions. Additionally, we provided specific examples of $E_{\mathrm{ac}}$ calculation for fluorine recoil energies $\mathcal{T}\in\{5.5,50.5,500.5\}\,\mathrm{keV}$ and alpha energies $E_{\alpha}\in\{5.6,6.1,7.8\}\,\mathrm{MeV}$. For all cases, a uniform energy deposition along the SRIM-derived mean ranges of F/He ions in C3F8 was assumed. In reality, however, each projected track length value is subject to certain variation $\sigma(l_{\mathrm{cyl}})$ due to collisions, which can be extracted from SRIM as longitudinal straggling Ziegler et al. (2010). This variation will inevitably translate into a spread in times $t_{\mathrm{s}}$ when the nucleated bubbles become spherical, as well as their radii $R_{\mathrm{s}}=R(t_{\mathrm{s}})$, according to Eqs. 5 and 6. Ultimately, we will obtain a distribution of radiated acoustic energies that will reflect the variations in the length of the energy deposition region $l_{\mathrm{cyl}}$.

For the nuclear recoil energies $\mathcal{T}$, we consider a range from 0.5 to $1000.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$. Below $0.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$, the probability of nucleating a stably growing bubble in C3F8 is zero regardless of the energy deposition geometry, which follows the Seitz threshold model Seitz (1958) and was confirmed both experimentally Amole et al. (2017, 2019a) and through the MD simulations for our LJ fluid. Taking into account recoil energies exceeding $1\,\mathrm{MeV}$ will not contribute much to the resulting distribution of acoustic energies, for the probability of such high $\mathcal{T}$ is very small. We verify this by invoking the kinematics of ${}^{241}\mathrm{Am}/^{9}\mathrm{Be}$ source neutrons, whose original spectra were measured in e.g. Marsh et al. (1995) with high precision. These neutrons enter the detector from the outside and propagate through the hydraulic fluid on their way to the target C3F8. By simulating this propagation in Geant4 Agostinelli et al. (2003) with the geometry of the PICO-40L detector, we obtain a shifted spectrum of ${}^{241}\mathrm{Am}/^{9}\mathrm{Be}$ neutrons, now directly incident on the target nuclei. We note that an ${}^{241}\mathrm{Am}/^{9}\mathrm{Be}$ source spectrum closely resembles that of a ${}^{252}\mathrm{Cf}$ source, which is why we omit a separate discussion of ${}^{252}\mathrm{Cf}$ neutron kinematics. Convolution of this spectrum with the cross sections of elastic neutron-nucleus scattering Brown et al. (2018) gives the distribution of the normalized probabilities of carbon and fluorine recoils shown in Fig. 10. We see that recoil energies above $1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ are less probable than the energies up to $100\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}$ by nearly 2 orders of magnitude, justifying our choice of the recoil energy range.

To generate an AP distribution resembling that obtained experimentally, we consider a sample of 11,000 neutron-induced recoils (among which 3,000 are those of 12C and 8,000 are those of 19F) and 5,000 energy depositions by each of the $\alpha$ particle types (II.4). The NR sample has a well-defined continuous distribution of energies according to Fig. 10. Given $N_{i}$ samples in the recoil energy range from $\mathcal{T}_{i}$ to $\mathcal{T}_{i}+\Delta\mathcal{T}$ ($\Delta\mathcal{T}=$5\text{,}\mathrm{k}\mathrm{e}\mathrm{V}$$), we can generate a Gaussian distribution of C/F ranges centered at the SRIM-derived mean values for the energy ${T}_{i}$ and having $\sigma(l_{\mathrm{cyl}})$ as a standard deviation. Eqs. 5 and 6 are then used to predict the time $t_{\mathrm{s}}$ when the nucleated bubbles will become spherical and their radii $R_{\mathrm{s}}$ at that point. From here, each NR instance in the distribution will have an associated bubble growth curve $R(t)$ generated from Eq. 1. In turn, $R(t)$ gives the acoustic power $\mathcal{P}(t)$ radiated at $t\geq t_{\mathrm{s}}$ (see Eq. 8). The acoustic pressure signal is then Fourier-transformed, and its normalized power spectrum is integrated over the frequency range of $1300\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ to obtain the acoustic energy $E_{\mathrm{ac}}$ emitted at $t\geq t_{\mathrm{s}}$ up to $0.01\text{\,}\mathrm{s}$. In this way, for each recoil energy $\mathcal{T}_{i}$ of carbon and fluorine ions, we are able to evaluate the acoustic energy emission at all stages of bubble growth following $t_{\mathrm{s}}$. This part of the calculation results in the low-$E_{\mathrm{ac}}$ nuclear recoil peak in the modeled distribution. The procedure of generating the distribution of acoustic energies radiated during expansion of the $\alpha$-induced bubbles is effectively identical, except only three populations of deposited energies are possible, which we assume to be equally probable.

Sec. IV presents the distribution of acoustic energies radiated by neutron-induced and $\alpha$-induced bubbles in a PICO-like bubble chamber. The uncertainties on $E_{\mathrm{ac}}$ values that are derived from the uncertainties on the model parameters ($A,B,\gamma,\nu_{\mathrm{l}},R_{\mathrm{s}}$, and $t_{\mathrm{s}}$) are not taken into account in generating such a distribution, as it would only result in a global shift of all AP values.

IV Results and discussion

Following renormalization of the simulated acoustic energies by the Gaussian mean of those radiated in nuclear recoil-induced bubble expansions, we obtain the distribution of AP values shown in Fig. 11a. We immediately observe that the derived AP pattern provides clear discrimination between neutron-induced nuclear recoils and $\alpha$ particles in $1300\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ frequency range. It can therefore be directly compared to the experimental AP distribution from the PICO-60 run Amole et al. (2017, 2019a) with $52\text{\,}\mathrm{k}\mathrm{g}$ of C3F8, which we reproduce in Fig. 11b.

Comparing the modeled and the experimentally obtained AP patterns, we report good qualitative agreement between the two, with larger AP values being a characteristic signature of $\alpha$ particle-induced bubbles in both cases. This outcome confirms that the spherical bubble growth model (1) is a plausible description of bubble dynamics for both NR- and $\alpha$-induced bubbles, as it correctly predicts the relative magnitudes of acoustic signals associated with the two bubble nucleation scenarios. The AP-based $\alpha$ particle discrimination used by the PICO collaboration thereby becomes fully motivated by a physical model for the first time.

Quantitatively, we find differences between the modeled and the simulated AP in both the means of the NR/alpha peaks and the spread of AP values about these means. The discrepancies between the exact positions of the corresponding peaks on Figs. 11a and 11b are not surprising. Indeed, in this study, we evaluated acoustic energy emitted in all directions during bubble expansion and converted it to APmod directly. In the experimental data, APexp is a measure of the acoustic signal magnitude as recorded by the piezoelectric sensors, the response of which is dependent not only on the total energy radiated in sound waves, but also on the exact sound wave propagation characteristics and the accompanying chamber wall vibrations. Several simplifications have been made prior to arriving at the acoustic emission model as well, such as describing C3F8 as a LJ fluid at the stage of molecular dynamics simulations. We therefore do not expect the APmod values to match the APexp ones exactly. The question of intrinsically higher NR peak resolution and clearer separation of the $5.6\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ and $6.1\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ alpha peaks is of a similar nature. In our model, only the energies deposited and the statistical variation in the respective ion ranges give rise to a continuum of radiated acoustic energies. In a more realistic scenario, the presence of noise in the recorded signal, limited transducer sensitivity, and the precision of bubble position reconstruction all contribute to the width of each observed peak, along with the aforementioned physical reasons. A full treatment of sound wave propagation in the PICO bubble chambers and instrumental responses to the original acoustic signal would be required to predict the exact smearing effects on the peaks from Fig. 11a. Such a study is beyond the scope of the present work and is in progress within the PICO collaboration.

V Conclusions

In this paper, we have attempted to physically motivate acoustic $\alpha$-particle background discrimination in C3F8-filled PICO bubble chambers. The model that we have constructed builds on predictions for bubble growth in the superheated C3F8-like Lennard-Jones fluid, which we simulated with the methods of molecular dynamics. The bubble evolution was found to depend heavily on initial conditions set by the magnitude and geometry of energy depositions by nuclear recoils and $\alpha$ particles but was otherwise driven by exactly the same physics in the two cases. We have shown that ${}^{222}\mathrm{Rn}\to\,^{210}\mathrm{Pb}$ chain $\alpha$ particles inevitably nucleate only a single bubble per each $\alpha$ track in the simulated fluid. This has given us the grounds to apply a universal approach in evaluating the acoustic energies radiated during bubble expansion by both $\alpha$-induced and nuclear recoil-induced bubbles, based on the volume change rate of a single bubble. With this simple idea, we were able to reproduce the main empirical finding of Aubin et al. (2008); namely, the stronger acoustic signature of all $\alpha$ events relative to the nuclear recoils. This feature was found to originate from the larger volume change rate of $\alpha$-induced bubbles at the times when most acoustic energy is emitted. The successful modeling of the acoustic spectrum validates our understanding of the bubble formation and acoustic emission based on the molecular dynamics simulations, applicable to both types of ionizing radiation considered. We thereby come to the first model-based justification for the approach used by PICO to discriminate $\alpha$ background from calibration neutron-induced nuclear recoils, which are kinematically similar to WIMP-induced recoils.

The present study can be further expanded by considering other types of background radiation relevant to underground bubble chambers. In particular, bubble nucleation by electron recoils has been of continued interest. The PICO collaboration has recently discovered that the Seitz model Seitz (1958) for bubble nucleation at moderate superheats only applies to nuclear recoils and alpha interactions Amole et al. (2019b). Electron recoils, due to the extraordinary suppression at the thermodynamic operation point of dark matter bubble chambers, follow a different model of threshold behavior than that predicted by the Seitz model. Therefore a separate set of molecular dynamics simulations would be required to investigate how the electron recoil bubble nucleation mechanism influences the timeline and magnitude of acoustic emission.

VI Acknowledgements

We wish to acknowledge an earlier study on molecular dynamics simulations of bubble nucleation in direct dark matter search detectors by Denzel et al. Denzel et al. (2016), which formed the foundation of the present work. In particular, we thank Philipp Denzel for useful clarifications on their publication and making the developed code publicly available via phdenzel.github.io. The computational resources and support were kindly provided by Compute Canada (www.computecanada.ca), SciNet, and WestGrid. We thank Gabriel Hanna, Bruce Sutherland, and John Beamish of the University of Alberta for helpful discussions on components of this study. We also wish to acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) for funding.

Appendix A Lennard-Jones potential parameters for C3F8

In molecular dynamics simulations, the Lennard-Jones (LJ) potential Jones (1924) is a common way to describe simple interatomic interactions pair-by-pair. The LJ potential is defined through the characteristic energy (potential well depth) $\epsilon$ and the characteristic length (effective atomic diameter) $\sigma$, which are the units used within MD simulations. Such a potential is typically cut off at a certain critical distance to reduce computational costs of a simulation. To minimize the error of such a truncation, the interaction force is shifted so that it goes to zero continuously, resulting in a truncated shifted-force Lennard-Jones (TSF-LJ) potential Toxvaerd and Dyre (2011). Following Denzel et al. (2016); Diemand et al. (2014), we choose the cutoff radius to be $2.5\sigma$, which ensures the deviation from the full-tail $u_{\mathrm{LJ}}$(r) to be smaller than $1.6\%$ Toxvaerd and Dyre (2011).

Since we aim to compare the results of the present study with the acoustic data from the PICO-60 run at \sim$$14\text{\,}\mathrm{\SIUnitSymbolCelsius} Amole et al. (2017), it is important to reproduce the major thermodynamic properties of C3F8 at this particular temperature, corresponding to $T_{0}=0.778\,\epsilon/k_{\mathrm{B}}$ in LJ units. For that reason, we repeated the standard liquid-vapor interface experiment (see e.g. Diemand et al. (2014)) for $0.778\,\epsilon/k_{\mathrm{B}}$ to extract equilibrium pressure $P_{\mathrm{eq}}$, liquid/vapor density $\rho_{\mathrm{l}}/\rho_{\mathrm{v}}$, and surface tension $\gamma$ in LJ units ($\epsilon\,\sigma^{-3}$, $m\,\sigma^{-3}$, and $\epsilon\,\sigma^{-2}$ respectively, where $m$ = $188.02\text{\,}\mathrm{u}$ for C3F8) at saturation. We then used the corresponding real values of these quantities from the REFPROP database Lemmon et al. (2018) to deduce the suitable value of $\sigma$, with $\epsilon$ fixed at $0.0318\text{\,}\mathrm{e}\mathrm{V}$ as discussed in Sec. II.1. We started the experiment by creating a simulation box with length and width $L_{x}=L_{y}=200\sigma$ and height $L_{z}=100\sigma$, to which we applied periodic boundary conditions, and placing a liquid slab in the middle of the box. The slab initially contained 1,240,000 “atoms” on a lattice, taking the whole $XY$-plane of the simulation box and $30\sigma$ in $z$ direction. The atoms were assigned Maxwell-distributed velocities corresponding to $T_{0}$. To preliminarily randomize their positions, we performed a microcanonical (NVE: constant particle number, volume, and energy) ensemble simulation for 100,000 time steps. A single time step was set equal to $0.0025\tau$, where $\tau=\sigma\sqrt{m/\epsilon}\approx$4.172\text{,}\mathrm{p}\mathrm{s}. In each run, we used a Langevin thermostat to control the temperature. After we obtained such stable liquid conditions, we allowed the system to evolve for 2,000,000 time steps as an *NVT* (constant particle number, volume, and temperature) ensemble. The purpose of the first 1,000,000 time steps was to obtain two well-equilibrated liquid-vapor interfaces, below and above the liquid slab, after which we were able to take the measurements $P_{\mathrm{eq}}$ and $\rho_{\mathrm{l}}/\rho_{\mathrm{v}}$ in the subsequent 1,000,000 time steps. The time averaging of these quantities was performed 20 times in 50,000-time-step-long intervals, which ultimately gave the $1\sigma$ scatter for the global weighted average Diemand et al. ([2014](#bib.bib24)). The equilibrium pressure was evaluated as the spatial average of the normal ($z$) component of the pressure tensor Walton et al. ([1983](#bib.bib25)). To obtain values for liquid and vapor densities, we produced a density profile of the system over the $z$-axis, with $\rho(z)$ spatially averaged over chunks of $0.5\sigma$ height. The hyperbolic tangent function Chapela et al. ([1977](#bib.bib26)) was then used to fit the time-averaged $\rho(z)$ dependence and extract the liquid and vapor densities as fit parameters. Finally, to compute surface tensions, we followed the Kirkwood-Buff approach Kirkwood and Buff ([1949](#bib.bib27)). [Table 3](#A1.T3) presents the results of the $P_{\mathrm{eq}}$, $\rho_{\mathrm{v}}$, $\rho_{\mathrm{l}}$, and $\gamma$ measurements in LJ units; also listed are the respective real values of these quantities extracted from REFPROP database Lemmon et al. ([2018](#bib.bib11)) in SI units. We observe that there is not a single value of $\sigma$ that would reconcile all of the measurements simultaneously. As discussed in [Sec. II.1](#S2.SS1), we pick $\sigma\approx 0.533\,$\mathrm{n}\mathrm{m} to approximate “intermolecular” interactions within C3F8 with the Lennard-Jones potential. This matches the simulation extracted liquid density at $14\text{\,}\mathrm{\SIUnitSymbolCelsius}$ with the corresponding REFPROP value and introduces $\lesssim 30\%$ errors on other quantities listed in Table 3.

Appendix B Heat spike simulation procedure

For the bubble nucleation simulation, we chose a cubic box extending from $-\frac{L_{i}}{2}$ to $\frac{L_{i}}{2}$ in each dimension ($i=\overline{x,y,z}$), with $L_{i}\approx 363.2\sigma$. This box was filled with 32,000,000 atoms on a face-centered cubic lattice, corresponding to $\rho_{\mathrm{l}}=0.668\frac{m}{\sigma^{3}}$. The latter value was chosen so as to reproduce the density of superheated C3F8 at $T_{0}=$ $14\text{\,}\mathrm{\SIUnitSymbolCelsius}$, $P_{\mathrm{l}}$ = $207\text{\,}\mathrm{k}\mathrm{P}\mathrm{a}$, which is equal to 1379 $\mathrm{\frac{kg}{m^{3}}}$ Lemmon et al. (2018). Periodic boundary conditions were applied in all three directions. Following the approach of Denzel et al. (2016) to bring an LJ fluid to a superheated metastable state, we started our system at a temperature exceeding $T_{0}$, namely $0.8\,\epsilon/k_{\mathrm{B}}$ ($22\text{\,}\mathrm{\SIUnitSymbolCelsius}$). This gave us stable liquid conditions, which were retained for 10,000 time steps in an NVE ensemble to equilibrate the system. The next step was to linearly lower the temperature to the target $0.778\,\epsilon/k_{\mathrm{B}}$ over the time interval of 15,000 steps, which was again done in an NVE ensemble. The fixed liquid density was now too low for $0.778\,\epsilon/k_{\mathrm{B}}$. Thus, the liquid became superheated. To equilibrate the system at these new conditions, it was evolved for another 30,000 time steps under NVE constraints, with temperature fixed for the first 15,000 steps.

Following the instant of a sufficiently energetic heat spike in the prepared superheated fluid, the approach of Denzel et al. Denzel et al. (2016) was to integrate the system in an NVE ensemble until the pressure in the tiny liquid volume ($\approx[$162\text{,}\mathrm{n}\mathrm{m}$]^{3}$) increased significantly. This pressure rise slowed the bubble growth down even on a ns scale, which would not have been a realistic description for the early stages of bubble expansion in any PICO-scale bubble chamber. One way to avoid this complication would be to simulate a much larger liquid volume, naturally resulting in increased computational costs. For this reason, we chose to deviate from the method described in Denzel et al. (2016) by switching from an NVE to an NPT ensemble for the post-spike part of the simulation. In NPT integration, the overall particle number, pressure, and temperature are preserved, which reproduces typical conditions of a large-scale bubble chamber at the stages of bubble expansion when most acoustic energy is emitted.

For $\alpha$-induced bubble nucleation simulations, we set up a simulation box with $L_{x}=1888.6\sigma\approx 1\,$\text{,}\mathrm{\SIUnitSymbolMicro m}$$ and $L_{z}=L_{y}=159.8\sigma\approx 85\,\mathrm{nm}$ in LAMMPS. As before, we required $\rho_{\mathrm{l}}=0.668\,m\,\sigma^{-3}=1379\,\mathrm{kg\,m^{-3}}$ and $T=0.778\,\epsilon/k_{\mathrm{B}}=14^{\circ}\mathrm{C}$. To superheat the liquid, we followed the same procedure as in the nuclear recoil case. However, several modifications to the energy deposition region geometry were made. We observed that if the cylinder radius $r_{\mathrm{cyl}}$ was set to $2\sigma$ as in the nuclear recoil simulation, the NPT integrators failed to maintain the requested $0.778\,\epsilon/k_{\mathrm{B}}$ temperature for the whole system due to overly high temperature within the cylinder. To overcome this issue, we enlarged the radius to $3\sigma$. The cylinder length was set equal to $L_{x}$.

Appendix C Bubble surface tracking

Throughout the MD simulations in LAMMPS, we made use of a grid of cubic cells of size $[4\sigma]^{3}$ and computed the numbers of atoms falling within these cells as averaged over $1000$ time steps ($\approx$0.01\text{,}\mathrm{n}\mathrm{s}$$) intervals. After each such interval, we saved a slice of the resulting spatial density distributions into a text file. With the energy deposition axis parallel to the $x$-axis and centered at $y=0$, $z=0$, we took the slice between $-2\sigma\leq z\leq 2\sigma$. This allowed us to record the 2D evolution of the post-spike densities in the $XY$ plane.

After that, we applied the measure.find_contours function of the scikit-image Python module Lorensen et al. (1987); Walt et al. (2014) to extract the bubble contour at each moment in time. This function acts upon the 2D density data and returns discrete sets of points ${x_{i},y_{i}}$ for all the contours where $\rho_{\mathrm{l}}(x_{i},y_{i})=c$ holds. With a suitable choice of $c$, one of the contours found with this method is that of the vapor-liquid boundary. Our tests showed that choosing the longest contour among those extracted is a perfectly robust method to pick up the bubble “surface” boundary if we require $c=0.4\leavevmode\nobreak\ m\,\sigma^{-3}$ $\approx$ 830\text{\,}\mathrm{k}\mathrm{g}$$\,\mathrm{m^{-3}}. This method of bubble surface tracking was developed as an alternative to recursive linking of vapor cells used in Denzel et al. (2016); Diemand et al. (2014). In this study, we are only interested in the evolution of a single bubble, which is why identification of vapor regions in the rest of the liquid volume as in Denzel et al. (2016); Diemand et al. (2014) would be redundant.

Appendix D Transition from nonspherical to spherical bubble shape and further dynamics

The differential form of Eq. 1 is universal for all spherical bubbles irrespective of how they were nucleated. Parameters $A,B,\gamma$, and $\nu_{\mathrm{l}}$ in Eq. 1 are also shared between the different bubble expansions in the same liquid. We may therefore fit a few of the simulated bubble growth histories $R(t)$ for these parameters and use the best-fitting values later on to predict the dynamics of any bubble expanding within our LJ fluid. For this purpose, we use the bubble evolution data from four simulations mimicking energy depositions by nuclear recoils. The summary of these simulations is given in Fig. D.1. We restrict the fit to those parts of the data where the bubbles are already spherical. To avoid overfitting, we find it advantageous to fix $A$ at $17.7\,\mathrm{m\,s^{-1}}$, which is the value we predict from Eq. 2 and Table 1 for C3F8. We let $B,\gamma$, and $\nu_{\mathrm{l}}$ be free parameters that are shared between the four data sets we described, and arrive at the results in Table 1.

While these parameters are universal for any bubble nucleated within the LJ fluid, the initial conditions for the spherical bubble growth equation (1) vary with deposited energy and ion track length. This is reflected in Figs. D.2 and D.3. The time $t_{\mathrm{s}}$ of transition to bubble sphericity has an apparent dependence only on the length of the energy deposition region. The bubble radius $R_{\mathrm{s}}$ at $t_{\mathrm{s}}$ is also affected by the energy deposited. We combine these observations in Eqs. 5 and 6.