1 Introduction↩︎

Sonoluminescence (SL) is the conversion of acoustic energy into light, occurring for example when a gas bubble suspended in liquid is driven by ultrasound to undergo periodic cavitation [1]. During collapse, the bubble wall accelerates to speeds approaching the local speed of sound in the surrounding liquid, compressing the enclosed gas to transient states of extreme temperature and pressure [2]. The phenomenon occurs naturally in biological systems, most notably in the cavitation produced by the snapping shrimp [3], and has attracted sustained interest for its potential medical applications, including in ultrasound-assisted lipectomy [4], acoustically activated drug delivery [5], and cavitation assessment in prosthetic heart valves [6].

Despite decades of theoretical and experimental investigation, the microphysics of the collapse in single-bubble sonoluminescence (SBSL) remains incompletely understood [2], [7]. The relevant events unfold over nanosecond time scales at microscopic spatial scales, making experimental characterization of the bubble’s interior thermodynamic state during collapse fundamentally difficult. Molecular dynamics (MD) simulation provides a feasible alternative for studying SBSL, with an ability to use minuscule time steps and simulate individual particle trajectories. Prior simulation studies of SL [8]–[14] have relied primarily on event-based hard-sphere MD implementations, which are inherently sequential and difficult to parallelize at large scale; to the authors’ knowledge, the largest such MD study reached \(10^7\) ensemble particles [10], which is deficient compared to the roughly \(10^{10}\) molecules involved in representative real-world SBSL experiments. The hard-sphere framework also provides no natural path toward including long-range Coulomb interactions between the charged carriers produced during ionization, an omission that becomes physically consequential at the extreme temperatures reached near the collapse minimum.

This work introduces the first time-based molecular dynamics (TBMD) formulation of SBSL. Rather than resolving particle interactions as instantaneous collision events, TBMD integrates continuous equations of motion at each timestep, allowing the interaction model to combine a soft-sphere Lennard–Jones potential with a damped shifted force (DSF) Coulomb interaction [15], [16] between ionized species. Implemented within the LAMMPS framework [17] with custom extensions for ionization, gas–wall energy exchange, and Keller–Miksis wall coupling, our proposed method scales to at least \(N_{\mathrm{ensem}} = 10^8\) ensemble particles, a tenfold increase over the prior state of the art. We apply the framework to a parameter sweep of the thermal accommodation coefficient \(\alpha_t\) and ionization model, quantifying how each modeling choice affects peak temperature, wall dynamics, and light emission, and we validate our results against the EBMD reference dataset of Schanz et al. [10]. Source code is released under the GPL-2.0 license.¹

The remainder of the paper is organized as follows. Section 2 reviews prior continuum and MD studies of SL. Section 3 describes the continuum Keller–Miksis framework used to drive the bubble wall in our simulations. Section 4 presents the TBMD model, including interaction potentials, ionization algorithm, and gas–wall coupling. Section 5 derives novel scaling rules for accurately setting physical quantities on ensemble particles. With all this in hand, Section 6 describes our simulation setup and diagnostics, and Section 7 presents a variety of numerical results. Section 8 examines convergence with ensemble size to assess how many simulation particles are needed for reliable simulation results. We conclude and discuss future directions in Section 9.

2 Related Work↩︎

Single-bubble sonoluminescence involves tightly coupled fluid dynamics, thermodynamics, and atomic-scale physics spanning at least twelve orders of magnitude in time. Prior numerical work has therefore developed along two largely parallel streams: continuum partial differential equation (PDE) models that resolve the macroscopic bubble dynamics and temperature field, and molecular dynamics simulations that track discrete particle dynamics near collapse. We review each of these research streams and then briefly discuss the MD cyberinfrastructure that enables the present work.

2.1 Continuum Modeling of Sonoluminescent Bubble Collapse↩︎

The foundation of all continuum SL modeling is an equation of motion for the spherical bubble wall. In its simplest form, this is the Rayleigh–Plesset (RP) equation, which balances inertial, viscous, and surface tension forces on an incompressible liquid shell [18]. An early compressible-liquid formulation was given by Gilmore [19], who applied the Kirkwood–Bethe hypothesis to reduce the flow equations around the bubble to a single ordinary differential equation for the wall velocity valid into the supersonic regime. At the extreme wall speeds reached during SL collapse, where the wall Mach number is not negligible, a compressible extension is required; the Keller–Miksis (KM) equation [20] incorporates acoustic radiation damping and remains accurate into the transonic regime. We adopt the KM equation as the bubble-motion model in this work. Broad reviews of bubble oscillation dynamics and the continuum SL modeling landscape are provided by Lauterborn and Kurz [21] and Brenner, Hilgenfeldt, and Lohse [22].

Coupling the bubble radius equation to the thermal and pressure state of the gas interior requires additional modeling. Kwak and Yang [23] and Kwak and Na [24] developed a self-contained hybrid PDE model that pairs the KM equation with a thermal boundary layer of thickness \(\delta(t)\) in the surrounding liquid and derives closed-form temperature profiles under both uniform and nonuniform pressure assumptions. This Kwak model forms the theoretical basis of the continuum reference solution used throughout our results; a description of the numerical procedure for solving the coupled ODE/PDE system is given in Section 3.5. A different continuum approach by Vuong and Szeri [25] solves the Navier–Stokes equations for the gas interior coupled to diffusive transport, predicting sharp shock-like features near the minimum radius. Hydrodynamic simulations of SL collapse, including the formation and propagation of compression waves inside the bubble, were carried out by Moss et al. [8], [11] using the full Euler equations, and by Popinet and Zaleski [12] for axisymmetric non-spherical bubble geometries. A pedagogical treatment of the single-bubble dynamics in the continuum framework, including a systematic comparison of ODE models of increasing complexity, is given by Vignoli et al. [13].

Water vapor trapped inside the bubble during collapse plays a physically important role that purely noble gas models omit. For instance, Storey and Szeri [26] and Toegel et al. [27] demonstrated that the endothermic dissociation of H\(_2\)O provides a significant internal energy sink that limits the peak temperature and has a measurable effect on the SL emission spectrum. Moss et al. [28] identified a new damping mechanism in strongly collapsing bubbles and proposed a modified RP equation that correctly predicts the rapid decay of post-collapse rebounds observed in SL experiments, resolving a longstanding deficiency of the classical RP equation, whose solutions persist until the next acoustic cycle. Sensitivity studies by Vazquez and Putterman [29] and Yasui [30] quantified how the ambient liquid pressure and temperature, respectively, shift the collapse intensity and emission brightness, helping to anchor continuum models to measurable experimental observables. More recent mesh-based simulations [31] of acoustic cavitation of bubbles in water have incorporated explicit water-vapor chemistry and spatial resolution of species and energy inside the bubble, further underscoring the importance of vapor dissociation in limiting peak temperatures and controlling radical production.

Continuum and hybrid radiative models have also been used to interrogate the ionization state and light emission mechanism of SL more directly. Moss et al. [32] computed optical emission from a partially ionized collapsing bubble including opacity effects; Hammer and Frommhold [33] extended the analysis to rare-gas bubbles containing water vapor; and An [34] incorporated chemistry, ionization, and competing radiative channels into a unified single-bubble SL emission model. These theoretical predictions are supported by direct experimental evidence: Flannigan and Suslick [35], [36] and Flannigan et al. [37] established transient plasma formation in argon bubbles spectroscopically; Chen et al. [38] resolved the time-dependent emission structure; and Khalid et al. [39] and Kappus et al. [40] characterized opacity effects and the energy budget of the light pulse, respectively. At the interface level, molecular simulation studies of vapor–liquid kinetic boundary conditions have shown that accommodation-like departures from simple Maxwell reflection arise naturally at nonequilibrium liquid interfaces [41], [42], underscoring the physical importance of the wall-coupling parameter \(\alpha_t\) explored in Section 7.2.

2.2 Molecular Dynamics Simulation of Sonoluminescence↩︎

Molecular dynamics provides a complementary perspective on SL collapse: discrete particle trajectories naturally resolve shock formation, local ionization, and statistical fluctuations without requiring closure assumptions for the equation of state or transport coefficients. The coupling between the MD gas interior and the surrounding liquid has been handled exclusively through a hybrid RP-MD (or KM-MD) approach, in which the MD simulation supplies the instantaneous gas pressure at the wall and the continuum equation advances the bubble radius.

The first MD study of SL was by Ruuth, Putterman, and Merriman [9], who modeled the gas interior as a hard-sphere noble gas driven by a Rayleigh–Plesset piston. Their event-based MD (EBMD) implementation, using fast tree-based collision-scheduling algorithms, reached up to \(10^6\) ensemble particles and included a collision-energy threshold ionization model—the predecessor of the algorithm adopted in our work. The study confirmed the formation of a central hot spot and estimated that the light emission duration scales with the ambient bubble radius. The same work noted that between \(10^9\) and \(10^{10}\) physical particles would be required to simulate a realistic SL bubble without ensemble-particle scaling, a gap that remains largely open to this day. Bass et al. [43], [44] subsequently developed a symmetry-reduced hard-sphere MD framework for spherically symmetric bubble implosion, enabling simulations of xenon bubbles with up to \(5\times10^7\) particles and revealing strong mass segregation during rapid collapse². These studies showed that algorithmic specialization can push SL bubble MD to much larger particle counts, but only by assuming exact spherical symmetry and remaining within the event-driven hard-sphere paradigm.

Kim, Kwak, and Kim [14] coupled a hard-sphere EBMD gas interior to the full Kwak–Na continuum model [24], enabling a direct comparison between MD trajectories and PDE predictions for a xenon bubble in sulfuric acid at \(10^6\) particles. Their study found good qualitative and quantitative agreement between the MD and continuum results for the velocity, temperature, and pressure evolution, validating the hybrid KM-MD framework. Among the most comprehensive of prior EBMD studies is Schanz et al. [10], who augmented the hard-sphere gas model with water vapor, vapor chemistry (dissociation, radical formation), and separate mass and energy transfer channels at the bubble wall. Their implementation reached up to \(10^7\) ensemble particles—though most production runs used \(10^6\)—and systematically explored the accommodation coefficient, driving amplitude, and noble-gas species, establishing a reference dataset that we use for direct comparison in Section 7.4. Ramsey [45] examined energetic cavitation collapse, including sonoluminescence, from the perspective of scalable simulation and identified GPU parallelism as a key enabling technology for the next generation of large-scale SL computations.

