Drift-kinetic PIC simulations of plasma flow and energy transport in the magnetic mirror configuration

Plasma flow in magnetic mirrors with converging–diverging magnetic fields is of interest for fusion, plasma processing, and electric propulsion applications. The magnetic mirror is one of the oldest schemes for magnetic confinement, aiming to create conditions for controlled fusion. The geometry of a highly diverging magnetic field is the basis for advanced divertor configurations designed to handle the heat flow to the walls of the fusion reactor.¹ Related physical phenomena occur in the magnetic nozzle² and magnetic mirror (converging–diverging) configurations, which are used to convert internal plasma energy into the directed kinetic energy of the supersonic plasma to create thrust in plasma propulsion devices. All these applications require a good understanding of the mechanisms of plasma flow acceleration and energy transport in the magnetic mirror configurations of various settings, e.g., see Refs. 3, 4, and references therein. A recent overview of the work in the context of electric propulsion can be found in Ref. 5.

With simple assumptions of quasineutrality and isothermal electrons, plasma flow in the converging–diverging magnetic mirror configurations is controlled by the robust profile of the electrostatic potential, which is uniquely determined by the magnetic field and the sonic condition at the mirror throat.^6–8 More generally, accelerating potential structures may be generated due to the effects of the magnetic field profile and mechanical apertures, e.g., at the interface of the plasma source with the expansion region,^9–12 as well as plasma pressure anisotropies in real and phase space, e.g., due to presence of several species with different temperatures. These accelerating structures can be quasineutral or involve space-charge effects.¹³ 

In weakly collisional plasmas typical for many applications, plasma expansion in the mirror magnetic field will naturally lead to the development of pressure anisotropy, even if the plasma source is initially isotropic. The pressure anisotropies are strongly affected by reflections due to the magnetic and electric fields. Overall, in weakly collisional regimes, electron, and ion dynamics in the magnetic mirror are non-local¹⁴ due to the global nature of the electric field formed by global ambipolarity and quasineutrality constraints. Such a non-local electric field may also be coupled to the boundary conditions such as sheaths at the boundaries of finite-length systems. All these nonlinear and kinetic phenomena typically require numerical simulations.

Previously, the effects of a finite ion temperature on plasma acceleration in the magnetic mirror were studied using a fluid two-pressure adiabatic model and Boltzmann electrons.⁷ Here, a full drift-kinetic model for ions and electrons is used to study the formation of the global electric field profile, energy transport, and resulting plasma acceleration.

The drift-kinetic model in the paraxial approximation used in our work is similar to the models used in Refs. 15, 16, 17, and 18. In Ref. 19, the stationary distribution functions for electrons and ions as functions of the integrals of motion were used to find the electrostatic potential by numerical iterations, imposing the quasineutrality and ambipolarity (current-free) conditions and assuming empirically that the region of doubly trapped electrons in the expanding part is fully populated. With the above assumptions, the electron and ion temperatures and heat fluxes were calculated as moments of the distribution function. A similar model was used to analyze parametric dependences, in particular, the effect of heat fluxes on the ion and electron energy balance.²⁰ The time-dependent Vlasov–Poisson model was used to study the effects of non-stationary electron trapping.¹⁶ The effects of collisions were studied in the time-dependent Boltzmann–Poisson model¹⁷ that included the Bhatnagar–Gross–Krook collision operator. The drift-kinetic paraxial nozzle was studied with an implicit, conservative particle-in-cell algorithm in Ref. 18, employing novel boundary conditions to properly describe the expansion-to-infinity situation typical for propulsion applications.

Similar to previous works, in this paper, we address the formation of the electron and ion pressure anisotropy, electron and ion energy fluxes, and their role in plasma flow acceleration with an emphasis on the energy balance and energy conversion between electrons and ions. The emphasis of our study is on the collisionless case, and, therefore, in the base case, the region of trapped electrons in the expander is disconnected from the source region. In this collisionless case, the trapped region can be filled only by transient (non-stationary) processes and as a result of the numerical noise inherently present in PIC simulations. As a result, we observe a small number of electrons in the trapped region. In this limit, we compare the results of our kinetic studies with the extended fluid model for ions²¹ that includes the collisionless heat fluxes beyond the two-pressure adiabatic model.²² Whenever possible, we compare our results with previous analytical results^6,7 and numerical simulations^15–17,19 using Vlasov–Poisson and Boltzmann–Poisson models.

Our study is focused on open mirror fusion applications,^23–27 where the energy flux from the mirror is expected to be absorbed on the wall. Therefore, we take into account the absorbing wall and investigate how the mirror effects (mirror forces and particle reflections) affect the potential profile (including the sheath at the absorbing wall) and energy transport for large magnetic field expansion proposed to reduce the overall energy losses.^1,28 Previous studies^1,28,29 emphasized the role of collisions (inside the plasma source) in forming the electron distribution function in the loss cone region. In our model, we consider a Maxwellian plasma source with a completely full loss cone—the worst case from the perspective of energy losses. Collisions also affect particles in the expander region, in particular, by providing a mechanism for electron trapping in that region. To evaluate these effects, we perform a separate study using the model collision operator for electron–neutral collisions.

Our simulation model is described in Sec. II. Section III presents the results of the simulations demonstrating the effect of the ion acceleration, formation of anisotropic distribution functions for electrons and ions, finite ion temperature effects, energy transport, and energy losses. Section IV describes the comparison of the results from WarpX³⁰ and EDIPIC³¹ and the results of implicit simulations with EDIPIC. The effects of collisions are considered in Sec. V. Section VI presents the summary and discussion.  Appendix A presents the Poisson equation in the paraxial model. Boundary conditions are discussed in  Appendix B.  Appendix C summarizes the extended hydrodynamic equations for ions.

