Modeling quantum processes in classical molecular dynamics simulations of dense plasmas

We present a method for treating quantum processes in a classical molecular dynamics (MD) simulation. The computational approach, called ‘Small Ball’ (SB), was originally introduced to model emission and absorption of free–free radiation. Here, we extend this approach to handle ionization/recombination reactions as well as nuclear fusion events. This method exploits the short-range nature of screened-particle interactions in a dense plasma to restrict consideration of quantum processes to a small region about a given ion, and carefully accounts for the effects of the plasma environment on two-particle interaction rates within that region. The use of a reduced set of atomic rates, corresponding to the bottleneck approximation, simplifies their implementation within an MD code. We validate the extended MD code against a collisional–radiative code for model systems under two scenarios: (i) solid-density carbon at conditions encountered in recent experiments, and (ii) high-density Xe-doped hydrogen relevant for laser fusion. We find good agreement for the time-dependent ionization evolution for both systems. We also simulate fast protons stopping in warm, dense carbon plasmas. Here, reasonable agreement with recent experimental data requires contributions from both bound electrons, as modeled by SB in the extended MD code, and free electrons; for the latter, use of the classical random phase approximation (RPA) formula instead of the MD prediction yields better agreement with the experiment, a result that can be attributed to the use of modified Coulomb potentials in MD simulations of electron–ion plasmas. Finally, we confirm that the fusion reaction rate obtained from an MD simulation agrees with analytical expressions for the reaction rate in a weakly screened plasma.

use of the classical random phase approximation (RPA) formula instead of the MD prediction yields better agreement with the experiment, a result that can be attributed to the use of modified Coulomb potentials in MD simulations of electron-ion plasmas. Finally, we confirm that the fusion reaction rate obtained from an MD simulation agrees with analytical expressions for the reaction rate in a weakly screened plasma.

Introduction
There now exists a considerable body of molecular dynamics (MD) research on dense, classical plasmas. Following seminal work by Hansen and colleagues [1][2][3], three decades of progress in computational hardware and software have enabled plasma MD to evolve from simple onecomponent models with ions moving in a uniform charge-neutralizing background, to fully dynamic mixtures of electrons and multiple ion species. For example, at the heart of the Cimarron Project (based at Lawrence Livermore National Laboratory) is a new, massively parallel MD code, 'ddcMD', designed to simulate high energy density plasmas undergoing thermonuclear burn [4][5][6][7]. With this capability one can now investigate quantitatively, and over a wide range of plasma conditions, basic properties such as dynamic structure factors and electrical and heat conductivities; time-dependent phenomena such as rates of temperature relaxation among electrons, ions and radiation; and stopping of fast ions.
It is well known that an isolated, classical system of ions and electrons interacting through pure Coulomb forces promptly collapses: with nothing to limit them, potential energies become increasingly negative, causing a corresponding, unphysical rise in particle kinetic energies. This run of events is similar to that experienced by an isolated N-body astrophysical system undergoing gravo-thermal collapse [8]. For the most part, past MD simulations have avoided this problem by construction (e.g. the use of modified Coulomb potentials), by choice of conditions (e.g. high temperatures) and/or by clever artifice (e.g. positrons and protons in a charge-neutralizing background). However, unlike stars and galaxies, ordinary plasmas are subject to the stabilizing effects of quantum mechanics. In particular, the existence of discrete bound states and stable ground states dramatically affects the dynamics of atomic systems.
There are two important reasons to introduce bound atomic states, and the quantum mechanical processes leading to their formation and destruction, into classical MD simulations of plasmas. The first is to avoid the non-physical Coulomb collapse experienced in purely classical MD. The second reason is to extend the capability of MD for studying dense, nonequilibrium plasmas of wide interest, including those generated by direct target interaction with intense beams of particles or radiation. In particular, MD simulations would gain an ability, similar to that of atomic kinetics codes, to make detailed connections with spectroscopic observations [9][10][11][12]. Our third motivation for this work is the desire to extend MD capability to encompass fusion reactions in burning plasmas.
In this paper we describe a method for dealing with quantized internal states of plasma ions in a classical MD simulation. This method is called 'Small Ball' (SB) and it was used initially to model emission and absorption of free-free radiation [13]. Here we extend it to include ionization and recombination processes that can occur when a projectile interacts with a target ion. Collision parameters are determined when the projectile enters a small spherical volume surrounding the target, its 'SB.' Probabilities for the various atomic processes are then constructed from accurate quantum cross sections that satisfy microscopic reversibility and, where needed, contain explicit three-body energetics. Additionally, there is an analogous nuclear version that allows one to treat nuclear fusion events.
The short-range nature of screened particle interactions in dense plasmas is key to the SB algorithm. Section 2 introduces the effective interaction potentials used in our MD simulations, and discusses plasma effects important for SB. The statistical physics underlying these ideas is reviewed in appendix A. Elaboration of the model is given in section 3. This is followed in section 4 by a description of the simplified atomic kinetics scheme we use to track charge states, ion by ion. The cross sections adopted for various processes are described in appendix B. After that, in sections 5-8, we describe specific implementation issues and test studies conducted to establish SB's credibility. Section 9 summarizes our main findings.

Interactions in a dense, equilibrium plasma
Important aspects of SB are based on the statistical mechanics of equilibrium systems. In particular, the use of quantum statistical potentials makes possible explicit simulation of electron and ion dynamics; the use of potentials of mean force makes it straightforward to allow for the average influence of neighboring charges on the dynamics of a colliding pair; and the adoption of equilibrium screening models makes possible a simple, approximate treatment of plasma effects on reaction cross sections. Since, in some cases, we wish to apply the SB model to nonequilibrium systems, the use of these equilibrium concepts requires some justification. First, we note that the equilibrium statistical quantities we use generally result in only modest corrections to collision energies and interaction potentials. Also, even in dense plasmas interacting with intense fluxes of particles or radiation, at any given time the bulk of the matter is in a nearequilibrium state with a reasonably well-defined temperature. To these physics points can be added a practical justification-similar (in some cases, identical) equilibrium statistical quantities are employed in atomic kinetics codes like the one we use in section 6 to establish the fidelity of SB. Since these atomic kinetics codes have a well-established history of reliable plasma modeling, we are mimicking that approach in order to validate the MD code in regimes in which both methods are reliable. On the other hand, Cimarron studies of fast-ion stopping in fully ionized plasmas have revealed that use of statistical potentials leads to an underestimate of energy transfer in collisions of these (non-equilibrium) ions with (equilibrium) plasma electrons [14]. We return to this point in section 7, which deals with SB calculations of stopping in partially ionized plasmas.

Quantum statistical pair potentials
Since the aforementioned simulations by Hansen and co-workers [1][2][3], it has been common to use Coulomb potentials modified at small inter-particle distances to avert rapid collapse of classical electron-ion plasmas. The physical basis of such modification is quantum mechanical, and arises when point particles are replaced by de Broglie waves whose scattering exhibits diffractive behavior. As reviewed by Jones and Murillo [15], if one approximates the potential energy in the quantum partition function for an N-body system as a sum of two-body terms-analogous to the familiar, classical approximation of pair-wise interactions-one can obtain an effective, temperature-dependent potential whose use yields correct equilibrium statistical properties of non-degenerate Fermi gases. Dunn and Broyles [16] first derived the form of the quantum statistical potential used in Hansen's work. For charges Z a and Z b in a gas at temperature T, diffractive effects lead to the modified interaction where ab = 2πh 2 /µ ab T is the thermal de Broglie wavelength of a particle whose mass µ ab is the reduced mass of the pair (a, b). Because of the large masses and the like signs of the charges, ion-ion interactions seldom experience quantum modifications. However, in dense but non-degenerate plasmas with temperatures T > few 10's of eV, the diffractive smearing represented by equation (1) is important for electron-electron and electron-ion interactions at short range. Pauli exclusion also modifies the Coulomb interaction between electron pairs (with like spin), and it too can be mimicked by a statistical potential such as the spin-averaged result of Deutsch et al [17], There exist a variety of published statistical potentials with similar behavior, and different ones have been adopted for particular simulations carried out as part of the Cimarron Project, but only equations (1) and (2) have been employed in the computations reported here. Whatever the choices, it should be noted that the efficacy of such potentials for studying dynamic properties of charged, possibly non-equilibrium, many-body systems is still an important and open subject of investigation as discussed in [5]. The SB model that we implement is largely independent of the choice of the statistical potentials because, as described in section 3, we minimize their effects on reaction probabilities for a colliding pair by construction.

