Molecular dynamics simulations of cavitation bubble collapse and sonoluminescence

The dynamics of the medium within a collapsing and rebounding cavitation bubble is investigated by means of molecular dynamics (MD) simulations adopting a hard sphere model for the species inside the bubble. The dynamics of the surrounding liquid (water) is modelled using a Rayleigh–Plesset (RP)-type equation coupled to the bubble interior by the gas pressure at the wall obtained from the MD calculations. Water vapour and vapour chemistry are included in the RP-MD model as well as mass and energy transfer through the bubble wall. The calculations reveal the evolution of temperature, density and pressure within a bubble at conditions typical of single-bubble sonoluminescence and predict how the particle numbers and densities of different vapour dissociation and reaction products in the bubble develop in space and time. Among the parameters varied are the sound pressure amplitude of a sonoluminescence bubble in water, the noble gas mixture in the bubble and the accommodation coefficients for mass and energy exchange through the bubble wall. Simulation particle numbers up to 10 million are used; most calculations, however, are performed with one million particles to save computer run time. Validation of the MD code was done by comparing MD results with solutions obtained by continuum mechanics calculations for the Euler equations.

Abstract. The dynamics of the medium within a collapsing and rebounding cavitation bubble is investigated by means of molecular dynamics (MD) simulations adopting a hard sphere model for the species inside the bubble. The dynamics of the surrounding liquid (water) is modelled using a Rayleigh-Plesset (RP)-type equation coupled to the bubble interior by the gas pressure at the wall obtained from the MD calculations. Water vapour and vapour chemistry are included in the RP-MD model as well as mass and energy transfer through the bubble wall. The calculations reveal the evolution of temperature, density and pressure within a bubble at conditions typical of single-bubble sonoluminescence and predict how the particle numbers and densities of different vapour dissociation and reaction products in the bubble develop in space and time. Among the parameters varied are the sound pressure amplitude of a sonoluminescence bubble in water, the noble gas mixture in the bubble and the accommodation coefficients for mass and energy exchange through the bubble wall. Simulation particle numbers up to 10 million are used; most calculations, however, are performed with one million particles to save computer run time. Validation of the MD code was done by comparing MD results with solutions obtained by continuum mechanics calculations for the Euler equations.