The motion of ions and electrons in our model is considered in the drift-kinetic approximation, assuming that the characteristic frequencies of all relevant processes are much smaller than the ion/electron cyclotron frequencies, $ω≪ωci≪ωce$⁠, and particles' Larmor radii are much smaller than the transverse length scale, so the motion of the guiding centers represents the particle dynamics.

A general setup of the simulations is represented in Fig. 1. In our model, we do not include ionization as the actual source of plasma. Our left boundary represents a junction of the simulated magnetic nozzle with the plasma source, where plasma is generated and heated. It is assumed that the residence time of particles in the plasma source is large compared to collision times, or equivalently, the mean free collision path is shorter than the length of the source. The electrons and ions are injected from the left boundary with the Maxwellian distribution and given temperatures assumed in the source. The electrons and ions are injected at equal rates. Any particles that are reflected back to the boundary (and crossing it from the right to the left) are reflected again back into the nozzle but resampled from the same Maxwellian distribution function assumed in the source for each species. Thus, the number of injected particles (per unit of time) in each species remains the same. In this setup, we may inject energy but not the charge so that the plasma injection is current-free (ambipolar). It remains current-free along the nozzle due to the flux conservation condition $nαVα∥/B=const$ inherent in our model for each species.

The electrons and ions are injected volumetrically in the narrow region near the left wall. Particle positions are sampled from the distribution $z/L=0.1×U1/4$⁠, where U is a uniform number from $[0;1)$ and L is the system length (1.2 m).

We do not employ any injection control nor impose any quasineutrality condition as it is done in alternative approaches in Refs. 17 and 18: the rates of injection for ions and electrons are fixed and kept equal so that the plasma flow is ambipolar (current-free). In this model, the only fixed parameters are the injection rate (equal for ions and electrons) and temperatures of the Maxwellian distributions from which injected particles are sampled. The resulting densities of electrons and ions throughout the nozzle (including the region near the left boundary) and actual energy distribution (the ion and electron effective temperatures) are established as a result of self-consistent particle and electric field dynamics including particle reflections from the magnetic mirror and the electric field. As noted above, a fraction of injected particles is reflected back to the wall by the magnetic field, and the electric field is formed self-consistently. Such particles are then reinjected with their velocities sampled from the original Maxwellian distribution of the same temperature. In what follows, this setup is referred to as the Maxwellian reflux model.³⁵ 

This simulation and injection model setup aims to represent plasma flow into the nozzle from the long plasma source with a given source of particles which is fixed by the chosen injection rate. The energy flux for each species is not fixed in this model: the particles gain energy during reflections back into the source and subsequent resampling. The energy flux into the nozzle per one ion–electron pair is a figure of merit of the mirror confinement systems and is one of the parameters of interest in our work.

In our work, we have used two open-access particle-in-cell (PIC) codes: the one-dimensional electrostatic direct implicit particle-in-cell (EDIPIC) code developed by D. Sydorenko,³¹ available at https://github.com/PrincetonUniversity/EDIPIC, and WarpX (Version 23.09), a massively parallel PIC code for kinetic plasma simulations,³⁰  https://github.com/ECP-WarpX/WarpX. We have modified both codes into the drift-kinetic versions. The EDIPIC is a momentum-conserving code, and for WarpX, we used the momentum-conserving option (it also has an energy-conserving option). The results from the two codes provide a useful benchmark comparison. Another important goal was to compare the results of the implicit EDIPIC code with the explicit WarpX, as elaborated in Sec. IV. The modified EDIPIC and WarpX codes can be found here: https://github.com/Maknagens/EDIPIC-Drift-Kinetic and https://github.com/Maknagens/WarpX-1D-Drift-kinetic.

In this form, it is a tri-diagonal matrix that can be efficiently handled using the Thomas algorithm. The magnetic field inhomogeneity effects in Eq. (14) were added to the Thomas algorithm used in the EDIPIC code. For WarpX, we directly solve Eq. (14) using the Thomas algorithm instead of using a prebuilt multigrid solver. This allows us to solve the Poisson equation exactly and implement the magnetic field effect. Additionally, the Maxwellian injection and reflux algorithms for the left wall were implemented in both codes. Finally, since WarpX lacked a default floating wall boundary condition, we included it analogously to the implementation in EDIPIC. All modifications to WarpX were made using the Python particle-in-cell modeling interface (PICMI), distributed with WarpX. We use the floating wall on the right wall, $z=L$⁠, and $ϕ=0$ on the left wall, $z=0$⁠.

As was mentioned in the previous subsection, the left wall has a fixed potential $ϕ=0$⁠. The right wall is absorbing and floating, which together with a global current-free condition guarantees the achievement of a stationary state. We have also tested alternative boundary conditions which are discussed in  Appendix B.

For the base case parameters, we consider the injection temperatures $T0e=300 eV$ for electrons, and $T0i=30 eV$ for hydrogen ions, with the injected flux (in units of the current) of $Iinj=10$ A/m². Some simulations are performed at different ion temperatures (30/300/600 eV), as specified at the top of the related graphs.

The base case magnetic field mirror and expansion ratios are $R=10$ and $K=50$⁠. The profile of the field from Eq. (15) is shown in Fig. 2(d). The size of the simulation domain is $L=1.2$ m. The PPC number (particle per cell) is 1000 on average. The Debye length is resolved in the whole domain—in the region of the largest concentration (at the left wall), we have 1.9 cells per Debye. There are 1179 cells distributed evenly across the domain, and the time step is $2.5×10−11$ s. Free streaming condition is $vT0eΔt/Δx=0.25$⁠, where $vT0e=2T0e/me$⁠. For the base case, we inject 7.6 particles (on average) at every time step. One particle is injected randomly when the generated random number is less than 0.6. There are no particles in the simulation domain at the start of the simulations.

