Particle-in-cell studies of fast-ion slowing-down rates in cool tenuous magnetized plasma

Ions with velocities greater than those of most particles in the surrounding plasma are found across all scales of plasma physics, from high energy cosmic rays passing through the interstellar medium¹ to heavy ions propagating into solid targets.² Understanding the fast-ion energy-transfer processes in fusion plasmas with an energetic component has also long been a topic of considerable interest. How rapidly fast ions slow down is important to the dynamics of and energy balance in their background plasmas. The research on magnetic fusion energy (MFE) has examined situations as diverse as intense circulating beams formed by the ionization of energetic neutral beams used for plasma heating³ to tenuous isotropic mono-energetic distributions of fusion-produced alphas.⁴ The hot core plasma of fusion devices has been particularly well-studied. There, the magnetic field has little effect on slowing down: fast-ion velocities are much less than the thermal velocity of the background electrons, and the electron gyroradius, ρ_e, is large compared to the Debye length, λ_De.

This work looks at a less well-studied (still weakly-coupled) regime where the fast-ion velocity, v_fi, is comparable to the electron thermal velocity, $vth,e$⁠, and $η≡λDe/ρe≲2$ (see Fig. 1). Both of these conditions might be expected to inhibit the electron response to the energetic ion, and hence reduce the efficiency of electron drag on the ions and the slowing-down rate. Such conditions are present in the edge plasma of many MFE devices, though are considered relatively unimportant in large devices whose plasma radii, a, are much greater than $ρi$⁠. We will describe why they can be important to power and ash exhaust in small MFE devices, where $a∼ρi$⁠. Note that we restrict attention to cases where the energetic particle density is much less than the background plasma density and does not include transit-time-type resonances, e.g., excluding streaming⁵ and GAE-type instabilities,⁶ respectively (either could increase the energy loss rate).

Rostoker,⁸ among others,⁹ derived a Fokker-Planck operator incorporating the effects of an external magnetic field; only the asymptotic cases of $η≪1$ and $η≫1$ were considered in detail. Much of the work building on these early results has tended toward the analysis of strongly-magnetized systems (⁠$η≫1$⁠), with either very large magnetic fields^7,10,11 or very low temperature plasmas.¹² The results obtained are not applicable to the case of one unmagnetized projectile species (with trajectories unperturbed by the magnetic field) scattering on a weakly magnetized field species, since the field particles are only approximately confined to move along the field lines. Analytical work has been done on slowing down rates in weakly magnetized systems, but explicit calculations of slowing down times have been limited to the subthermal (⁠$v≪vth,e$⁠) regime, e.g., Hassan and Watson.¹³ Hassan and Watson found a correction to the field-parallel slowing down rate of −5.2%, and an asymmetry in the cross-field versus field-parallel slowing down rates of −2% (slower cross-field slowing down) for η = 2 and $v=1.34vth,e$ (see Fig. 2). One way forward in examining how best to extend these various analytical theories to our regime of interest where $η∼1$ was to study the issue via direct simulation. Recent work on direct simulations of slowing down in a molecular dynamics system¹⁴ and others^15,16 suggest that the problem may now be computationally feasible using particle-in-cell (PIC) codes.

The classical slowing down theory was first formulated by Rosenbluth et al., in 1957,¹⁷ followed by Rostoker's extension of the theory to a homogenous background magnetic field in 1960,⁸ using the BBGKY theory approach. Here, we briefly discuss nonmagnetized classical slowing down theory in order to emphasize a few of its scaling properties in the context of PIC codes and formulate an appropriate model for comparison with subsequent simulations.

To avoid confusion when comparing simulation results with theory, we use the following definition of the (energy) slowing-down time:

We start from the Fokker-Planck equation, following Bellan's treatment¹⁸ of Rosenbluth. We use the first velocity moment and assume that the fast ions (the test particles, subscript T) are mono-energetic and the background (subscript F for field particles, both electrons and ions) particles are Maxwellian, obtaining (MKS units)