Introduction
In a strong spherical collapse of a cavitation bubble, the bubble medium is highly compressed and heated. This is demonstrated, for example, by the phenomenon of cavitation bubble luminescence, i.e. by light emitted by a bubble in the final stage of collapse (see [1] and, for reviews, [2][3][4]). Furthermore, chemical reactions, e.g. dissociation and radical formation, take place in the bubble as it turns into a chemical microreactor in these conditions [5][6][7][8][9][10]. Moreover, in a strong bubble collapse the state of the bubble medium may become inhomogeneous with compression and possibly shock waves arising [11,12], see, however, [13,14]. Thus hydrodynamical, thermodynamical and chemical evolution have to be considered simultaneously and, moreover, they have to be resolved on small space (of the order of µm and even nm) and short time scales (of the order of ns and even ps) for correct modelling.
Collapsing cavitation and sonoluminescence (SL) bubbles have been approached theoretically by constructing models of increasing sophistication. There are two very different lines of approach, the first one treating the liquid and bubble interior as a continuum that can be split down to arbitrarily small parts and the other one observing the discreteness of matter that may enter into the problem on small scales, representing the bubble and its surrounding liquid by interacting particles in molecular dynamics (MD) simulations. A third, hybrid approach treats only the bubble interior by MD simulations and treats the surrounding liquid, on the other hand, as continuum.
The conventional approach is via continuum models. An equation for the bubble radius (or surface) as a function of time is formulated, leading in its simplest forms to ordinary differential equations (see the books [15,16] and the review [17] as a starting point) and in more elaborate forms to a set of partial differential equations of motion for the liquid and bubble interior together with partial differential equations for heat and mass diffusion and separate equations for the equation of state and the evolution of chemical species with their reactions and spatial distribution (see, in particular, [6][7][8][9][10][11][12][13][14] and also [18][19][20][21][22][23][24][25][26][27][28][29][30] for a variety of models in various states of approximation). The final continuum approach including all effects has not yet been done, as too many phenomena are involved that are not yet properly solved separately. Among these are evaporation and condensation problems in the presence of different species and reaction rates at extreme pressures and temperatures as well as plasma formation.
The second line of numerical modelling-MD simulation-resorts to discretization, where the matter is made up of interacting particles. The bubble interior is represented by its molecules, modelled by interacting and reacting particles, instead of by a homogeneous medium with an appropriate equation of state and chemical reaction rate equations. The liquid confining the bubble may also be modelled by particles, for instance as a Lennard-Jones (LJ) fluid, i.e. by particles interacting via the (continuous) LJ pair potential, or else in a hybrid model, as a continuum with prescribed liquid motion or, interactively, with the MD calculations coupled to a corresponding continuum equation.
The first significant approach to bubble collapse and bubble luminescence with an (MD) particle model known to the authors starts with 725 hard sphere particles in a spherical volume (the bubble) made shrinking at a prescribed constant velocity with Mach numbers of 0.5, 1.0, 1.5 and 15 with respect to the sound velocity in the undisturbed gas [31]. The particles move according to Newton's laws. Reflecting boundary conditions are adopted at the wall. Twelve independent computations were averaged for smoothing the results (temperature, density and pressure versus distance from the wall) that nevertheless show substantial fluctuations because of the low number of particles. No realistic conclusions can be drawn from this start-up model, for instance on temperatures, but some trends can be predicted. The temperature inside the bubble increases 'dramatically' [31] with increasing Mach number. Heavy, non-reacting gases (particles) gain more energy from the wall collisions than lighter ones. This means that heavier particles collapsed at the same Mach number reach higher temperatures. This is as found experimentally for sonoluminescing bubbles when going from He to Xe as gas in the bubble [2,3] and confirmed later by continuum mechanics calculations [13]. Up to a Mach number of 1.5, the bubble interior is found to be almost homogeneous. At Mach number 15, a strong increase in temperature develops towards the centre of the bubble. Similar calculations with more than a million hard sphere particles in a bubble collapsed at a Mach number of about 1.5 confirm the strong increase in temperature towards the centre of the collapsed bubble in this kind of gas ( [32], figure 21). Moreover, a comparison with a continuum mechanics calculation based on Euler's equations showed reasonable agreement [33][34][35][36]. The differences (a steep shock for the Euler equations and the singularity at the bubble centre) can be explained by the dissipation missing in the Euler equations but being naturally included in the MD calculations. Twodimensional (2D) MD simulations with hard spheres and again a constantly contracting cavity have been compared with continuum solutions based on Euler's equations and Navier-Stokes equations [37]. The Navier-Stokes calculations are found to be in excellent agreement with the MD calculations validating the MD approach and vice versa. The results can be transferred to 3D bubbles [37].
The work with prescribed constant bubble wall motion has been extended to a prescribed time-dependent wall motion by taking the Rayleigh equation with a van der Waals gas inside the bubble as liquid motion input to the MD calculations [38,39]. For the MD calculations, particles simulating pure noble gas atoms are taken. For consistent results between inner motion and outer liquid motion, the course of bubble collapse has to be estimated beforehand as prescribed motion, because there is no backreaction from the MD calculations on the liquid flow. Ionization of the particles is introduced as an energy loss mechanism. Reflecting and heat bath boundaries are applied, and besides hard spheres, variable soft spheres [40,41] are also used. Although no realistic temperature values can be given, as, for instance, water vapour is not included, calculations confirm the trend in temperature rise from light to heavy particles (He to Xe) [13,31], the build-up of strong inhomogeneities in the bubble upon strong collapse [32,35,36] and mixture segregation [25,28,36]. From the duration of high temperature inside the bubble in dependence on the number of particles (three cases: 10 4 , 10 5 and 10 6 ) a proportionality of bubble luminescence pulse width to the bubble radius is derived as found experimentally ( [32], figure 38; see also [42], figure 10).
Besides these hybrid starts with prescribed wall motion, there have also been pure MD approaches for a collapsing bubble, where soft sphere particles (interaction potential ∝ r −12 , r being the centre-of-mass distance) in the bubble interact with their likes and with a surrounding LJ fluid [43]; that is, full-scale MD simulations with long-range interaction forces have been done. In the investigation, the liquid is made up of 62 500 particles in a box, in which a spherical void is suddenly produced and filled by insertion of different numbers (0, 500 and 1372) of non-condensable soft sphere particles besides the vapour particles that evaporate into the void during a period of equilibration. A sudden compression lets the system develop. The highest temperatures arise when there are no non-condensable particles inside, i.e. there is only water vapour. Condensation has been found to be significant. When the bubble shrinks by 1/5 of its original size, about half of the vapour molecules condense at the wall. In a similar means of investigation, an LJ liquid has been taken, where a cavity has been created as a 'hole' in the liquid, followed by a rapid filling of the cavity by particles evaporating from the liquid [44]. An oscillatory filling is observed in this numerical MD experiment, reminding us of the empty Rayleigh bubble collapsing under constant pressure in the liquid [45]. In an extended investigation [46] of the same problem with about 10 7 LJ particles (with truncated potential) of pure liquid argon in a box with a central void, a comparison is made with the empty Rayleigh bubble collapse. As MD calculations naturally include the compressibility of the liquid, a tension wave is observed running into the liquid. The filling of the void with argon gas and the subsequent compression lead to high temperatures in the centre. Thus a different, albeit similar, dynamics to that of the empty void in an incompressible liquid is observed.
The hybrid method was extended by Metten et al [33][34][35][36][47][48][49] through the coupling of the bubble interior, which is modelled by MD with hard spheres, with the continuum Rayleigh-Plesset (RP) equation for the surrounding liquid. The coupling of the two fluid parts is effected by determining the MD gas pressure inside the bubble and applying the boundary pressure as input gas pressure to the RP equation (RP-MD model). Furthermore, there must be stated collision rules with the (moving) wall. Reflecting as well as thermal boundary conditions are introduced. Also water vapour with evaporation and condensation at the bubble wall is considered. Extensive calculations are done (mostly with about one million particles, this being found to be necessary for sufficient resolution) for a typical RP sonoluminescence bubble under various conditions (adiabatic, isothermal, quadratic thermal boundary layer, different gas mixtures, evaporation and condensation, and chemical reactions). Also aspherical bubbles are considered for constantly collapsing and RP-MD coupled collapsing bubbles. The main results are that, for a typical sonoluminescence bubble, there are strong compression waves running to the centre upon collapse even when water vapour, its dissociation and chemical reactions are included; that water vapour is trapped at the centre of the bubble and counteracts the temperature increase [9]; that mixture segregation occurs as predicted by continuum mechanics calculations [25]; and that the energy focusing mechanism is quite robust to perturbations of the collapsing surface, as might be expected from experiments [50].
MD calculations with 10 6 ensemble particles and coupling to the Keller-Miksis equations have been done by Kim et al [55], mainly for a xenon bubble in sulfuric acid with heat transfer at the bubble wall, but neglecting mass transfer (evaporation and condensation) and chemical reactions. The results are compared with a Navier-Stokes model that gives results in between the extremes of adiabatic and heat bath boundaries for the MD calculations, confirming gross features of the theoretic-numerical approaches in either direction. The hard sphere model represents a realistic description of particle interaction in highly inhomogeneous systems (avoiding unacceptably small time steps as with a continuous interaction potential) and/or in high-temperature environments with high particle velocities, where the attractive part of the LJ potential is negligible [56]. The collapse of a bubble, in particular the most interesting final stage, is characterized by strong density waves and by centre temperatures in excess of 10 000 K [2,3]. In this scenario the hard sphere model is a reasonable approach.
New results may be expected from MD simulations of collapsing bubbles making use of the intrinsic, favourable MD properties of the approach, for example the generic inclusion of interior mass diffusion and heat conduction. Therefore, this work presents the results of the influence of various parameters in the MD simulation approach on the hydrodynamic and thermodynamic quantities inside a collapsing cavitation bubble in water. In particular, water vapour and its chemistry are included in the investigation.