The simulations reach a steady state at $∼$ 100 μs for $R=10$ and $K=50$ regime. The steady state is reached after many ions bounce between the left wall and reflections in the Yushmanov potential, so the injected ion/electron fluxes equilibrate at the right wall. The steady-state energy balance is reached approximately at the same time. The number of reflections required to reach steady-state increases with R since the injection rate is the same and the density in the source region increases due to improved confinement.

The profiles of the potential, ion velocity, ion concentration, etc., shown in what follows are averaged with the 3 μs time frame with 100 snapshots to suppress noise.

Our main results in this paper are based on full kinetic descriptions for ions and electrons. It is of interest to compare these with the quasineutral approach, in which the condition $ni=ne$ is used together with the Boltzmann relation $n=n0 exp(eϕ/Te)$ for the electron density, where $n0$ represents the ion concentration at the left wall, and the electron temperature $Te$ is assumed constant. Using the ion density obtained from our kinetic solution and assuming isothermal electrons, one can extract the potential as $eϕ=Te ln (n/n0)$ and subsequently move the ions according to the electric field from this expression. This approach shall be referred to below as the hybrid model and was used to characterize basic features of plasma acceleration in the magnetic nozzle in Ref. 8.

To highlight the effects of the mirror field on plasma acceleration, we first consider the plasma flow in the uniform magnetic field for the same injection model, boundary conditions, and $T0e=300$ eV and $T0i=30$ eV. This comparison, shown in Fig. 2, also provides a useful check of our results with respect to the previous work in Ref. 35, where the plasma flow under a similar injection model was considered in one-dimensional problem without magnetic field. In our model, the latter corresponds to the case of the uniform magnetic field. A notable feature of plasma flow in the uniform magnetic field is the formation of the source sheath at the left wall. The source sheath occurs as a result of the higher electron thermal velocity (respectively, the higher particle flux) when the electrons and ions are injected at equal rates. The positive electric field is then generated near the injection region retarding the electrons. As a result, the ions are accelerated supersonically in the narrow (of the order of a few Debye length) transition layer of the source sheath between the plasma source and quasineutral plasma region as it is visible on the density, ion velocity, and potential profiles in Figs. 2(a)–2(c). The potential drop in the source sheath depends on the ratio of the ion to electron injection temperatures and will decrease for larger ion temperatures³⁵ and even may become positive, changing the electric field to negative, e.g., as it was shown in Ref. 8. The rapid and discontinuous transition at the boundary of the plasma source and expansion region was also detected in experiments, e.g., Ref. 9.

Plasma streaming toward the absorbing material wall has been considered in numerous conditions for fusion and low-temperature plasma applications, including the classical sheath problem at the plasma-wall interface. In the simplest model, a plasma boundary sheath is formed to equilibrate ion and electron fluxes to the floating wall, leading to the Bohm condition $Vi>cs$ at the sheath entrance. In the standard model, the ions are assumed to become supersonic due to the acceleration in the quasineutral pre-sheath region. In the case of a uniform magnetic field, similar to Ref.35, the collector sheath is formed near the right (absorbing) wall, even though the ions are accelerated supersonically in the source sheath, as seen in Fig. 2(b); the electric field is zero across the quasineutral region, and there is no pre-sheath, see Fig. 2(c).

The source sheath practically disappears for the converging–diverging mirror configuration, but the collector sheath remains, Figs. 2(c) and 2(e). The case of the converging–diverging magnetic field demonstrates improved plasma confinement due to the mirror effect: one can observe a plasma density increase in the source region, Fig. 2(a). Also, one can note an increased value of the total potential drop across the mirror, pointing to the larger number of reflected electrons and thus to the lower electron energy losses, Fig. 2(c). Improved confinement is inherently coupled to plasma acceleration, which has a classical feature of the acceleration in the magnetic de Laval nozzle, where plasma flow becomes supersonic at the sonic point located at the nozzle throat. The plasma flow velocity is monotonically increasing along the nozzle with the sonic point close to the location of the maximum of the magnetic field. This is the point where the plasma flow is equal to the local ion sound velocity, as shown in Fig. 2(b), where Mach velocity $M≡V∥i/cs$ and potential are normalized with the electron temperature of the injection, $T0e$⁠, and $cs=T0e/mi$⁠. The sonic point condition regularizes the singularity that may occur at this point in the quasineutral plasma. As it was shown in our previous work,⁸ our simulations here further demonstrate that the transonic accelerating velocity profile is a robust attractor for plasma flow in the magnetic mirror.

It is worth noting that the exact expression for the ion sound velocity depends on the used plasma model. It is easy to define the ion sound velocity for isothermal electrons and cold ions⁶ and for finite ion temperature with CGL model.⁷ However, additional effects such as heat fluxes, ionization, and charge-exchange modify the definition of the ion-sound velocity and the location of the sonic point condition.^7,36

Nevertheless, one can see from Fig. 2(b) that the sonic point condition $M=1$ at the maximum magnetic field with $∂B/∂z=0$ is well satisfied in our base case when the contribution of the ion temperature to the ion sound velocity can be neglected. One notes a small jump in the velocity gradient that occurs at the sonic point, Fig. 2(b). It was shown previously that the velocity gradient at the sonic point is determined by the second derivative of the magnetic field at the maximum, Eq. (14) in Ref. 6. In our model, the second derivative of the magnetic field in Eq. (15) is not continuous at the mirror throat. This explains the jump in the velocity in Fig. 2(b) at $z′=0.5$⁠, as predicted by the analytical theory.⁶ 

