Scholarly article on topic 'Modeling quantum processes in classical molecular dynamics simulations of dense plasmas'

Modeling quantum processes in classical molecular dynamics simulations of dense plasmas Academic research paper on "Physical sciences"

Share paper
Academic journal
New Journal of Physics
OECD Field of science

Academic research paper on topic "Modeling quantum processes in classical molecular dynamics simulations of dense plasmas"



Home Search Collections Journals About Contact us My IOPscience

Modeling quantum processes in classical molecular dynamics simulations of dense plasmas

This content has been downloaded from IOPscience. Please scroll down to see the full text. View the table of contents for this issue, or go to the journal homepage for more

Download details: IP Address:

This content was downloaded on 28/09/2014 at 09:26 Please note that terms and conditions apply.

New Journal of Physics

The open access journal for physics

Modeling quantum processes in classical molecular dynamics simulations of dense plasmas

S P Hau-Riege1'3, J Weisheit2, J I Castor1, R A London1, H Scott1 and D F Richards1

1 Lawrence Livermore National Laboratory, Livermore, CA 94550, USA

2 Department of Physics & Astronomy, University of Pittsburgh, Pittsburgh, PA 15260, USA


New Journal of Physics 15 (2013) 015011 (39pp)

Received 29 May 2012 Published 22 January 2013 Online at


Abstract. 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,

3 Author to whom any correspondence should be addressed.

Content from this work may be used under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

New Journal of Physics 15 (2013) 015011 1367-2630/13/015011+39$33.00

© IOP Publishing Ltd and Deutsche Physikalische Gesellschaft

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.


1. Introduction 2

2. Interactions in a dense, equilibrium plasma 3

2.1. Quantum statistical pair potentials................................................4

2.2. Pair distributions and potentials of mean force....................................5

2.3. Hypernetted chain examples......................................................6

3. Small ball (SB) model of quantum events 9

3.1. Preliminary considerations........................................................9

3.2. Reaction probabilities..............................................................10

3.3. The role of plasma screening......................................................12

4. Atomic kinetics for MD simulations 14

5. SB implementation 17

6. Atomic kinetics validation 19

6.1. X-ray photoionization and Auger decay..........................................19

6.2. Ionization evolution................................................................20

7. Proton stopping in dense plasmas 22

7.1. Model comparisons................................................................22

7.2. Experimental validation: proton stopping in warm graphite......................25

8. Fusion-relevant statistics in ddcMD—initial results 26

9. Summary 28 Acknowledgments 29 Appendix A. Plasma screening constants 29 Appendix B. Atomic physics for SB 32 References 38

1. Introduction

There now exists a considerable body of molecular dynamics (MD) research on dense, classical plasmas. Following seminal work by Hansen and colleagues [1-3], three decades of progress in computational hardware and software have enabled plasma MD to evolve from simple one-component 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-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, non-equilibrium 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-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.

2. 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 non-equilibrium 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 near-equilibrium 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.

2.1. Quantum statistical pair potentials

Since the aforementioned simulations by Hansen and co-workers [1-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 Za and Zb in a gas at temperature T, diffractive effects lead to the modified interaction

UDb (r) = (Za Zbe2 / r )[1 - exp(-2n r/Aab)], (1)

where Aab = ^2nh2/^abT is the thermal de Broglie wavelength of a particle whose mass /xab 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],

UePe(r) = T(ln2) exp{-(4nr2)/[(ln2)A2e]}. (2)

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.

2.2. 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/P and contained in volume Q, is

Here, UN is the system's total interaction energy, and {rN} denotes the coordinates of all the particles. If one assumes that UN is the sum of all pair-wise interactions Uab, 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 Na + Nb + • • • = N, positions {rai; sbj; ...} = {rN}, and mean densities pa = Na/ Q, etc. Then, the local density at r2 of species 'b' surrounding a particle of species 'a' at r1 can be written in terms of another integral involving the coordinates of all particles (see [18]),

^ ^ (Nb — 8ab)Q Cv _ _ _ _