Hybrid Rayleigh-Plesset molecular dynamics model
The basis of the model used is a hard sphere MD code suggested by Rapaport [56]. Every possible action considered relevant is represented as an event, occurring at a precise time in the future. Events can be particle specific (collisions, wall collisions and cell crossings) or system specific (snapshot of current state, end of the run, etc). The chronological order of these events is managed by a binary structure, the event tree. The event tree gives the time and a set of determining data of the event closest to the current time. The system time is advanced to this event time. As hard sphere particles have no long-range interaction, only the directly participating particles have to be considered. Their position vector is updated and the consequences of the event are applied. This way, the simulated particles 'jump' from event to event; a costly calculation of full trajectories is avoided.
The particles have the following properties: mass, position, velocity, diameter, rotation and chemical status. The rotational degrees of freedom are modelled as spin of rough spheres with an angular velocity ω that can be transferred upon particle collisions.
Typical bubbles considered in this work contain a very large number of particles. A sonoluminescence bubble with an equilibrium radius of R e = 4.5 µm holds several 10 11 particles. Obviously, this number is too large for current hardware to do a one-to-one simulation of a bubble of this type. While still achieving reasonable run times, at present it is only possible to use several millions or several tens of millions of simulated particles. To bridge the gap, it is necessary to introduce ensemble particles that represent a certain number of real particles. The necessary transformations are chosen in such a way that the mass (N r m r = N m), energy (N r T r = N T ) and volume (N r d 3 r = N d 3 ) scale correctly, where N is the number of particles, m is the mass of a particle, T is the temperature, d is the diameter of the particle and the subscript r indicates the real particles. As a side effect, the mean free path l = V /( √ 2πd 2 N ) is larger than in reality by the factor (N r /N ) 1/3 . Therefore diffusive properties, such as thermal conduction, viscosity and particle diffusion are also amplified by this factor.
There have been approaches in parallelizing code for hard sphere simulations [57]. This is an interesting aspect to achieve larger particle numbers; however, due to the complexity of our model and the necessary rollback protocol, a parallelization system was not implemented.
The particles are confined to a spherical region (a bubble), whose boundary movement is determined by the MD-enhanced RP equation: P l being the pressure in the liquid at the bubble wall: and p ∞ being the far field pressure: R is the bubble radius,Ṙ is the velocity of the bubble wall,R is the acceleration of the bubble wall and a superscript dot means derivation with respect to time, t. ρ l denotes the density of the liquid, c l the speed of sound in the liquid, σ the surface tension, ν the kinematic viscosity and P(R, t) the gas pressure at the bubble wall. p 0 is the static pressure, p ac a time-varying pressure, here taken as a sinusoidal sound wave with sound pressure amplitude p a and frequency ν a . The viscosity of the liquid is represented by the last term in (2) and radiation of sound by the last term in (1). The gas pressure at the bubble wall, P(R, t), and its time derivative,Ṗ(R, t), are taken directly from the hard sphere simulation via the energy transferred by a given number of particle-wall collisions. These values are used in (1) for some time interval t int . The solutions obtained are thereupon used in the wall collision prediction during the next interval. The interval duration t int is calculated so that the number of wall collisions matches a specified number. This value should be chosen as high as possible without sacrificing accuracy. Therewith, the interior of the bubble is coupled to the bubble wall. This hybrid model is called the RP-MD model.

Energy transfer
When particles collide with the collapsing bubble wall, heat is transferred to the bubble medium. In the model, this is accounted for by a change of the particle velocity and an according change in the bubble wall temperature. The heat quantity P Q released by N colliding particles during a time t is calculated by where T i is the temperature of the incoming particles, T W is the wall temperature and k B is the Boltzmann constant. The temperature of the particles after reflection at the wall is given by By the choice of the thermal accommodation coefficient, α t , the amount of heat conduction can be influenced. By varying α t , any condition from an isothermal boundary (α t = 1) to a completely adiabatic boundary (α t = 0) can be realized.
The energy transferred will heat the water around the bubble. For an estimate of the bubble wall temperature, a thermal boundary layer in the liquid is used [58,59]: T l (r, t) is the liquid temperature for r > R, and T ∞ is the liquid temperature far away from the bubble. Equation (6) describes a quadratic temperature profile of thickness δ.
In continuum models with heat transfer, boundary layers are also used at the inside (gas side) wall of the bubble [6,19,60]. This is not necessary here, as in MD calculations this is accounted for by virtue of the method.

Water vapour
During the comparatively long period of bubble growth vast numbers of molecules evaporate into the bubble. Most of these will condense back into the liquid through the bubble wall in the collapse phase. The initial vapour content for the simulation is calculated by the ideal gas law and the assumption that the partial pressure of the vapour equals the saturation vapour pressure, valid for a bubble not being in the final collapse phase [9,60]. Equation (2) has to be extended by the saturation vapour pressure, p * v (T ∞ ). For the molecular simulation, an expression for the evaporation or condensation probability is needed. Following Yasui [60], a phase change probability P pc for water is calculated: Here, α v is the accommodation coefficient for evaporation and condensation, also called the mass accommodation coefficient or sticking coefficient, as it is equal to the ratio of vapour molecules sticking to the phase interface to those impinging on it [19]. p v (R) is the current vapour pressure at the bubble wall and C a correction factor near one. If P pc is positive, a vapour particle colliding with the bubble wall condensates with probability P pc . The negative case (P pc < 0) describes evaporation: a new particle will be put into the bubble with probability −P pc at the opposite wall of the bubble and with velocity opposite to the colliding particle. The choice of α v is of prominent importance for the bubble conditions, as will be shown later. Equation (7) is valid only up to the critical temperature of water, T cr = 647 K. For higher temperatures, phase change is not considered. This limitation of the model will induce only minor effects on the results, as in the final stages of the collapse most of the vapour is trapped in the bubble centre and condensation is negligible. Furthermore, the time scale of the bubble wall motion is much faster than that of diffusion [9].
For some calculations, a temperature-dependent mass accommodation coefficient, α v (T W ), has been used as in [7]: A similar functional dependence was investigated recently for its influence on bubble oscillations in [30,61].
The temperature-dependent mass accommodation coefficient function, α v (T W ), is scaled by a factor, A v , for stretching or compressing the curve for simulation applications: When this function is used in the calculation, the mass accommodation scaling factor A v is given.

Vapour chemistry
Due to the high temperatures reached in the collapse phase, much of the vapour remaining in the bubble will be dissociated and chemical reaction will take place [62].  [63]. These reactions can be divided into two main groups, which have to be treated differently: bimolecular and catalytic reactions. The former draws the required energy solely from the kinetic energy of the particles; for the latter, a third particle is needed, acting as a reaction catalyst. For bimolecular reactions of the kind A + B C + D the following thoughts are relevant. According to classic collision theory, particles need to collide in order to react. Furthermore, the collision energy, E c , has to exceed the reaction activation energy, E a , and the relative orientation of the molecules has to be suitable for an effective collision. All molecules are treated as spheres; thus orientation considerations are not applicable to the model. Whether a reaction takes place is determined by a reaction probability, which is given as the ratio of the total reaction cross section, σ tr , to the total impact cross section, σ ti = π(d 1 + d 2 ) 2 /4. Note that σ tr can be determined as with m AB , v AB and δ AB being the reduced mass of the educts, the modulus of the vectorial difference in velocity and the Kronecker symbol, respectively. is the gamma function, and and η are constants in the rate coefficients k(T ) = T η exp (−E a /(k B T )) of the macroscopic rate equations of the reactions [64]. Equation (10) gives σ tr in dependence on the particles involved and their excess energy. If the probability evaluation shows that the reaction takes place, the educts are changed to the products and are given new velocities, preserving the mass, momentum and energy.
As there are no three-body collisions incorporated into the model, catalytic reactions of the kind A + B + M C (+D) + M require the use of a two-step method. In the first step, an activated molecule (AB) * is formed when two potentially reacting particles A and B collide. This complex is treated as a single particle with a certain lifetime t comp . If (AB) * collides with a third particle M, it reacts into the product C (+D). If no such collision occurs within t comp , Table 1. Chemical reactions for water molecules and their constituents as provided by GRI-Mech 3.0 [63] and by [6,7]. is given in m 3 mol −1 , the activation energy E a /k B in K and the reaction enthalpy E D in kJ mol −1 . η is dimensionless. M denotes a catalytic particle.