In the mirror field, there are two complementary effects leading to plasma acceleration. One is the apparent acceleration due to the reflection by the mirror force on the left side from the magnetic field maximum. This is simply the filtering effect: particles with low values of parallel (axial) velocity are reflected by the mirror, so the passing particles, on average, have larger net axial velocity. Such passing particles are further accelerated by the mirror force on the right side from the magnetic field maximum. Thus, apparent acceleration is simply a ballistic effect in the mirror field and does not involve the electric field. True acceleration of ions is related to the electric field, which occurs due to the plasma expansion and associated density gradient. The generated electric field accelerates ions and reflects electrons. In this process, the electron (thermal) energy is converted into the kinetic energy of accelerated ions. Note that for the base case parameters here, $T0e=300$ eV and $T0i=30$ eV, the ion temperature is much lower than that for the electrons, and the ion thermal effects are not significant.

The conversion of the electron thermal energy into the ion kinetic energy (facilitated by the electric field) can be seen from the electron temperature profiles in Fig. 3. The perpendicular electron temperature drops drastically in the expander region for $z′>0.5$⁠: the electron mirror force acceleration in this region occurs at the expense of the electron perpendicular energy. Since the amount of the converted energy depends on the energy of passing electrons, and due to the asymmetry of the magnetic field (mirror region— $R=10$ and expansion region— $K=50$⁠), it leads to an increase in parallel temperature to $∼1.1T0e$⁠.

The behavior of the ion temperature is different from that of electrons: in Fig. 3(d) for the evolution of the ion perpendicular energy, one can see a weak trend suggested by the conservation of the perpendicular adiabatic moment for individual particles. The increase in $T⊥i$ in the converging part of the mirror is much weaker than expected from the perpendicular CGL (Chew–Goldberger–Low)²² adiabatic constant, $S⊥=p⊥/nB$⁠, $p⊥=(m/2)∫v⊥2fd3v$ so that $S⊥$ is not conserved along the length of the mirror. The ion parallel temperature drops dramatically in the expander region, as seen in Fig. 3. One can associate this with the strong density and magnetic field dependence in the parallel CGL adiabatic constant, $S∥=p∥B2/n3$⁠, $p∥=m∫v∥′2fd3v$⁠, $v∥′$ is random velocity $v∥=V∥+v∥′$⁠). In general, however, both two-pressure adiabatic invariants $S⊥$ and $S∥$ are not conserved due to the presence of the parallel heat fluxes. The role of the ion heat fluxes is discussed in further detail in Sec. III E.

The development of anisotropy in electron and ion pressures is a natural result of collisionless plasma expansion in the magnetic field. Anisotropy of the ion and electron distribution can be well characterized by the two-dimensional maps in $v⊥,v∥$ space, as shown in Figs. 4 and 5 for different axial locations in the mirror with $R=10,K=50$⁠. The distribution function is averaged over the interval $[z′,z′+Δz′]$⁠, $Δz′=0.02$ for each location.

The ion distribution function at $z′=0$ is an isotropic Maxwellian distribution with an empty loss cone for $v∥<0$⁠. Except for the particles with low energy, in the vicinity of the origin $v∥≃v⊥≃0$⁠, the shape of the ion loss cone at $z′=0$ is mostly determined by the reflections from the magnetic mirror with $R=10$⁠. Further down along the axial direction, the ions gradually form the beam distribution due to the acceleration in the electric field, and the magnetic field effects due to the conversion of $v⊥$ into $v∥$ following the conservation of $μ$⁠.

It is worth noting that our hybrid quasineutral model with $T0i=30 eV$ may experience numerical oscillations in the expansion region, leading to ion heating. The results in Figs. 6 and 7 were obtained with an increased number of particles per cell (PPC) which eliminated the noise. The mesh resolution was taken $0.01 m$⁠, the time step— $2.6×10−10 s$⁠, and the average PPC number was $104$⁠.

The fluid theory with the two-pressure adiabatic model for ions predicts⁷ that the finite ion pressure enhances plasma acceleration in the nozzle. This effect is related to the mirror force on ions with a finite perpendicular energy. This prediction is confirmed in our full kinetic simulations for two values of the ion temperature in the plasma source, 30 and 600 eV, Fig. 6(b).

In our simulations, for $T0e=300$ eV and $T0i=600$ eV, at the location of the maximum magnetic field, we observe $T∥i≃0.2T0i$ and $T∥e≃0.95T0e$⁠, from Figs. 7(a)–7(c). Then, in the units of the injection temperature for electrons, Eq. (19) gives $vs=2.15T0e/mi$⁠, which is fairly close to the ion velocity observed in simulations, see Fig. 6(b). One can note the jump in the velocity gradient at the sonic point, Fig. 8, due to the discontinuity of the second derivative of the magnetic field at the mirror throat.

One also has to note that we do not expect the full agreement of the simulations with Eq. (19), which was derived with the assumption of the constant and uniform (and isotropic) electron temperature and neglect of the ion heat fluxes. As one can see from Fig. 7, the electron temperature is not constant and not isotropic. The parallel electron temperature is slowly decreasing in the converging part of the mirror (to the left of the mirror throat). There is a significant drop in the perpendicular electron temperature in the expander, Fig. 7(b). The parallel ion temperature for $T0i=600 eV$ rises in the expander, exhibiting a temperature increase similar in nature to that observed for electrons. Additionally, it should be noted that in full kinetic theory, taking heat fluxes into account, the definition of ion-sound velocity changes, and the location of the sonic point may differ slightly from that of the magnetic throat.

One can observe a reduction of plasma density at higher ion temperatures, Fig. 6. This can also be explained by the increase in plasma exhaust velocity, including the region near the plasma source (to the left of the magnetic throat). In our simulations (for both values of the ion temperature), the injected particle flux is the same. Therefore, the increased ion velocity results in a lower plasma density.