Pbgba (S' S) = --- [8(ri — Sa 1 — Sb2)Jexp( — P Un{Sn})d{^N}. (4)

This equation serves to define the two-particle distribution function gba(r2, Si). 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 r2 depends only on the distance r = In — S21, so gab(n, S2) = gab(r) = gba(r).

From the barometric equation for the density of particles under the influence of an external field of force, viz., n(r) = n0 exp[—PU(r)], we see that the distribution function gab(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 Hab (r)that represents the average screening of this interaction by other charges, viz.,

Limiting values, Hab (0), which we refer to as plasma screening constants, are of particular importance to our SB model (cf section 3.3).

UMF(r) = -(1/j8) ln[gab(r)].

UMF (r ) = Uab (r ) - Hab (r ).

2.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-Huckel (DH) expression

UTC(r) ^ Uf (r) = ZaZb e2 exp(-r/D)/r, (7)


(4n e2 / T Pi Z2

is the Debye length involving all plasma species 'i.' This model is appropriate only when there is a large number of particles ND 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 Pion of all plasma ions, is the typical ion-ion separation Rion = (3/4npion)1/3. The second is Re, similarly defined in terms of the mean electron density pe = (Z)pion, 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 -Ze: Rz = (3Z/4npe)1/3 (cf equation (A.9) and [19]). The relation among these radii is

Rion =(Z)1/3 Re = ((Z)/Z)1/3 Rz. (9)

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 HeZ (0) . Because HNC theory cannot treat bound electrons explicitly, we approximate the influence of Nbnd such electrons on a free-electron-ion interaction via a sum of two quantum statistical terms, UeZ(r) = [UeIZ(r) + NbndUeP(r)]. (Effectively, bound electrons are 'pasted' onto their parent nucleus.)

The first plasma is warm, partially ionized, solid-density carbon, pion = 1023 cm-3, at temperatures T = 25 or 50 eV. At the lower temperature, the plasma is taken to be 2/3 C1+ and 1 /3 C2+; at the higher temperature, the ion fractions are reversed. (The values of ND 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 gee 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 Rion, 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).

....... ......... ......... ' : —(c):

. 25eV, C++,e 50eV, C++,e ------ " 25eV, C+,e ...... 50eV, C+,e ........ ll ........1 1

Figure 1. HNC results for warm, partially ionized, solid-density carbon. (a) Radial dependence of all six pair distributions gab (r/Rion) at 25 eV, (b) electron distributions geZ and (c) screening functions HeZ(r/RZ).

Figure 1(b) shows just the electron distributions geZ surrounding C1+ and C2+, but for both temperatures; here, each pair distribution is plotted as a function of r/RZ, where RZ is that particular charge state's ion sphere radius. Note that the free electron density is nearly equal to its mean value pe in the outer half of each ion sphere's volume. All four curves reveal essentially complete screening beyond r = RZ.

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 pion = 1026cm-3, and the

Figure 2. HNC results for hot, partially ionized, high-density H-Xe mixture. (a) Radial dependence of all ten pair distributions gab (r/Rion), (b) electron distributions geZ and (c) screening functions HeZ(r/RZ).

D/T/Xe mass fractions are 0.2/0.3/0.5. The first panel of figure 2(a) shows the radial dependence of all the pair distributions gab (r/Rion) 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 XeZ + ion sphere contains about Z/2 DT ions, mostly in the sphere's outer region. Figure 2(b) shows the electron distribution geZ surrounding Xe44+ and, to illustrate high-Z dependences, the electron distribution surrounding Xe48+, obtained from an otherwise

identical HNC calculation; here again, each pair distribution is plotted as a function of r/RZ, where RZ 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 = RZ.

Figures 1(c) and 2(c) show the screening functions HeZ (r/RZ) for Cl+ and C2+ (figure 1(c)) and Xe44+ and Xe48+ (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-Huckel 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.

3. Small ball (SB) model of quantum events

3.1. 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, RB(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 PY of a particular inelastic reaction 'y '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 UMF (r) are negligible for separations r > RZ, which implies that electron encounters with individual ion targets effectively begin at r = RZ. Therefore, we specify that RB ^ RZ. We assume that a pair's 'asymptotic' collision energy s0 is established at RZ, and we also assume that, once inside RZ, a colliding pair has fixed total energy. Hence, the relative kinetic energy is s(r) = [s0 — UMF(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 RB, 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 s0 via the potential of mean force, becomes more likely. Finally, one must recognize that each ion's RZ-value changes when it undergoes an ionization or recombination, so a single, fixed SB radius is inappropriate for a simulation with changing ionization.

3.1.2. 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, Za

and Zb, as long as the SB radius RB used for such an event is smaller than Rion. How much smaller can be set by considering the two kinds of ion-ion reactions presently of interest, thermonuclear fusion of light ions (e0 ^ 10keVamu-1), and energy loss of fast, non-thermal

ions (£0 ^ 10MeVamu-1) via ionization of heavy elements. For these energies, the classical distance of closest approach, r0 ~ ZaZbe2/£0, is also much less than the mean inter-ion spacing Rion. To ensure that candidate events for these kinds of reactions are not missed, one needs ion-ion SB radii larger than r0. Taking the target to be the more massive ion, with charge Za, values RBb > 0.01 Zaa0, where a0 is the Bohr radius, satisfy this condition as long as the collision energy exceeds several keV. (Section 8 contains further discussion of these points.)

3.2. Reaction probabilities

We develop the formula for the reaction probability PY 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 RB larger than the range of interaction between the pair (a, b), but still smaller than the mean separation of targets. At this value of RB, incident particles have a density equal to their mean value, pb = Nb/ft, and their flux has a Maxwellian distribution of energies £0: f (£0) = (4£0/nT3)1/2exp(-£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

dvb = P (4nRB) [4w(£0)] f (£0) d£0,

where w(£0) is the (relative) velocity. Equating the number of 'y time, PY times dvb, to that predicted by kinetic theory, we obtain

Py dvb = P [w(£0)M^0)] f (£0) d(£0), where aY is the ordinary reaction cross section. This yields

Py[ideal gas] = ay(£0)/nRB, (12)

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-Huckel interaction, viz Rion < D. The distribution of relative energies £B at RB < Rion 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 — UDH (RB). Moreover, the density of particles (b) at RB is not the mean density pb; rather, from the discussion in section 2.2, we have

Pb(Rb) = Pbgba(Rb) = Pb exp[-UaDbH(Rb)/T].

reactions per SB per unit (11)

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

PY [dilute plasma] =

ay [£B + UDH ( Rb )]

£b + UDH ( Rb ) £b

^y(£q )

( Rb )

where the first line is written in terms of the relative energy sB 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 RB due to an attractive (repulsive) interaction. (Note that the energy intervals are identical: d£0 at £0, and dsB = d£0 at £B = £0 — Uab(RB).) 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 UMF(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 UDbH(RB) is simply replaced by the corresponding potential of mean force, UMF(RB). 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 dsB = d£0 at sB = s0 — UMF(RB). 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.

3.2.1. Electron-ion events. A value of £0 determined at r < RZ 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 s0 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 RB = RZ. Then, UMF (RB) = 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.

3.2.2. Ion-ion events. Equations (14) and (15) exhibit a peculiarity of our classical model for repulsive, ion-ion collisions. Only collision energies s0 > Uab (RB) > 0 will result in energies £B > 0 at RB. This means that the contribution of lower-energy events to the total rate of ion-ion reactions 'y' will be excluded if Uab (RB) exceeds that reaction's threshold, Asy. 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 > 102 Asy ([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.

[£b + UaMF(Rb)] [SB + UMbF(Rb)!

Also, if accurate potentials of mean force for ion-ion pairs are unavailable, we believe it is best to keep the SB radius at RB = 0.01 Zaa0, 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).

3.3. The role of plasma screening

As we noted earlier, the presence of UabF(RB) in equation (15) relates to the environment's influence on the flux of incident particles at RB, 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 Hab (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 Rab (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

Kab(£0, r) dr

where the reaction-dependent quantity SY is only a weak function of energy, r0 is the classical distance of closest approach and the pair's effective wavenumber inside the barrier Uab (r) = ZaZbe2/r is

ay(£o) = (Sy/£o) exP \-2 Í L Jr

/ X ,2^abVTT Uab(r) - £0

Kab (£0, r) = J -—^[Uab (r) - £o] = —J---. (17)

h Aab V n T

This nuclear cross section involves the true Coulomb interaction, without modification.

To find the dominant influence of plasma screening, one assumes that SY is constant and that the effect of screening on r0 is ignorable. Then, with the replacement

Uab (r) ^ UOf ~ Uab (r) - Hab (0), (18)

it becomes evident that, with screening, the tunneling factor is computed at the effective energy £ = £0 + Hab (0). (Here and below, a tilde denotes a screening-related quantity.) This leads directly to the result that the screened fusion cross section aY satisfies the simple scaling relation

fusion: £0 ó Y (£0) = £ óy (£). (19)

Hence, the plasma-screened version of equation (15) for fusion is

PY [screened (a, b) fusion] =

Óy [£0 + Hab (0)]

YL V ' ' J nRB

' £0 + Hab (0)

£0 - UMbF ( Rb )

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 RB due to the pair's repulsive

interaction. This reaction probability leads to Salpeter's fusion rate enhancement factor of exp[p Hab (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 HeZ (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 < Rz/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 geZ 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 Hab (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 UeZ (r) — HeZ (0). Then, by inspection of the Schrodinger equation, it follows that if an orbital (nI) has unscreened energy Eirif—11), its perturbed, screened energy

~ ( z_i) ( Z_1)

is E(f ) = Ení ) — HeZ (0). Alternatively, one can state that screening has reduced its binding energy ¡nZ—1) by the amount AI(Z—11) = | HeZ (0) |:

ti—1) - nz—1) = I(Z—1) — AI(Z—1). (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: oy (s0; Inl) = oy (s0; Inl). (22)

We use this same formula to model screening effects on the cross section for ionization by fast ions, £0 ^ fewMeVamu-11, 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,

Na Y exp[—PEj)] = NaGa(P)exp[—pE0a)]. (23)

CL shifts the eigenvalues and truncates the summation, but leaves the numerical value of Ga 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 p z and p z — of an ion's adjacent charge states involves the same CL quantity AI( z-1) :

Pe P Z

_P Z-1_

Saha screened

= exp[+ß AI( Z-1) ]

Pe P Z

_P 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-Htickel 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.

4. 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 pZ, j = (NZ, j /ft) 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 pe 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, pZ,gnd - PZ = Estates j Pz, j and pZ+1,gnd - PZ+1 = Estates i Pz+1,i. In the second step, all the ground-state populations are evolved in time according to a collection of effective rates CZ for ionization and recombination, viz.

dpZ,gnd dpZ y^recomb . _ y^ioniz ioniz . y^recomb\

—dt— - "dT = pz+1 Cz+1 + pZ-1 Cz-1 - pZ\CZ + CZ ). (25)

In general, the effective rate CZ?niz contains contributions due to photoionization and two-body ionizing collisions, the latter being proportional to the density of incident particles. The

effective rate CZecomb generally contains contributions due to radiative, dielectronic and three-body recombinations, the first two being proportional to pe and the last to p2 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 cze+omb 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í^comb (j) that an electron in that state reaches the ground state (j = gnd) before being ionized. Similarly, effective ionization rates C^niz are composed of a sum of terms representing direct ionizations, plus excitations from the ground state to other states multiplied by the probability piZoniz( j) that an electron in the excited state is then ionized before undergoing any other inelastic process. Thus, one has the expressions

C*+omb = czcomb^ — gnd) + ^ C™^(c — j)p^(j),