No.
Reaction (AB) * decays into the educts. This is the 'Lindemann mechanism' for gas phase reactions. Both steps of the mechanism are treated as a bimolecular reaction with E a = 0 and σ tr /σ ti = 1. The reactions used in the model are given in table 1.

Light emission
In order to simulate the process of light emission, a model based on Bremsstrahlung is used [11]. Discussing the details of the light emission processes is outside the scope of this paper [65,66]. The occurrence of Bremsstrahlung requires high temperatures ionizing the gas. When the-now charged-gas particles are decelerated or accelerated, light is emitted. It is assumed that the system is locally in thermal equilibrium and that only first-order ionization occurs.
Using the Saha equation the total Bremsstrahlung output per volume integrated over all frequencies, P Br , can be calculated as Here, q denotes the degree of ionization, n is the number density of atoms, T is the temperature and χ 1 is the first ionization energy (He: 24.59 eV; Ar: 15.77 eV; Xe: 12.13 eV). An integration over only the part of the spectrum visible in water (λ > 200 nm) yields the results used in this work. It has to be noted that ionization is not done on a particle level, nor does the emission have any impact on the hydrodynamics of the model. The temperature and the number of the atoms in a cell of given size are averaged to achieve a spatially resolved light emission profile.

Simulation conditions
The initial distribution of the particles in the bubble (full spherical volume) is done by placing the particles on a regular lattice. Each component of the velocity of the particles is assigned a random value according to the Maxwell-Boltzmann distribution: The simulation is started either at the maximum radius of the bubble cycle (evaluated by a simple polytropic model for a van der Waals gas) or at a point where the bubble wall is already moving inward. In the latter case, the initial particle velocities have to be modified considering the movement. Distributing the particles randomly proves to be unnecessary owing to fast mixing of the system. To simplify the calculations, the following values are set to unity: the diameter of the largest species, the mean mass of all particles in the initial bubble, the Boltzmann constant k B and the initial temperature T 0 .
In order to obtain the local conditions in the bubble for visualization, the simulation region (full sphere) is divided into a number of measurement cells. Utilizing the spherical symmetry of the system, spherical shells of equal volume are used. The properties of all particles residing in a given cell are averaged to derive localized values for hydrodynamic and thermodynamic quantities.
For comparison with data in the literature, a bubble in water driven at a sound field frequency of ν a = 26.5 kHz for different sound pressure amplitudes (1.2 p a 1.45 bar) and (equilibrium) bubble radii, R e , (3.85 µm R e 7.72 µm) is investigated. Water is used at T 0 = 300 K, p 0 = 100 kPa with the parameters density ρ l = 996.6 kg m −3 , sound velocity c l = 1481 m s −1 , surface tension σ = 0.072 75 N m −1 and kinematic viscosity ν = 8.569 × 10 −7 m 2 s −1 .

Validation
As only one MD simulation [31] existed for spherical symmetry at the beginning of this project, the validation of the code developed had to be ensured. This was done first by extending the MD calculations with a constantly contracting sphere to particle numbers of 10 6 ([32], figure 21). Next, a comparison was made with the same problem in the (continuum) formulation of the Euler equations [35,36]. For the continuum mechanics calculations a finite-volume Godunov scheme for the spherical Euler equations in Lagrange coordinates for a van der Waals gas has been developed and compared with the corresponding MD calculations (figure 1). The steep shock developing early in the continuum solution is due to the inviscid nature of the Euler equations, whereas the smooth, but steep, temperature wave in the MD calculation is a result of the naturally included diffusion processes.

Results
Introductory results are given for a bubble with typical sonoluminescence properties: the equilibrium radius, R e , is 4.5 µm, the acoustic driving frequency, ν a , is taken as 26.5 kHz with  Figure 2 shows the density and temperature of the bubble for 700 ps during collapse and rebound in its space-time evolution. During the collapse phase the majority of the particles gather near the bubble wall. The wall reaches its maximum velocity (approximately 1400 m s −1 ) about 300 ps before the bubble reaches its minimum radius. As the bubble shrinks rapidly, the density in the outer shells grows and simultaneously many particles are driven to the bubble centre. A fast rise in density can be seen (from 50 to 500 kg m −3 within 500 ps), although the wave never steepens into a shock. The velocity of the particles thermalizes quickly in the bubble centre and temperatures up to 17 500 K are obtained. At the time of maximum core temperature, the bubble density is roughly homogeneous around a value of 525 kg m −3 (about half the density of normal water), whereas the temperature follows a bell-shaped curve from the centre to the wall (see also below for temperature dependences). Subsequently, a density wave runs back to the wall, where most of its energy is absorbed. The moment this wave starts reaching the bubble wall marks the minimum bubble radius. The fast particles, colliding with the wall, induce a high wall pressure, which causes the bubble wall to halt and reverse its movement. At this moment, a temperature of nearly 9000 K of the particles near the wall is found. The bubble centre cools down rapidly. The peak width at half-height of the hot spot is about 400 ps.

Simulation accuracy
Any MD calculation with ensemble particles needs a minimum number of these particles for modelling with sufficient accuracy. To find this number that is simultaneously the number of As particles can diffuse, condensate and evaporate, N i denotes the initial particle number. Already a 'low' particle number of N i = 10 5 is able to predict averaged properties reasonably well (see figure 3). Increasing the particle number by a factor of 100 changes the hydrodynamic and thermodynamic values by less than 10%. Only the vapour content partly shows larger variations ( figure 3(d)). The more particles are used for simulation, the more (scaled, therefore 'real') H 2 O particles remain in the bubble. This effect is due to overestimated particle diffusion, which scales with (N r /N i ) 1/3 . The additional water vapour (figure 3(d)) leads to a larger minimum radius, R min , (figure 3(a)) and reduces average bubble temperatures at R min ( figure 3(b)). Nevertheless, the average density at R min is lower ( figure 3(c)).
For a spatially resolved analysis of the bubble interior, larger numbers of particles give smoother results. Figure 4 shows the temperature and the density both at the centre and near the wall. Some definite trend with (ensemble) particle number can only be seen for the temperature in the bubble near the wall ( figure 4(c)). The other curves are more or less just smoothed versions of the curves calculated for lower particle number. For larger particle numbers, the temperature and density gradients around the time of minimum radius grow (not shown here, see [35]). While for low particle numbers the density distribution is nearly uniform, a gradient towards the bubble wall is developing for higher particle numbers, caused by the slower diffusion. Figure 5 shows a space-time plot of the temperature evolution inside the bubble for the typical parameters chosen, calculated with 10 7 initial ensemble particles. The maximum temperature reached in this case is 16   10 6 ensemble particles already approximate the temperature evolution in space and time. The maximum temperature gets lower with increasing particle number, as the amount of water vapour trapped is higher (see section 3.2). The real particle number in a bubble as considered here is of the order of N r = 10 11 . Judging from the scaling calculations, it can be expected that a simulation with N i = N r yields results that differ from those presented by only about 10%. A simulation with N i = 10 6 ensemble particles therefore should reproduce a direct real particle simulation reasonably well and predict the trends correctly. For this reason, the majority of the simulations are done with about 10 6 (initial) ensemble particles. The introduction of a particle number-dependent mass accommodation