It is interesting to look at the effect of the expansion ratio on the plasma flow in the mirror, keeping the same length of the mirror, as shown in Fig. 8. As expected, an increase in the expansion ratio results in the reduction of plasma density near the absorbing wall. Respectively, the sheath width is increasing. A notable effect is that the total potential drop across the whole mirror remains roughly constant, with the potential drop across the sheath and quasineutral regions also being very similar, as shown in Figs. 8(c) and 8(e). Figures 8(d) and 8(f) show how the potential profiles depend on the magnetic field. One may note that in the expander region, the potential does not really follow the magnetic field even for the quasineutral region, i.e., the potential is different for different magnetic field profiles even though the values of the expansion ratio are the same, while the total potential drop across the quasineutral region stays roughly the same, 8f. In these simulations with the constant mirror length, the increasing sheath width for large expansion becomes a large fraction of the total length of the expander as it is seen in Fig. 8(e) for the deviation from quasineutrality and in the plasma velocity profiles, Fig. 8(b). To further explore this behavior, we performed simulations by increasing the length of the expander, using Eq. (15) for the magnetic field and extending the expander length to larger values of $z′>1$⁠, Fig. 9. Figures 9(c) and 9(f) show the potential profiles for different expander lengths $Ltot$⁠. Now, the length of the quasineutral region and non-quasineutral sheath are much better separated. One can observe from Figs. 9(a), 9(c), 9(e), and 9(f) that plasma density continues to drop and plasma expands quasineutrally without any further acceleration and keeping the same total potential drop across the whole quasineutral region. The same effect was reported in Ref. 16. The total potential drop (and related plasma acceleration) ceases to increase since the perpendicular particle energy (both in ions and electrons) is fully “used up” and there is no more free energy reservoir for the acceleration. The plasma velocity no longer increases, and plasma density continues to decrease following $n∼B$⁠, corresponding to the expansion at constant velocity. This expansion occurs quasineutrally, as shown in Fig. 9(e).

It is worth reminding that due to multiple reflections of electrons back to the source, energy taken from the source can be much higher than the electron temperature in the source, $η>1$⁠. It is reasonable to assume that the bulk of the electrons in the source have the Maxwellian distribution due to the energy confinement time (residence time in the source) being much larger than the collisional time. It is important to note however that the distribution function coming from the source in the loss region of the phase space (loss cone with a cutoff due to the electric field) is expected to be different from the assumed Maxwellian. The distribution function in this region is established due to the balance of the collisional diffusion in the velocity space and the axial streaming of the thermal electrons away. Since in practical applications, the collisional electron–electron frequencies are many orders of magnitude smaller than the bounce frequency, the distribution function in this region is strongly depleted compared to that of the Maxwellian distribution, thus reducing the electron energy losses.^29,32 This is an important effect that has to be considered in the calculation of the actual axial losses from the magnetic mirror systems.^1,28 A full analysis of the effects of collisions is beyond the scope of the present paper. Our reflux model for the plasma source makes an assumption that the loss cone is full: the electrons in the loss regions are randomly replenished from the original Maxwellian at every return to the source region. From the perspective of the energy losses from the confinement volume, it is the worst possible situation. Therefore, our collisionless model with the Maxwellian reflux provides an upper bound for energy losses. The analysis of the electron and ion energy transport also confirms the energy conservation in the simulations and establishes the energy conservation properties for the collisionless base case.

Equation (27) shows that the total energy flux consists of the convective flux of the thermal (random) energy, the convective flux of the kinetic energy of the directed flow, and the heat fluxes. The ion and electron energy fluxes are shown in Figs. 10(a) and 10(c). The ions get their energy from the initial injection and from the accelerating electric field. The electron energy flux injected at the source gradually decreases due to the reflections of electrons by the mirror and electric forces. At the right boundary, the electron energy flux is determined by a small fraction of the electrons that overcome the total potential barrier in the Yushmanov potential and reach the wall.

The electrons lose their energy in reflections by transferring it to the electric field, which subsequently accelerates ions. A substantial exchange between electrons and ions mediated by the electric field is illustrated in Figs. 10(a) and 10(c). The total energy exchange between plasma and the electrostatic field is described by the last term in Eq. (30). For the current-less situation $J∥=0$⁠, the total energy in the electric field remains constant so that the sum of the electron and ion energy fluxes remains constant. Note that a momentum-conserving scheme is used in these simulations. We observe a weak variation of the total energy flux in the converging part of the mirror where plasma density is high. This numerical error is further discussed in Sec. V.

For the case of higher ion temperature, $T0i=600$ eV, the total energy loss for an ion–electron pair is around $9T0e$⁠, where the order of $3T0e$ (⁠ $≈3T0i/2$⁠) ions get from the injection, $5T0e$ is transferred from the electrons, and $T0e$ is carried by the over-barrier electrons, Fig. 10(c). The latter value is about the same energy as for the case of the lower ion energy. Since, in the current-less case, the total energy flux is constant in the axial direction, the net energy flux to the wall is a sum of energy fluxes in electrons and ions at $z′=0$⁠. These fluxes can be calculated from the anisotropic distribution functions at $z′=0$⁠, as shown in Figs. 5 and 4. The shape of the electron and ion distribution at this location is determined by the loss cone which, for each species, is determined by the combination of the magnetic mirror and electric field reflections. It is worth noting that the injected flux Maxwellian distributions carry very large energy fluxes, e.g., of the order of $Q∥≃nvtαTα$ for each species. However, a large fraction of injected particles are confined, reflected back to the wall, and then replaced by random particles from the distribution with an original injection temperature, so the net energy flux of particles inside the loss cone is much lower.

For all our simulations, the energy loss factor can be directly related to the total potential drop as $η≈(−eϕtot+T0e+3T0i/2)/T0e$⁠, e.g., Fig. 11 for the same parameters as in Fig. 8 shows that $η$ for the K values between $K=10$ and $K=1000$ are roughly the same.