j >gnd

_ (26) CZoniz = CZoniz(gnd — c) +J2 CZxck(gnd — j)pZoniz(j).

j >gnd

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., pZ?niz(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!^comb(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, two-body 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 ^ nCL = [Z^Ih/AI(Z—1) — 1/2] - few, with Ih = 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 excited-state 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-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 {pZ} in the MD simulations. To motivate the bottleneck picture, we recall that radiative lifetimes of excited states vary with principal quantum number roughly as nx, 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 CZ, defined in equations (26), reduce to the bottleneck coefficients Bz,

c recomb CZ+1

_x mecomb _ \ A ^»recomb/ _

^ BZ+1 — CZ+1 (c ^ j

j =gnd

CZT ^ BTZ — CTZ (gnd ^ c),

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],

(1/rnrad) — [2 x 1010Z4/n4-5] s-1,

(1/C11) — {[3 x 10-8Pe(cm-3)] (Ih/T)1/2n4/Z2} s-1

Together these two timescales yield

7 x 1017 cm

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.,

precomb[bottleneck] = V PZecomb(c ^ j), Z j=gnd Z j (30)

P^ [bottleneck] = PZ°niz (gnd).

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 nCL < n*, we let plasma screening terminate the summation in equation (27); otherwise, the bottleneck approximation itself introduces an effective reduction in ionization potentials, AI(Z-1) = [ Z2/(n * + 1/2)2] IH. A consequence of this simplification is that populations in levels between n* and nCL 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.

5. 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 AI(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 pe and temperature T, which are calculated by averaging over the full simulation volume.

We verified the statistical assumptions underlying the SB model by extracting time-averaged statistical properties of electrons entering the SB from MD calculations. Specifically, we followed the evolution of 50 eV carbon plasmas of atomic density 1023 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 A.

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 gei(RB) from the SB statistics. We also have calculated gei directly from the distribution of electron-ion distances. Table 1 shows that the two ways of deriving gei 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

's 'E 300

§ 250

n i— 200

LU < 150

¡¡r' 100

LU 150

!Q i— 150

Kinetic Energy E (eV)

Small Ball i i

MBx 1.10 -

Kinetic Energy E (eV)

Small Ball i i

MBx0.94 -

MB .........

Kinetic Energy E (eV)

Figure 3. Histograms of the distribution of the kinetic energy of electrons entering the SB for (a) C1+, (b) C4+ and (c) C6+. Also shown are the 50 eV Maxwell-Boltzmann distribution and the same distribution scaled by a factor of 1.15, 1.10 and 0.94, respectively; these are the electron-ion pair distribution values given in table 1.

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.

Table 1. Electron-ion pair correlation function in carbon at 50 eV with an atomic density of 1023 cm-3, extracted from MD calculations and calculated with HNC theory. For C6+, the HNC algorithm did not converge for this temperature.

gei (r = 1 A) from gei (r = 1 A) from pair- gei (r = 1 A)

Ionization SB dynamics distance statistics calculated with

state within MD within MD HNC theory

C1+ 1.15 ± 0.01 1.14 ± 0.01 1.1499 ± 0.001

C4+ 1.10 ± 0.008 1.09 ± 0.01 1.1020 ± 0.001

C6+ 0.94 ± 0.006 0.94 ± 0.003 —

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.

6. 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].

6.1. X-ray photoionization and Auger decay

We simulated the evolution of solid-density carbon exposed to high-intensity 2keV x-ray radiation with an x-ray flux of 8.5 x 104 photons A-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 three-body recombination are important [42], with the ionization process dominating and leading to higher charge states than shown in this test case.)

o °-6

it 0.4

0 x i-

0 2 4 6 8 10

Time (fs)

Figure 4. Ionization dynamics in a solid-density neutral carbon system that is exposed to 4keV x-ray radiation with a flux of 8.5 x 104 photons A-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%.

6.2. Ionization evolution

We also simulated the evolution of two model systems, solid-density carbon (1023 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 16000 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 1025 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 5keV 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. Figure 6(a) shows the evolution of the average ionization of Xe. The systems equilibrate within

1 1 i i

MD Cretin

Z=0 Z=0 _

-Z=1 Z=1

-Z-2 ----1=2

Z=3 Z=3

« l\ l\ \\ _ - Z=4 - Z=4

'/ \ '/ \ V ¡I Si'J^C,

Time (fs)

Time (fs)

Figure 5. Equilibration dynamics in a 50 eV carbon plasma with a density of 1023 atoms cm-3 for initially (a) singly-ionized carbon and (b) four-times-ionized carbon. Shown are the MD results and the results from the collisional-radiative model Cretin.

tens of femtoseconds. The predictions from the kinetics code Cretin are overlaid on figure 6(a). Both the rate of change and the asymptotic long-term values agree very well.

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

X f 48

X f 48

MD Cretin

H—I I I I I I I-1-1-1 I I I I I I-1-1-1 I I I I I I


H—I I I I I I I-1-1-1 I I I I I I-1-1-1 I I I I I I


Time (fs)

Figure 6. Equilibration dynamics in a fully-stripped 5keV H (plus 2% Xe) plasma with a density of 1025 atoms cm-3 (a) without radiative processes, (b) with radiative recombination (optically-thin plasma) and (c) with radiative recombination and photoionization (optically-thick plasma). Shown are the MD results and the results from the collisional-radiative model Cretin.

and downward radiative rates cancel. In all the cases shown in figure 6 the MD predictions for the rate of change and the asymptotic values agree very well with the predictions from Cretin.

7. Proton stopping in dense plasmas

7.1. Model comparisons

For fast protons with energy E0 = 1M0V02 and charge Z0 = 1, we first compare SB's binary encounter (BE) stopping model (described in appendix B) with other results. The discussion is simplified if we deal with a target ion's dimensionless stopping number nZ, defined such that

1 d E o 2 / Z o e2 \2

---0 = 8n a01A n z (Vo) (31)

pz ds 0 H\h V0 7 Z ( 0) V }

c 'a. a.

90 80 70 60 50 40 30 20 10 0


j_i_i_i_i_■ ■ ■

Proton energy (MeV)

Figure 7. Stopping numbers for neutral C, Ne, Ar and Kr atoms; solid curves represent predictions of the BE model used in SB (appendix B), while dashed curves represent SRIM program fits to experimental data [43].

is the stopping per unit density of those ions in the target medium. As derived in appendix B, the BE model gives

nBE(Vo) = £ NnlS(Vo/Wnl), (32)


where, for each subshell, Nnl is its number of electrons, wnl = V2Inl/me is their mean speed and the BE function S 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 MeVamu-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 E0-values [44-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, pC = 1023 atoms cm-3, but (artificially) having various uniform charge states Z (neutral to fully-ionized). These plasmas are irradiated by protons having an initial kinetic energy of 1 MeV. First, figure 8(a) shows the mean proton energy loss (averaged over ten simulations), as a function of time and as determined by the ddcMD code without the stopping contribution of bound electrons. Overlaid in this figure are stopping predictions based on the

i— CD C LU

4 6 Time (fs)

Figure 8. (a) Kinetic energy loss of a 1 MeV proton (averaged over ten protons) in 50 eV solid-density carbon plasmas at different, fixed ionization stages; impact ionization is not included. Overlaid in solid lines are losses to plasma electrons, as given by the random phase approximation, equations (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 pC = 1023 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

nRPA = Z ln(AEmax/h^e). (33)

Here, AEmax = 2 me V02, E0 = 1/2M0 V02 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., nZ = nRPA + nlE, 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.

7.2. 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 A-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 C2+ 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 16000 carbon atoms and 96000 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



15 --------


JLF i-1-1

SRIM[bnd+free] ............

MD[bnd]+RPA[free] .........

MD[bnd+free] -

RPA[free] --------

MD[free] ........

5 ----------------------!

0.01 0.1 1

Proton energy (MeV)

Figure 9. Proton stopping power in a 14 eV solid-density carbon plasma, both without and with bound electrons, labeled MD[free] and MD[bnd + free], respectively. Overlaid are SRIM fits to proton stopping data for cold carbon [43], labeled SRIM[bnd + free], results from the classical random-phase approximation (RPA[free]), as well as a recent experimental data point obtained at JLF [5]. Finally, with the curve labeled 'MD[bnd] + RPA[free]', we show the sum of the RPA results and the bound-electron portion of the ddcMD results, taken as the difference of the MD results with and without proton-impact ionization.

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.

8. 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 RB. 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 RB sufficiently small

(0.01 A) that the screening function Hpp(RB) is the same as the screening constant Hpp(0). (See appendix B.) We perform the simulation for a time t, and whenever the separation of a pair crosses RB 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 RB: an approximate trajectory is used to predict the time of the inward crossing of r = RB. If this time exists and lies in the interval [t, t + At), the tallies are made of the kinematic quantities interpolated to the crossing time.

The number of events is expected to be nRB(v)p(RB)t, where (v) is the mean relative speed of the particles at this state and p(RB) is the corresponding density. It has been assumed that the relative velocity is isotropically distributed in the incoming hemisphere, which we can

5 3000

6 2500

1 2000

^ 1500

j§ 1000 E

^ 500 0

0 0.002 0.004 0.006 0.008 0.01 Impact parameter b (A)

Figure 10. Frequency of the impact parameter b for particle pairs arriving at the audit radius. The value of b is RB sin 0. The distribution of b as obtained in the ddcMD run is shown by the red line. The green line is the relation expected if the angular distribution of 0 is isotropic weighted by a factor cos 0, i.e., f (b) a b, as described in the text.

verify. The local density p(RB) is in principle related to the global particle density by the expression equation (13) involving the pair correlation function gpp(RB), which, in turn, can be expressed in terms of the potential of mean force as exp(-UppF(RB)/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 030912 particles, half protons and half electrons, in a cubic box 63 A on a side. This corresponds to a proton density of 1025 particles per cm3. 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 x 10-7 fs. The statistical potentials suggested by Hansen were used, as described earlier.

The number of RB-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 x2/d.f. is 1.178; there is a 23% probability that x2 would be this large by chance, an insignificant result for identifying a deviation from the expected distribution.) The histogram of sin 0, where 0 is the angle between the relative velocity vector and the line of centers is consistent with the expected distribution that is cos 0 times an isotropic function, as shown in figure 10. (The probability distribution function (pdf) of 0 is proportional to sin 0 for an isotropic velocity distribution, but for the pairs crossing the surface r = RB in a specified time the pdf of 0 is proportional to sin 0 cos 0, which means that the pdf of b = RB sin 0 is proportional to sin 0, 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 gpp(RB) is 77 266. The value of gpp (RB) 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 RB 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 (r = 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 — a 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.

9. 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 Hij (r) and the potential of mean force UM (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 Hij (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 inter-particle distances r are greater than the SB radius RB, while two-body quantum mechanics is applied for r < RB. The radius RB 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 RB. These probabilities are expressed in terms of isolated-atom cross sections as modified by the screening function Hab (RB). 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 RB. The SB picture is also applied to the calculation of atomic ionization and recombination processes, for which it is important that RB 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.


We thank our Cimarron colleagues for advice throughout the course of this research; R M More and S B Hansen were particularly helpful during its early stages. Work of the LLNL authors was performed under the auspices of the US Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344; that of JW was supported by LLNL contracts to the University of Pittsburgh.

Appendix A. Plasma screening constants

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

F[Nnt][N, T, A] = -T ln(Qn/An). (A.1)

Suppose that the N-body system includes Na particles of type a, Nb of type b and Nab 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 Hab (0) is related to the change in this part of the system's free energy resulting from the formation of one more composite; specifically,

Hab(0) =- ^ (3FNm/9N^ ^AN, (A.2)

i =a ,b,ab

where the reactants'A N 's are negative and the composite's A N is positive [49, 50].

Given an analytical model of the plasma's electrostatic energy of interaction, UN[{Na}, T, A], the evaluation of equation (A.2) is easily accomplished after the free energy term is determined from a coupling constant integration [51],

F[Nm =\ UN(ne2) dn. (A.3)

The integration represents a slow 'turn-on' of the Coulomb interaction, with scaled strength ne2, throughout the system. Below, we use this technique to reproduce Debye-HUckel (weak screening limit) and ion sphere (strong screening limit) expressions for the screening constants Hab (0).

The electrostatic potential of a point charge Za, weakly screened in a medium containing numerous mobile charges, has the well-known Debye-Huckel form,

0DH(r) = (Zae/r) exp(-r/D), (A.4)

where the screening length D is given by equation (8). At small r-values, $DH(r) reduces to the Coulomb potential of particle 'a' plus a constant term, 0DH = —Zae/D, representing the potential due to all other charges. Thus, the total energy of interaction for this N-body model is

2 / \ 3/2 UDH = -2eD £ (NaZa2) = —e3 (n/Qr)1/2 £NaZa2 . (A.5)

all species \ a J

When this expression is used in equation (A.3), the coupling constant integration gives the corresponding free energy as

FDH = 3 UDH. (A.6)

With this result we readily obtain the partial derivatives needed to calculate screening constants, namely,

(d F^nt]/9 Na) = —(Zae)2/2 D. (A.7)


Now consider the fusion of two nuclei, with charges Za and Zb. The screening constant for this case is easily determined to be

HDbH(0) = [(Za + Zb)2 — Za2 — Zb] = ^Dt. (A.8)

This formula also applies to 'recombination,' when an electron (charge Zb = — 1) unites with an ion Za.

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 Ra yields overall charge neutrality, (4n R3aNe/3Q) = Za. Inside this sphere the electrostatic potential is

oa (r) =

i — — 3—r-

2 Ra\ R2

and outside it is zero. At small r-values, oa (r) reduces to the Coulomb potential of the nucleus itself plus the constant potential ^aiS) = —3Zae/2Ra due to the Za electrons.

In this model the plasma's total interaction energy UN is just the sum of internal interaction energies of individual ion spheres. For each ion sphere, the electrostatic energy U1^ has two parts—that of the nucleus interacting with the electrons, and that of the sphere's electrons interacting with each other,

IS 3(Zae)2 3(Zae)2 9(Zae)2

uas = —+ V a } =--^^. (A.10)

a 2 Ra 5 Ra 10 Ra

By inserting UN = ions(NaUf) (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,

FN = UN = - —J] [Za2 Na / Ra ]

£[ Z^ Na ]. (A.11)

In the second line above we use the mean inter-electron separation Re, 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 Za and Zb. Use of equations (A.2) and (A.11) yields Salpeter's constant for plasmas with strong screening,

fusion: HIS (0) = — [(Za + Zb )5/3 - Z5J3 — Zf]. (A.12)

In the ion sphere model, recombination results in one less central ion of a particular charge Za, and one more of charge Za — 1. The consequence of this is a screening constant different from equation (A.12), namely

9e2 9(Z e)2

recombination: H™(0) =—— [(Za — 1)5/3 — Z5J3} = ( a )

10RZ a J 10R;

3 Zae2

1 — {1 — ¿)"

2 Rz (A13)

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 Za ^ 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 Za, arising from the other plasma charges:

recombination: H™(0) = and HjS(0) ~ e^f. (A.14)

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-Huckel 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

OSP(r) in the small- and large-r regions about a central ion of charge Za. When these inner and

outer solutions, and their first derivatives, are matched at an intermediate r-value, the results

yield the potential at the charge Za due to the other charges. Its value gives the electron-ion

screening constant

SP SP T i[3(Z* + 1)Ka + 1]2/3 — 1}

HesaP (0) = e0f =--^----J-. (A.15)

) 2(Z* + 1)

Here, Ka = —e^DH/T = Zae2/DT, and Z* is defined in terms of the plasma's various ion charges as Z* = (Z2^/(Z)ions.

Table A.1. A comparison of CL predictions for two plasmas discussed in the text.

CL (eV)


C2+ + * 18.4 24.6 32.1 44.7

Xe2+ + * 689 2386 5621 7670

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 UN [{Na}, T, Q] 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 approximates actual screening constants Hea (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 HeZ (r/RZ) for C1+ and C2+ 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 HeZ (0) is larger for the higher carbon charge state. Figure 2(c) plots Xe ion screening functions for two different 5keV, highly compressed DT plasmas; in each case there is about 2% xenon (all Xe44+ or all Xe48+), by number. Again, limiting values of the screening functions are apparent as r ^ 0, but the cross-over near 10—2RZ 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 (C2+ + e) and for (Xe44+ + 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.

Appendix B. Atomic physics for SB

In this appendix we describe the specific cross sections and rates associated with the atomic processes used in the simulations. We use the symbol Z to represent the ionic charge seen by an active electron; thus, when bound, that electron belongs to an ion of net charge Z — 1. Binding energies for inner and valence subshells of neutral atoms are taken from Hartree-Fock-Slater

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 an(RR)(£0) for radiative recombination of a free electron, energy e0, into state (n£) bound by energy Ini, and the cross section for the inverse photoionization (PI) event:

„(rRR) [(nl)Nnl-1 £o ^ (nl)Nnl] = (4l + 3 - Nnl)

2IH ihœ 4£o \Ih

(œ). (B.1)

Here, 1 ^ Nni ^ 2(2i + 1) is the number of equivalent electrons in the recombined ion, hœ = £0 + Ini is the energy of the emitted photon, a is the fine structure constant and an(pI)(œ) is the photoionization cross section per (ni)-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

26n ala n ' 1 x 3

(œ) =

3V3 Z2 Vhœ

B.1.2. Bremsstrahlung. We use the standard expression [37, 56]


16na3 a2

Z 2 Ih\ 1

[1+ n(œ)]gff (œ),

3^3 \ ^0

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

n(^) = (n 2 c3/h v3 )W M (B.4)

is the mean occupation number for photons of energy h v in a radiation field with spectral energy density W (v). The third factor, gff, is the quantum free-free Gaunt factor; one has gff ~ 1 unless the relative collision energy is high, £0 > Z2IH, 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 £i = £0 — hv, we start with the corresponding absorption cross section, given by the microscopic reversibility principle as [37]

affabs) (£1, œ ^ £0 ) =

gff (œ).

This quantity, which has units (length4 time), must be multiplied by the flux F(v) dv = [cW(v)/hv] dv of photons in (v, v + dv) to obtain the effective, two-body cross section for (electron-ion) scattering with absorption of radiation in that v-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).

B.1.3. Collisional ionization and three-body recombination. Treatment of these processes in an MD simulation requires ionization cross sections that are differential in the energy of the ejected electron (£2), so that the distribution of energies of the scattered electron (£1) is correct. We use the simple Mott expression, based on a classical BE picture, for differential cross section of an electron with orbital quantum numbers (ni) and binding energy Ini (see [57]),

daMCott[(nl)£0 ^ £1 £2] 4na2

Inl + £2

Inl + £2

£0 — £2

£0 - £2

Energy conservation gives £0 — = Inl + £2, so this formula is symmetric with respect to £1 and £2. By convention, the ejected electron is identified as the one with lesser energy (viz., £2 ^ £1). The total ionization cross section for a target with Nni equivalent (ni)-electrons is defined as

vS[(nl)Nill£o ^ (nl)Nnl-1 £1 £2]

4n ag IH Nni

£0 Inl

= Nn i f Jo

(£0 - Ini)/2

£0 - I

daM0tt[(nl)£0 ^ £1 £2]

,lo4 -

£0 + Inü \ Inl,

For comparison, the familiar Lotz formula [58], which has a simpler, Born form,

aLCtz(£0) - 2.8na^(IH/£0Ini)Nn i [log(£0/4i)],

is about a factor of two larger in most of the important, near-threshold regime, 1 ^ (£0/Ini) ^ 3. However, as discussed by Salzmann [59], the Lotz formula is known to somewhat overestimate experimental values in this energy region, and the discrepancy is larger for higher target charges. We believe, therefore, that cross sections determined from equation (B.6) are sufficiently accurate for use in SB calculations.

Invoking microscopic reversibility, the cross section for the inverse process of three-body recombination can be written using equation (B.6), as [60]

a(3R)[(nl)Nnl-1 £1 £2 ^ (nl)Nnl£0] = (4l + 3 - Nni)

n2h3 £0 da(CI)[(nl)£0 ^ £1 £2] d£2

me £1 £2

As in the case of free-free absorption, this three-body reaction cross section has units (length4 time); when averaged with respect to the flux of non-recombining electrons it yields an effective, two-body cross section for the 'collisional' recombination, £2 ^ ni. Using the Mott formula, and assuming a Maxwellian flux of electrons whose mean number density is pe, we obtain

<Rt [(nl)Nnl-1 £2 ^ (nl)Nnl] = Pe (w(£1 ) f (£1 )aIM30Rt)[(nl)Nnl-1 £1 £2 ^ (nl)Nnl£0])

=Pe(4i + 3 — Nn*)(2n«0^[vj Ht' Y) '

where, with y = £2 / T and Z = h 1/ T, the function Y ( y, Z) is

Y ( y ,Z) =

L(y + Z)2 Z J

( y + Z)

eZ E1(Z)



and E1 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 E0 = 1 /2M0 V02 ~ 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 E0 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., dE0/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-64]. To define a differential cross section for atomic ionization by a fast charge (Z0, M0, E0), we start with the BE result of Garcia et al [64] for the total ionization cross section of an electron bound by energy Inl (/n^, if CL is included) and having a mean orbital speed wnl = V2Ini/me:

(E o) = 4n —

Z 0 e2 hVo


The quantity S is an algebraic function of the velocity ratio £ = (V0/wnl) = VE0me/InlM0. When £ < (V2 — 1)/2 = 0.2071, S(£) = 0; when 0.2071 < £ < 1.2071,

(f) = 3

and when £ > 1.2071,

2(3 + 2^2) 8(f + 1)


8(£ + 1) 8(£ - 1)

4? - 1

8?2(£ - 1)2 4£(£ - 1)J


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 wnl for each occupied subshell.

We obtain a simple differential cross section, per electron, from the fact that the dominant dependence on energy transfer AE0 = qInl, q ^ 1, from projectile to an (nl)-electron, in a BE is proportional to (1/AE0)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(AE0) = Inl dq at q:

Pt fHq) dq = 2dq/{q3 [1 - (lnl / Eo)2 ]}, 1 < q ^ (Eo / In t).


This gives the differential cross section:

da„(lFI)( E0 )

= 8n a;

HjVt) c^ t1 - (Inl/E0)21-1 s(^-l).

21-1 r,

d(A E o )

Here, we exhibit the subshell dependence of the velocity ratio, viz., (V0/wni) ^ Çni. target as a whole, contributions from all occupied subshells must be summed,

da(FI)( E 0 ) d(A E 0 )



( E 0 )

d(A E 0 )

(B.16) For the


An integration of this expression over the range of energy transfers indicated in equation (B.15) gives a(FI)(E0), 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 (ni)-electron is

2 Ini E 0

(A E 0 (nl)> =


Ini + E 0

This prediction is consistent with data on energy loss in noble gases [45]. A cross-section-weighted average of this result with respect to subshells yields the mean energy lost per ionization in collisions at energy E0 with a particular ion species,

(A E0 ( E0 )> =

/^(FI)( E 0 ).


1 d E(

= (AE0(E0)> a(FI)(E0) = 8na021

£ Nnl3(£nl)


£ Nnian(FI)(E0,£ni) (AE0(ni)>


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 pz and a total of Nbnd bound electrons each, can be written in the form

z0 e2\2 hV0)


This is the same as the well known Bethe expression, useful for projectile energies, E0, much greater than 1 MeV, except that the BE quantity in square brackets has replaced the so-called Bethe stopping number, n!ethe = [Nbndln(AEmax/Ibnd)], where AEmax = 2 meV02 is the maximum energy transfer in a collision with an electron at rest, and Ibnd is a mean transition energy for all of the ion's bound electrons (weighted by oscillator strengths) [44-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 E0 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 E0 their stopping number curves are shifted to the left.

0.001 0.01 0.1 1 Proton energy (MeV)

Figure B.1. BE stopping number, per electron, for electron binding energies of 3, 10, 30, 100, 300, 1000 and 3000 eV (left to right).

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(to) dto = [cW(to)/hto] dto be the photon flux in (to, to + dto) associated with the radiation's spectral energy density W(to). The rate of photoionization associated with the cross section an(pI)(to) is

r™ = / an<PI)(to)F(to) dto. (B.21)

Therefore, the probability that an ion experienced this event during a simulation time-step At is

P™ = {1 - exp[-At r™ ]}. (B.22)

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 (At < 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(pl) ~ At rpI)is invariably valid.

After the creation of a vacancy in an inner subshell, A = (n0£0), of a light element, it is most likely that the ion undergoes a radiationless Auger relaxation involving two electrons, from higher subshells B = (n1) and C = (n2£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 Aabc . The same argument as above leads to the following Auger probability for our MD simulations:

pAACger) = {1 - exp[-At Aabc]}. (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 PA"Buger) — AtAABC 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 aA^c . To obtain approximate Auger rates AABC for atomic ions of net charge Z > 1,

(1) 'AB

[30,41], viz.

we scale AA'BC according to the number of electrons N(z)present in the ion's relevant subshells

A(Z) _ A(1) AABC = aABC


NC1} - 8bc


© US Government.


[1] Hansen J-P, McDonald IR and Pollock E L 1975 Phys. Rev. A 11 1025

[2] Hansen J-P and McDonald I R 1981 Phys. Rev. A 23 2041

[3] Frenkel D and Smit B 2002 Understanding Molecular Simulation 2nd edn (San Diego, CA: Academic)

[4] Richards D F et al 2009 SC '09: Proc. Conf. on High Performance Computing Networking, Storage and

Analysis (New York)

[5] Graziani F R et al 2012 High Energy Density Phys. 8 105

[6] Glosli J N, Graziani F R, More R M, Murillo M S, Streitz F H, Surh M P, Benedict L X, Hau-Riege S,

Langdon A B and London R A 2008 Phys. Rev. E 78 025401

[7] Benedict L X, Glosli J N, Richards D F, Streitz F H, Hau-Riege S P, London R A, Graziani F R, Murillo M S

and Benage J F 2009 Phys. Rev. Lett. 102 205004

[8] Binney J and Tremaine S 2008 Galactic Dynamics 2nd edn (Princeton, NJ: Princeton University Press)

section 7.3

[9] Chung H-K, Chen M H, Morgan W L, Ralchenko Y and Lee R W 2005 High Energy Density Phys. 1 3

[10] Scott H A and Hansen S H 2010 High Energy Density Phys. 6 39

[11] Calisti A and Talin B 2011 Contrib. Plasma Phys. 51 524

[12] Calisti A, del Rio Gaztelurrutia T and Talin B 2007 High Energy Density Phys. 3 52

[13] Glosli J, Graziani F, More R, Murillo M, Streitz F and Surh M 2009 J. Phys. A: Math. Theor. 42 214030

[14] Surh M P, Ellis I N, Glosli J N, Graziani F R, Krauss W D, Murillo M S, Richards D F and Streitz F H 2011

IEEE Trans. Plasma Sci. 39 2620-1

[15] Jones C S and Murillo M S 2007 High Energy Density Phys. 3 379

[16] Dunn T and Broyles A A 1967 Phys. Rev. 157 156

[17] Deutsch C, Gombert M M and Minoo H 1978 Phys. Lett. A 66 381 Deutsch C, Gombert M M and Minoo H 1979 Phys. Lett. A 72 481

[18] Hansen J-P and McDonald I R 2006 Theory of Simple Liquids 3rd edn (Amsterdam: Academic)

[19] Salpeter E E 1954 Aust. J. Phys. 7 373

[20] More R 2012 private communication

[21] Walling R S and Weisheit J C 1988 Phys. Rep. 162 1

[22] Cox J P 1984 Principles of Stellar Structure vol I (New York: Gordon and Breach) chapter 17

[23] Seidel J, Arndt S and Kraeft W D 1995 Phys. Rev. E 52 5387

[24] Fisher D V and Maron Y 2003 J. Quant. Spectrosc. Radiat. Transfer 81 147

[25] Stewart J C and Pyatt K D 1966 Astrophys. J. 144 1203

[26] Weisheit J and Murillo M S 2005 Handbook of Atomic, Molecular and Optical Physics 2nd edn (New York:


[27] Hansen S B 2012 private communication

[28] Bates D R, Kingston A E and McWhirter R P 1962 Proc. R. Soc. Lond. A 267 297

[29] Fontes C J, Abdallah J Jr, Bowen C, Lee R W and Ralchenko Yu 2009 High Energy Density Phys. 5 15

[30] Weisheit J C 1974 Astrophys. J. 190 735

[31] Hau-Riege S P, London R A and Szoke A 2004 Phys. Rev. E 69 051906

[32] Hulse R A 1983 Nucl. Technol./Fusion 3 259

[33] Byron S, Stabler R C and Bortz P I 1962 Phys. Rev. Lett. 8 376

[34] Wanless D 1971 J. Phys. B: At. Mol. Phys. 4 522

[35] Kastner S O 1980 J. Quant. Spectrosc. Radiat. Transfer 23 327

[36] Bethe H A and Salpeter E E 1957 Quantum Mechanics of One- and Two-Electron Atoms (New York:

Academic) p 269

[37] Sobelman 11 1992 Atomic Spectra and Radiative Transitions 2nd edn (Berlin: Springer) chapter 9

[38] Fisher VI et al 1997 Phys. Rev. A 55 329

[39] Scott H A 2001 J. Quant. Spectrosc. Radiat. Transfer 71 689-701

[40] Landau L D 1937 Sov. Phys.—JETP 7 203

Lifshitz E M and Pitaevskii L P 1981 Physical Kinetics (Oxford: Pergamon) Spitzer L 1967 Physics of Fully Ionized Gases (New York: Interscience)

[41] Emma P et al 2010 Nature Photon. 4 641

[42] Hau-Riege SP et al 2012 Ultrafast transitions from solid to liquid and plasma states of graphite induced by

x-ray free-electron laser pulses Phys. Rev. Lett. 108 217402

[43] Ziegler J F 2004 Nucl. Instrum. Methods B 219-220 1027

[44] Atzeni S and Meyer-ter-Vehn J 2004 The Physics oflnertial Fusion (Oxford: Clarendon) chapter 11

[45] Dalgarno A 1962 Atomic and Molecular Processes ed D R Bates (New York: Academic) chapter 15

[46] Landau L D and Lifshitz E M 1977 Quantum Mechanics 3rd edn (Oxford: Pergamon) section 149

[47] Benedict L X et al 2012 Molecular dynamics simulations and generalized Lenard-Balescu calculations of

electron-ion temperature equilibration in plasmas unpublished

[48] Mao D, Mussack K and Dappen W 2009 Astrophys. J. 701 1204

[49] Widom B 1963 J. Chem. Phys. 39 2808

[50] Jancovici B 1977 J. Stat. Phys. 17 357

[51] Fetter A L and Walecka J D 1971 Quantum Theory of Many-Particle Systems (New York: McGraw-Hill)

section 7

[52] Feynman R P, Metropolis N and Teller E 1949 Phys. Rev. 75 1561

[53] Liberman D A 1979 Phys. Rev. B 20 4981

[54] Yeh J J and Lindau I 1985 At. Data Nucl. Data Tables 32 1

[55] Faussurier G, Blancard C and Renaudin P 2008 High Energy Density Phys. 4 114

[56] Landau L D and Lifshitz E M 1987 Classical Theory of Fields 4th edn (Oxford: Pergamon) section 70

[57] Kim Y-K and Rudd M E 1994 Phys. Rev. A 50 3954

[58] Lotz W 1967 Astrophys. J. Suppl. 14 207

[59] Salzmann D 1998 Physics in Hot Plasmas (Oxford: Oxford University Press)

[60] Sobelman 11, Vainshtein L A and Yukov E A 1981 Excitation of Atoms and Broadening of Spectral Lines

(Berlin: Springer) chapter 3

[61] Chandrasekhar S 2005 Principles of Stellar Dynamics (New York: Dover Phoenix) reprint of original 1942


[62] Gryzinski M 1965 Phys. Rev. 138 336

[63] Gerjuoy E 1966 Phys. Rev. 148 54

[64] Garcia J D, Gerjuoy E and Welker J E 1968 Phys. Rev. 165 66

[65] McGuire E J 1969 Phys. Rev. 185 1

[66] McGuire E J 1971 Phys. Rev. A 3 587