Vapour trapping
Bubble dynamics is sensitive to the amount of water vapour remaining in the bubble upon collapse. Due to its smaller weight (when considering argon as noble gas), the vapour trapped in the bubble tends to accumulate in the bubble centre, where it counteracts the development of high temperatures. At the final stages of the collapse, the bubble can be seen as a quasi-adiabatic system; therefore the additional mass of the water vapour leads to reduced compression heating. Furthermore, thermal energy can be taken up by rotational degrees of freedom of the vapour molecules, even without considering endothermal chemical reactions. Therefore, the choice of the accommodation coefficient for evaporation and condensation of H 2 O, α v , is crucial for the results obtained and their conformity to experimental data. The literature provides values for α v varying from 0.01 to 1.0 [67][68][69]. Fujikawa and Akamatsu [19] give results for a set of values: α v = 0, 0.01, 0.04, 0.1 and 1. Yasui [7] used a temperaturedependent value for α v , ranging from 0.35 (T W < 350 K) to 0.05 (T W > 500 K), reproduced in (8), while Storey and Szeri [9] chose a constant, α v = 0.4.
To clarify the role of α v in the RP-MD model, values of α v between 0.05 and 0.7 have been tested to show its relevance to results and to obtain a reasonable value for the current calculations. No chemical reactions are included in these simulations in order to isolate the effect of α v . Figure 6 shows the development of centre temperature (figure 6(a)) and water vapour density at the centre of the bubble (figure 6(c)) in a short time span before and after the for higher α v . This shift can be explained by the decreasing minimum radius: as fewer particles remain in the bubble, the bubble collapses further in radius and it takes a shorter time for the reflected density wave to hit the bubble wall and stop the collapse. Compared to the results of Storey and Szeri [9], the water vapour distribution in the bubble is much more uniform (see figure 6(d)), regardless of α v . For high values of α v (α v > 0.25), the distribution is virtually homogeneous. Figure 7(a) shows the total number of water vapour molecules in the bubble for different α v and figure 7(b) the corresponding radius-time curve. It can be seen that the differences in water vapour content are most distinct in the main collapse and the first afterbounce (an order of magnitude when going from α v = 0.1 to α v = 0.4). The magnitude of the first afterbounce is affected markedly by the differences in water vapour content. It is noteworthy that the lower the condensation probability the later the time of minimum water vapour content is reached, as the equilibration process is slower. Beyond the afterbounces, the bubble motion becomes slow enough to allow the vapour to equilibrate with the saturation vapour pressure, regardless of α v (t > 5 µs, not shown).
A picture of vapour trapping (with argon as noble gas) is given in figure 8 in a space-time plot for a typical sonoluminescence bubble. It is seen that during final collapse the argon molecules mainly gather at the wall, whereas the water vapour molecules gather in the bubble centre. Table 2 summarizes the findings. The MD calculation with p a = 1.2 bar is included for a comparison with the results of Storey and Szeri [9] (for their case II with α v = 0.4 as used throughout). When using the temperature-dependent accommodation coefficient, α v (T W ), (8) with A v = 0.4 (9), the results are close, except for the maximum centre temperature.
In MD simulations with ensemble particles, diffusion does not scale linearly with the particle number (due to a larger mean free path); therefore a lower accommodation coefficient is needed to obtain similar results as compared to a continuum model (see also the section on simulation accuracy, section 3.1). Experimental [2,70] and continuum numerical [71] results indicate that the temperatures inside the bubble at collapse are in the range of 10 000-20 000 K, assuming a Bremsstrahlung mechanism as the origin of the light emission. This points to a low α v for the MD calculations ( figure 6(a)). Therefore, in the following simulations, a value of α v = 0.1 is used. This value is also near the value of α v = 0.12 [69], recently measured for stationary evaporation and in normal conditions, however. Colussi and Hoffmann [72] also arrive at α v = 0.1 from mean free path arguments and the kinetic theory of gases. They predict spatially inhomogeneous conditions inside the bubble, as found here in MD calculations. In the context of laser-induced spherical bubbles and their collapse and rebound, a value of α v = 0.075 has been found to approximate measurable quantities-such as minimum bubble radius and the first maximum radius after rebound-best [73].
The case of a temperature-dependent mass accommodation coefficient, α v (T W ), (8), with scaling factor, A v , according to (9), has been investigated for its fit with experimental data. The best fit is found for a scaling factor of A v = 0.3-0.4 (see also table 2). This is in accordance with Table 2. Effect of variations in α v and comparison of the RP-MD model with a continuum model (CM) of [9] with respect to minimum bubble radius, maximum temperature, maximum average temperature and water vapour content. Argon and water vapour. T max is the overall highest temperature in the bubble; T av,max is the highest average temperature. The amount of water vapour molecules still in the bubble at minimum radius is given both in absolute numbers and as molar fraction. R e = 4.5 µm, ν a = 26.5 kHz, α t = 0.3.