An assumption of isothermal electrons is often used in simple theories of plasma flow in the magnetic mirror. It is of interest to compare the results of the full kinetic calculations with the hybrid quasineutral model with isothermal electrons. This comparison is shown in Figs. 12–14 for $T0e=300$ eV and $T0i=600$ eV. Obviously, the quasineutral model does not reproduce the sheath physics. The plasma potential in the full kinetic model closely follows the potential from the hybrid model in the converging part of the mirror. In the expanding part, the potential and ion velocity profiles start to diverge strongly, especially for large expansion ratios. In part, the difference occurs because of the influence of the sheath structure. Another, probably more important, reason is the electron cooling and the anisotropy of the electron pressure in the expander region. In the converging part, the perpendicular electron energy remains pretty much constant, Fig. 13(b), due to a large number of electrons reflected by the mirror. The reduction of the parallel electron energy in the converging part is not significant, Fig. 13(b), so the electron temperature in this region may be considered approximately isotropic and isothermal, similarly as in the hybrid model. In Ref. 19, the authors observed the same behavior. As a result, all plasma parameters (density, ion velocity, and potential) in the converging part of the mirror in the hybrid and kinetic models are very similar, as seen in Figs. 12(a)–12(c). In the expander region, the electron perpendicular temperature drops sharply, Fig. 13(b), since many of the electrons with large values of the perpendicular energies are reflected before the mirror maximum. The parallel electron temperature also decreases, Fig. 13(a). Therefore, for large expansion ratios, electron cooling that is fully included in the kinetic model is a dominant mechanism of the difference between the kinetic and hybrid models. One can see from Figs. 12(b), 12(c), and 12(e) that the potential and ion velocity profiles start to diverge well before the deviations from the quasineutrality become noticeable. A similar effect is also observed in Figs. 6(c) and 6(e).

The main differences between the hybrid and full kinetic models are in the potential and ion velocity profiles, Figs. 12(b) and 12(c). The total potential drop and ion velocity in the hybrid model increase with the expansion ratio, $δϕ≃Te ln K$⁠. In the hybrid model with the isothermal electrons, the energy transfer from electrons to ions remains large (formally infinite) in the expander region, resulting in the higher kinetic energy of the accelerated ions, diverging logarithmically with K, Fig. 12(b). In the full kinetic model, the changes in the parallel and perpendicular electron pressure are related to the compression/expansion work on electrons, the work of the electric fields, and heat fluxes, as seen in Eq. (27) for the total energy flux. In our injection and reflux model, the energy transferred from the electrons to ions is fixed by the amount of energy injected into the electron component. As a result, the total potential drop (and the final ion energy) quickly saturates for large K and becomes independent of K, see Figs. 18(c) and 18(b).

The ion temperature and ion heat fluxes (both perpendicular and parallel) are very similar in the hybrid and kinetic models, as seen in Figs. 13(c), 13(d), 14(a), and 14(b). The ion temperature effects on plasma flow in the mirror configurations were studied previously with two-pressure anisotropic ion pressure models in Refs. 7, 39, and 40. The simplest model⁷ is based on the two-pressure adiabatic CGL model,²² which assumes that the heat fluxes are zero, which is, in general, not justified. In weakly collisional plasmas relevant to many applications, the closures for the heat fluxes are difficult.⁴¹ Different closures for the ion heat flux due to collisional, ad hoc free streaming corrections, wave-particle interactions, and anomalous transport contributions were proposed^41–43 and implemented to model plasma flow in the scrape-off layer (SOL) and flux-expanding divertors.^39,40,44,45 The role of passing and trapped particles and ambipolar potential on the heat fluxes along open magnetic field lines was discussed in Refs. 46 and 47.

We have previously used an extended hydrodynamic model for ions to model collisionless plasma flow in the mirror field, taking into account finite ion temperature effects.²¹ In this model, the heat fluxes $q∥i$ and $q⊥i$ are included in the evolution of the ion parallel and perpendicular pressures. The heat fluxes themselves are calculated from the time-dependent evolution equations for $q∥i$ and $q⊥i$⁠. These equations are closed by the assumption that all fourth-order moments are calculated by using a two-temperature Maxwellian distribution. For completeness, the full system of extended hydrodynamic equations for ions is given in  Appendix C. In Ref. 21, characteristics of plasma acceleration were obtained using extended fluid equations for ions, Boltzmann approximation for electron density, and the quasineutrality assumption, as in the hybrid model. In Figs. 15 and 16, we compare the results from Ref. 21 with our hybrid model, which uses the same Boltzmann approximation for electrons and hydrodynamic solution for ions with $T0i=T0e=300$ eV and symmetric mirror $R=K$⁠. It is interesting to note that for a high mirror ratio $R=100$⁠, the general behavior of the ion temperature and heat flux are similar in the extended fluid model and hybrid calculations, Figs. 15 and 16. The parallel ion temperature and parallel heat fluxes are in good agreement in both models, Figs. 16(a) and 16(c). The perpendicular ion temperature and perpendicular heat fluxes are somewhat different. The notable difference is observed in the ion velocity profile, Fig. 15(b), while the plasma density and potential are fairly close, Figs. 15(a) and 15(c). It appears that the difference in the ion velocity originates from the large perpendicular ion temperature (in the region near the maximum magnetic field) in the fluid model, Fig. 15(b). The mirror force due to the large perpendicular ion velocity (temperature) provides additional ion acceleration in the expander region, resulting in higher final ion velocity.