Pair distributions and potentials of mean force
Determining the plasma's influence on atomic and nuclear processes requires accurate knowledge of screening on length scales of the order of the mean inter-particle spacing, and less. Particle distributions in an MD simulation contain this information, but as r-values decrease the statistics get poorer. For SB, we chose to obtain the requisite plasma screening properties from theoretical models discussed in appendix A and in section 2.3, below. Here, we introduce the key quantities.
In the canonical ensemble formalism, the configuration integral for an N-body system, in thermal equilibrium at temperature T = 1/β and contained in volume , is Here, U N is the system's total interaction energy, and { s N } denotes the coordinates of all the particles. If one assumes that U N is the sum of all pair-wise interactions U ab , it is straightforward to express formally the distribution of the other N − 1 particles with respect to one particle. Suppose the plasma has multiple species, with particle numbers N a + N b + · · · = N , positions { s ai ; s bj ; . . .} = { s N }, and mean densities ρ a = N a / , etc. Then, the local density at r 2 of species 'b' surrounding a particle of species 'a' at r 1 can be written in terms of another integral involving the coordinates of all particles (see [18]), This equation serves to define the two-particle distribution function g ba ( r 2 , r 1 ). By inspection it is evident that such distribution functions are symmetrical in their indices. Also, for a homogeneous system of particles interacting via central potentials, this local density at r 2 depends only on the distance r = | r 1 − r 2 |, so g ab ( r 1 , r 2 ) = g ab (r ) = g ba (r ). From the barometric equation for the density of particles under the influence of an external field of force, viz., n(r ) = n 0 exp[−βU (r )], we see that the distribution function g ab (r ) can be used to define an effective energy of interaction for the pair (a, b), in the presence of all the other particles and thermally averaged with respect to their possible configurations, This energy, called the potential of mean force (MF), is usually written as the pair's direct interaction minus a quantity H ab (r )that represents the average screening of this interaction by other charges, viz., Limiting values, H ab (0), which we refer to as plasma screening constants, are of particular importance to our SB model (cf section 3.3).

Hypernetted chain examples
Analytical screening models, based on limiting high-or low-temperature approximations, are often used to estimate dense plasma effects and to provide approximate screening constants. Their development also serves to highlight the underlying statistical physics concepts. An overview is presented in appendix A. However, to obtain accurate screening information over a wide range of plasma conditions, one needs more sophisticated methods to compute pair distributions. Here, we use the hypernetted chain (HNC) method [18], originally applied to the structure of liquids, to illustrate general points needed below. Our HNC calculations employ the modified Coulomb interactions, equations (1) and (2). In the weak coupling (i.e. low density, high temperature) limit, the HNC potential of mean force tends to the familiar Debye-Hückel (DH) expression where is the Debye length involving all plasma species 'i.' This model is appropriate only when there is a large number of particles N D within the Debye sphere; in other cases, the HNC and the Debye results will differ significantly. Subsequent discussion utilizes several other, interrelated radii. The first of these, which involves the mean density ρ ion of all plasma ions, is the typical ion-ion separation R ion = (3/4πρ ion ) 1/3 . The second is R e , similarly defined in terms of the mean electron density ρ e = Z ρ ion , where Z is the plasma's average ion charge. Third, for each ion of charge Z there is an ion sphere radius defined such that the net charge of free electrons within the ion sphere is -Z e: R z = (3Z /4πρ e ) 1/3 (cf equation (A.9) and [19]). The relation among these radii is Our two HNC examples, below, involve conditions representative of the plasmas studied in sections 6-8. These HNC calculations are used again in appendix A to assess different analytical approximations for electron-ion screening constants H eZ (0) . Because HNC theory cannot treat bound electrons explicitly, we approximate the influence of N bnd such electrons on a free-electron-ion interaction via a sum of two quantum statistical terms, U eZ (r ) = [U D eZ (r ) + N bnd U P ee (r )]. (Effectively, bound electrons are 'pasted' onto their parent nucleus.) The first plasma is warm, partially ionized, solid-density carbon, ρ ion = 10 23 cm −3 , at temperatures T = 25 or 50 eV. At the lower temperature, the plasma is taken to be 2/3 C 1+ and 1/3 C 2+ ; at the higher temperature, the ion fractions are reversed. (The values of N D for these plasmas are 0.56 and 0.99, respectively.) The first panel of figure 1(a) shows the radial dependence of all six pair distributions for T = 25 eV. As expected, at small r-values the higher ion charge gives rise to higher electron density and lower ion density. The sharp drop in all three ion-ion pair distributions indicates that every ion strongly excludes other ions from its own ion sphere. The decrease in g ee at small r-values is a consequence of the Pauli repulsion, equation (2). The important point to recognize is that, because the pair distributions all tend to unity beyond R ion , the potentials of mean force-the effective, screened interactions between particle pairs-tend to vanish rapidly beyond that distance. This follows directly from equation (5).   shows just the electron distributions g eZ surrounding C 1+ and C 2+ , but for both temperatures; here, each pair distribution is plotted as a function of r/R Z , where R Z is that particular charge state's ion sphere radius. Note that the free electron density is nearly equal to its mean value ρ e in the outer half of each ion sphere's volume. All four curves reveal essentially complete screening beyond r = R Z .
Our second HNC example is a hot, T = 5 keV, compressed DT plasma with a small admixture of highly stripped xenon ions; the total ion density is ρ ion = 10 26 cm −3 , and the   figure 2(a) shows the radial dependence of all the pair distributions g ab (r/R ion ) for the case when the specified Xe charge is 44+. As expected, at small r-values this very large Xe charge gives rise to electron (ion) densities that are much higher (lower) than those surrounding D and T ions. Every Xe ion strongly repels other ions but, even so, each Xe Z + ion sphere contains about Z/2 DT ions, mostly in the sphere's outer region. Figure 2(b) shows the electron distribution g eZ surrounding Xe 44+ and, to illustrate high-Z dependences, the electron distribution surrounding Xe 48+ , obtained from an otherwise identical HNC calculation; here again, each pair distribution is plotted as a function of r/R Z , where R Z is that particular charge state's ion sphere radius. Even in this plasma with high-and low-Z components, the electron density is nearly uniform in the outer portion of each Xe ion sphere, and, again, there is essentially complete screening beyond r = R Z . Figures 1(c) and 2(c) show the screening functions H eZ (r/R Z ) for C 1+ and C 2+ (figure 1(c)) and Xe 44+ and Xe 48+ (figure 2(c)) for the corresponding conditions. These will be discussed further in appendix A. The screening behavior seen in these two examples, which differs markedly from the Debye-Hückel prediction [equation (7)], is typical of dense plasmas and supports our SB picture of incident charges (electrons or ions) interacting with one target ion at a time.

Preliminary considerations
Our approach represents a probabilistic treatment of various inelastic quantum events that can occur in the close encounter of a pair of particles, and the SB construct is used to evaluate these probabilities for MD simulations. The SB itself is a sphere of fixed radius that is centered on each target ion. The size of the SB, R B (Z ) is constrained by considerations that pertain to the establishment of a collision's initial conditions, as described below. For a dense plasma it is of the order of the ion-sphere radius. The probability P γ of a particular inelastic reaction 'γ ' is computed once per encounter, at the time when the incident particle enters the target's SB. Elastic scattering corresponds to the probability that 'nothing happens', in which case the particles complete their encounter on trajectories set by the MD forces.
3.1.1. Electron-ion events. As our HNC results illustrate, in dense plasmas the electron-ion potentials of mean force U MF eZ (r ) are negligible for separations r > R Z , which implies that electron encounters with individual ion targets effectively begin at r = R Z . Therefore, we specify that R B R Z . We assume that a pair's 'asymptotic' collision energy ε 0 is established at R Z , and we also assume that, once inside R Z , a colliding pair has fixed total energy. Hence, the relative kinetic energy is ε(r ) = [ε 0 − U MF eZ (r )]. These simplifications, which limit the role of all other particles to just their establishment of the potential of mean force, enable us to determine ε 0 from an 'audit' of the collision at any point within the ion sphere. However, typically there are Z − 1 other free electrons, and perhaps some other ions, already inside the target's ion sphere when a new collision begins. As one decreases R B , multi-center scattering from more than one ion becomes less likely, but a close encounter with some spectator electron, which would affect our calculation of ε 0 via the potential of mean force, becomes more likely. Finally, one must recognize that each ion's R Z -value changes when it undergoes an ionization or recombination, so a single, fixed SB radius is inappropriate for a simulation with changing ionization.

Ion-ion events.
The HNC calculations also show that ions in dense plasmas tend to avoid each other's ion spheres, the behavior being more pronounced for larger charges. This suggests that multi-center scattering effects can be avoided in the collision of two ions, Z a and Z b , as long as the SB radius R ab B used for such an event is smaller than R ion . How much smaller can be set by considering the two kinds of ion-ion reactions presently of interest, thermonuclear fusion of light ions (ε 0 10 keV amu −1 ), and energy loss of fast, non-thermal ions (ε 0 10 MeV amu −1 ) via ionization of heavy elements. For these energies, the classical distance of closest approach, r 0 Z a Z b e 2 /ε 0 , is also much less than the mean inter-ion spacing R ion . To ensure that candidate events for these kinds of reactions are not missed, one needs ion-ion SB radii larger than r 0 . Taking the target to be the more massive ion, with charge Z a , values R ab B > 0.01Z a a 0 , where a 0 is the Bohr radius, satisfy this condition as long as the collision energy exceeds several keV. (Section 8 contains further discussion of these points.)