Energy transfer
During bubble oscillation, energy is exchanged between the bubble interior and liquid. By varying the thermal accommodation coefficient, α t (see section 2.1), the amount of energy transferred between the particles in the bubble and the surrounding medium (water) upon collisions is influenced. The literature provides values for α t between 0.5 and 1 [69]. Yasui [74], in his investigation of thermal conductivity on bubble dynamics, used α t = 1 (and α v = 0.04). As with α v , the conditions of applicability of a specific value for α t are widely unknown. Therefore the importance of this parameter, α t , has been tested in the range between 0.1 (near adiabatic boundary condition) and 1 (isothermal or heat bath boundary condition). Figure 9 shows the evolution in time of bubble radius, average temperature, wall velocity and average pressure for α t = 0.1, α t = 0.3 and α t = 1 for the near collapse and rebound phase. For the same values, figure 10 shows temperature and pressure at the bubble centre and near the bubble wall, allowing insight into the dynamics of the inhomogeneous interior of the bubble.
For higher values of α t , the average temperatures ( figure 9(b)) decrease, while the conditions in the bubble become more inhomogeneous (figure 10(a) (centre) in comparison with figure 10(c) (wall)). The wall velocity (figure 9(c)) is considerably higher in the last ns of the collapse, as the particles near the wall move with wall velocity and cause a comparatively low wall pressure. Due to the high collapse speed at α t = 1, a small minimum radius of R min = 0.64 µm is reached, compared to R min = 0.75 and R min = 0.93 µm for α t = 0.3 and α t = 0.1, respectively. At the time of minimum radius (t = 0), a huge pressure build-up to values of ≈1.5 GPa averaged can be seen ( figure 9(d)). The cases with lower values of α t show less extreme variations in pressure, due to the less violent collapse. The radius curve (figure 9(a))  shows that for α t = 1, as soon as the outgoing pressure wave hits the bubble wall it is thrown back. A high negative wall velocity and a 'wavy' radius curve occur. When using α t = 0.3 this effect is less pronounced; for α t = 0.1 it is completely absent. For α t = 1 (heat bath boundary condition), the particles transfer all their kinetic energy to the bubble wall, which leads to the formation of a high-density shell at the wall during collapse ( figure 10(d)). In the final stage of the collapse, many of these particles are driven to the bubble centre ( figure 10(b)), where they focus and cause a steep rise in density and temperature. Almost 25 000 K is reached in the bubble centre; a huge temperature gradient is present ( figure 10(a)). The situation for α t = 1 is depicted in a space-time plot for the density and the temperature in the bubble in figure 11, making the inhomogeneities directly visible. A comparison with figure 2 (α t = 0.3) shows the more pronounced spatial inhomogeneity developing for α t = 1.
Contrariwise, α t = 0.1 (near adiabatic boundary condition) leads to a much more uniform density and temperature distribution, even in the final collapse phase. While the average temperatures are higher (as almost no energy is transferred out of the bubble medium), the peak temperatures are lower than in the α t = 1 case. About 21 000 K is obtained, with strong fluctuations, however. At α t = 0.3, an intermediate case between isothermal and adiabatic boundary conditions, the heat-up of the centre is much slower compared to the case of α t = 1, but much faster than in the case of α t = 0.1. About 20 100 K is obtained. A value of more than 21 000 K would have been expected, but this value is within the fluctuations that still prevail even with 10 6 initial particles. The temperature gradient is less pronounced.
For selecting an appropriate value of α t for the calculations to come, the centre temperatures are not a good choice as they are only known approximately from measurements [4].
The light emission pulse widths, however, have been determined to good accuracy experimentally [75][76][77][78]. For typical sonoluminescence bubbles of the kind investigated here the light pulse widths cover a range of 100-250 ps. When comparing the calculated pulse widths in the present case (1010 ps for α t = 0; 320 ps for α t = 0.1; 125 ps for α t = 0.3; 45 ps for α t = 1), the best fit is α t = 0.3. Moreover, the collapse velocity has been measured by Mie scattering [79]. The maximum collapse velocity found experimentally ranges between 1200 and 1600 m s −1 . Again, α t = 0.3 gives the best fit with its maximum collapse velocity of 1530 m s −1 , the other values being 1090 m s −1 for α t = 0.1 and 1880 m s −1 for α t = 1.
For definiteness and for giving the most probable results, the calculations presented in the following are mostly done taking the intermediate value of α t = 0.3. From the variations in α t presented, conclusions can be drawn cum grano salis for the extreme cases.

Light emission
An example of light emission in a space-time plot for helium, argon and xenon is given in figure 12, together with the temperature development in the bubbles. In all three cases, the temperature forms a more or less confined 'hot spot' at the centre of the bubble with strong temperature increase and an increase in extent from helium to xenon. The same is found for the light emission calculated according to (12). For purposes of presentation, different scales are used for the power of the light emission due to the huge differences.
In order to analyse the possibilities of enhancing the collapse, as, for instance, seen in the increase in light emission, a series of simulations has been carried out with variation of the driving pressure, p a . For a given driving pressure, only certain equilibrium radii, R e , result in diffusive equilibrium, allowing a stable driving of the bubble. These values are also highly sensitive to the gas pressure in the liquid far away from the bubble [80]. An implementation of the diffusive equilibrium model discussed in [81] was used to derive R e for various p a (see also [82,83]). A gas concentration far away from the bubble, c ∞ , being a fraction of the saturation gas concentration, c sat , has been used: c ∞ /c sat = 0.001. The parameters and results for six cases are shown in table 3. Figure 13 presents the light emission in its time evolution and spatial extent for four different sound pressure amplitudes ( p a = 1.33, 1.35, 1.4 and 1.45 bar).
It is evident from table 3 that an increase in temperatures is not possible by simply increasing p a . Rather the opposite effect can be seen. As p a is increased, the growth of R min (resulting from the larger number of noble gas atoms because of the increase in the equilibrium radius, R e ) is not fully compensated for by an appropriate growth in R max , leading to a smaller compression ratio, C r = R max /R min . In order to achieve a comparable C r the driving frequency, ν a , would have to be decreased to give the bubble enough time to grow. Furthermore, much more water vapour evaporates into the bubble for larger p a (because of the larger R e ), leading to a greatly increased heat capacity. The trapping mechanism scales linearly with increasing p a , as the percentage of vapour in the bubble at the end of the collapse does not change appreciably. The additional water vapour inhibits a rise in temperature. As C r is lowered further, an increase in p a actually leads to reduced temperatures in the bubble centre (case 1 ( p a = 1.33 bar): 26 745 K; case 6 ( p a = 1.45 bar): 21 996 K).
The total light emission, P SL , shows a different behaviour. As expected, the peak light emission at the bubble centre is much higher for case 1, compared to cases 3, 5 and 6 (see figures 13(a) and (b)), but, as a result of the small centre volume, this is of little relevance for P SL . The outer shells of the light-emitting region have a much greater impact, as can be seen in The increase in p a leads to an enlargement of the emitting region; the outer shells of this region contribute the main part to P SL , due to their larger volume. At the same time, the pulse width t SL increases (figure 13(c)): for case 1, an FWHM of 70 ps is observed. This value grows continuously with rising p a ; 170 ps are seen in case 6. The pulse widths and trends correspond well to the experimental results in [75][76][77][78][79].
Overall, the enlargement of the emitting region and pulse width leads to a fivefold increase of P SL in the parameter domain studied (case 1: 0.096 pJ; case 6: 0.478 pJ). Sonochemical yields are scaled up comparably, as indicated by the development of OH production. This finding is supported by [84], where an increase of OH yield in dependence on acoustic power in multi-bubble sonoluminescence was found experimentally. Toegel et al [85] tried to enhance sonoluminescence by reducing the driving frequency to produce larger bubbles and also found Table 3. Acoustic pressure ( p a ) variations, associated equilibrium radius R e and simulation results for six cases including vapour chemistry. 'C r ' is the compression ratio (R max /R min ), 'OH' denotes the maximum OH yield in real particle number and 'H 2 O' is the real particle number of water vapour trapped in the bubble just before dissociation occurs. 'H 2 O (%)' gives the molar percentage of water vapour at this point in time. ' t SL' is the FWHM (full-width at half-maximum) of the light emission and 'P SL ' the total emitted light energy.  the additional water vapour to hinder stronger collapses and higher temperatures. Indeed, stronger collapses and higher temperatures must be sought at higher frequencies and, moreover, at larger static pressures [17,82,83]. Table 4. Species segregation including vapour chemistry. Simulation results for different noble gases (NG) and noble gas mixtures (1 : 1) of helium (He), argon (Ar) and xenon (Xe). T max is the absolute maximum temperature in the bubble; T av,max is the maximum average temperature. R min denotes the minimum bubble radius, v W,max is the maximum wall velocity, p W,c is the wall pressure during collapse at the inside of the bubble at a time 200 ps before the minimum radius and p W,max is the maximum wall pressure (typically 30-50 ps after the time of the minimum radius). 'O 2 ' denotes the O 2 yield in real particles and 'H 2 O' the real particle number of water vapour trapped in the bubble just before dissociation occurs. R e = 4.5 µm, ν a = 26.5 kHz, p a = 1.3 bar,