In this work, we use WarpX and EDIPIC (with modifications as described above in Sec. II C) to compare their performance and determine the practical limits for the simulations with both codes. Explicit PIC simulations of realistic dimensions (a few meters in length) and densities of the order of $1018$ $m−3$ and larger are resource-intensive and may become impractical. While the current release of WarpX is explicit, EDIPIC³¹ employs the implicit algorithm.⁴⁸ Therefore, it is of interest to extend the study into the high-density regimes, exploring the advantages of the implicit EDIPIC code.

To compare the WarpX and EDIPIC results, we have performed the simulation of the symmetric mirror $R=K$ with different mirror ratios for $T0e=300 eV$⁠, $T0i=600 eV$⁠. These simulations are performed with $Δt=5×10−11 s−1$ and $Δz=2×10−3 m$⁠, which are two times larger compared to those for the base case. The simulation results for $R=2$ and $R=10$ perfectly agree as seen in Fig. 17. However, the results for $R=100$ for WarpX (not shown in Fig. 17) diverge due to the violation of the condition $Δt<0.2ωpe−1≃2×10−11 s$⁠. The latter condition is a standard requirement for the explicit PIC numerical stability.⁴⁹ To achieve the agreement between EDIPIC and WarpX for $R=100$⁠, as shown in Fig. 17, the time and space resolution in WarpX simulations were increased by a factor of four. The number of macroparticles was kept the same for a faster simulation. At the same time, the parameters for EDIPIC simulations remained the same for all values of R.

To further test the implicit algorithm of EDIPIC, we simulate different injection rates with the same spatial and time resolution as in base case $Δt=2.5×10−11 s−1$ and $Δz=10−3 m$⁠, Fig. 18. At the low injection flux of $Iinj=10$ $A/m2$⁠, the maximum density in the mirror (before the throat) is $n≃4.2×1015$ $m−3$⁠, thus giving $Δz/λDe≃1.9$ and $ωpeΔt≃0.09$⁠. The plasma density scales linearly with the injection flux, so at the largest injection flux we tested, $Iinj=5×104$ $A/m2$⁠, the plasma density reaches $n≃2.1×1019$ $m−3$⁠, Fig. 18(a). At this density, our simulation is performed with $Δz/λDe≃37$ and $ωpeΔt≃6.5$⁠.

One can see from Figs. 18(a) and 18(b), that with the increase in the injection flux, the density amplitude rescales while global profiles of plasma density, potential, and ion velocity remain unchanged. The sheath size at the right boundary (visible in the potential and ion velocity profiles) decreases, corresponding to the rescaled densities. However, as it is seen in Fig. 18(c), at the largest tested injection rate of $Iinj=5×104 A/m2$⁠, the irregular fluctuations in the potential profile appear. They are not smoothed out by our standard averaging procedure done over 300 snapshots uniformly spaced over 9 μs (the averaging window for the implicit case was increased for this set of simulations). Large fluctuations at $Iinj=5×104 A/m2$ are apparent (under accepted averaging parameters) in the potential profile but are not visible in the ion density nor in the ion velocity. Less apparent oscillations exist at the lower values of the current as well. The oscillations occur in the region of largest density, as Fig. 19(a) shows oscillations of the electric field for $Iinj=2×104 A/m2$ and Fig. 19 for $Iinj=5×104 A/m2$⁠. Figure 19 shows five consecutive snapshots of the electric field. The snapshots are averaged with a moving window of the width $100Δz$ to reduce noise. Figures 19(a) and 19(b) show that the amplitude of the oscillations (between the snapshots) and the width of the region where the oscillations are localized increases with the increase in plasma density. Figure 18(d) shows that the energy conservation starts to deteriorate for the currents larger than $Iinj=103 A/m2$⁠, corresponding to the density larger than $n≃4.2×1017$ $m−3$⁠.

Therefore, the EDIPIC implicit scheme allows us to overstep the Debye length and the electron plasma frequency conditions roughly by a factor of 10 but retain good energy conservation. These results are in agreement with general recommendations for the direct implicit scheme (with $D1$ spatial smoothing⁵⁰) that conserve energy as long as $ωpeΔt≃30$ for $Te/meΔt≃0.3Δz$⁠.⁵⁰ In our case, we have $Te/meΔt≃0.18Δz$⁠. It is interesting to note that even for considerably larger time steps and spatial grid sizes the global density, velocity, and potential profiles do not change and retain the sheath physics while the energy conservation may deteriorate, as seen in Fig. 18(d). Additionally, the behavior of the numerical heating/cooling is non-monotonical with changes in time and spatial steps. As one can see in Fig. 18(d), the energy conservation is the best for an intermediate density corresponding to the current $Iinj=103 A/m2$⁠, the green line in Fig. 18(d). This behavior is reminiscent of that described in Ref. 50. Simulations of this scale typically require 2–3 days on the Beluga (The Digital Research Alliance of Canada) cluster with 160 CPUs.

The important role of collisions on energy losses in open mirror systems has long been recognized and studied theoretically, through simulations, and experiments.^{1,28,29,38,50–52} A recent work⁵³ has considered the roles of intricate details of the collision operator on plasma confinement in mirror geometry. The comprehensive study of such effects requires a detailed analysis of the role of collisions both in the plasma source and the expander region, which is beyond the scope of our paper. Here, we address the role of collisions and the related effects of trapped electrons in the expander region on the plasma potential and electron pressure anisotropy. For our base case, with $T0i=600 eV$⁠, we include electron–neutral (e–n) elastic collisions with a uniform background of Hydrogen atoms across the entire length of the system. The cross section for the process was taken from Ref. 54. The collisions are modeled using the Monte Carlo null-collision method as perfectly elastic collisions. In simulations, the collision frequency is calculated based on the energy dependence of the cross section, thus becoming a function of the electron energy varying along the nozzle. To characterize different collisional regimes, we use $0.7T0e$ for the total electron temperature as an effective value which is equal to the average value in the fully collisional regime in the expander. The same $0.7T0e$ value is used to estimate an average bounce frequency, giving $νb=vthe/L=7.2×106 s−1$ for $L=1.2$ m. These parameters are summarized in Table I.