where $ξF=mF2κTFu$ is the ratio between the fast ion velocity u and the field particle thermal velocity, while m_f and T_f are the mass and the temperature of the field particles, erf is the error function, and μ_F is the reduced mass. One of the approximations made by a PIC code is that each macroparticle in the simulation represents some number of real particles, denoted by the clumping factor ζ. The clumping factor modifies the properties of a macroparticle in the following way:

Note that all velocities and relevant frequencies are preserved by Eq. (3), as well as the Debye length. Incorporating these transformations into Eq. (2), we have

For the plasma conditions, we are interested in namely, $vT≫vth,i$⁠; the contribution of the background ions to the slowing down force will be negligible, so we restrict Eq. (4) to electrons only (field particles will henceforth be denoted by subscript e). We now consider the quantity $(1+ζTmTζeme)$ and note that if $ζT≳ζe$ and $mT∼mp$⁠, then $(1+ζTmTζeme)≈ζTmTζeme$⁠, simplifying Eq. (4) further

Note that Eq. (5) is independent of the electron ζ. Using the asymptotic approximations of $∂∂ξe(ξe−1erf(ξe))$ for large and small ξ_e, we have

Converting to energy, where $W=12ζTmTu2$⁠, and substituting into Eq. (6), we obtain the following two expressions:

For the $ξe≪1$ case, we have

As our particular simulations involve projectile particles with the mass of ⁴He, we finally specialize Eq. (8) to fast ions with $mT=4mp$ and $qT=Zfie$

This simple model will be used as the basis for comparison with simulation results (both the subthermal benchmarks and the non-subthermal scaling studies). In particular, it is the $Zfi2$ factor that will allow for orders-of-magnitude acceleration in slowing-down time simulations, when Z_fi is scaled up.

Baldwin and Watson derived “magnetized Rosenbluth potentials” for the case of a magnetized plasma where the ion Larmor radius is much larger than the Debye length (so straight-line ion trajectories are assumed), but the electron Larmor radius can be comparable to or less than the Debye length. When evaluated for the case of fast ions interacting with a Maxwellian background when $ξe≲1$⁠,¹³ this modified $h̃$ potential has the following form:

[K₀ and K₁ are modified Bessel functions of the second kind, with arguments of s², $Φ$ is the error function, r is $v||/vth,e$⁠, s is $v⊥/vth,e2$⁠, and $Γie$ is $4π(Zfie2/4πϵ0mfi)2$].

Mathematica was used to evaluate the friction coefficients using Eq. (10); the difference between the cross-field (⁠$v||=0$⁠) and field-parallel (⁠$v⊥=0$⁠) slowing down rates is plotted as a function of ζ_T and η in Fig. 2. For η = 2 and $ξT=1.34$ (the case of a 1 MeV triton), the cross-field slowing down rate (⁠$∼1/τ⊥$⁠) is 2%–3% smaller than the field-parallel slowing down rate (⁠$∼1/τ||$⁠), a small correction to the unmagnetized case. This is consistent with the intuition that as the system is only weakly magnetized, electron mobility near fast ions is not severely impacted.

To see this, consider a 1 MeV triton moving across a 5 T background magnetic field and a counter-propagating 100 eV electron, with some perpendicular separation $d<λDe$⁠. If $d=λDe/2$⁠, then the electron feels the effect of the triton for $∼0.075$ ns, and over the course of the interaction moves $∼0.03$ nm toward the triton (a very small fraction of its gyroradius) and gains $∼3×10−29$ J. For a background plasma density of 10¹⁴/cc, a stream of such counter-propagating electrons could remove 1 MeV from the fast ion in $∼3$ ms.