Species segregation
The homogeneity or inhomogeneity of the bubble interior is of importance to thermodynamic properties, as the speed of sound is dependent on the density of the medium. A spatial change in sound speed would influence the build-up of compression waves in the bubble. While for most of the time the mixture composition can be considered uniform, this assumption becomes invalid for the final stages of the collapse, where pronounced pressure and temperature gradients occur. Thermal diffusion (the Soret effect; lighter particles are driven to high temperature regions) as well as pressure diffusion (lighter particles are driven to low pressure regions) counteract conventional diffusion, leading to a segregation of species [25].
In order to examine the effects of mixture segregation in a luminescing bubble, computations with different noble gases (helium, argon and xenon) as well as mixtures of these gases were conducted. Table 4 gives numerical results for all six cases examined. Figure 14 shows the results for bubbles containing one noble gas only (besides water vapour and its decomposition and reaction products), figure 15 gives a space-time plot for the densities of a helium/water-vapour bubble and a xenon/water-vapour bubble and figure 16 covers bubbles with a mixture of two noble gases (molar ratio 1:1) and water vapour. Figures 14(a) and (c) illustrate the particle distribution within the bubble at the time of minimum radius for bubbles with a single noble gas (helium, argon or xenon). A distinct segregation of the bubble content, dependent on molecular weight, is observed. For argon and xenon (both of which are heavier than water vapour) a build-up of water vapour molecules in the bubble centre can be seen. The effect is more pronounced with the heavier xenon. When helium is used, water vapour remains mainly in the outer shells of the bubble, as the very light helium accumulates in the centre. The profiles show a roughly linear increase or decrease within the bubble (outside the very centre) at this moment. The maximum average temperature decreases from ≈14 300 K for xenon via ≈12 300 K for argon to ≈10 100 K for helium. These results are well within expectations, considering that the speed of sound decreases from helium via argon to xenon. However, another effect is to be considered: the heavier the particles, the more energy can be transferred from the bubble wall motion to the inside of the bubble (see table 4 and compare [13,31]). The wall velocity v W shrinks as molecular weight is increased, because the heavier particles induce a higher pressure at the bubble wall (example given at t = −200 ns). Only at the very last moment of the collapse, when density approaches its highest values, does the wall pressure p W become larger for lighter molecules due to the higher wall velocity (see p W,max ). The average kinetic energy of a xenon atom is about twice as high as for a helium atom in the last stages of the collapse. When the particle motion is thermalized as the particles focus in the bubble centre, this additional amount of energy leads to much higher central temperatures. The maximum temperature of a xenon bubble is ≈24 000 K, while a helium bubble does not reach 14 000 K (compare the space-time plot in figure 12). The temperature gradient is much steeper for the xenon bubble. It is noteworthy that the more pronounced species segregation in the xenon bubble, compared to the argon bubble, does not prevent the much higher temperatures, indicating that species segregation plays a minor role in the hydrodynamic evolution of the bubble. It might be more important in bubble chemistry, though.
It can be seen that the chemical activity is significantly higher in the xenon bubble. The O 2 yield, for instance, is more than double the O 2 yield of a helium bubble (see table 4). This can be partly attributed to the higher temperatures, but is also due to the particle distribution. As stated above, the light helium tends to accumulate in the (hot) bubble centre, restricting water vapour to the (colder) outer regions. The reduction in vapour concentration in the areas of highest temperatures leads to further reduced dissociation. In the xenon and argon cases, the opposite effect leads to a high water vapour concentration in the hot bubble centre. Figures 14(b) and (d) show the differences in temperature and water vapour molar fraction between the bubble centre and the bubble wall during collapse and rebound. It is evident that species segregation starts long before the temperature gradient becomes relevant. In particular in the xenon case, segregation is already pronounced 45 ns before R min is reached, while a notable temperature gradient develops just ≈2 ns before this moment. The pressure gradient (not shown) shows a similar timescale. Segregation is strongest at the time of minimum radius and quickly relaxes, because the outgoing pressure wave drives the lighter particles from the centre back to the wall. During expansion, the gradient becomes slightly negative. By looking at these results it can be concluded that species segregation is not mainly induced by the extreme gradients reached at the time of minimum radius, but rather develops during the whole collapse phase. As a concentration gradient is already present, the inward-travelling density wave intensifies the effect only by a small amount. The reflected pressure wave weakens the concentration gradient and leads to a fast relaxation of the segregation. The pressure diffusion seems to dominate the thermal diffusion at this point. These findings are in agreement with the results of Storey and Szeri [25], where a bubble with a He-Ar mixture (no water vapour) was simulated using a Navier-Stokes model. While the collapse phase shows quite similar results, the relaxation of segregation takes much longer in [25] (≈50 ns). A slight reversal of segregation during bubble growth is present in both models.
The results of calculations with two noble gases (cases IV-VI in table 4 and figure 16) generally confirm the considerations for cases I-III: particles segregate according to molecular weight-the greater the difference in molecular weight, the more pronounced the segregation. Temperatures achieved and chemical activity fall in between the results seen with only one of the two noble gases. A linear dependence of molar fraction from the centre to the wall cannot be seen in the computations with three components. This means that a third component has a nonlinear influence on the results. For all cases, a very comparable number of H 2 O molecules remain in the bubble. This shows that the trapping mechanism is influenced only marginally by species segregation. When a relevant difference in water vapour at the bubble wall is reached, most of the vapour molecules have already left the bubble.