The main effect of collisions is the isotropization of the electron distribution function and the related trapping of electrons in the expander. In Fig. 20, electron distribution functions are shown for two locations and different collision frequencies. As expected, the number of trapped particles increases with the collision frequency, and the trapped region shrinks as a result of the concomitant modification of the potential profile.

The profiles of the parallel and perpendicular electron and ion temperatures for different collisionalities are shown in Fig. 22. The results in Fig. 21 show that collisions, despite being present throughout the whole system, do not affect the region before the throat of the magnetic nozzle; they are apparent only in the expansion part. The collisions modify the shape of the ion velocity (b) and potential profile (c) but not the total potential drop or the plasma density profile (a). The plasma potential becomes more Boltzmann-like (see Fig. 12). However, collisions nearly equalize the parallel and perpendicular temperatures for highly collisional regimes, $νen/νb$ $>0.1$⁠, which continuously decrease with the expansion and the decrease in the magnetic field, see Figs. 22(a) and 12(b). For large collisionalities, the potential variations shift away from the absorbing plate toward the maximum of the magnetic field (without changes in the total potential drop across the nozzle) with simultaneous isotropization and a decrease in the temperature. These results are generally in agreement with the results of Ref. 17.

The electron cooling is a result of plasma expansion, the related conversion of thermal energy into the kinetic energy of ions (via the electric field), and electron heat conductivity. Often, electron cooling is characterized by the polytropic constant $γ$⁠, for the isotropic pressure $p≃nγ$⁠. Such characterization is theoretically limited. Its validity depends on the details of the electron distribution function, e.g., degree of pressure isotropization, i.e., isotropic $p=(2p⊥+p∥)/3$ vs anisotropic $p⊥≠p∥$ pressure, the magnitude of heat fluxes and other possible factors, e.g., electric and magnetic fields.^46,55 The polytropic relation is one of the simplest assumptions of the equation of state for a general moments closure problem in weakly collisional regimes. Nevertheless, it is often used to characterize electron cooling in experiments and modeling.^17,56–60 In our simulations, the polytropic index varies along the nozzle as is shown in Fig. 23 and is not too sensitive to the collisionality (varies between $γ=1.2$ for the collisionless case and $γ=1.26$ for the highly collisional case near the nozzle exit). Experiments under various conditions report values from almost isothermal $γ≃1$ to adiabatic $γ=5/3$ and even larger, as discussed and summarized in Ref. 61. We would like to emphasize that the notion of the polytropic index is limited in scope; for example, the effects of pressure anisotropy are outside the polytropic equation of state.

We have presented the effective drift-kinetic PIC model to study plasma flow, acceleration, and energy transport in the geometry of the magnetic mirror. Despite of its limitations related to paraxial approximation, this model provides useful insights into the physics of plasma acceleration and axial energy losses in the magnetic mirror and magnetic nozzle.

We have employed implicit EDIPIC code to study the regimes approaching realistic plasma parameters with high density, realistic lengths, and electron/ion mass ratio. The results from EDIPIC and WarpX codes well agree in the explicit regime with small time steps and small grid size. The EDIPIC implicit algorithm allows the simulations to overstep the Debye length and electron plasma frequency conditions. We find that the energy is well conserved up to $Δz/λDe≃10$ and $ωpeΔt=6.5$⁠. The simulations with even larger time and spatial size remain stable and show the same “universal” profiles of plasma density, velocity, and potential, thus opening the possibilities for practical simulations of realistic dimensions and plasma densities. We have found that practical limits of the implicit algorithm exist roughly at $Δz/λDe≃37$ and $ωpeΔt≃6.5$ when the profile of averaged potential starts to depart from the “universal” profile.

We have shown directly in our kinetic model how the finite ion temperature in the plasma source increases plasma acceleration out of the mirror and axial energy losses. Full kinetic results, in particular, the heat flux, are compared with the results obtained in the hybrid formulation with kinetic ions, isothermal Boltzmann electrons, and quasineutrality. Furthermore, the ion kinetic results are compared with the results of the fluid closure model, where the collisionless heat fluxes are calculated from the extended (higher order, Grad type) hydrodynamic equations for the time evolution of the heat fluxes. It is shown that the kinetic results reasonably agree with the results from the extended hydrodynamic closure model. It is suggested that such a fluid closure model can also be useful to characterize the thermodynamic properties of the electron component.

We have studied the relative contributions of the heat fluxes to the energy transport and energy conversion in each component. Also, the formation and development of the anisotropies in the electron and ion distribution function, and related electron and ion pressure anisotropies are demonstrated at different locations along the mirror. The heat fluxes for the parallel and perpendicular energy (both for electron and ion components), responsible for the development of the pressure anisotropies, are explicitly determined here. These results show the role of the heat fluxes at different locations and can be used to describe thermodynamic properties^59,61,62 of the compressing and expanding flow of plasma in the mirror, i.e., the transitions between limiting cases of the isothermal and two-pressure adiabatic regimes that may exist at different regions in the mirror. We find that the heat flux magnitudes are largest in the region of the maximum magnetic field, Fig. 14.

The energy balance is investigated in detail, revealing the energy conversion between electrons and ions, ultimately defining the total axial energy losses in the mirror. The global integral relations are derived that relate the energy input/loss in each species to the total potential drop along the mirror, Eqs. (31) and (32). It is shown that the total potential drop across the mirror (including the sheath) saturates at large expansion and remains largely constant for large K. This behavior is different from the logarithmic dependence of the total potential drop in the hybrid model with isothermal electrons.