The PIC code used for these simulations is LSP,¹⁹ as has been used in previous efforts in this project, including initial magnetized slowing down simulations in 2D,²⁰ detailed studies of the effects of macroparticle clumping and collision models in LSP on slowing down,²¹ and two-stream instability studies in 2D showing enhanced slowing down for beams.²² All runs were performed using an energy-conserving explicit particle pusher¹⁹ and an explicit field solver; in the absence of fast ions, LSP conserves energy to better than 1 part in 10⁶ over the course of a 30 ns simulation. As a consequence of the use of explicit algorithms, an appropriate cell size was chosen to resolve λ_De (as well as ρ_e when including an external magnetic field), typically by a factor of 5–7. Note that although LSP includes a number of different subgrid collision models, we disabled them to avoid double counting the collisions, a point we shall return to later in this paper. A collisional model is not appropriate for these simulations because ionization and fluid effects are unimportant at the temperatures and densities of these weakly coupled plasmas, and the number of particles in the well-resolved Debye sphere is large. Thus, the frictional drag that the fast ions experience is entirely due to collective effects and communicated via the electric field grid.

One consequence of discretizing space is that the simulation Coulomb logarithm takes on a different value than it would as calculated from the plasma parameters in use; specifically, the lower bound on the integral is the cell size rather than the interparticle spacing. This simulation Coulomb logarithm was computed as

Here, the factor of 2 results from the use of higher-order particle shapes in LSP (to reduce particle self-force). The effects of this change were investigated by varying the cell size.

All frequencies (notably ω_pe and Ω_e) were resolved to at least 1 part in 100, and at least 60 periods of ω_pe were included in each simulation; the combination of these conditions ensured that the background electron distribution stayed Maxwellian. When a background magnetic field was used, it was uniform and oriented in the +Z direction.

A characteristic simulation consisted of a box in 3D Cartesian space with a side length of about $∼5λDe$ and periodic boundary conditions. Each cell size was no larger than one-seventh the Debye length, and with 64 particles per cell of each background species; there were at least 10⁵ electron macroparticles in each Debye sphere. The electron density in these simulations was chosen to be 10¹²/cc (unclumped), in order to relax some of the computational constraints, while the applied magnetic field was chosen to give a range of η between 0 and $≲2$⁠. This placed the simulation in the same dimensionless parameter space (of $η≲2$⁠) as 100× higher density and 10× higher magnetic field. Finally, both background species were created uniform in space, but Maxwellian in energy, with the electrons at 100 eV and the ions at 1 eV.

In order to maintain charge neutrality and avoid spurious charge buildup, the fast ions were injected with a contingent of neutralizing electrons with the same temperature as those in the background plasma. Most commonly, the fast ions were injected with isotropically directed velocities and homogenous in space, with monoenergetic starting energies. However, non-isotropic distributions were also used to explicitly study field-parallel versus cross-field slowing down rates. In order to not perturb the background plasma, the density of fast ions was kept to a very small fraction of the background, approximately 10⁶/cc, or about 1 fast ion per Debye sphere, a percentage that is comparable to MFE-device scrape-off layers. Also, in order to avoid nonlinear effects due to high-Z projectiles,²³ the fast ion charge was limited to $Z=2000e$⁠, giving any fast ion at most 1% of the charge in a given Debye sphere.

The majority of the simulations were run on NERSC clusters (Hopper, Edison, and Cori), with some performed on PPPL's Dawson cluster. Typical job sizes were 192 or 768 processors, with run times from 6 to 24 h; the restart functionality provided by LSP was used to extend simulations where needed. Computational efficiency was found to be around 2 ps of simulation time per hour of charged time, with slightly increased efficiency for the 192-processor jobs. Overall, a typical simulation consumed 20 000–50 000 CPU hours. However, the computational efficiency declined significantly as additional regions were used. Of the $∼6$ million CPU hours represented in Fig. 5, $∼1.1$ million were consumed by the 384-cell runs, and $∼3.3$ million were consumed by the 576-cell runs.

We benchmarked the PIC simulations for fast ions with velocities well below the electron thermal velocity and zero magnetic field, to explore the effects of varying the macroparticle clumping factor and charge. The subthermal regime was chosen to avoid wake effects from superthermal particles, which would require a much larger simulation volume. The fast ions were injected isotropically and homogeneously with the mass of an alpha particle and either 50 keV or 100 keV of energy (before clumping) each.