Reaction probabilities
We develop the formula for the reaction probability P γ in steps, beginning with a dilute, weakly interacting gas where ordinary kinetic theory is valid, and ending with a dense plasma. Reaction rate expressions used below involve Maxwellian particle fluxes and a target at rest, but the probability definitions we extract are general, and not limited to these special circumstances. The rationale is analogous to the use of detailed balance under thermal equilibrium conditions to obtain the unique relation between, say, Einstein's A and B coefficients.
In the first step, we consider an essentially ideal gas for which one can choose an SB radius R B larger than the range of interaction between the pair (a, b), but still smaller than the mean separation of targets. At this value of R B , incident particles have a density equal to their mean value, ρ b = N b / , and their flux has a Maxwellian distribution of energies ε 0 : f (ε 0 ) = (4ε 0 /π T 3 ) 1/2 exp(−ε 0 /T ). For a target at rest, these energies correspond to incident collision energies. For the target's SB, the inward radial flow of particles of type (b) per unit time having energies in (ε 0 , ε 0 + dε 0 ), is where w(ε 0 ) is the (relative) velocity. Equating the number of 'γ ' reactions per SB per unit time, P γ times dv b , to that predicted by kinetic theory, we obtain where σ γ is the ordinary reaction cross section. This yields which is just the result our intuitive notion of the cross section would suggest. Next, consider the situation in a dilute plasma where the mean distance between ions is less than the screening length in the Debye-Hückel interaction, viz R ion < D. The distribution of relative energies ε B at R B < R ion is still Maxwellian, but these are not asymptotic collision energies; rather, from the aforementioned assumption of energy conservation for the colliding pair we have ε B = ε 0 − U DH ab (R B ). Moreover, the density of particles (b) at R B is not the mean density ρ b ; rather, from the discussion in section 2.2, we have The analogous calculation of the number of reactions per sphere per unit time with energies in (ε B , ε B + dε B ) now leads to the result where the first line is written in terms of the relative energy ε B as the projectile enters the SB, and the second, in terms of the corresponding asymptotic collision energy ε 0 . In both versions, the bracketed factor corrects for the focusing (defocusing) of flux at R B due to an attractive (repulsive) interaction. (Note that the energy intervals are identical: dε 0 at ε 0 , and dε B = dε 0 at .) The general form of equation (14) also follows from a consideration of the invariance of trajectories in phase space [20].
As described earlier, we assume that the many-body potential U MF ab (r ) represents the average gain or loss of kinetic energy, relative to ε 0 , by a pair (a,b) when their separation in a dense plasma decreases from the ion sphere radius (where their mutual interaction is nil and, hence, their pair distribution function is essentially unity). This makes for a straightforward extension of our SB probability scheme to dense plasmas: in both equations (13) and (14) the potential U DH ab (R B ) is simply replaced by the corresponding potential of mean force, U MF ab (R B ). It follows that With this formula, MD with SB predicts the same average number of reactions per target per second as kinetic theory for the same plasma conditions and for the corresponding energy intervals, dε 0 at ε 0 , and . Hence, at this stage, SB has no plasma 'effects.' Section 3.3, below, describes how we introduce a simple version of plasma screening into these reaction probabilities.

Electron-ion events.
A value of ε 0 determined at r < R Z can be less than zero if the event's local environment has particle densities and, hence, many-body potentials far from their statistical averages. If these relatively rare events are ignored, or if ε 0 is somehow adjusted, SB no longer yields the same reaction rates as kinetic theory (though the differences are very small). Alternatively, one can avoid this issue altogether by choosing R B = R Z . Then, U MF ab (R B ) = 0 and the audited energy is ε 0 ; additionally, the dense plasma reaction probability reduces to the ideal gas form, equation (12). The price paid for this considerable simplification is a greater likelihood of undetermined multi-center effects involving other ions near the edge of the target's ion sphere. (14) and (15) exhibit a peculiarity of our classical model for repulsive, ion-ion collisions. Only collision energies ε 0 > U ab (R B ) > 0 will result in energies ε B > 0 at R B . This means that the contribution of lower-energy events to the total rate of ion-ion reactions 'γ ' will be excluded if U ab (R B ) exceeds that reaction's threshold, ε γ . Although true, this fact is unimportant for the following reasons. (i) For fusion reactions there is no actual threshold, just increasingly smaller reaction cross sections as ε 0 decreases. Consequently, the candidate events excluded from consideration are unimportant. (ii) For atomic reactions, it is well known that incident ions cannot transfer large fractions of their energy to bound electrons; indeed, detailed numerical studies of many different ion-impact reactions show that inelastic cross sections are very small unless ε 0 > 10 2 ε γ ( [21], and discussion in section 4.3). Since the ion-induced transitions of interest in our atomic kinetics studies are ground-state ionizations (see discussion in section 4.1), any excluded low-energy events contribute negligibly to ion-impact ionization rates.

Ion-ion events. Equations
Also, if accurate potentials of mean force for ion-ion pairs are unavailable, we believe it is best to keep the SB radius at R ab B = 0.01Z a a 0 , and to ignore the many-body correction (which is less important at the high ε 0 -values of interest); then, equation (15) reduces to the dilute plasma form, equation (14).