Despite their respective contributions, prior SL MD studies share a set of structural limitations that motivate the present work. Interaction model: Most prior large-scale SL MD simulations use hard spheres, and even the short-range model of Kim et al. [14] provides no natural route to including long-range Coulomb forces between the charged carriers produced by ionization. At the extreme temperatures reached during SL collapse, the ionization fraction is non-negligible [9], [10], and the Coulomb interaction between charge carriers exerts a physically meaningful restoring pressure that ought to be incorporated. Scalability: EBMD requires maintaining a global priority queue of future collision events, an operation that is inherently sequential in time and difficult to parallelize across distributed-memory nodes or GPUs. None of the prior EBMD implementations exploits GPU acceleration. Particle count: The maximum ensemble size achieved in prior work (\(10^7\) in Schanz et al. [10]) falls two to three orders of magnitude short of particle counts seen in real-world experiments; while our present results (\(N_\mathrm{ensem} = 10^8\), Section 7.1) do not fully close this gap, we obtain an order-of-magnitude increase over state of the art—despite paying substantial additional costs for incorporating additional modeling terms.

The present work addresses each of these limitations by adopting a time-based MD (TBMD) formulation with a soft-sphere Lennard–Jones potential supplemented by a damped-shifted-force Coulomb interaction [15], full multi-stage ionization, and a KM-coupled continuum boundary, implemented in a scalable framework that reaches \(10^8\) ensemble particles. To our knowledge, this is the first application of TBMD to single-bubble sonoluminescence.

2.3 Large-Scale Molecular Dynamics Infrastructure↩︎

The computational demands of large-scale particle simulations have driven sustained development of general-purpose MD packages and parallelization strategies. The LAMMPS package [17], [46] provides a flexible, domain-decomposed framework for particle-based simulation that scales efficiently across thousands of CPU cores and has well-optimized GPU force kernels. GPU acceleration for general MD was studied systematically by Glaser et al. [47], who demonstrated near-linear strong scaling across GPU clusters for systems of tens of millions of particles. At the exascale frontier, the EXAALT project [48] targets atomistic simulations across macroscopic length scales, ambitiously targeting billion-particle MD. Large-scale MD has also been applied directly to bubble collapse: Vedadi et al. [49] simulated shock-induced nanobubble collapse in water using reactive MD; Chen et al. [50] systematically studied system-size effects in bubble collapse MD with up to 4.5 billion water molecules; and Asano [51] recently demonstrated 100-billion-atom acoustic cavitation MD on Fugaku, marking the current frontier of large-scale bubble simulation. At the particle-method level, the use of simulated or ensemble particles to represent many physical particles is standard in particle-in-cell and DSMC methodologies [52], [53], and weighted-particle extensions have been used to control fluctuations for rare ionized species in weakly ionized flows [54]. Our ensemble-particle formulation is closest in spirit to this broader macro-particle tradition, although the specific scaling laws are tailored here to dense, collapsing LJ–Coulomb gas dynamics. The present work is based on a version of LAMMPS with significant extensions for SBSL simulation, including ionization, wall coupling, and ensemble-particle scaling. Our SL framework therefore inherits all the parallelization, scalability, and modular aspects of LAMMPS.

3 Theory↩︎

The dynamics of a sonoluminescing bubble couple the radial motion of the bubble wall, the energy exchange between the hot compressed gas and the surrounding liquid, and the spatiotemporal distribution of temperature and pressure inside the bubble. We adopt the theoretical model developed by Kwak and collaborators [23], [24] as the continuum reference framework, which resolves each of these components through a coupled system of ordinary and partial differential equations. The model serves two roles in this work: it provides the initial and boundary conditions used to initialize the molecular dynamics simulation (Section 6), and it supplies the continuum baseline against which MD results are compared throughout Section 7. No single prior publication has collected a complete description of the numerical solution procedure for this model, so we present the governing equations and our implementation in the following subsections.

3.1 Bubble motion↩︎

The radial dynamics of the bubble wall are governed by the Keller–Miksis equation [20], which accounts for liquid compressibility and acoustic radiation damping: \[\frac{d R_b}{dt}=U_b , \label{eq:PDE:Rb}\tag{1}\] \[\left(1-\frac{U_b}{C}\right)R_b \frac{dU_b}{dt} +\frac{3}{2}U_b^2\left(1-\frac{U_b}{3C}\right) =\frac{1}{\rho_{\infty}}\left(1+\frac{U_b}{C} +\frac{R_b}{C}\frac{d}{dt}\left[P_B-P_s\left(t+\frac{R_b}{C}\right)-P_{\infty}\right]\right) , \label{eq:PDE:ddRb}\tag{2}\] where \(R_b\) is the bubble radius, \(U_b\) is the bubble wall velocity, \(C\) is the sound speed in the liquid, and \(\rho_{\infty}\) is the ambient liquid density. The pressure due to the driving acoustic wave is \(P_s(t)=-P_A \sin\omega t\), where \(P_A\) is the driving pressure amplitude and \(\omega=2\pi f\) with \(f\) the driving frequency. The liquid pressure outside the bubble wall \(P_B(t)\) is related by the pressure inside the bubble wall \(P_b(t)\): \[P_B=P_b(t)-\frac{2\sigma_s}{R_b}-\frac{4 \mu U_b}{R_b} , \label{eq:out95press}\tag{3}\] where \(\sigma_s\) and \(\mu\) are the surface tension and dynamic viscosity of the liquid medium.

3.2 Liquid dynamics↩︎

Although mass transfer through the bubble interface (evaporation and condensation of liquid molecules) is neglected, energy exchange between the gas and the surrounding liquid must be retained because it heats the liquid shell and thereby affects the predicted temperature at the bubble wall. Heat transfer is assumed to occur through a thermal boundary liquid shell of thickness \(\delta(t)\), governed by [55]: \[\left[1+\frac{\delta}{R_b}+\frac{3}{10}\left(\frac{\delta}{R_b}\right)^2\right] =\frac{6 \alpha}{\delta} -\left[2\frac{\delta}{R_b}+\frac{1}{2}\left(\frac{\delta}{R_b}\right)^2\right]\frac{dR_b}{dt} -\delta\left[1+\frac{1}{2}\frac{\delta}{R_b}+\frac{1}{10}\left(\frac{\delta}{R_b}\right)^2\right]\frac{1}{T_{bl}-T_{\infty}}\frac{dT_{bl}}{dt} , \label{eq:PDE:delta}\tag{4}\] where \(\alpha\) is the thermal diffusivity of the liquid, \(T_{bl}\) is the liquid temperature at the bubble wall, and \(T_{\infty}\) is the liquid temperature far away from the bubble.

3.3 Temperature distribution in the liquid thermal boundary layer↩︎

The temperature profile within the liquid boundary layer is assumed to be quadratic [56]: \[\frac{T-T_{\infty}}{T_{bl}-T_{\infty}}=(1-\xi)^2 ,\] where \(\xi=(r-R_b)/\delta\) and \(T_{bl}\) is again the temperature at the bubble wall. This quadratic profile satisfies the boundary conditions: \[T(R_b,t)=T_{bl}, \quad T(R_b+\delta,t)=T_{\infty}, \quad \left(\frac{\partial T}{\partial r}\right)_{r=R_b+\delta}=0 . \label{eq:temperature32boundary}\tag{5}\] The temperature profile given in Eq. 5 is valid provided the characteristic time scale of bubble evolution is much longer than the translational relaxation time of the molecules [57] yet comparable to or longer than the vibrational relaxation time [58].

3.4 Temperature and pressure distribution inside the bubble↩︎

The gas interior can be described in two ways: via a baseline uniform model, or via a nonuniform correction of said model. We first solve for a spatially uniform baseline state, in which the density and pressure are taken to be uniform inside the bubble—an approximation that is reasonable except near the bubble’s minimum radius. Spatial nonuniformity can then be reconstructed as a secondary correction to this baseline solution, accounting for the nonuniform density and pressure induced by bubble wall acceleration, a transient effect that becomes significant on nanosecond time scales near collapse.

3.4.0.1 Uniform baseline