Figure 3 shows a typical time history for a high-Z simulation, along with an exponential fit; while the time histories of individual particles may diverge substantially from each other, the average exhibits good agreement with an exponential fit. Although the simulation is just under 8 ns long, the fast ions experience a drop in the energy of approximately half an e-folding (factor of $∼ 0.85$⁠).

In order to match the subthermal model with the simulation data and verify that we were using the correct calculation for the simulation Coulomb logarithm, a scan of the cell size was carried out, resulting in Figs. 4 and 5. In these runs, the simulation volume was held constant, while the number of cells per linear dimension was scanned.

Finally, the overall trend in τ_s with ζ_fi and Z_fi is shown in Fig. 6; the subthermal prediction includes the simulation $lnΛ$ factor, as determined from Fig. 5.

The simulations in this section encompass a range of magnetic field strengths up to 1 T, and initial fast ion energies of up to 3.6 MeV. Figure 7 presents a scan of energy at three magnetic field strengths. The slowing down time increases with initial fast ion energy in all cases, although slowing down times with and without a background magnetic field are comparable.

Figure 8 presents a higher-resolution scan in the magnetic field, with Z_fi = 2000 to decrease the required wall time. Three simulations with different starting particle microstates, for both the background particles and the fast ions, were performed for each data point. The slowing down time is 50% longer for a field of 1 T as compared to the 0.3 T case, while the 0.5 T case is approximately the same as the 0.3 T case. Figure 9 presents aggregate energy time histories for the fast ions at seven of the magnetic field strengths in Fig. 8. Finally, Fig. 10 presents a decomposition of the field-parallel and cross-field slowing down rates for the magnetic field scan simulations.

The subthermal benchmarking simulations follow the scaling with a fast ion charge, Z_fi, and a fast ion clumping factor, ζ_fi, as predicted by Eq. (9). For these simulations, the addition of a small background magnetic field (⁠$ρe/λDe=0.64$ for 5 kG and 0.32 for 10 kG) did not result in a systematic difference with the zero magnetic field case across a range of charges and fast ion energies.

The scaling shown in Fig. 6, an extension of previous work,²¹ allows for an effective speed-up on the order of 10⁶, vastly reducing the computational resources required for each simulation. This effect was harnessed to speed up the energy and magnetic field scaling studies involving non-subthermal fast ions.

The excellent agreement between the cell size scan data and the subthermal model when evaluated with the appropriate computational Coulomb logarithm in Fig. 5 verifies the formula in Eq. (11). Also, the predictions made by the classical slowing down theory are usually only guaranteed to order $lnΛ−1$⁠; however, the slowing-down times in Fig. 5 match the classical prediction to within a few percent when the effective computational $lnΛs∼1$⁠.

The magnetic field scan at Z_fi = 2000 (Fig. 8) reveals an unexpected structure for intermediate field strengths around 0.5 T. This was not observed to the same extent in simulations with smaller values for the fast ion charge. Possible explanations include nonlinear effects due to high-Z, and/or coupling between the fast ions and waves²² in the background plasma. The oscillations seen in the fast ion energy time histories from these data (Fig. 9) appear to lend some support to the wave hypothesis. Notably, the frequency of the oscillation is approximately equal to the lower hybrid frequency, pointing at the possible heating of the background plasma by the fast ions via waves. Additional simulations concerning this point are on-going.

The cross-field and field-parallel slowing down times in Fig. 10 were determined from projections of fast ion trajectories onto the cross-field and field-parallel directions. No consistent anisotropy (with respect to B_z) is present in these data, and more detailed simulations would be required to reveal the 2%–3% differences predicted in Fig. 2.

Although the slowing-down time increases with increasing fast ion energy, the slowing-down time for 1 MeV particles is at most a factor of 2 longer than for 100 keV particles. The addition of a weak magnetic field adds another factor of 1.5 (at η = 2). Thus, the overall increase in the slowing-down time for 1 MeV fast ions is a factor of 3 over the subthermal prediction.