The role of plasma screening
As we noted earlier, the presence of U MF ab (R B ) in equation (15) relates to the environment's influence on the flux of incident particles at R B , and in particular it is not an indication of a dense plasma effect on the reaction's cross section. In SB such an effect involves the screening constants H ab (0), and what follows is an elaboration of our model to include screening, in an average way.
We begin with the well-studied process of thermonuclear fusion. With it, one can confirm the SB treatment of plasma screening without extra complications involving bound electrons. Cox [22] gives a detailed description of the standard kinetic theory treatment of fusion reactions, with and without screening. Because essentially all collisions in multi-keV plasmas are at energies ε 0 too small to surmount the peak of the Coulomb barrier at the contact distance R ab (the sum of the two nuclear radii), the dominant energy dependence of the fusion cross section involves the ions' probability of tunneling through that barrier, the so-called Gamow penetration factor. Only s-wave scattering, which adds no centrifugal term to the barrier, is important. This unscreened fusion cross section for a pair (a, b) has the well known form where the reaction-dependent quantity S γ is only a weak function of energy, r 0 is the classical distance of closest approach and the pair's effective wavenumber inside the barrier U ab (r ) = Z a Z b e 2 /r is This nuclear cross section involves the true Coulomb interaction, without modification. To find the dominant influence of plasma screening, one assumes that S γ is constant and that the effect of screening on r 0 is ignorable. Then, with the replacement it becomes evident that, with screening, the tunneling factor is computed at the effective energyε = ε 0 + H ab (0). (Here and below, a tilde denotes a screening-related quantity.) This leads directly to the result that the screened fusion cross sectionσ γ satisfies the simple scaling relation Hence, the plasma-screened version of equation (15) for fusion is In the numerator, the energy change from ε 0 relates to screening's enhancement of tunneling, while in the denominator it relates to the flux correction at R B due to the pair's repulsive interaction. This reaction probability leads to Salpeter's fusion rate enhancement factor of exp[β H ab (0)] when the usual simplifications are made in evaluation of the integral of equation (11) [19,22]. We now turn to the effects of plasma screening on cross sections for atomic ionization and recombination processes. The relevant transition matrix elements involve bound and free electron wave functions, whose form in a dense plasma environment is not well understood (see, e.g., [23,24]). We therefore limit our treatment of screening to the simple model of continuum lowering (CL), which involves only the electron-ion screening constants H eZ (0) [25,26]. Calculations of these screening quantities involve free electrons and point-charge ions, leaving the influence of any bound electrons unaccounted for.
Density functional computations of electron distributions in dense, partially ionized matter show that bound electrons tend to be localized within the inner portion of the ion sphere, say r < R z /2 [27], so only free electrons near the origin sense the partially screened nuclear charge of the central ion. Figures 1(b) and 2(b), which illustrate the sensitivity of g eZ to different point Z-values, indicate the magnitude of changes in the pair distribution-at small r-values-likely to be caused by bound electron screening. However, most free electrons (and ions) are well outside the core of bound electrons, and they sense only the central ion's net charge. Using figures 1(b) and 2(b) again, we see that electron radial distributions become relatively less sensitive to Z as r increases toward the ion sphere radius. Further, the derivations in appendix A make it clear that screening constants H ab (0) are determined mainly by the distribution of numerous but more distant charges. We therefore expect bound electron effects on screening constants to be small, and ignorable for present purposes.
To arrive at the CL picture of plasma screening, one assumes that the effective interaction between an active, bound electron in an ion of net charge Z − 1 and the ionic core (of net charge Z) is approximately U eZ (r ) − H eZ (0). Then, by inspection of the Schrödinger equation, it follows that if an orbital (n ) has unscreened energy E (Z −1) n , its perturbed, screened energy isẼ (Z −1) . Alternatively, one can state that screening has reduced its binding energy I (Z −1) n by the amount I (Z −1) = |H eZ (0)|: (21) Because bound orbitals are unchanged by the addition of the constant term, as are energy differences between bound levels, cross sections for bound-bound processes also are unchanged. This means that, in the CL scheme, SB probabilities for bound-bound transitions (collisional or radiative) would be unaltered by screening. The situation is different for ionizations and recombinations. For those processes, an electron-ion screening constant has the effect of lowering reaction thresholds, but otherwise leaving cross section unchanged, viz., plasma screened ionization and recombination:σ γ (ε 0 ; We use this same formula to model screening effects on the cross section for ionization by fast ions, ε 0 few MeV amu −1 , because the theory we employ involves energy transfer in binary (unscreened) Coulomb collisions (see appendix B). The free energy of partially ionized plasma includes, in addition to translational (ideal gas) and interaction terms, an electronic (internal energy) contribution from each ionic species 'a', which may be written as a sum over states, CL shifts the eigenvalues and truncates the summation, but leaves the numerical value of G a (β) little changed in situations of practical interest. Therefore, when the Saha equilibrium ionization formula is derived from a free energy minimization procedure, the plasma-screened ratio of densities ρ Z and ρ Z −1 of an ion's adjacent charge states involves the same CL quantity I (Z −1) : The screening enhancement factor in this density ratio appears to be the same as that in Salpeter's fusion reaction rates. In cases of weak screening (the Debye-Hückel limit), it actually is the same quantity, but in cases of strong screening (the ion-sphere limit) it is different, as discussed in appendix A in connection with equations (A.8), (A.12) and (A.13).
To summarize, our treatment of plasma screening produces changes in SB fusion probabilities that lead to the familiar Salpeter enhancement. For ionization and recombination processes, our treatment produces changes in ionization thresholds (CL) but not in the analytic form of the reaction cross sections. SB event probabilities defined by equation (15) can readily include CL by substituting the screened cross sections, equation (22), for the cross sections described in appendix B. In equilibrium situations, these CL modifications of electron-impact processes yield screened Saha density ratios.

Atomic kinetics for MD simulations
Most treatments of atomic kinetics in non-equilibrium plasmas are modern versions of the 'collisional-radiative' model of Scott and Hansen [10], Bates et al [28] and Fontes et al [29]. In this approach, one solves a set of coupled rate equations for the average densities ρ Z , j = (N Z , j / ) of ions in various quantum states 'j.' Choices are made for the ions and states to be included, the quantum numbers needed to label the states, the role of the radiation field and surrounding plasma (e.g. optical depth effects and CL), as well as the kinds of collisional and radiative processes important for populating and de-populating the states under consideration. The electron distribution is usually taken to be a Maxwellian at a given temperature. The temperature may be specified, or change according to an energy balance equation, while the local electron density ρ e is assumed to be an average value, which is related to the ion densities by charge neutrality.
Because inelastic bound-bound transitions usually are faster than transitions that cause ionization or recombination, and because excited-state populations usually are much smaller than ground-state populations, one approach to the collisional-radiative model breaks the calculation of populations into two steps: in the first, steady-state populations of all excited states of each ion Z are determined for specified values of the much larger, adjacent ground-state populations, ρ Z ,gnd ρ Z = states j ρ Z , j and ρ Z +1,gnd ρ Z +1 = states i ρ Z +1,i . In the second step, all the ground-state populations are evolved in time according to a collection of effective rates C Z for ionization and recombination, viz.
In general, the effective rate C ioniz Z contains contributions due to photoionization and twobody ionizing collisions, the latter being proportional to the density of incident particles. The effective rate C recomb Z generally contains contributions due to radiative, dielectronic and threebody recombinations, the first two being proportional to ρ e and the last to ρ 2 e . If Auger emissions following inner-shell ionization are included, more charge states are directly coupled and the right hand side of equation (25) contains additional terms [30,31]. This approach works well at low densities, and at high densities when plasma constituents are ionized at least into the M-shell. Moreover, as the discussion below will make clear, it is also relatively straightforward to implement with the SB model.
The effective recombination rates C recomb Z +1 defined by Scott and Hansen [10], Bates et al [28] and Fontes et al [29] are composed of a sum of terms, each of which corresponds to a transition from the continuum 'c' directly into a particular state 'j' of the recombined ion Z, multiplied by the probability p recomb Z ( j) that an electron in that state reaches the ground state ( j = gnd) before being ionized. Similarly, effective ionization rates C ioniz Z are composed of a sum of terms representing direct ionizations, plus excitations from the ground state to other states multiplied by the probability p ioniz Z ( j) that an electron in the excited state is then ionized before undergoing any other inelastic process. Thus, one has the expressions Extensive atomic data are required for accurate collisional-radiative calculations of plasmas having numerous bound states and several species, since the probabilities in equation (26) are determined from the complete matrix of transition rate coefficients among each ion's manifold of bound states. Opacity effects require additional quantities relating to lineshapes. For the SB model, data needs are exacerbated because-as we describe in appendix B-cross sections differential in energy, not just rate coefficients, are required to determine how energy is shared in three-body, collisional ionization and recombination events. Even so, in principle SB is a workable model if one tracks the states occupied by each ion's bound electrons and allows for spontaneous radiative decays between collision events.
When the plasma density is well below that of normal solids, there is a much simpler kinetics picture, the coronal model, wherein there are no effects of plasma screening or opacity. One assumes that excited-state populations are negligible, which means that ionization events arise only from ground states, i.e., p ioniz Z ( j > gnd) = 0. Additionally, recombinations to excited states are assumed to proceed to the ground state via radiative cascade without collisional interruption, i.e., p recomb Z ( j > gnd) = 1. Thus, only the total radiative and dielectronic recombination rate coefficients are needed, three-body recombinations being ignorable. With all these simplifications, the degree of ionization of any element in a plasma with Maxwellian distributions becomes a function of temperature only (see, e.g., [32]). As long as total, twobody recombination cross sections for radiative and dielectronic recombination, summed over all states of the recombined ion, are available, the SB scheme is easily implemented within the coronal model. For coronal plasma conditions, CL is unimportant and SB reaction probabilities are those of dilute plasma, equation (14).
When the plasma's density is near or above that of normal solids, CL severely limits the bound states in an ion of net charge Z − 1 to those having principal quantum number no larger than some small value, n n CL = [Z I H / I (Z −1) − 1/2] ∼ few, with I H = 13.6 eV. Consequently, at the high densities of most interest in the Cimarron Project, atomic data sets needed for a collisional-radiative kinetics calculation, or for the use of the SB model in an MD simulation, are much more modest in size. In either approach one keeps track of some excitedstate populations and accounts for radiative decays occurring between inelastic collisional events.
Fortunately, there is an even simpler description of the kinetics-the so-called 'bottleneck approximation' [33][34][35]-that is particularly well suited to MD with the SB model. By adopting it, SB can ignore all bound-bound collisional and radiative transitions and track only distinct charge-state populations {ρ Z } in the MD simulations. To motivate the bottleneck picture, we recall that radiative lifetimes of excited states vary with principal quantum number roughly as n x , with 4 < x < 5 [36,37] whereas lifetimes against collisional destruction (most frequent are n → n + 1 transitions induced by electron impact) vary roughly as n −y , with 3 < y < 4 [38]. These strong, opposing n-dependences lead to sharp variations in the probabilities of equation (26) in the neighborhood of a particular principal quantum number n * , which is the bottleneck n-value. At n = n * the total rate of depopulation is smallest. Above n * the recombination probabilities are very small, while below n * these probabilities are almost unity. In contrast, below n * the excited-state ionization probabilities are very small, while above n * the ionization probabilities are almost unity. In essence, rapid collisional mixing among bound states above n * causes them to act as a pseudo-continuum. Taking account of this information, the collisional-radiative coefficients C Z , defined in equations (26), reduce to the bottleneck coefficients B z , with the effective continuum beginning above level n * Note that these bottleneck coefficients are similar to those of the coronal model, except that this truncated sum of recombination terms includes three-body events; also, of course, the appropriate reaction probabilities are those of a dense plasma, equation (15).
To determine a numerical value for n * , we use hydrogenic expressions for radiative and collisional lifetimes of states in an ion of net charge Z − 1 [28], Together these two timescales yield n * = Z 6 7 × 10 17 cm −3 ρ e T I H 1/2 1/8.5 (29) for the principal quantum number nearest the bottleneck. As examples, consider the two plasmas discussed at the end of section 2. For highly stripped Xe ions in a hot, compressed D/T/Xe plasma, we obtain n * 3, so only one or two excited levels contribute to the recombination sums in equation (27). For the warm, solid carbon plasma, we obtain n * = 1, which means no excited states are included in the recombination sums for any carbon ion. As a practical matter, we round n * down to the next integer, but require it to be no smaller than the principal quantum number of the ion's ground state.
It is straightforward to incorporate the bottleneck approximation into our SB model. From the form of equation (11), we see that is only necessary to replace the coefficients defined in equation (27) by their corresponding reaction probabilities, viz., These reaction probabilities involve recombination and ionization cross sections for various radiative and collisional events that we present in appendix B, and they can be modified (by means of equation (22)) to include CL effects. If n CL < n * , we let plasma screening terminate the summation in equation (27); otherwise, the bottleneck approximation itself introduces an effective reduction in ionization potentials, I (Z −1) = [Z 2 /(n * + 1/2) 2 ]I H . A consequence of this simplification is that populations in levels between n * and n CL are effectively treated as part of the next higher charge state. Therefore, bottleneck values of Z will tend to overestimate those of a kinetics calculation that actually tracks all excited-state populations. Results presented in section 6 show the accuracy of the SB/bottleneck scheme under various dense plasma conditions.

SB implementation
We have implemented the quantum processes described above in our simulation code ddcMD [6,7] using the SB model. Appendix B contains a description of the cross sections and binding energies we adopted. We use the Stewart-Pyatt [25] model to describe the CL (but see appendix A for concerns about the accuracy of this widely-used prescription). When the CL I (Z −1) exceeds the binding energy of an atomic subshell, we exclude atomic processes involving that specific subshell. Three-body recombination is implemented using an effective two-body cross section given by equation (B.10). This requires an electron density ρ e and temperature T, which are calculated by averaging over the full simulation volume.
We verified the statistical assumptions underlying the SB model by extracting timeaveraged statistical properties of electrons entering the SB from MD calculations. Specifically, we followed the evolution of 50 eV carbon plasmas of atomic density 10 23 cm −3 . We considered three plasmas, each homogeneous and containing carbon that was one-, four-, or six-times ionized. We did not allow atomic or nuclear processes to occur. We monitored the electrons crossing the SB radius, chosen to be 1 Å.
The distributions of the kinetic energies are shown in figures 3(a)-(c). In all three cases we found that the general shapes of the distributions are consistent with a Maxwellian of temperature 50 eV, but they are scaled by a constant factor. According to discussion involving equation (4), this scale factor is the electron-ion pair correlation function, so these electron statistics allow us to extract g ei (R B ) from the SB statistics. We also have calculated g ei directly from the distribution of electron-ion distances. Table 1 shows that the two ways of deriving g ei from the MD agree with the HNC theory. These results confirm our statistical expectations for SB events.
The radiation field is represented by an isotropic homogeneous photon spectrum. During the simulation, we track the photon energy distribution in about 100 bins that are spaced logarithmically in photon energy. The radiation cross sections are averaged over the photon bin width. For laser-line radiation, we introduce additional narrow-energy bins if required. For photoionization processes, the photon fluxes associated with each bin are given by the product of the velocity of light and the photon specific number density, integrated over the bin's energy range. We account for all electrons and atomic nuclei in our simulation. Each electron is either bound or free. We simulate microcanonical and canonical ensembles, so that during each atomic physics reaction, an electron can change its state between bound and free, but it is not destroyed or created. Instead, bound electrons are assigned atomic subshell parameters (n, l) and required to follow the trajectory of their host nucleus. In this way, ions with bound electrons still can be treated as point charges with effective interaction potentials given by equations (1) and (2): the force on the ion is just the sum of the forces on the nucleus and its co-located bound electrons. This approach helps the symplectic integration algorithm used in ddcMD [5] by maintaining energy conservation because electrons can simply be released from the nucleus upon ionization, thereby avoiding sudden changes in the force field.

Atomic kinetics validation
We validated our SB MD model by comparing it to a non-local-thermodynamic-equilibrium collisional-radiative atomic kinetics code (Cretin) [39]. In the Cretin model, the electrons and ions are each assumed to have Maxwellian distributions but with (possibly) different temperatures. Coupling between the electrons and ions is calculated with a rate equation based on the classical theory of Landau and Spitzer [40].

X-ray photoionization and Auger decay
We simulated the evolution of solid-density carbon exposed to high-intensity 2 keV x-ray radiation with an x-ray flux of 8.5 × 10 4 photons Å −2 fs. These conditions are representative of high-peak-brightness experiments performed at free-electron-laser facilities such as the Linac Coherent Light Source [41]. In fact, an electron-ion equilibration experiment was recently performed to measure the evolution of electron and ion temperatures of graphite exposed to such a high-intensity x-ray beam [42], and one of our goals has been to simulate the evolution of the solid-density plasma using ddcMD. In these simulations, we consider only photoionization and relaxation by Auger decay. X-rays are absorbed predominantly through inner-shell photoionization of the K shell. Both Auger decay and outer-shell photoionization lead to depletion of the 2s and 2p subshells. Figure 4 shows the MD prediction of the ionization dynamics, and overlaid are the corresponding atomic kinetics results. The agreement here is very good. (For the proper simulation of the experiment, both electron impact ionization and threebody recombination are important [42], with the ionization process dominating and leading to higher charge states than shown in this test case.) Figure 4. Ionization dynamics in a solid-density neutral carbon system that is exposed to 4 keV x-ray radiation with a flux of 8.5 × 10 4 photons Å −2 fs. Shown are the MD results and the results from the collisional-radiative model Cretin. The number triples (a, b, c) indicate the number of electrons in the 1s, 2s and 2p subshells. We only show the populations that are larger than 5%.

Ionization evolution
We also simulated the evolution of two model systems, solid-density carbon (10 23 atoms cm −3 ) at conditions encountered in recent electron-ion equilibration and proton-stopping experiments, and high-density Xe-doped hydrogen relevant for laser fusion [5].
We first discuss the solid-density carbon system. We prepared two different carbon plasmas, one that is singly ionized and another that is four-times ionized. Each had 16 000 carbon atoms. We allowed both systems to relax without atomic processes turned on, using a Langevin thermostat to constrain the electron and ion temperatures to 50 eV. Maintaining the thermostat, we then instantly switched on collisional ionization and three-body recombination and let the system evolve to steady state. Figures 5(a) and (b) show the evolution of the different ionization states. It can be seen that the system equilibrates within a few femtoseconds since collisional processes are very effective at solid density. Also shown in figures 5(a) and (b) are the predictions from Cretin which was configured so that the atomic physics processes closely matched the MD model. It can be seen that both the rate of change of the populations of the different ionization states and the asymptotic long-term values agree fairly well. The remaining discrepancies are due to minor differences in interaction cross sections (quantities that are not well known for dense plasmas) and the aforementioned differences in charge-state counting by our bottleneck approximation.
The second system we considered was a hydrogen-xenon mixture with 2% Xe by number and a total density of 10 25 atoms cm −3 . Two such plasmas were prepared, one that contains 52 times ionized ground-state Xe and one that contains 44 times ionized ground-state Xe. The hydrogen was assumed to be fully ionized throughout the simulation. We equilibrated the H/Xe mixtures at a temperature of 5 keV in the absence of atomic or nuclear processes. Similar to the carbon case, we then instantly turned on collisional ionization and three-body recombination and let the system evolve to a steady state while maintaining the thermostat.  We have found that radiative processes can play a significant role in the equilibration dynamics of hydrogen-xenon mixtures. In figure 6(b) we included radiative recombination in addition to collisional ionization and three-body recombination, representing an optically-thin plasma in which the radiation is allowed to escape. As expected, the additional recombination channel leads to a lower Xe ionization state. In figure 6(c) we show the case of equilibration in an optically-thick plasma. Here, we assume an infinite medium which accumulates a radiation field in time as a result of plasma emission. This provides photoionization in addition to the processes accounted for in the previous simulations. In this case, the long-term value for the average ionization is similar to the case without radiation. This is to be expected since with radiation trapping the situation soon tends to one of radiative detailed balance, in which upward  [43].
is the stopping per unit density of those ions in the target medium. As derived in appendix B, the BE model gives where, for each subshell, N nl is its number of electrons, w nl = √ 2I nl /m e is their mean speed and the BE function is given by equations (B.13) and (B.14). Figure 7 shows fast proton stopping numbers for several neutral atoms, as predicted by the BE model and as inferred from stopping data (compiled and fitted by the SRIM program, [43]). Noble gases offer the cleanest comparisons, because results for these targets are not complicated by condensed matter effects. However, any consequence of shell structure, such as is evident in the BE curves, cannot be captured by the functional form the SRIM program uses to fit the numerous experimental results for each element. More importantly, the range of validity of the BE stopping model seems limited to projectile energies below a few MeV amu −1 . Even so, it is still a useful complement to the Bethe approximation, which is reliable at high projectile energies but which overestimates stopping at low E 0 -values [44][45][46].
Of course, the applications focus of our MD simulations is dense plasmas, and we are interested mainly in the contribution to fast-ion stopping made by electrons in partially ionized atoms. For these situations there exists neither data nor the atomic physics information needed to evaluate Bethe's stopping formula. Further, even if that information were available, it would not include the differential cross sections required by MD to treat individual ionization events. These circumstances lead us back to the BE model as the best choice, at present, for SB.
To gain insight into plasma ionization effects on stopping, we again consider 50 eV carbon plasmas of solid density, ρ C = 10 23 atoms cm −3 , but (artificially) having various uniform charge states Z (neutral to fully-ionized).  (31) and (33). Because data are plotted every time step, there are spikes corresponding to proton-carbon collisions, in progress, that were not sufficiently close to result in a significant energy transfer. (b) Stopping numbers for different carbon charge states, in a uniformly ionized plasma of ion density ρ C = 10 23 cm −3 and temperature T = 50 eV. In each case, the stopping numbers represent the sum of bound and free electron contributions, as given by equations (31)- (33). Dashed curves show the effect of CL under these plasma conditions. RPA dielectric function of a free electron gas [45], which yields an effective stopping number per carbon nucleus of Here, E max = 2 m e V 2 0 , E 0 = 1/2M 0 V 2 0 is the maximum energy transfer in a collision with an electron at rest, and ω e is the electron plasma frequency, which in these runs increases as √ Z . We see that, without a contribution from bound electrons, the simulation correctly shows negligible stopping in neutral carbon. For ionized carbon, the RPA dielectric function consistently predicts greater stopping than that occurring in the MD simulations. The discrepancy, which grows with Z, can be traced to the fact that our simulations use the (weaker than Coulomb) statistical potential, equation (1), to mediate energy exchange between plasma electrons and a fast proton. For any collision at small impact parameter, equation (1) results in a center-of-mass deflection that is smaller than the Rutherford expression, and hence (with reference to a fixed reference frame) a smaller energy transfer. A direct consequence is an underestimation of the true ion stopping. This explanation is supported by the fact that simulations involving pure Coulomb interactions between free charges of like sign yield stopping results that are much closer to the RPA prediction [5].
In figure 8(b), bound electron contributions are added to the stopping produced by free electrons to produce a total effective stopping number, viz., η Z = η RPA Z + η BE Z , for a carbon atom in various charge states; shell structure effects are still distinct for the lower Z-values. Differences among the curves illustrate the diminished effectiveness of bound electrons for fast ion stopping. Also plotted are curves showing the consequence of CL under these plasma conditions. The stopping enhancement due to this effect can be understood qualitatively by reference to figure B.1 (appendix B): CL reduces electron binding energies and thus, for a fixed projectile energy, increases each electron's contribution to its ion's stopping number.

Experimental validation: proton stopping in warm graphite
Experimental validation of our code is very important but is complicated by the limited phasespace volume that is accessible to both simulations and experiments. As part of the Cimarron Project, the stopping power of protons in warm dense graphite was measured [5]. Proton beams were generated by irradiating thin foils with intense, ultra-short optical laser pulses at the Jupiter Laser Facility (JLF) to volumetrically heat solid-density graphite foils that were several micrometers thick to a temperature of about 14 eV. Subsequently, the change in the velocity spectrum of a second, independent proton beam traversing the target was measured. It was found that the stopping power of 0.5 MeV protons is about (8.4 ± 2) eV Å −1 . The kinetics code Cretin estimated that the degree of carbon ionization in this experiment was approximately Z = 2.2.
In [5] it was pointed out that the simulated proton stopping power of a 20 eV carbon plasma with structureless C 2+ and two free electrons, but no bound electrons, is about three times smaller than the measured value, and it was suggested that the contribution of bound electrons to the stopping at warm-dense-matter conditions was an important, missing process. Our SB model allows us to investigate this point using MD, with the caveat that under these warm-dense-matter conditions all free electrons will still be treated classically. Figure 9 shows the calculated dependence of the proton stopping power on kinetic energy in a 14 eV carbon plasma. In these simulations, we equilibrated the solid-density plasma with 16 000 carbon atoms and 96 000 electrons at a temperature of 14 eV, resulting in an average ionization of Z = 2.1. Subsequently, while maintaining the thermostat, we launched ten protons of equal speed into this plasma, and monitored their drop in kinetic energy. The calculations were then repeated for different proton speeds. Our simulation results show that inclusion of stopping by bound electrons improves the agreement between MD and the experimental results at 0.5 MeV, but that a discrepancy of about 50% remains. However, if we keep the MD bound-electron contribution but replace the free-electron stopping calculated from MD with the classical RPA results (labeled MD[bnd] + RPA[free] in figure 9), we obtain good agreement with the experimental results. This is in concert with previous Cimarron studies of fast-ion stopping, as discussed just above, and reinforces the need for caution when using statistical potentials in electron-ion MD simulations of non-equilibrium phenomena. Ongoing Cimarron Project studies of electron-ion temperature equilibration [47] offer another example in this regard.

Fusion-relevant statistics in ddcMD-initial results
We have made ddcMD calculations of a proton and electron plasma in order to collect statistics of proton-proton collisions that reach a small distance of closest approach R B . The intent is to assess the validity of expressions such as equation (15) for the probability of a fusion event by directly counting the number of candidate events. We set the value of R B sufficiently small (0.01 Å) that the screening function H pp (R B ) is the same as the screening constant H pp (0). (See appendix B.) We perform the simulation for a time t, and whenever the separation of a pair crosses R B in the inward direction during the time step, we record kinematic information about the particles. Care is taken to ensure that the crossing events are counted once and only once when a pair approaches R B : an approximate trajectory is used to predict the time of the inward crossing of r = R B . If this time exists and lies in the interval [t, t + t), the tallies are made of the kinematic quantities interpolated to the crossing time.
The number of events is expected to be π R 2 B v ρ(R B )t, where v is the mean relative speed of the particles at this state and ρ(R B ) is the corresponding density. It has been assumed that the relative velocity is isotropically distributed in the incoming hemisphere, which we can verify. The local density ρ(R B ) is in principle related to the global particle density by the expression equation (13) involving the pair correlation function g pp (R B ), which, in turn, can be expressed in terms of the potential of mean force as exp(−U MF pp (R B )/T ). The mean speed v is the mean of the Maxwellian distribution of the relative particle fluxes, which is v times the Maxwellian distribution of the particle densities. We seek to verify these relations.
For our simulation we had 5 030 912 particles, half protons and half electrons, in a cubic box 63 Å on a side. This corresponds to a proton density of 10 25 particles per cm 3 . The temperature was chosen to be 5 keV. In order to reduce the amount of computation and thereby allow the desired number of events to be accumulated, the proton mass was reduced by a factor 100 compared with the physical proton mass. The duration of the simulation was 0.166 fs. The time step was 5 × 10 -7 fs. The statistical potentials suggested by Hansen were used, as described earlier.
The number of R B -crossing events, as described above, was 77 115 in this simulation. The histogram of the relative speeds recorded for these events agreed very well with v times the Maxwellian distribution of the relative velocity at temperature T; a chi-squared test indicated no significant difference. (The value of χ 2 /d.f. is 1.178; there is a 23% probability that χ 2 would be this large by chance, an insignificant result for identifying a deviation from the expected distribution.) The histogram of sin θ , where θ is the angle between the relative velocity vector and the line of centers is consistent with the expected distribution that is cos θ times an isotropic function, as shown in figure 10. (The probability distribution function (pdf) of θ is proportional to sin θ for an isotropic velocity distribution, but for the pairs crossing the surface r = R B in a specified time the pdf of θ is proportional to sin θ cos θ , which means that the pdf of b = R B sin θ is proportional to sin θ, i.e., to b.) These results, while unremarkable, do support the hypothesis that the velocity distribution of particle pairs in a close encounter is the same as the velocity distribution in the bulk plasma, as expected from the canonical phase-space density function, but in contradiction to some suggestions [48].
The predicted number of events if the HNC method is used to compute g pp (R B ) is 77 266. The value of g pp (R B ) itself is 0.7523. (This takes into account a small drift of T during the run of about 10 eV.) Thus the potential of mean force at R B is about 0.28 times T. The screening function given by the HNC theory is very small, 0.00238 times T, reflecting the fact that this plasma is very weakly coupled ( = 0.01). The screening constant with the Debye model is 0.00246. The number of events in the ddcMD run was 77 115, with a 1 − σ statistical error of 0.36% = 278. Thus both the HNC model and the Debye model are well within the error bar of the simulation. A more strongly coupled plasma than this one would be needed to be able to distinguish the models.

Summary
In sections 2-4 we laid the theoretical foundation for the SB method of treating atomic processes with a MD calculation. First, section 2 elucidated the connection between the plasma screening function H i j (r ) and the potential of mean force U MF i j (r ). In appendix A this is extended to a discussion of the alternative models of CL, since the magnitude of the CL can be identified with H i j (0). The potential of mean force and the screening function can be obtained directly from MD simulations and also from approximate methods such as HNC; examples of these were given. The ambiguity associated with finding the right screening function was also brought out in appendix A.
The essence of the SB method is that many-body classical dynamics is applied when interparticle distances r are greater than the SB radius R B , while two-body quantum mechanics is applied for r < R B . The radius R B is chosen to be larger than the relevant de Broglie wavelengths, to justify neglecting QM outside the SB, and such that the screening function is well defined within the SB, since the effect of the plasma environment is conveyed entirely by it. These ideas were developed in section 3, and in particular we gave a framework for expressing rates of quantum processes in terms of probabilities of events for particle pairs that arrive at R B . These probabilities are expressed in terms of isolated-atom cross sections as modified by the screening function H ab (R B ). The modifications consist of an energy shift by the potential of mean force, and also a focusing correction. One established instance of this procedure is the calculation of the thermonuclear fusion probability given that a pair of reactants has arrived at R B . The SB picture is also applied to the calculation of atomic ionization and recombination processes, for which it is important that R B is larger than the bound atomic orbitals in question. The application of the SB method to a full atomic kinetics calculation was treated in section 4. This is formally the same as the isolated-atom kinetics method (collisional-radiative model) except that the rates of the atomic processes are modified for the SB as described in section 3, and in more detail in appendix B. Considerable simplification results when the full set of collisional-radiative equations is replaced by the particular, restricted set that corresponds to the 'bottleneck' approximation; in these, one can ignore all bound-bound transitions.
Results from our SB implementation within the ddcMD code were described in sections 5 and 6. Section 5 concentrated on the statistical properties of the plasma and how the MD screening function compared with the alternative HNC model. Satisfactory agreement was obtained. In section 6 we compared the time-dependent kinetics computed with SB in ddcMD to the isolated-atom kinetics code Cretin, in the case that the atomic rates were made to correspond as closely as possible. Good agreement was found, and use of the 'bottleneck' scheme in dense plasmas was justified. Section 7 discussed stopping in partially ionized matter and then applied the SB method to an experimental problem: stopping of fast protons by warm carbon, as measured in an experiment on the JLF. In this case the MD results did not compare well with the experiment unless the free-electron part of the stopping was replaced with classical RPA theory. And, last, section 8 validated our expression for the fusion reaction rate in a dense plasma, given by equation (15), by directly comparing to counting statistics obtained from an MD simulation of a proton and electron plasma.
In this appendix, we summarize derivations of known formulae for plasma screening constants appropriate in the weak and strong screening limits, i.e., when the number of particles inside a Debye sphere is large or small, respectively. Additionally, we summarize key aspects of the model of Stewart and Pyatt [25], which interpolates smoothly between these limiting cases. Any of these expressions can be used by the SB model to evaluate screened reaction probabilities. Finally, for the two plasma examples discussed in section 2.3, we compare the CL predicted by these simple models with numerical results obtained from HNC calculations.
The contribution to the Helmholtz free energy arising from particle interactions plays the central role in concepts discussed here. In terms of the configuration integral defined in equation (3), this free energy term is Suppose that the N-body system includes N a particles of type a, N b of type b and N ab particles that are (a, b) composites, i.e., pairs whose separation is nil-for example fused nuclei or electrons recombined with ions. The screening constant H ab (0) is related to the change in this part of the system's free energy resulting from the formation of one more composite; specifically, where the reactants' N 's are negative and the composite's N is positive [49,50]. Given an analytical model of the plasma's electrostatic energy of interaction, U N [{N a }, T, ], the evaluation of equation (A.2) is easily accomplished after the free energy term is determined from a coupling constant integration [51], The integration represents a slow 'turn-on' of the Coulomb interaction, with scaled strength ηe 2 , throughout the system. Below, we use this technique to reproduce Debye-Hückel (weak screening limit) and ion sphere (strong screening limit) expressions for the screening constants H ab (0). The electrostatic potential of a point charge Z a , weakly screened in a medium containing numerous mobile charges, has the well-known Debye-Hückel form, where the screening length D is given by equation (8). At small r-values, DH a (r ) reduces to the Coulomb potential of particle 'a' plus a constant term, φ DH a = −Z a e/D, representing the potential due to all other charges. Thus, the total energy of interaction for this N-body model is When this expression is used in equation (A.3), the coupling constant integration gives the corresponding free energy as With this result we readily obtain the partial derivatives needed to calculate screening constants, namely, Now consider the fusion of two nuclei, with charges Z a and Z b . The screening constant for this case is easily determined to be This formula also applies to 'recombination,' when an electron (charge Z b = −1) unites with an ion Z a . Equation (A.8) was used by Salpeter [19] in his pioneering investigation of screening effects on fusion rates in stars. Therein, he also considered the ion sphere model as an approximation of the strong screening limit. In this model, each plasma ion is surrounded and completed shielded by electrons uniformly filling a sphere whose radius R a yields overall charge neutrality, (4π R 3 a N e /3 ) = Z a . Inside this sphere the electrostatic potential is and outside it is zero. At small r-values, IS a (r ) reduces to the Coulomb potential of the nucleus itself plus the constant potential φ (IS) a = −3Z a e/2R a due to the Z a electrons. In this model the plasma's total interaction energy U IS N is just the sum of internal interaction energies of individual ion spheres. For each ion sphere, the electrostatic energy U IS a has two parts-that of the nucleus interacting with the electrons, and that of the sphere's electrons interacting with each other, (A.10) By inserting U IS N = ions (N a U IS a ) (this sum over species does not include electrons) into equation (A.3) and performing the integration, we obtain the free energy of interaction for the ion sphere model of the plasma, In the second line above we use the mean inter-electron separation R e , defined in connection with equation (9); it is independent of any particular central ion's charge. Now, consider again the fusion of two nuclei of charges Z a and Z b . Use of equations (A.2) and (A.11) yields Salpeter's constant for plasmas with strong screening, In the ion sphere model, recombination results in one less central ion of a particular charge Z a , and one more of charge Z a − 1. The consequence of this is a screening constant different from equation (A.12), namely The ion sphere is a primitive, Thomas-Fermi-like model, and as such becomes more realistic as the central charge increases; hence, there is frequent use of the last expression above, which is valid when Z a 1. In both the weak and the strong screening limits, equations (A.8) and (A.13), respectively, the screening constants obtained for recombination are equal to minus the energy an electron has in the potential, at an ion Z a , arising from the other plasma charges: This fact motivates the third electrostatics model of interest here. Stewart and Pyatt [25] developed their model to provide a screening approximation that varies smoothly between the limiting Debye-Hückel and ion-sphere results. It is based on finite-temperature Thomas-Fermi theory [52,53], extended to include all plasma charges, not just those within an ion sphere. After several simplifications of Poisson's equation, including a restriction to non-degenerate electrons, Stewart and Pyatt obtain approximate solutions for SP a (r ) in the small-and large-r regions about a central ion of charge Z a . When these inner and outer solutions, and their first derivatives, are matched at an intermediate r-value, the results yield the potential φ SP a at the charge Z a due to the other charges. Its value gives the electron-ion screening constant .
Here, K a = −eφ DH a /T = Z a e 2 /DT , and Z * is defined in terms of the plasma's various ion charges as Z * = Z 2 ions / Z ions . The simplicity of this model, the fact that it reproduces known results for the limiting cases of weak and strong screening, and the fact that it can deal with situations far from Saha ionization equilibrium (an ability that HNC theory does not have) have led to wide use of Stewart-Pyatt screening constants to estimate plasma CL, as described in section 3.3 and as used in this work. However, it should be noted that the Stewart-Pyatt model does not provide the potential produced by the plasma at the site of an electron, so it cannot determine the total interaction energy U N [{N a }, T, ] that is needed to compute ion-ion screening constants for fusion events. Additionally, even though equation (A.15) exhibits the correct limiting behavior, how well eφ SP a approximates actual screening constants H ea (0) in plasmas with moderate screening is a topic that merits careful study. Here, we use HNC results, obtained for the plasma examples discussed in section 2.3, to provide two instructive comparisons. Figure 1(c) plots the screening functions H eZ (r/R Z ) for C 1+ and C 2+ ions in solid-density carbon at temperatures T = 25 eV or 50 eV. Under these warm dense matter conditions, the T-dependence is much weaker than the dependence on ion charge. The flat behavior of the screening functions in the limit r → 0 is what one expects on the basis of established results for one-component plasmas [43], and at each temperature the absolute magnitude of H eZ (0) is larger for the higher carbon charge state. Figure 2(c) plots Xe ion screening functions for two different 5 keV, highly compressed DT plasmas; in each case there is about 2% xenon (all Xe 44+ or all Xe 48+ ), by number. Again, limiting values of the screening functions are apparent as r → 0, but the cross-over near 10 −2 R Z and the ripples at larger separations are puzzling details not foreseen in the related HNC data plotted in figures 1(a) and (b); these features may indicate that conditions in these D/T/Xe plasmas are almost too extreme for convergence of the HNC numerical algorithm. Table A.1 collects CL values (the absolute magnitude of the screening constants for recombination), as determined by HNC theory and by the DH, ion-sphere (IS) and Stewart-Pyatt (SP) formulae, for (C 2+ + e) and for (Xe 44+ + e) in the 25 eV and 5 keV plasmas described above. For both of these realistic examples, SP is closer to IS than to DH, while HNC is closer to DH than to IS. For the carbon ion the ratio of HNC to SP CL values is 1.3, while for the xenon ion the agreement is even less satisfactory, this same ratio being 2.3. We have too few HNC results to adopt them for our SB simulations and so, for consistency, we have used only SP values. Further investigation of this problem is underway. calculations [54]; for all other cases, we use values from a screened-hydrogenic model [55]. CL can be included in cross sections dependent on binding energies by means of the substitution indicated in equation (21).

B.1. Cross sections for electron impact events
B.1.1. Radiative recombination. The Milne relation [37] connects the cross section σ (RR) n (ε 0 ) for radiative recombination of a free electron, energy ε 0 , into state (n ) bound by energy I n , and the cross section for the inverse photoionization (PI) event: Here, 1 N n 2(2 + 1) is the number of equivalent electrons in the recombined ion,hω = ε 0 + I n is the energy of the emitted photon, α is the fine structure constant and σ (PI) n (ω) is the photoionization cross section per (n )-electron. For photoionization of neutral atoms in their ground state, we use Yeh and Lindau's Hartree-Slater subshell cross sections [54]. For all other cases, we use the semi-classical Kramers formula We use the standard expression [37,56] to determine the probability of radiation emission in the angular frequency interval (ω, ω + dω) when an electron of initial relative energy ε 0 is scattered by an ion of charge Z. The quantity in the first square brackets is Kramers' differential free-free cross section. The factor [1 + n(ω)] accounts for stimulated emission, where is the mean occupation number for photons of energyhω in a radiation field with spectral energy density W (ω). The third factor, g ff , is the quantum free-free Gaunt factor; one has g ff 1 unless the relative collision energy is high, ε 0 > Z 2 I H , and/or most of that energy is lost to radiation.
To evaluate the probability of absorption in the time reversed event, where the 'incoming' electron has relative energy ε 1 = ε 0 −hω, we start with the corresponding absorption cross section, given by the microscopic reversibility principle as [37] σ (abs) ff (ε 1 , ω → ε 0 ) = 2 7 π 3 a 5 0 3 √ 3c This quantity, which has units (length 4 time), must be multiplied by the flux F(ω) dω = [cW (ω)/hω] dω of photons in (ω, ω + dω) to obtain the effective, two-body cross section for (electron-ion) scattering with absorption of radiation in that ω-interval. Equations (B.3) and (B.5) are appropriate when electrons and ions are weakly coupled and plasma screening has little effect on trajectories. Strong electron-ion coupling tends to suppress bremsstrahlung, but this complication is not treated here (see [26] for a brief discussion). and E 1 is the exponential integral function. Equation (B.10) is used to compute the SB probability of collisional recombination into a particular state.

B.2. Cross sections and related quantities for fast ion impact ionization
Interest here centers on projectiles with non-relativistic energies in the range E 0 = 1/2M 0 V 2 0 ∼ 0.01 − 1 MeV amu −1 . These energies, which encompass those of fusion products slowing down in a burning plasma, and also those of most light ion beams used to produce and probe warm dense matter [44], represent something of a quantum theory 'no man's land,' because generally E 0 is too high for coupled-channel scattering calculations to be practical, yet too low for the Born or Bethe approximations to be reliable. Moreover, experiments tend to focus on the important problem of charged particle stopping, viz., dE 0 /ds, but they cannot provide the detailed cross section information needed for an MD simulation. (Of course, such data do provide essential benchmarks for validation of simulation results.) Therefore, at present, one must rely on some classical description of the inelastic scattering of fast ions by partially ionized atoms.
Chandrasekhar's work on gravitational scattering in stellar systems [61] provides the best framework for classical treatments of charged-particle collisions [52][53][54][55][56][57][58][59][60][61][62][63][64]. To define a differential cross section for atomic ionization by a fast charge (Z 0 , M 0 , E 0 ), we start with the BE result of Garcia et al [64] for the total ionization cross section of an electron bound by energy I n (Ĩ n , if CL is included) and having a mean orbital speed w nl = √ 2I n /m e : and when ξ > 1.2071, . (B.14) It is important to note that, in the BE model, the charge of the target ion is irrelevant; its nucleus serves only to anchor bound electrons in space with some velocity distribution, here chosen to be a delta function at w n for each occupied subshell.
We obtain a simple differential cross section, per electron, from the fact that the dominant dependence on energy transfer E 0 = q I n , q 1, from projectile to an (n )-electron, in a BE is proportional to (1/ E 0 ) 3 [63]. With this, it is straightforward to determine an approximate differential cross section for ion impact ionization by multiplying the total cross section, equation (B.12), by the probability that the amount of energy transferred to the bound electron is in the interval d( E 0 ) = I n dq at q: P (FI) n (q) dq = 2 dq/ q 3  Here, we exhibit the subshell dependence of the velocity ratio, viz., (V 0 /w n ) → ξ n . For the target as a whole, contributions from all occupied subshells must be summed, An integration of this expression over the range of energy transfers indicated in equation (B.15) gives σ (FI) (E 0 ), the total cross section for ionization by a fast ion.
With the above results, one can compute quantities needed to determine the contribution of ionization to the slowing of fast ions in a plasma. In an ionization event, the mean energy transfer to a bound (n )-electron is This prediction is consistent with data on energy loss in noble gases [45]. A cross-sectionweighted average of this result with respect to subshells yields the mean energy lost per ionization in collisions at energy E 0 with a particular ion species, Given this energy loss quantity, the slowing of the fast ion (its rate of energy loss per unit path length) due to a species of plasma ions having number density ρ z and a total of N bnd bound electrons each, can be written in the form This is the same as the well known Bethe expression, useful for projectile energies, E 0 , much greater than 1 MeV, except that the BE quantity in square brackets has replaced the socalled Bethe stopping number, η Bethe Z = [N bnd ln( E max /I bnd )], where E max = 2 m e V 2 0 is the maximum energy transfer in a collision with an electron at rest, and I bnd is a mean transition energy for all of the ion's bound electrons (weighted by oscillator strengths) [44][45][46]. Note that, unlike the Bethe formula, the effective BE stopping number readily accommodates CL that arises in dense plasmas.
In figure B.1 we plot individual electron contributions to the BE stopping number (the bracketed quantity in equation (B.20)), for a wide range of electron binding energies and projectile kinetic energies. It is clear from these curves that, at high projectile energies, every bound electron contributes equally to the stopping predicted by the BE model. However, as E 0 decreases, loosely bound electrons become relatively more important for fast ion stopping. These outer electrons are the ones most affected by CL, and so their removal by this phenomenon can significantly reduce the bound electron contribution to stopping (but add to the free electron contribution). However, it is also clear from the sequence of curves in figure B.1 that binding energy reductions can increase the stopping contributions of tightly bound electrons since, effectively, for a fixed E 0 their stopping number curves are shifted to the left.

B.3. Photoionization and Auger decay probabilities
These two processes do not use the SB model because they do not involve particle collisions. Instead, probabilities of photoionization and Auger decay (following the creation of an innersubshell vacancy via x-ray absorption) are treated with a simple Monte-Carlo-type prescription.
As introduced in connection with the bremsstrahlung process, let F(ω) dω = [cW (ω)/hω] dω be the photon flux in (ω, ω + dω) associated with the radiation's spectral energy density W (ω). The rate of photoionization associated with the cross section σ (PI) n (ω) is After each time-step, and for each ion, a random number in [0, 1] is compared with this probability to determine whether the ion has undergone photoionization. Our simulation timesteps for electron-ion plasmas are so small ( t < 10 −18 s) that, even in the presence of intense radiation fluxes, the probability for an ion to undergo more than one event is negligible and the linearized expression P (PI) n t (PI) n is invariably valid. After the creation of a vacancy in an inner subshell, A = (n 0 0 ), of a light element, it is most likely that the ion undergoes a radiationless Auger relaxation involving two electrons, from higher subshells B = (n 1 1 ) and C = (n 2 2 ) that may be identical. One of these electrons fills the vacancy and the other is ejected into the continuum. (In heavy elements, there can be a series of such events before the ion is stabilized.) In the usual notation, this Auger decay rate is written as A ABC . The same argument as above leads to the following Auger probability for our MD simulations: (B.23) Again, after each time-step, and for each ion with an 'A' subshell vacancy, a random number in [0,1] is compared with this probability to determine whether the ion has undergone an Auger decay; here too the linearized probability expression P (Auger) ABC tA ABC is invariably valid because our time-steps are so small.
For singly-ionized atoms, we use McGuire's Auger decay rates [65,66]; these we denote A (1) ABC . To obtain approximate Auger rates A (Z ) ABC for atomic ions of net charge Z > 1, we scale A (1) ABC according to the number of electrons N (Z ) present in the ion's relevant subshells [30,41], viz.