The thermal conductivity of the gas inside the bubble is assumed to vary linearly with temperature: \[k_g=AT+B . \label{eq:conductivity}\tag{6}\] A linear fit to the conductivity data of Boulos et al. [59] provides a good approximation in the range \(200\,\mathrm{K} < T < 3000\,\mathrm{K}\) [18]. In practice, the linear form is applied up to \(25{,}000\,\mathrm{K}\), above which \(k_g\) is capped at \(k_g' = 5\,\mathrm{W\,m^{-1}\,K^{-1}}\). Applying Fourier’s law with this conductivity model yields the following radial temperature profile under the uniform assumption [58]: \[T_b(r)=\frac{B}{A}\left(-1+\sqrt{\left (1+\frac{A}{B}T_{b0}\right)^2-2\eta\left(\frac{A}{B}\right)(T_{bl}-T_{\infty})\left(\frac{r}{R_b}\right)^2 }\right) .\] The wall temperature \(T_{bl}\) is obtained by applying the boundary conditions of Eq. 5 : \[T_{bl}=-\frac{B}{A}(1+\eta)+\left( \frac{B}{A}\right) \sqrt{(1+\eta)^2+2\left(\frac{A}{B}\right)\left(T_{b0}+\frac{A}{2B}T_{b0}^2+\eta T_{\infty}\right)} ,\] \[\eta=\frac{R_b}{\delta}\frac{k_1}{B} . \label{eq:eta}\tag{7}\] Applying energy conservation to the uniformly distributed gas yields evolution equations for the center temperature and pressure: \[\frac{dT_{b0}}{dt}=-\frac{3(\gamma-1)T_{b0}}{R_b}\frac{dR_b}{dt}-\frac{6(\gamma-1)k_1(T_{bl}-T_{\infty})T_{b0}}{\delta R_b P_{b0}} , \label{eq:PDE:Tb0}\tag{8}\] \[\frac{dP_{b0}}{dt}=-\frac{3\gamma P_{b0}}{R_b}\frac{dR_b}{dt}-\frac{6(\gamma-1)k_1(T_{bl}-T_{\infty})}{\delta R_b} , \label{eq:PDE:Pb0}\tag{9}\] where \(T_{b0}\) and \(P_{b0}\) are the gas temperature and pressure at the bubble center, respectively, \(\gamma\) is the specific heat ratio, and \(k_1\) is the thermal conductivity of the liquid.

3.4.0.2 Nonuniform correction

The abrupt temperature rise and subsequent rapid quenching driven by bubble wall acceleration occur on a time scale distinct from the conduction-dominated evolution captured by the uniform baseline. This motivates incorporating a nonuniform correction. Building on that baseline rather than treating the interior pressure as spatially constant, this correction decomposes the gas density into a spatially uniform part \(\rho_0\) and a radially dependent correction \(\rho_r\): \[\begin{align} \begin{aligned} \rho_g &=\rho_0+\rho_r \\ u_g &=\frac{\dot{R_b}}{R_b}r \\ P_b&=P_{b0}(t)-\frac{1}{2}(\rho_0+\frac{1}{2}\rho_r)\frac{\ddot{R_b}}{R_b}r^2 , \label{eq:density} \end{aligned} \end{align}\tag{10}\] where \(\rho_0 R_b^3=\mathrm{const}\), \(\rho_r=ar^2/R_b^5\), and a dot denotes a time derivative. The constant \(a\) is related to the total gas mass \(m\) inside the bubble by \[\frac{a}{m}=\frac{5}{4\pi}(1-N_{BC}) ,\] where \(N_{BC}=(P_{b0}R_b^3/T_{b0})/(P_{\infty}R_0^3/T_{\infty})\) is a constant obtained from the uniform baseline and \(R_0\) is the equilibrium radius. The nonuniform density and pressure distributions yield an inertial correction to the temperature profile: \[\begin{align} \begin{aligned} T(r)&=T_b(r)+T'_b(r) \\ T'_b(r)&=-\frac{1}{40(\gamma-1)k_g'}\left(\rho_0+\frac{5}{21}\rho_r \right)\left[(3\gamma-2)\frac{\dot{R_b}\ddot{R_b}}{R_b^2}+\frac{\dddot{R_b}}{R_b}\right]r^4+C(t) \\ C(t)&=\frac{1}{20(\gamma-1)}[(3 \gamma-2)\dot{R_b}\ddot{R_b}R_b+\dddot{R_b}R_b^2]\times \left[\frac{\delta'}{k_l}\left(\rho_0+\frac{5}{14}\rho(R_b)\right)+\frac{R_b}{2 k_g'}\left(\rho_0+\frac{5}{21}\rho(R_b)\right)\right] . \end{aligned} \label{eq:nonuniform} \end{align}\tag{11}\] The coefficient \(C(t)\) is determined from the interfacial heat–flux condition \(k_g' \,\partial_r T_b = k_l \,\partial_r T_l\) at the bubble wall, where \(T_l\) denotes the temperature inside the quadratic liquid thermal boundary layer of thickness \(\delta'\). In the continuum nonuniform model of Kwak & Na [60], \(\delta'\) is treated as a semi–empirical gas–side thermal layer thickness, and Eqs. 11 yield a perturbative inertial correction \(T_b'(r,t)\) that is intended to remain small relative to the conduction–controlled background temperature \(T_b(r,t)\).

However, when evaluated in the strongly driven SL regime considered here, the authors found that \(T_b'(r,t)\) becomes comparable to or larger than \(T_b(r,t)\) and may even render the predicted gas temperature negative near the interface at certain collapse stages, indicating that the linearized correction lies outside its validity range. For this reason, we retain \(T_b(r,t)\) as the quantitative theoretical temperature profile used for comparison with MD results, and regard further investigation into nonuniform correction models as future work.

We note that we verified our implementation against the continuum model of Kwak & Na by reproducing their air–bubble predictions under the same acoustic and thermophysical conditions reported in Kwak & Yang [23], obtaining excellent agreement for both the bubble radius evolution and the associated temperature field (Fig. 1 (a)). However, when the same nonuniform formulation, including the inertial correction \(T_b'(r,t)\) of Eqs. 11 , is applied to the strongly driven SL configuration considered in this work, the model produces unphysical results: at certain collapse stages it predicts negative gas temperatures near the bubble wall and substantially distorted thermal profiles (Fig. 1 (b)). We conjecture that these pathologies arise because the quantity driving the correction term, \((3\gamma-2)\dot{R}_b\ddot R_b/R_b^2 + \dddot R_b/R_b\), becomes sufficiently large that \(T_b'(r,t)\) is no longer a small perturbation to the conduction-controlled background temperature \(T_b(r,t)\), violating the formal assumptions of the perturbative derivation in Kwak & Na [60]. For this reason, in our simulations, we retain the full continuum dynamics for the mean bubble evolution while restricting ourselves to using the uniform model \(T_b(r,t)\) as the quantitative theoretical temperature profile. This choice captures the dominant nonuniform conduction structure and avoids introducing the non-physical temperature artifacts that arise when \(T_b'(r,t)\) is applied outside its regime of validity.

3.5 Numerical solver↩︎

Our numerical solver integrates the five-variable uniform baseline system (Eqs. 1 –9 ) using MATLAB’s ode45 solver. To ensure consistency with our molecular dynamics experiments (detailed in subsequent sections), we initialize the system with the following conditions: \(R_b = R_0 = 4.5\,\mu\)m, \(U_b = 0\) m/s, \(P_b = P_{\infty} = 101325\) Pa, \(T_{bl} = T_{b0} = T_{\infty} = 300.0\) K, and \(\delta = 0.3 R_0\). The remaining simulation parameters are set as follows: \(f = 26.5\) kHz, \(\gamma = \frac{5}{3}\), \(P_A = 1.3 P_{\infty}\), \(C = 1481\) m/s, \(\mu = 0.001\) Pa\(\cdot\)s, \(\sigma_s = 0.072\) N/m, \(\rho_{Ar} = 1.603\) kg/m\(^3\), \(A = 2.682 \times 10^{-5}\) W/mK\(^2\), \(B = 1.346 \times 10^{-2}\) W/mK, \(k_1 = 0.61\) W/mK, and \(N_{BC} = 1.316\). The nonuniform inertial correction \(T_b'(r,t)\) of Eqs. 10 and 11 is evaluated post hoc from the baseline solution solely to generate Fig. 1 (b); as discussed above, it is not used as a quantitative prediction.

To maintain numerical stability and accuracy, we employ a two-stage adaptive time-stepping scheme. First, we perform a preliminary simulation with a fixed timestep of \(10^{-10}\) s to obtain a coarse estimate of the time at which the bubble reaches its minimum radius. Based on this estimate, we then rerun a refined simulation with a timestep of \(10^{-15}\) s, starting approximately 5 ns before the coarse minimum-radius time and running for a total of 10 ns. Throughout this work, \(t_{\min}\) denotes the minimum-radius time identified by this fine-run simulation. Because the timestep decreases by five orders of magnitude between the two stages, the refined \(t_{\min}\) does not in general coincide with the midpoint of the stage-2 window; in our configuration it falls at approximately 4.6 ns from the start of that window rather than at the expected 5 ns. This offset is reflected in the simulation time axes of the results section.

4 Molecular Dynamics Simulations↩︎

This section describes the molecular dynamics component of our hybrid continuum–MD framework. The principal challenge in simulating sonoluminescence is the need to track large particle ensembles through extreme thermodynamic states on nanosecond time scales. We address this using the LAMMPS molecular dynamics package [17] as the computational backbone, extended with custom modules for ionization, gas–wall energy exchange, and Keller–Miksis wall coupling.

4.1 Hard-sphere and soft-sphere models↩︎

In the hard-sphere model, each atom is represented as a rigid sphere that interacts only through instantaneous, perfectly elastic binary collisions; no interaction force acts between particles at any other time. The soft-sphere model, by contrast, assigns each particle a continuous, spherically symmetric potential (e.g., Lennard–Jones) with a finite cutoff. Although the hard-sphere model reproduces noble gas viscosity accurately [53] and shows little difference from soft-sphere results in dilute conditions [61], its validity in a dense, strongly compressed system near the minimum bubble radius is unclear. Furthermore, a continuous potential provides a natural route to include long-range Coulomb forces between charged carriers produced during ionization, enabling a more physically complete simulation.

4.2 Time- and event-based molecular dynamics↩︎

Time-Based Molecular Dynamics (TBMD) relies on integrating the equations of motion using discrete time steps, such as through the Verlet or leapfrog algorithms. At each time step, the forces on the particles are computed, and their positions and velocities are updated accordingly. Event-Based Molecular Dynamics (EBMD), on the other hand, takes a different approach by focusing only on events of interest, such as collisions between particles or changes in potential energy states. Instead of updating particle positions at every small time step, the algorithm predicts the time of the next event and skips directly to it. In the SL context, the timestep and initial conditions must be chosen carefully to achieve a stable and accurate simulation. There is a one-to-one correspondence with the preceding subsection: EBMD is based on the hard-sphere model, whereas TBMD is based on the soft-sphere model. The remainder of this section assumes TBMD and soft spheres for the particles.

We note that in the context of SBSL, a potentially drastic downside of the EBMD approach is that as millions or even billions of particles are compressed into a microscopic bubble, the number of collisions increases exponentially, incurring great computational cost to perform searches for the next collision. In this highly compressed regime, TBMD may also have a drawback in that two or more particles may non-physically pass through each other if the time step is too large (i.e., collisions are missed). However, we take the view that a small number of missed collisions is an acceptable price to pay for being able to take a relatively large time step, and our numerical results demonstrate that our TBMD approach produces results that align with the literature while enabling substantially larger particle counts.

4.3 Gas dynamics inside the bubble↩︎

In contrast to the hard-sphere model, where collisions are resolved only at discrete events based on energy and momentum conservation, the present approach requires computing the pairwise interaction potential at each timestep. The total pair interaction energy combines a short-range Lennard–Jones term and a Coulombic contribution computed via the damped shifted force (DSF) method [15]: \[E_{\text{pair}}(r) = E_{\text{LJ}}(r) + E_{\text{Coul-DSF}}(r),\] where each contribution is computed with its own cutoff radius.

4.3.0.1 Lennard–Jones potential

Short-range molecular interactions are governed by the Lennard–Jones 12-6 potential, which accounts for van der Waals forces:

\[E_{\text{LJ}}(r) = \begin{cases} 4 \epsilon \left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6} \right], & r < r_{\text{LJ}} \\ 0, & r \geq r_{\text{LJ}} \end{cases}\] where \(r\) is the interatomic separation, \(\epsilon\) is the depth of the potential well, and \(\sigma\) is the characteristic length scale at which the potential energy is zero. The interactions are truncated at a cutoff distance \(r_{\text{LJ}}\), beyond which the potential is set to zero. For argon we use \(\sigma = 3.401\,\text{\AA}\) and \(\varepsilon/k_B = 116.81\,\text{K}\); all force-field and cutoff parameters are collected in Table 3.

4.3.0.2 Coulombic interaction with the DSF method

For charged particles, the electrostatic interaction is modeled using the damped shifted force (DSF) potential [15], which provides a computationally efficient alternative to long-range summation methods such as Ewald [62] or particle–particle particle–mesh (PPPM) [52] while ensuring smooth truncation and local charge neutrality within a cutoff radius \(r_{\text{Coul}}\):

\[E_{\text{Coul-DSF}}(r) = \begin{cases} q_i q_j \left[ \frac{\operatorname{erfc}(\alpha r)}{r} - \frac{\operatorname{erfc}(\alpha r_{\text{Coul}})}{r_{\text{Coul}}} + \left(\frac{\operatorname{erfc}(\alpha r_{\text{Coul}})}{r_{\text{Coul}}^2} + \frac{2\alpha}{\sqrt{\pi}} \frac{e^{-\alpha^2 r_{\text{Coul}}^2}}{r_{\text{Coul}}} \right)(r - r_{\text{Coul}}) \right], & r < r_{\text{Coul}} \\ 0, & r \geq r_{\text{Coul}} \end{cases} ,\] where \(q_i\) and \(q_j\) are the particle charges, \(\alpha\) is the damping parameter, and \(\operatorname{erfc}(\cdot)\) is the complementary error function. We set the damping parameter to \(\alpha = 3/r_{\text{Coul}}\), which places the Gaussian screening width at one-third of the cutoff so that the residual interaction at the cutoff, \(\operatorname{erfc}(\alpha r_{\text{Coul}}) = \operatorname{erfc}(3) \approx 2.2\times10^{-5}\), is a negligible fraction of the bare Coulomb value. In the production runs, both cutoffs are set equal: \(r_{\text{Coul}} = r_{\text{LJ}} = 2\,d_{\text{ensem}} \approx 2.15\,\sigma_{\text{ensem}}\), ensuring Coulomb interactions are resolved across the full LJ range. The ionization reactive cutoff is \(r_{\text{ion}} = g^{1/2}d\), where \(d = 3.66\,\text{\AA}\) is the argon atomic diameter and \(g = N_{\text{real}}/N_{\text{ensem}}\).

4.4 Gas ionization inside the bubble↩︎

Near the minimum bubble radius, the extreme pressure and temperature are sufficient to ionize a significant fraction of the gas. Ionization removes kinetic energy from the gas and thereby limits the peak temperature, making its proper treatment essential for accurate predictions. We follow a strategy analogous to that of Ruuth et al. [9]: whenever the combined kinetic energy \(E = E_i + E_j\) of a candidate pair exceeds the ionization potential \(\chi\) of the less-ionized atom, ionization is assumed to occur with certainty. Both particles’ velocities are rescaled by the same factor \(s = \sqrt{2(E-\chi)/(3E)}\), giving each particle \(k \in \{i,j\}\) a post-ionization kinetic energy \[E_{k,\mathrm{new}} = \frac{2(E - \chi)}{3E}\,E_k , \label{eq:ionize}\tag{12}\] where \(E_k\) is the pre-event kinetic energy of particle \(k\). The total energy extracted from the pair is \(\chi + \tfrac{1}{3}(E-\chi)\): the first term is the ionization potential, and the second accounts for the freed electron taking one-third of the remaining kinetic energy as it thermalizes with the gas—the two-loss mechanism identified by Ruuth et al. [9]. The model tracks how many electrons each particle loses, enabling the use of successive ionization energies for future ionization events. This also helps determine local ionization levels and properly account for Coulomb forces. A check for potential ionization events is conducted at every timestep. Specifically, if two atoms, \(i\) and \(j\), are within the ionization reactive cutoff \(r_\mathrm{ion}\) of each other and their combined kinetic energy exceeds the ionization energy of at least one of the atoms, the pair becomes a candidate for ionization. The reactive cutoff \(r_\mathrm{ion}\) is an independent parameter, distinct from the Lennard–Jones and Coulomb cutoffs, that sets the range within which a collision can trigger ionization. During the compression phase of sonoluminescence, each atom may encounter multiple neighboring atoms eligible for ionization. To resolve this, each atom examines its list of potential ionization partners and selects the one with the shortest distance as its “sole” ionization candidate.

Once this labeling is complete, a pair \((i, j)\) is deemed an “eligible” ionization pair if and only if \(i\) labels \(j\) as its sole ionization partner and vice versa. If the combined kinetic energy is sufficient to ionize either \(i\) or \(j\), the model prioritizes ionizing the atom with the lower ionization level, thereby favoring the atom that requires less energy to ionize further. This approach does not aim to strictly replicate the continuous nature of ionization in the physical world; instead, it provides a discrete approximation, with finer timesteps yielding more accurate representations of ionization dynamics.

4.5 Gas-wall dynamics and energy transfer↩︎

In modeling liquid–gas interactions, the surrounding liquid is treated as a continuum, and gas–liquid coupling arises from collisions of gas molecules with the bubble wall. The gas–liquid interface is represented by a repulsive Lennard–Jones wall potential of the same 12-6 form described in Section 4.3. Because the liquid is treated as a continuum, heat exchange between gas and liquid occurs at the boundary through molecular–wall collisions. Each collision leads to an adjustment of particle velocities, thereby capturing thermal effects at the boundary.

The heat energy \(P_Q\) transferred by \(\Delta N\) particles colliding with the wall within a timestep \(\Delta t\) is given by: \[P_Q = \frac{3}{2} k_B \frac{\Delta N}{\Delta t} (T_i - T_r) = \frac{3}{2} k_B \alpha_t \frac{\Delta N}{\Delta t} (T_i - T_w), \label{eq:heatflux}\tag{13}\] where \(T_i\) is the temperature of the incoming particles, \(T_w\) is the wall temperature, and \(k_B\) is the Boltzmann constant. The thermal accommodation coefficient \(\alpha_t\) governs the degree of heat conduction between the gas and liquid. The temperature of particles reflected from the wall, \(T_r\), is then \[T_r = (1 - \alpha_t)\,T_i + \alpha_t\,T_w.\]

By selecting an appropriate \(\alpha_t\), one can control the thermal boundary condition: \(\alpha_t = 1\) corresponds to an isothermal boundary, whereas \(\alpha_t = 0\) corresponds to a completely adiabatic boundary. Because the wall potential is continuous, a single physical collision spans multiple timesteps. To avoid applying the thermal energy transfer multiple times within one contact, we track for each molecule \(p\) the timestep \(t^p\) of its most recent wall contact, updated whenever the molecule lies within the gas–wall interaction cutoff. Energy transfer is applied at the onset of each new contact event, identified as the first timestep for which \(t^p\) advances by more than one step from its previous value.

4.6 Energy redistribution within the liquid shell↩︎

The total energy stored in the liquid shell is updated at each timestep to account for three contributions: \[E_{\mathrm{new}} = E_{\mathrm{old}} + P_Q\,\Delta t + E_v + E_{\mathrm{acoustic}}, \label{eq:shell95energy}\tag{14}\] where \(P_Q\) is the gas-to-liquid heat flux at the bubble wall (Eq. 13 ), \(E_v\) captures the energy change due to advection as the liquid shell expands or contracts, and \(E_{\mathrm{acoustic}}\) is the work performed by the external acoustic field on the shell.

The advection term is \[E_v = C_l \rho_\infty T_\infty (V_{\mathrm{new}} - V_{\mathrm{old}}), \label{eq:Ev}\tag{15}\] where \(C_l\) is the specific heat capacity of the liquid, \(\rho_\infty\) is the ambient liquid density, \(T_\infty\) is the ambient temperature, and \(V_{\mathrm{new}} - V_{\mathrm{old}}\) is the change in volume of the liquid shell. The acoustic work term is \[E_{\mathrm{acoustic}} = P_s(t)\cdot 4\pi R_b^2\,\dot{R}_b\,\Delta t, \label{eq:acoustic95work}\tag{16}\] where \(P_s(t) = -P_A\sin\omega t\) is the acoustic driving pressure. This term ensures that the energy bookkeeping is thermodynamically consistent across the full acoustic cycle, regardless of the duration of the MD simulation window.

Given the updated total energy in the liquid shell, the bubble wall temperature \(T_{bl}\) is recovered through the quadratic temperature distribution of Eq. 5 . At each timestep, the MD-computed gas pressure at the bubble wall \(P_b\) and the heat flux \(P_Q\) are returned to the continuum solver to update \(R_b\), \(U_b\), and \(T_{bl}\), closing the hybrid loop.

4.7 End-to-end algorithm for one timestep↩︎

The components described above—velocity-Verlet integration, the Lennard–Jones/DSF pair interaction, the discrete ionization model, the gas–wall heat bath, and the Keller–Miksis continuum solver—are advanced together once per MD timestep within a single LAMMPS run loop. Figure 2 shows how the data flows between the MD particle system and the continuum bubble-wall state, and Algorithm 3 gives the complete step-by-step procedure for one timestep. The two solvers are coupled in a staggered manner: the gas pressure \(P_b\) and wall temperature \(T_{bl}\) measured by the MD system at the current step are consumed by the continuum update, whose new wall radius \(R_b\) and velocity \(U_b\) define the moving boundary seen by the MD particles at the next step.

5 Ensemble-Particle Representation and Scaling↩︎

Even with our TBMD formulation, simulating the \(N_\mathrm{real}\sim 10^{10}\) gas atoms inside a sonoluminescing bubble at fully atomistic resolution is computationally challenging. We therefore represent the gas by \(N_\mathrm{ensem} \ll N_\mathrm{real}\) ensemble particles (EPs), each standing in for \(g = N_\mathrm{real}/N_\mathrm{ensem}\) microscopic atoms that occupy similar positions and velocities in phase space. The goal of the mapping is to preserve the dominant mesoscopic observables—temperature, short-range collisional hardness, pressure, and ionization statistics—while reducing the particle count to a more computationally tractable level.

Our proposed design principle is as follows. Replacing \(N_\mathrm{real}\) atoms with \(N_\mathrm{ensem} = N_\mathrm{real}/g\) EPs reduces the number density by a factor of \(g\). To maintain the same local packing and interaction frequency, all characteristic length scales are enlarged by \(g^{1/3}\), restoring the volume per particle to its physical value. All energies are scaled by \(g\) so that the energy per microscopic atom, and hence the temperature, is unchanged. The following subsections apply this principle to each interaction in the model.

5.1 Mass and kinetic energy↩︎

Each EP aggregates the mass of \(g\) atoms while leaving the velocity unchanged: \[m' = g\,m, \qquad v' = v. \label{eq:ep95mass}\tag{17}\] The per-EP kinetic energy is therefore \(K' = \tfrac{1}{2}m'v^2 = g\,(\tfrac{1}{2}mv^2)\); since \(v' = v\) and the temperature formula uses the microscopic mass \(m\), the temperature is unchanged. The Maxwellian velocity distribution and all characteristic thermal velocities and flight times are preserved exactly.

5.2 Interaction potential scaling↩︎

5.2.0.1 Lennard–Jones potential

Reducing the number density by \(1/g\) would soften short-range encounters and alter the virial pressure contribution. To maintain the effective collisional hardness, the LJ parameters are rescaled as \[\sigma' = g^{1/3}\sigma, \qquad \epsilon' = g\,\epsilon, \qquad r_c' = g^{1/3} r_c. \label{eq:ep95lj}\tag{18}\] This enlarges the interaction range to the EP length scale and increases the well depth so that EP–EP collisions reproduce the local stress contribution of \(g\) microscopic pairs.

5.2.0.2 Coulombic DSF interaction

The length and damping parameter of the DSF potential follow the same geometric scaling, \[r_c' = g^{1/3} r_c, \qquad \alpha' = g^{-1/3}\alpha, \label{eq:ep95dsf95cutoff}\tag{19}\] which preserves the shape of the screened Coulomb field at the EP scale. Charges are left unscaled (\(q_i' = q_i\)). Each EP has mass \(m' = gm\) but carries a single ionic charge \(q\), so without scaling its Coulomb acceleration would be \(1/g\) that of the physical ion. Setting \[s_{\mathrm{pair}} = g, \label{eq:ep95dsf95scale}\tag{20}\] restores the correct ionic acceleration. This leaves the Coulomb virial undercounted by \(g^{1/3}\); since the ionized fraction reaches at most 6–12%, the impact on the total pressure is small.

5.3 Ionization scaling↩︎

Ionization encounters are screened by the reduced EP number density. The reactive radius is enlarged as \[r_{\mathrm{ion}}' = g^{1/2} r_{\mathrm{ion}}, \label{eq:ep95rion}\tag{21}\] which compensates for the lower encounter rate and yields an approximately invariant ionization frequency per EP. The ionization energy threshold is scaled in proportion to the energy budget of the \(g\) microscopic constituents, \[E'_\mathrm{ion} = g\,E_\mathrm{ion}, \label{eq:ep95eion}\tag{22}\] thereby preserving the ratio of relative kinetic energy to ionization cost.

5.4 Summary↩︎

Table 1 collects the complete EP parameter mapping. Together, these transformations preserve the local forces, collisional hardness, ionization statistics, and charging dynamics that govern the macroscopic pressure, temperature, and light emission of the sonoluminescing gas, while reducing the simulated particle count by up to four orders of magnitude relative to the physical system.

6 Simulation Procedures↩︎

6.1 Physical setup↩︎

All simulations in this work use argon; Table 2 also lists helium and xenon for reference. The successive ionization potentials are drawn from the NIST Atomic Spectra Database [63] and are applied directly as the ionization thresholds in the collision-energy criterion described in Section 4.

The physical molecule count inside a sonoluminescing argon bubble is obtained from the ideal gas law: \[P_{\infty} \left(\frac{4 \pi}{3} R_0^3\right) = N_{\mathrm{real}}\, k_B T_{\infty}.\] For the setup used throughout this paper (\(R_0 = 4.5\,\mu\mathrm{m}\), \(T_{\infty} = 300\,\mathrm{K}\), \(P_{\infty} = 1\,\mathrm{atm}\)), this gives \(N_{\mathrm{real}} \approx 9.3 \times 10^9\) molecules, which is computationally challenging to simulate at full atomistic resolution. We therefore employ the ensemble-particle (EP) representation described in Section 5, where each EP stands in for \(g = N_{\mathrm{real}}/N_{\mathrm{ensem}}\) microscopic atoms. Production simulations use ensemble sizes \(N_{\mathrm{ensem}} \in \{10^6, 10^7, 10^8\}\), corresponding to EP diameters of \(21.1\,d\), \(9.8\,d\), and \(4.5\,d\) respectively (with \(d = 3.66\,\text{\AA}\) for argon), and behavior with respect to \(N_{\mathrm{ensem}}\) is examined in Section 7.1.

6.2 Simulation time window and initialization↩︎

Prior EBMD studies of sonoluminescence initialize the bubble at its maximum radius and simulate the entire inward collapse, which is computationally feasible in event-based codes because computation is proportional to the number of collision events rather than wall-clock time. TBMD, by contrast, advances every timestep regardless of activity; simulating a full acoustic half-cycle at the sub-femtosecond resolution required near collapse would be prohibitive. We therefore initialize the simulation with the bubble already contracting and restrict the simulation to the 10 ns fine-resolution window identified in Section 3.5. The continuum solution provides the bubble radius, wall velocity, and boundary temperature at the window start time, which are used directly to initialize the MD system. This starting point is justified by the prior finding that both uniform and nonuniform continuum theories remain accurate through the pre-collapse approach, so that the MD need only resolve the final compression and rebound where the two descriptions diverge.

Particle positions are initialized to reproduce the radial density profile \(\rho(r)\) from the continuum solution. Particles are placed on a cubic lattice with spacing \(2d_{\mathrm{ensem}}\) spanning the bubble volume, where \(d_{\mathrm{ensem}} = g^{1/3} d\) is the ensemble-particle diameter, \(d = 3.66\,\text{\AA}\) is the argon atomic diameter, and \(g = N_{\mathrm{real}}/N_{\mathrm{ensem}}\); sites are accepted or rejected in parallel according to the local value of \(\rho(r)\), yielding a collision-free initial configuration for all ensemble sizes without pairwise overlap checks. Initial velocities are drawn from the Maxwell–Boltzmann distribution at the local temperature \(T(r)\) predicted by the continuum model, with isotropically randomized directions; the bulk inward velocity field \(u_g = (\dot{R}_b/R_b)r\) is not superimposed at \(t_0\). The inward-moving wall naturally imparts this bulk momentum to the gas through particle–wall collisions during the pre-collapse window, so no explicit correction is needed at initialization.

6.3 Diagnostics and post-processing↩︎

Macroscopic flow quantities are obtained by averaging particle data within spherical shells of equal volume. The radial profiles of bulk velocity, temperature, density, and pressure are computed as \[\begin{align} u_g(r) &= \frac{1}{N_r} \sum_{i=1}^{N_r} \frac{\boldsymbol{v}_i \cdot \boldsymbol{r}_i}{|\boldsymbol{r}_i|}, \notag\\ T_b(r) &= \frac{m}{3 k_B N_r} \sum_{i=1}^{N_r} \left|\boldsymbol{v}_i - u_g(r)\hat{\boldsymbol{r}}_i\right|^2, \notag\\ \rho_b(r) &= \frac{N_r}{V_r}, \notag\\ P_b(r) &= \rho_b(r) k_B T_b(r) + \frac{1}{3 V_r} \sum_{i=1}^{N_r} \boldsymbol{r}_i \cdot \boldsymbol{f}_i, \label{eq:MD95properties} \end{align}\tag{23}\] where \(N_r\) is the number of particles in the shell between \(r-\delta r\) and \(r+\delta r\), \(V_r\) is the shell volume, \(m\) is the microscopic atomic mass (not the EP mass \(m'=gm\)), and \(\boldsymbol{f}_i\) is the total force on particle \(i\) (the last term is the virial pressure contribution). The instantaneous bubble radius and wall velocity are advanced by the Keller–Miksis equation (Eq. 2 ), driven at each timestep by the gas pressure \(P_b\) computed from the MD particle data; the liquid boundary-layer thickness \(\delta(t)\) is prescribed by the continuum solution (Eq. 4 ), and all other quantities evolve self-consistently from the particle dynamics.

Raw per-bin diagnostics are subject to Poisson shot noise whenever a shell contains few particles, most severely near the bubble center at early times (small shell volume) and near the wall after rebound (gas rarefaction). To suppress this noise, we merge consecutive bins adaptively from the outermost shell inward until each merged group contains at least \(N_{\mathrm{min}} = 100\) ensemble particles, and assign each group the group-averaged temperature, density, and pressure. All particle-count and density fields are then rescaled by the factor \(g = N_{\mathrm{real}}/N_{\mathrm{ensem}}\) to convert from simulated to physical particle numbers; this rescaling applies to the number-density term in \(P_b\) only, as the virial term is already physically correct through the \(g\)-scaled EP forces. The scalar peak temperature \(T_{\mathrm{max}}\) reported in the results tables is evaluated only over bins whose rescaled count satisfies \(N_{\mathrm{count}} \geq 10\), excluding wall-adjacent bins that contain too few physical particles to yield a statistically meaningful estimate. Space–time temperature heatmaps use a robust display cap at the \(99.9\)th percentile of the non-empty-bin temperature distribution, which suppresses extreme outlier cells while preserving the interior radial structure.

7 Results↩︎

We structure the results as follows. Section 7.1 considers behavior with respect to ensemble particle number \(N_\mathrm{ensem}\) and validates the ensemble-particle (EP) scaling derived in Section 5. Section 7.2 investigates the sensitivity to the thermal accommodation coefficient \(\alpha_t\), which governs the gas–liquid energy exchange at the bubble wall. Section 7.3 examines the role of pair-interaction physics, comparing simulations with no ionization, short-range ionization only, and full-range ionization with Coulomb forces. Section 7.4 places our results in the context of prior simulation work through a comprehensive parameter-sweep table.

7.1 Trends with increasing ensemble particle number↩︎

A key question motivating our study is: how many particles is sufficient to effectively have a “fully-resolved” SBSL simulation? In other words, do we observe convergent or similar behavior as the number of EPs is increased? From another perspective, we wish to validate our EP model and scaling factors by determining whether macroscopic observables behave similarly as \(N_\mathrm{ensem}\) increases. We vary \(N_\mathrm{ensem} \in \{10^6,\,10^7,\,10^8\}\) while holding all other parameters fixed (\(\alpha_t = 1.0\), full-range ionization with Coulomb forces), spanning two orders of magnitude and reaching \(10^8\) particles—to our knowledge the largest sonoluminescence MD simulation performed to date.

Figure 4 shows the time evolution of the bubble radius, volume-averaged temperature, and volume-averaged pressure near the collapse point for the three particle counts. The three MD radius trajectories are mutually consistent, confirming that the ensemble-particle scaling captures the macroscopic collapse kinematics at all three resolutions. All three share the qualitative collapse shape of the Keller–Miksis prediction, but diverge substantially near the collapse minimum: the MD minimum radius (\(R_\mathrm{min} \approx 0.488\)–\(0.498\,\mu\text{m}\)) is roughly half the continuum value (\(0.883\,\mu\text{m}\)). Near the collapse minimum, several core assumptions of the Keller–Miksis continuum framework are known to break down simultaneously, including the local thermodynamic equilibrium required for continuum transport, an ideal-gas equation of state with constant \(\gamma\), and a thermal conductivity calibrated for \(T < 3{,}000\,\mathrm{K}\) [9], [14], [23]. The continuum model therefore captures only the gross features of the collapse dynamics and is retained as a qualitative reference; Section 7.4 provides quantitative validation of our results against prior MD work. The time of minimum radius \(t_\mathrm{min}\) shifts slightly later as \(N_\mathrm{ensem}\) increases (from \(4.548\,\text{ns}\) at \(10^6\) to \(4.622\,\text{ns}\) at \(10^8\)), and correspondingly the minimum radius \(R_\mathrm{min}\) increases modestly from \(0.488\,\mu\text{m}\) to \(0.498\,\mu\text{m}\). Both effects are consistent with finer resolution reducing spurious numerical stiffness during the most compressed phase. The volume-averaged temperature is nearly consistent across the three runs: the peak values are \(\langle T\rangle_\mathrm{max} = 10{,}967\,\text{K}\), \(10{,}870\,\text{K}\), and \(10{,}852\,\text{K}\) for \(10^6\), \(10^7\), and \(10^8\) particles respectively, a variation of less than \(1.1\%\) over two orders of magnitude of \(N_\mathrm{ensem}\). The Kwak continuum reference, which neglects ionization, peaks at \(17{,}494\,\text{K}\) (approximately \(61\%\) above the ionizing MD result), directly attributing the temperature suppression to ionization rather than to ensemble-particle resolution. The average pressure shows a more pronounced but still systematic trend, with the peak decreasing monotonically from \(5.70\,\text{GPa}\) at \(10^6\) to \(4.95\,\text{GPa}\) at \(10^8\) (\(\sim\)​13% reduction), indicating that lower-resolution runs slightly overestimate the sharpness of the compressive collapse. The Kwak continuum volume-averaged pressure peaks at \(1.55\,\text{GPa}\) (distinct from the wall pressure \(p_{W,\mathrm{max}}=1.488\,\text{GPa}\) listed in Table 4), substantially below the MD values; this gap reflects the shallower Kwak collapse (\(R_\mathrm{min} = 0.883\,\mu\text{m}\) versus \(0.49\)–\(0.50\,\mu\text{m}\) in MD) rather than a difference in thermodynamic modeling.

Figure 5 shows the radial profiles of temperature, pressure, and mass density as a function of the fractional radius \(r/R_b\). The temperature profile (panel a) is evaluated at the time of maximum center temperature \(t(T_{\mathrm{center,max}})\), where the center is defined as the region \(r \leq R_b/10\); this snapshot captures the sharpest hot-core structure, which occurs slightly before minimum radius. The pressure and density profiles (panels b–c) are evaluated at minimum bubble radius \(t_\mathrm{min}\), where the compression is greatest. All three profiles share the same qualitative structure across resolutions: a hot, dense, high-pressure interior that decreases toward the bubble wall. The temperature profiles display a pronounced hot core: the center temperature (averaged over \(r \leq R_b/10\) at \(t(T_{\mathrm{center,max}})\)) reaches \(15{,}536\,\text{K}\), \(15{,}354\,\text{K}\), and \(15{,}555\,\text{K}\) for \(N_\mathrm{ensem} = 10^6\), \(10^7\), and \(10^8\) respectively, occurring at \(0.266\,\text{ns}\), \(0.346\,\text{ns}\), and \(0.431\,\text{ns}\) before minimum radius. The center temperature does not converge monotonically with \(N_\mathrm{ensem}\), reflecting the inherently higher statistical noise in the small inner bins. The radial pressure and density profiles are more stable across resolutions: the overall shape and interior magnitude are consistent, with the \(10^8\) case yielding slightly lower and smoother profiles owing to better statistical sampling of each radial shell. Taken together, the behavior in Figs. 4–5 supports the conclusion that the ensemble-particle scaling developed in Section 5 successfully preserves the dominant mesoscopic physics at all three resolutions, and that \(N_\mathrm{ensem} = 10^8\) provides the best-resolved reference for quantitative comparisons.

7.2 Effect of thermal boundary condition↩︎

The thermal accommodation coefficient \(\alpha_t\) governs the rate of energy exchange between gas molecules and the bubble wall upon collision: \(\alpha_t = 0\) corresponds to a fully adiabatic wall (specular reflection, no energy transfer), while \(\alpha_t = 1\) corresponds to a fully isothermal wall (complete thermalization to the local wall temperature). Physical estimates for \(\alpha_t\) in SL systems are not directly accessible from experiment; Schanz et al. [10] adopt \(\alpha_t = 0.3\) as a representative value, while Kim et al. [14] consider both limiting cases. We compare \(\alpha_t \in \{0.0,\,0.5,\,1.0\}\) at fixed \(N_\mathrm{ensem} = 10^6\) with full-range ionization, to quantify how strongly the choice of \(\alpha_t\) shifts the predicted observables.

Figure 6 shows the time evolution of the volume-averaged temperature, center temperature, and wall temperature for the three accommodation coefficients. The volume-averaged temperature, shown in panel (a), decreases monotonically with increasing \(\alpha_t\): the peak values are \(14{,}257\,\text{K}\), \(12{,}404\,\text{K}\), and \(10{,}967\,\text{K}\) for \(\alpha_t = 0.0\), \(0.5\), and \(1.0\) respectively. This monotone reduction reflects the increasing rate at which energy is transferred from the gas to the liquid as the accommodation coefficient grows. More strikingly, panels (b) and (c) reveal a physically distinct behavior at the bubble center and wall. The center temperature exhibits the opposite ordering: \(T_{\mathrm{center,max}} = 14{,}380\,\text{K}\), \(14{,}989\,\text{K}\), and \(15{,}536\,\text{K}\) for \(\alpha_t = 0.0\), \(0.5\), and \(1.0\). That is, stronger wall coupling produces a hotter center despite a lower global average. The wall temperature follows the same trend as the center, increasing from \(2{,}429\,\text{K}\) to \(3{,}050\,\text{K}\) as \(\alpha_t\) increases. This behavior is physically consistent with the gas–wall energy exchange model of Section 4: a larger \(\alpha_t\) increases the heat flux \(P_Q\) transferred from the hot gas to the liquid shell at each wall collision, which raises the liquid-side boundary temperature \(T_{bl}\) through the shell energy update. Since gas molecules adjacent to the wall are thermalized toward this higher \(T_{bl}\) upon collision, the measured gas-side wall temperature rises in tandem. A further consequence of increasing \(\alpha_t\) is a systematic advance of the collapse: the time of minimum radius shifts from \(t_\mathrm{min} = 4.720\,\text{ns}\) at \(\alpha_t = 0\) to \(4.548\,\text{ns}\) at \(\alpha_t = 1.0\), and the thermal peak moves progressively closer in time to the radius minimum. At \(\alpha_t = 0.0\) the average temperature peaks approximately \(0.55\,\text{ns}\) before \(t_\mathrm{min}\); at \(\alpha_t = 1.0\) this interval narrows to only \(0.17\,\text{ns}\), indicating that stronger wall thermalization sharpens and delays the thermal concentration toward the final collapse stage.

The apparently paradoxical inversion between the average and center temperature trends is resolved by examining the spatial temperature structure. Figure 7 shows the space–time evolution of the gas temperature inside the bubble for each boundary condition, with the instantaneous bubble radius overlaid. At \(\alpha_t = 0.0\) (panel a), the temperature field is relatively broad and spatially uniform: energy is retained throughout the gas, the average is high, but the center does not develop a strong concentration relative to the bulk. As \(\alpha_t\) increases, the wall continuously extracts energy from the outer gas layers, steepening the radial thermal gradient. At \(\alpha_t = 1.0\) (panel c), the result is a pronounced hot core coexisting with a relatively cool near-wall region: energy is concentrated inward rather than distributed uniformly. The total thermal energy budget of the gas is reduced (lower average temperature) but the spatial focusing is stronger, driving the center to higher absolute temperatures than in the adiabatic case. The intermediate case \(\alpha_t = 0.5\) (panel b) bridges these extremes cleanly, displaying an intermediate degree of radial stratification. Since the primary radiative source in SL is the hot central plasma, these results indicate that the isothermal boundary condition, despite lowering the bulk temperature, produces the most localized and intense light-emitting region. The choice of \(\alpha_t\) therefore has qualitative as well as quantitative consequences for the predicted light emission, and experimental calibration of this parameter is important for reliable emission estimates.

7.3 Effect of pair interactions and ionization↩︎

Prior MD studies of sonoluminescence have modeled gas particles as either hard spheres [9], [10] or with a purely short-range Lennard–Jones potential [14], neglecting both the energy cost of ionization events and the long-range Coulomb interaction between the resulting charge carriers. At the extreme temperatures reached during SL collapse, however, significant ionization is expected: for argon under the conditions studied here, the Saha equilibrium predicts appreciable ionization fractions above approximately \(15{,}000\,\text{K}\). We compare three levels of pair-interaction physics at fixed \(N_\mathrm{ensem} = 10^6\) and \(\alpha_t = 1.0\): no ionization (LJ potential only), short-range ionization, and full-range ionization. Both ionized configurations use the identical Lennard–Jones and DSF Coulomb interactions and differ only in the ionization reactive cutoff \(r_\mathrm{ion}\) (see Table 3): the short-range case sets \(r_\mathrm{ion} = r_\mathrm{LJ} = 2\,g^{1/3}d\), confining ionization events to within the LJ interaction range, while the full-range case extends it to \(r_\mathrm{ion} = g^{1/2}d\), allowing ionization over a wider separation. We note that ionization events in the simulation are treated as irreversible; the light-emission estimates below are therefore obtained as a post-processing diagnostic by applying the Saha–LTE equation [64] to the local MD temperature and number-density fields, following the approach of Schanz et al. [10]. In the strongly-coupled plasma regime reached near collapse, ionization-potential lowering would increase the true ionization fraction beyond these Saha estimates [39]; the emission values below should therefore be interpreted as order-of-magnitude lower bounds rather than quantitative predictions.

Figure 8 compares the volume-averaged temperature, center temperature, and total simulation charge as functions of time for the three interaction models. The volume-averaged temperature (panel a) decreases monotonically with the level of ionization physics included: the peak values are \(14{,}975\,\text{K}\) (no ionization), \(11{,}557\,\text{K}\) (short-range), and \(10{,}967\,\text{K}\) (full-range), representing reductions of \(23\%\) and \(27\%\) relative to the no-ionization baseline. The center temperature (panel b) shows a far more dramatic effect: the no-ionization run reaches \(T_\mathrm{center,max} \approx 33{,}000\,\text{K}\), while both ionized runs remain near \(15{,}500\)–\(16{,}400\,\text{K}\) — a suppression of roughly a factor of two. This large discrepancy arises because, in the absence of any ionization energy sink, kinetic energy at the bubble center accumulates without bound during the final stages of compression, producing unrealistically high local temperatures. Once ionization is active, the energy consumed by successive electron-removal events acts as a distributed internal energy sink that prevents runaway heating of the hottest gas. The difference between the two ionized configurations is modest by comparison (\(\Delta T_\mathrm{center} \approx 870\,\text{K}\) relative to the short-range case): the larger ionization reactive cutoff \(r_\mathrm{ion}\) of the full-range case admits ionization events over a wider separation near the hot core, deepening the ionization energy sink and further lowering the peak center temperature. Panel (c) confirms that the net simulation charge remains identically zero in the no-ionization run, while both ionized models accumulate charge during the collapse, reaching peaks of \(5.78\%\) (short-range) and \(6.74\%\) (full-range) of the ensemble at the time of most intense ionization—which occurs slightly after the thermal maximum, as the charge carriers persist into the post-collapse phase.

The practical consequence for the predicted light output is illustrated in Fig. 9, which shows the space–time distribution of a Saha–LTE emission estimate (bolometric, band-limited to \(\lambda \geq 200\,\text{nm}\)) inferred from the local MD temperature and density fields. The no-ionization run produces a peak emission approximately three orders of magnitude larger than either ionized case. This dramatic difference arises from two compounding mechanisms. First, the ionization algorithm described in Section 4 removes kinetic energy from colliding particles at each ionization event, directly lowering the gas temperature that enters the Saha calculation: the center temperature in the ionized runs (\(\sim15{,}500\)–\(16{,}400\,\text{K}\)) is roughly half that of the no-ionization run (\(\sim33{,}000\,\text{K}\)). Second, the Saha–LTE ionization fraction scales exponentially with inverse temperature, \(q \propto \exp(-\chi/k_BT)\), and the emission proxy scales as \(q^2 n^2 T^{1/2}\); a factor-of-two reduction in \(T\) therefore produces an orders-of-magnitude reduction in \(q\) that enters the emission as \(q^2\). The three-order-of-magnitude drop in emission is thus an amplification of the factor-of-two temperature suppression through the Saha exponential, not a direct proportionality. To preserve the spatial structure within each panel, the heatmaps use independent per-panel display scales; absolute cross-panel comparison should be made using the scalar totals in Table 5. In all three cases the emission is sharply concentrated near the bubble center in a brief temporal window around the collapse. These results demonstrate that neglecting ionization not only introduces a quantitative bias in the temperature prediction, but leads to emission estimates that are qualitatively inconsistent with physically expected magnitudes. Extending the ionization reactive cutoff provides a further, physically motivated refinement beyond the short-range ionization treatment alone.

7.4 Comparison with prior work↩︎

Table 4 presents a three-way cross-model comparison of peak observables among this work, the Kwak continuum model, and the Schanz et al.EBMD reference, all under identical driving conditions. Table 5 collects the full eleven-configuration TBMD parameter sweep. We first discuss scalar peak observables in this subsection; the qualitative structure of the collapse—the relative timing of the velocity and temperature maxima, the radial profiles, and the shape of the emission pulse—is taken up separately in Section 7.5.

The reported \(T_\mathrm{max}\) in both tables is evaluated over radial bins with at least ten particle counts (after the spatial smoothing described in Section 6.3), to avoid domination by statistically underpopulated wall-adjacent bins. The \(f_\mathrm{ion}\) column reveals two trends: it decreases with increasing \(\alpha_t\) because stronger wall thermalization lowers the collapse temperature available for ionization, and it increases with \(N_\mathrm{ensem}\) at fixed physics parameters because finer EP resolution captures the hot-core ionization events more completely.

The Kwak continuum (PDE) row in Table 4 uses our implementation of the Kwak & Na model [24] solved under the same acoustic and thermophysical conditions as the MD runs. This model predicts \(T_\mathrm{max} = 34{,}927\,\text{K}\) and \(R_\mathrm{min} = 0.883\,\mu\text{m}\). The elevated temperature, arising because the continuum formulation contains no explicit ionization energy sink, is quantitatively on the same order as our no-ionization TBMD runs (\(T_\mathrm{max} = 39{,}617\)–\(45{,}820\,\text{K}\)), confirming that the Kwak PDE and our MD approaches predict qualitatively similar collapse extrema when ionization is absent.

Several trends are immediately apparent within the TBMD parameter sweep. First, our no-ionization runs produce \(T_\mathrm{max}\) values substantially above the Schanz reference value of \(18{,}870\,\text{K}\), confirming that ionization is an essential physical process for obtaining realistic temperature predictions. The absence of such extreme temperatures in the Schanz et al.results, despite their model not incorporating an explicit ionization energy sink during the MD itself, is attributable to two features of their modeling approach. Their event-based hard-sphere formulation yields a different compressibility and thermal equation of state than our soft-sphere LJ model. More significantly, their model includes water vapor chemistry: a substantial quantity of \(\mathrm{H_2O}\) molecules is trapped in the bubble during collapse and subsequently dissociates, providing an internal energy buffer that prevents extreme center-temperature excursions [10]. In our argon-only TBMD model, the explicit collision-based ionization algorithm [9] plays the analogous energy-limiting role, and our ionized configurations consequently bracket the Schanz temperature range closely. Second, for configurations that include ionization, \(T_\mathrm{max}\) falls in the range \(17{,}272\)–\(21{,}265\,\text{K}\), bracketing the Schanz value comfortably. Our best-resolved result, \(N_\mathrm{ensem} = 10^8\) with full-range ionization and \(\alpha_t = 1.0\), yields \(T_\mathrm{max} = 17{,}796\,\text{K}\) and \(T_\mathrm{av,max} = 10{,}852\,\text{K}\). The \(\alpha_t = 0.5\) full-ionization run (\(T_\mathrm{av,max} = 12{,}404\,\text{K}\)) offers a closer match to the Schanz average temperature, consistent with their adopted value of \(\alpha_t = 0.3\) lying between our \(0.0\) and \(0.5\) cases. Third, our minimum bubble radii (\(R_\mathrm{min} = 0.488\)–\(0.733\,\mu\text{m}\)) are systematically smaller than the Schanz value of \(0.750\,\mu\text{m}\). We attribute this in part to the absence of water vapor in our argon-only model: in Schanz et al., the trapped vapor exerts additional internal pressure that resists the final compression and arrests the collapse at a larger radius. The maximum wall velocity in our simulations (\(936\)–\(1{,}271\,\text{m\,s}^{-1}\)) is correspondingly lower than the Schanz value (\(1{,}510\,\text{m\,s}^{-1}\)), consistent with this softer rebound. Wall pressure values are broadly comparable across models. Overall, the agreement between our ionization-inclusive TBMD results and the EBMD reference is physically reasonable given the different interaction models and the absence of vapor chemistry in the present work.

7.5 Qualitative comparisons↩︎

Section 7.4 compared scalar peak values; an equally important test of a model is whether it reproduces the shape of the collapse. Despite differences in interaction model, gas species, and liquid host, all MD simulations of single-bubble sonoluminescence share a set of universal structural features: the wall velocity and average temperature both reach their maxima before the bubble arrives at its minimum radius, a pronounced hot core develops in the radial thermal structure, and the emission pulse is narrow and concentrated near the centre. We examine these features in turn through Figures 10, 5, and 11, using our best-resolved TBMD run (\(N_\mathrm{ensem}=10^8\), full-range ionization, \(\alpha_t=1.0\)) and comparing against the prior MD results of Schanz et al. [10] and Kim et al. [14]. We do not revisit the scalar comparison of Table 4; instead we reuse the standalone Kwak continuum solution introduced there as a shared baseline. Because Kim et al.adopt the same Kwak & Na model [24] as their theoretical reference, this baseline links all three MD frameworks, even though the differing gas–liquid system of Kim et al.permits only the qualitative comparison made here.

7.5.0.1 Collapse dynamics

Figure 10 shows the bubble radius, wall velocity, and volume-averaged temperature over the final \(1.5\,\text{ns}\) before and after minimum radius, for both the TBMD run and the standalone Kwak PDE solution (the Keller–Miksis wall driven by the Kwak continuum gas model rather than by the MD gas pressure), plotted relative to each model’s own \(t_\mathrm{min}\). The radius trajectories agree closely, confirming that the KM-coupled continuum model reproduces the large-scale collapse kinematics faithfully. The wall velocity (panel b) peaks inward before \(t_\mathrm{min}\) and reverses sign at rebound; the corresponding peak inward speeds (TBMD \(1{,}133\,\text{m\,s}^{-1}\), Kwak PDE \(728\,\text{m\,s}^{-1}\)) are the values listed in Table 4. The volume-averaged temperature (panel c) also peaks before \(t_\mathrm{min}\) in both the MD and PDE solutions. This pre-minimum timing of both the velocity and temperature maxima is a universal feature of all MD SL simulations: it reflects the fact that the gas heats most strongly when the inward wall speed is greatest, which occurs during the final inward rush rather than at the moment of maximum compression [10], [14].

7.5.0.2 Radial structure at collapse

The radial profiles already shown in Fig. 5 also illustrate this universal collapse structure. All three runs share the same center-peaked temperature structure (panel a): the hot core reaches approximately \(15{,}000\)–\(16{,}000\,\text{K}\) and decreases to roughly \(4{,}500\)–\(6{,}000\,\text{K}\) at the outer edge, with the \(N_\mathrm{ensem}=10^8\) profile the smoothest and lower-resolution runs exhibiting statistical oscillations in the high-density inner core. This bell-shaped radial temperature distribution is consistent with the profiles reported by Schanz et al. [10] and the center-to-wall gradient described by Kim et al. [14] for different gas–liquid systems, confirming it as a consistent feature of MD sonoluminescent collapse. The density profile (panel c) shows a nearly uniform elevated interior with a density enhancement near the bubble wall: during the final collapse stage the inward wall velocity substantially exceeds the thermal diffusion speed of the outer gas shell, so particles accumulate at the moving boundary faster than they can redistribute, producing a wall-adjacent density pile-up consistent with the advection-dominated compression regime. This wall density enhancement is reported in prior MD studies [9], [14] and is more pronounced for \(N_\mathrm{ensem}=10^6\) because the larger EP diameter amplifies the per-particle contribution to the outer bins, a coarsening effect that diminishes as \(N_\mathrm{ensem}\) increases. The pressure profile (panel b) decreases monotonically from center to wall at \(t_\mathrm{min}\), with the interior pressure approximately twice the wall-region value for the best-resolved run; the gradient is more pronounced at higher \(N_\mathrm{ensem}\), where reduced statistical noise reveals finer radial pressure structure. This center-to-wall gradient is consistent with the non-uniform pressure profile predicted by the Kwak model (Eq. 10 ), in which the inertial correction during collapse introduces a quadratic \(r^2\) dependence that reduces the pressure toward the bubble wall.

7.5.0.3 Emission pulse

Figure 11 shows the total spatially integrated Saha–LTE emission \(L_\mathrm{tot}(t)\) for all three ensemble sizes at fixed full-range ionization and \(\alpha_t=1.0\), alongside a normalised comparison and the radial emission distribution at peak emission time. In all three cases the emission pulse peaks before minimum radius: the peak occurs \(0.221\,\text{ns}\), \(0.315\,\text{ns}\), and \(0.440\,\text{ns}\) before \(t_\mathrm{min}\) for \(N_\mathrm{ensem} = 10^6\), \(10^7\), and \(10^8\) respectively. The peak shifts further before \(t_\mathrm{min}\) as \(N_\mathrm{ensem}\) increases, consistent with the finer spatial resolution of the inner core resolving the early onset of hot-core formation more precisely. The absolute peak emission increases with \(N_\mathrm{ensem}\): \(3.45\times10^{-4}\,\text{mW}\), \(4.06\times10^{-4}\,\text{mW}\), and \(5.28\times10^{-4}\,\text{mW}\) for \(10^6\), \(10^7\), and \(10^8\). The normalised panel (b) shows that the pulse shape is consistent across resolutions and the pulse width is narrow, on the order of \(0.1\,\text{ns}\), consistent with the \(0.070\)–\(0.170\,\text{ns}\) full-width at half-maximum reported by Schanz et al. [10] for similar driving conditions. Panel (c) shows the radial emission profile at peak emission time: emission is strongly concentrated near the bubble centre, consistent with the central hot-core structure seen in panel (a) of Fig. 5 and with the spatially resolved emission profiles of Schanz et al. (their Fig. 13).

Kim et al. [14] study a xenon bubble in sulfuric acid with equilibrium radius \(R_0 = 15\,\mu\text{m}\) and driving frequency \(f = 37.8\,\text{kHz}\)—a different gas species, larger equilibrium radius, and acidic solvent. Even so, the same qualitative signature appears in their results: temperature peaks before minimum radius, a hot core is present in the radial profile, and the emission pulse is short relative to the acoustic period. These commonalities across hard-sphere EBMD [10], [14], soft-sphere TBMD (this work), and the shared Kwak continuum reference confirm that the essential mesoscopic physics of sonoluminescent collapse is robust to the choice of molecular interaction model and simulation methodology.

8 Convergence and Computational Performance↩︎

Two aspects of the numerical implementation warrant examination beyond the physics results: the accuracy of the time integrator as a function of LJ cutoff radius, and the computational cost of the production ensemble sizes.

8.1 Integrator convergence↩︎

LAMMPS employs the velocity–Verlet time integration scheme [65], [66], which achieves global second-order accuracy in time when the pair potential is at least \(C^1\) continuous [67]. When a truncated Lennard–Jones potential is used with a short cutoff, however, the force exhibits a discontinuity at \(r = r_c\), breaking this continuity and degrading the formal convergence order of the integrator.

To quantify this effect, we designed a controlled 12-particle test system in which particles are placed symmetrically on the coordinate planes and directed toward the origin with equal speed components. This configuration isolates pairwise LJ forces from boundary and many-body artifacts, providing a clean measurement of integrator accuracy. Figure 12 shows the measured convergence order as the timestep is refined for several LJ cutoff values. With cutoffs in the range \(10\)–\(20\,\text{\AA}\), the scheme recovers the expected second-order scaling; with shorter cutoffs the rate degrades substantially, confirming that loss of force continuity directly reduces integrator accuracy. This behavior is consistent with the potential shapes shown in panel (b): at short cutoffs the interaction energy drops abruptly at \(r_c\), whereas at larger cutoffs the tail is nearly flat and the truncation is effectively smooth.

These results justify the cutoff choice used in the large-scale SL simulations. Although computational constraints preclude arbitrarily large cutoffs at \(N_\mathrm{ensem} = 10^8\), an excessively short cutoff would bias both thermodynamic observables and integrator accuracy. The selected value \(r_c \approx 2.15\,\sigma\) is large enough that the truncated LJ interaction is smooth in the dynamically relevant separation range, preserving near-second-order integrator behavior, while keeping the per-timestep force evaluation tractable.

8.2 Computational performance↩︎

The three production ensemble sizes span nearly three orders of magnitude in CPU-hours: \(N_\mathrm{ensem} = 10^6\) requires 983 CPU-hours on 128 cores, while the \(10^7\) and \(10^8\) runs each use 2,048 cores and require 55,600 and 975,000 CPU-hours, respectively. Figure 13 shows how runtime is distributed across LAMMPS components at all three scales. The distribution is strikingly stable: pair-force evaluation (LJ + DSF Coulomb) and ionization together consume just over half the total compute time at every scale, with the ionization share holding at 21–22% across three orders of magnitude in \(N_\mathrm{ensem}\). This invariance reflects two competing effects: as \(N_\mathrm{ensem}\) decreases, the enlarged EP reactive radius \(r'_\mathrm{ion} = g^{1/2}r_\mathrm{ion}\) increases ionization neighbor-list work per particle, while the reduced absolute particle count offsets this cost in equal measure.

The temporal distribution of computational work reveals a further structure. Figure 14 plots the normalized step throughput and adaptive timestep over the 10-ns simulation window for all three ensemble sizes. Throughput is approximately 15-fold higher during the pre-collapse phase than at the collapse minimum near \(t \approx 4.2\)–\(4.7\,\mathrm{ns}\): as the bubble radius decreases, particles concentrate near the center and sharply increase both pair-force and ionization neighbor-list work per timestep. The symmetric recovery after the rebound confirms that the bottleneck is geometric in origin, tied directly to the collapse dynamics rather than to any cumulative state change such as a growing ionization charge distribution.

9 Conclusions and Future Work↩︎

This work has developed a time-based molecular dynamics framework for simulating single-bubble sonoluminescence within a hybrid continuum–MD formulation. Unlike prior event-based codes, which are restricted to hard-sphere interactions and event-driven time advancement, the present approach integrates a continuous Lennard–Jones plus damped-shifted-force Coulomb potential at each timestep, enabling self-consistent treatment of ionization state and long-range electrostatics throughout the collapse. To make such simulations tractable at physical particle counts, an ensemble-particle (EP) scaling formalism maps \(N_\mathrm{real}\sim 10^{10}\) physical atoms onto up to \(N_\mathrm{ensem} = 10^8\) simulated EPs while preserving temperature, pressure, and ionization statistics, representing a tenfold increase over the largest ensembles reported in prior EBMD studies. Gas dynamics couple bidirectionally to a continuum Keller–Miksis solver at each timestep, and the resulting predictions are benchmarked against the EBMD reference dataset of Schanz et al. [10] under identical driving conditions.

The parameter sweep over ionization model and thermal accommodation yields several clear conclusions. Ionization is by far the dominant regulator of peak temperature: removing the ionization energy sink raises \(T_\mathrm{max}\) from the physically plausible range of \(17{,}272\)–\(21{,}265\,\mathrm{K}\), which brackets the Schanz et al.reference of \(18{,}870\,\mathrm{K}\), to \(39{,}617\)–\(45{,}820\,\mathrm{K}\), a factor-of-two inflation that persists regardless of the electrostatic model chosen. The thermal accommodation coefficient \(\alpha_t\), by contrast, primarily governs the spatially averaged temperature at collapse with limited influence on the peak local compression temperature, a finding consistent with the thermal boundary layer being thin relative to the bubble radius at the collapse minimum. Scalar observables show consistent, systematic behavior as \(N_\mathrm{ensem}\) increases from \(10^6\) to \(10^8\), validating the EP scaling formalism and confirming that the \(10^8\) configuration provides the most reliable reference for all reported quantities.

Several directions warrant future investigation. The most direct is to eliminate the EP approximation altogether: with the scaling formalism validated at \(N_\mathrm{ensem} = 10^8\), a direct \(g = 1\) simulation at the full physical particle count (\({\sim}10^{10}\)) would eliminate modeling and approximation errors due to the use of EPs and provide a numerical ground truth against which the EP results can be assessed. A second direction is the extension to reactive gas mixtures, since the present model is restricted to monatomic argon; incorporating water vapor, dissolved noble-gas impurities, or nitrogen would require a reactive force field or a QM/MM description of bond-breaking chemistry, and would enable direct comparison with spectroscopic observations of molecular emission during the collapse pulse. Another interesting direction would be to consider other fluids and fluid models surrounding the sonoluminescing bubble, such as compressible fluids modeled with more complex equations of state [68]. Finally, replacing the empirical Lennard–Jones potential with a machine-learning interatomic potential (MLIP) trained on ab initio data at SL-relevant conditions would remove the assumption that LJ parameters calibrated at ambient conditions remain valid at the extreme temperatures and pressures reached during the collapse.

Acknowledgments↩︎

This work was funded in part by an ORAU Ralph E.Powe Junior Faculty Enhancement Award. The authors were supported by NSF Grant No.. This work used Anvil at Purdue University through allocation PHY250136 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program [69], which is supported by U.S.National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. D.H.thanks Seth Putterman at UCLA for very helpful conversation about sonoluminescence at the beginning of this project. S.C.thanks Ho-Young Kwak at Chung-Ang University for many generous and insightful emails about sonoluminescence.

10 Supplementary Space–Time Heatmaps↩︎

The following nine panels provide supplementary space–time views supplementing the results of Sections 7.1–7.2.

11 Radial Profiles Across Parameter Sweeps↩︎

The following nine panels provide supplementary radial-profile views across all three parameter sweeps in a unified format.

12 Scalar Parameter Sweep Summary↩︎

The following six panels visualise all peak scalar observables across the full eleven-configuration parameter sweep (Table 5).

References↩︎

---

1. https://github.com/csh-apprentice/lammps-tbsl↩︎

2. Note, the actual particle counts in these simulations are on the order of \(5.5 \times 10^4\)–\(1.5\times 10^6\); they simulate one cone of a bubble, which by symmetry they argue gives results equivalent to simulating a bubble with \(5\times10^7\) particles, but since they only simulate a cone of the bubble their total particle counts are never that large.↩︎