Reactor-scale aneutronic (utilizing D-³He) FRC designs have been discussed previously.^24,25 Unlike many MFE reactor designs, the FRC designs referenced here do not rely on any core heating by fusion-produced alphas; all heating power is delivered via RF. Consequentially, the energy transport focus shifts to the efficient extraction of those alphas and other fusion products, both to remove ash from the core and recover the produced energy.

Computational studies of fusion product orbits in a reactor-scale FRC using the RMF code²⁶ have indicated that many fusion products will be born in orbits that traverse the separatrix.²⁷ In the cool scrape-off layer (SOL) plasma outside the separatrix, the fusion product experiences orders of magnitude more slowing down than inside. As a result, the fusion product loses energy to the SOL plasma; this “airbraking” effect causes the orbit to move out of the core and change character from a betatron to a cyclotron. The transition from an orbit traversing the separatrix to a cyclotron orbit on only one side of the separatrix is most likely to occur in the SOL, at which point the fusion product is then exhausted along the open field lines of the SOL. This process represents a convenient method for efficiently transporting ash out of the core,²⁵ facilitating both energy extraction (via conversion of the heated SOL plasma to electricity) and neutron production mitigation (rapid removal of T produced by D-D reactions).

Using Eq. (8), $τs=0.001$ s for T in the SOL (for the conditions of $ne=1014$/cc, $Eth,e=100$ eV, $lnΛ=14.7$⁠). Assuming that fusion products spend approximately 10% of their time in the SOL and 6 e-folding times (which should reduce the triton's energy from 1 MeV to $∼20$ keV), the slowing-down time is t_s = 0.060 s. As the average lifetime of T is approximately 20 s before fusion, this indicates that the vast majority of D-D-produced T should slow down and be lost along the SOL open field lines rather than fuse. Meanwhile, the combination of the scaling studies for the fast ion energy and the background magnetic field suggest an overall increase in slowing down times of no more than a factor of 2 as compared to the unmagnetized subthermal prediction. Thus, comparing the triton slowing-down to the triton fusion time, we estimate that the neutron production rate should be at least two orders of magnitude lower than in a conventional D-T tokamak, substantially reducing neutron wall loads and simplifying the engineering requirements of the reactor.

We have performed PIC simulations of fast ions slowing down on a cool weakly magnetized plasma for various macroparticles and charge scaling factors. Our simulations preserve the scalings of those factors predicted by the classical unmagnetized slowing down theory in a background plasma with unmagnetized and weakly magnetized electrons. This result indicates that using fast ion macroparticles with hundreds of times more charge than a corresponding real particle can allow relatively short simulations to access slowing down physics on an effectively longer timescale. Specifically, the simulated slowing down rate can be increased by a factor of over 10⁶ using fast ion macroparticles with greatly increased charge, allowing millisecond-scale slowing down physics to be accessed by a nanosecond-scale simulation. Simulations of MeV-scale fast ions (with $ξfi≡vfi/vth,e>1$⁠) indicate that the unmagnetized subthermal slowing down prediction should be increased by no more than a factor of 3 to account for the combined effects of superthermal fast ions and a weak magnetic field (⁠$1<η≡Ωce/ωpe<2$⁠). This result indicates that fusion products in a small aneutronic FRC reactor can be quickly exhausted via slowing down in the SOL, for efficient power and ash extraction as well as neutron production mitigation.

We thank Dr. G. Hammett and C. Swanson for useful suggestions, and the Program in Plasma Science and Technology for support under DOE Contract No. DE-AC02-09CH11466. Early contributions to this study were made by A. Creely, E. Paul, and E. Kolmes. Most simulations were performed on the Cori, Edison, and Hopper clusters at NERSC. Most plots were produced with Python3/Matplotlib.²⁸ 

The digital data for this paper can be found at http://arks.princeton.edu/ark:/88435/dsp01x920g025r.