Chemical evolution
Chemical reactions occur in the final stage of bubble collapse, when temperatures are high enough to dissociate the water vapour trapped. To illustrate the influences of chemical reactions on bubble temperature and other quantities, a comparison with calculations neglecting chemical reactions is presented in figures 17(a) and (c). Two versions are considered when not regarding chemical reactions: one simulation incorporates the excitation of rotational degrees of freedom for the water molecules (curve 3) and one neglects this effect (curve 1). Curve 2 incorporates chemical reactions, but no rotations, of the water vapour according to table 1. The reason for not including rotation is that it is not known what consequences a chemical reaction of a rotating molecule has for the rotational energy of the product(s). Therefore, the rotation model was not Table 5. Comparison of five different collapse scenarios defined by considering rotation or vapour chemistry with different α v,OH , the mass accommodation coefficient for the OH molecule. T av,max is the maximum average temperature, R min denotes the minimum bubble radius, 'OH cond ' is the number of real particles condensing at the wall during collapse and for about 30 µs afterwards, 'OH prod ' denotes the OH yield in real particles in the first collapse phase. R e = 4.5 µm, ν a = 26.5 kHz, p a = 1.3 bar,   The bubble collapses more strongly, i.e. to a smaller R min , when water molecule rotation or vapour chemistry is considered. This is as expected, because in both cases energy is consumed, leading to slower pressure build-up inside the bubble to stop the collapse. Despite the stronger collapse, the temperatures (both maximum and average) are lower for the simulations with rotation as well as for the simulation with vapour chemistry. This again is as expected, as in the case of rotation kinetic energy is lost for the particles upon collision and in the case of vapour chemistry the number of particles sharing the energy is increased by dissociation of the water vapour molecules.
In figures 17(b) and (d), the OH production and decay is given, as the OH radical is of interest for sonochemistry (see, e.g., [8]). Figure 17(b) gives the total OH content in real particle numbers in the bubble for three different values of the mass accommodation coefficient α v,OH : 0, 0.001 and 0.01. No large differences are to be seen for the total OH content with respect to the different values of the OH accommodation coefficient. The condensation of the OH molecules at the bubble wall is given in figure 17(d), again in real particle numbers (note the different scale as compared to figure 17(b)). The condensation has been followed for about 30 µs beyond the collapse. The number of condensed OH molecules during this time span is given in table 5. It strongly depends on the OH mass accommodation coefficient, α v,OH . Because this coefficient is not known, some estimates are given as a prediction. As for the trend, a factor of ten increase from α v,OH = 0.001 to 0.01 leads to about a factor of ten increase in condensed particles. This factor has also been found in the continuum model calculations of [9], albeit the absolute numbers there are lower at p a = 1.2 bar driving.
In figure 18, more details of the vapour chemistry inside the bubble are presented. Figure 18(a) shows the (real) number of particles for the H 2 O molecules (the fast decay near the final collapse is due to dissociation because of the high temperatures obtained) and the production and decay of the species OH, O 2 and H 2 . On a much smaller scale, the production and decay of the further species H, H 2 O 2 , HO 2 and O take place ( figure 18(b)). In figure 18(c), the temperatures inside the bubble are given for the centre and at the wall, together with the average temperature throughout the bubble. It is seen that the maximum temperature reaches about 18 900 K, whereas the maximum average temperature is about 11 500 K (see table 5) and the maximum wall temperature is only about 7 500 K.
The density of the OH species is given in figure 18(d) in its time dependence in the centre and at the wall of the bubble. Again, the distribution is not uniform. A strong concentration in the centre is observed. The decay is fast; the molecules are almost gone after 1 ns. Compared to the continuum model [9], the spatial OH distribution is much more homogeneous here. This outcome may be attributed to the faster diffusion caused by the 10 6 ensemble particles only.
For better visualization, figure 19 shows the space-time evolution of the temperature ( figure 19(a)) and of the five densities for H 2 O, OH, H, O 2 and H 2 . It is seen that the dissociation products OH and H appear first according to the temperature rise; the formation of the reaction products O 2 and H 2 , as expected, sets in afterwards.

Conclusions
Cavitation bubble collapse is a violent process difficult to approach experimentally as well as theoretically or numerically. A new, self-consistent, hybrid approach has been investigated, the RP-MD model, where the bubble interior is modelled by hard spheres in MD simulations and the surrounding liquid by the usual continuum RP equation, bubble interior and liquid coupled by the MD pressure at the wall as input to the RP equations. A typical sonoluminescence bubble in water at low driving frequency (26.5 kHz) has been investigated and compared with continuum calculations and experiments. The investigations explore a new pathway to describe strong bubble collapse and the associated phenomena, including light emission. The point of departure is the molecules filling the bubble. As a first approach, they are taken as hard spheres, a good choice, when long-range attractive forces are negligible. This is the case at high velocities and high temperatures as encountered in strong bubble collapse. At present, the number of particles that can be handled easily in a computer is limited to about 10 7 . This makes for small bubbles only, below about a micrometer. A partial remedy is via ensemble particles as used here, only partial because not all physical properties can be scaled with the same factor simultaneously (see section 2). The mean free path is larger; thus particles diffuse faster, but the relation ordering is retained. Faster diffusion smoothes steep gradients, as is indeed observed with increasing the number of real particles combined to one ensemble particle (see section 3.1). Thus, the inhomogeneity inside the collapsing bubble is (slightly) underestimated in the present investigation (see sections 3.4 on light emission, section 3.5 on species segregation and section 3.6 on chemical evolution).
A comparison has been made with continuum model results [9]. Close agreement is found for water trapping (see table 2) and mixture segregation trends. Concerning vapour chemistry, the at times small number of product particles may give only coarse, fluctuating results (see, for instance, figure 19(d), the coarse distribution of H particles in the bubble), not yet well suited for detailed comparison. The general trends, however, are confirmed by the MD calculations presented.
Unfortunately, the RP-MD approach does not remove the problem of mass and energy transfer at the bubble wall, where MD and continuum dynamics meet. At this border, the continuum is in need of diffusion constants exemplified by appropriate mass and thermal accommodation coefficients, α v and α t . These coefficients can either be determined by MD calculations on both sides of the bubble wall or by experiments. Both ways are difficult at present and are unable to supply the necessary data in dependence on temperature and pressure. Relational orders, however, seem not to be affected, i.e. trends are possible to predict. This includes the inhomogeneity in temperature and the associated inhomogeneity in chemical reaction product distribution. The transition of particles into the liquid, however, cannot be properly predicted without these coefficients, because the particles then leave the range of their MD existence. Thus the MD approach to the interior improves the description of interior dynamics, but still suffers from the continuum boundary that also delegates its unknowns into the interior via the boundary condition.
As a solution, a three-region model is suggested. A thin shell of MD fluid (with 10 6 or 10 7 particles, say) may be inserted at the wall, shifting the gas phase-liquid boundary into the liquid. The MD liquid-continuum liquid boundary should be easier to handle and produce smaller errors than the present boundary, as no evaporation and condensation have to be considered. When the liquid layer is large enough, also the heat flux may not need special consideration. A new problem may arise because of the shock waves radiated that are now included in the RP model. For a liquid, the long-range attractive forces play a role and a continuous potential LJ fluid as in [43,44,46] may be used. This could solve the deadlock concerning the diffusion properties of collapsing bubbles.