MCTrans++: a 0-D model for centrifugal mirrors

Abstract

The centrifugal mirror confinement scheme incorporates supersonic rotation of a plasma into a magnetic mirror device. This concept has been shown experimentally to drastically decrease parallel losses and increase plasma stability as compared with prior axisymmetric mirrors. MCTrans++ is a dimensionless (0-D) scoping tool which rapidly models experimental operating points in the Centrifugal Mirror Fusion Experiment (CMFX) at the University of Maryland. In the low-collisionality regime, parallel losses can be modelled analytically. A confining potential is set up that is partially ambipolar and partially centrifugal. Due to the stabilizing effects of flow shear, the perpendicular losses can be modelled as classical. Radiation losses such as bremsstrahlung and cyclotron emission are taken into account. A neutrals model is included, and, in some circumstances, charge-exchange losses are found to exceed all other loss mechanisms. We use the SUNDIALS ARKODE library to solve the underlying equations of this model; the resulting software is suitable for scanning large parameter spaces, and can also be used to model time-dependent phenomena such as a capacitive discharge. MCTrans++ has been used to verify results from prior centrifugal mirrors, create an experimental plan for CMFX and find configurations for future reactor-scale fusion devices.

1 Introduction

Axisymmetric mirror machines were at one time attractive as thermonuclear fusion devices because of their engineering simplicity, high beta and steady-state operating capability. However, these devices are plagued by poor axial confinement and the interchange instability (Post 1987). The centrifugal mirror confinement scheme, in contrast, incorporates supersonic rotation of a plasma into a conventional axisymmetric magnetic mirror device. Centrifugal confinement greatly reduces loss rates and is beneficial for removing impurities (Lehnert 1971). Other concepts for improved confinement and stabilization in a magnetic mirror (or ‘open trap’) are outlined in Ryutov et al. (2011).

Supersonic rotation has been demonstrated to improve axial confinement (Teodorescu et al. 2010) and stability (Huang & Hassam 2001), thereby reducing perpendicular losses too. This concept was first demonstrated with the Ixion device at Los Alamos National Lab (Baker & Hammel 1961; Baker, Hammel & Ribe 1961), and variations have been constructed at the Institute of Nuclear Physics in Russia (PSP-2, Volosov (2009)) and the University of Maryland (MCX, Ellis et al. (2005)), but none of these are currently in operation.

Most recently, construction has been completed and first plasma achieved at the Centrifugal Mirror Fusion Experiment (CMFX) at the University of Maryland (Romero-Talamas et al. 2022). This paper discusses the details and results from MCTrans++, a dimensionless (0-D) scoping tool which is primarily used to model experimental operating conditions in CMFX. Additionally, MCTrans++ can be used to predict the performance of reactor-scale centrifugal mirrors, as well as verify results from centrifugal mirror experiments mentioned previously.

There are significant engineering challenges to developing a fusion power plant based on the centrifugal mirror, but the aim of this paper is to demonstrate what is physically achievable if experimental concerns can be overcome. However, it is worth mentioning that two major challenges are impurity ions sputtered off plasma-facing surfaces and avoiding electrical breakdown from the necessary high voltages. Other research to tackle these issues is ongoing in parallel to this theoretical work.

The rest of this paper is structured as follows. In § 2, we survey the basis for the underlying physics, including formulae for parallel and perpendicular loss rates. In § 3, optional features of MCTrans++ are discussed, including time dependence, neutrals and radiation models and alpha heating. Results and discussion are covered in § 4 and conclusions in § 5. Appendix A is a literature review of shear flow stabilization, both in centrifugal mirrors and tokamaks. Appendix B presents a list of symbols used in this work.

2 Physics model

A diagram of CMFX presents the most important features, namely the magnetic field, central conductor, insulators and grounded chamber and ring electrodes (figure 1). The central conductor is biased by a large voltage from a capacitor bank, and a large radial electric field generates supersonic $\boldsymbol {E} \times \boldsymbol {B}$ flows. The potential drop from the centre electrode to the grounded ring electrodes is held fixed, and it has been shown previously that the velocity and temperature profiles across the radial width of the plasma are approximately parabolic (Huang & Hassam 2001; Romero-Talamás et al. 2012). We take the transverse length scale $a$ to be half the width, whereas the parallel length scale $L$ is the axial extent of the plasma.

Figure 1. Diagram of CMFX. The axial magnetic field is generated by superconducting magnets, and the electric field by directly biasing the central electrode and grounding the vacuum chamber. Supersonic rotation is achieved by the $\boldsymbol {E} \times \boldsymbol {B}$ velocity. The magnetic field terminates on electrically insulating surfaces so that the voltage can vary across the field lines.

The grounded ring electrodes are an improvement over the designs in other centrifugal mirrors, where the outer radius of the plasma was limited by the chamber wall (as in MCX and Ixion) or a liner coincident with the vacuum field (as in PSP-2). The ring electrodes present a smaller surface area for plasma–material interactions, thereby decreasing sputtered impurities into the core. The innermost and outermost flux surfaces are those which intersect the central conductor and ring electrode, respectively. We can therefore calculate the width of the plasma $2a$ as the distance between these flux surfaces at the midplane. Additionally, the simplest, and most transparent, geometric approximation of the magnetic field is the ‘square well’ model. In this model we assume the field lines to be straight, directed in the $z$-direction, with a step-function shape (see inset in figure 1).

Prior centrifugal mirrors, like PSP-2 (Abdrashitov et al. 1991), and current experiments, like GDT (Beklemishev et al. 2010) and TAE's C-2W (Binderbauer et al. 2015), use edge-biasing with concentric ring electrodes to establish a radial electric field profile in the plasma. Ryutov et al. (2011) provides a good theoretical understanding for stabilization via the combination of edge-biasing and ‘line-tying’. This is different from the CMFX approach, which has insulating end plates and a central conductor. Because the total potential drop and the boundary conditions are fixed, diffusive transport smooths the voltage into a singly peaked profile (Huang & Hassam 2001; Romero-Talamás et al. 2012). While experiments like GDT and C-2W do use some amount of biasing to create shear-flow stabilization, centrifugal mirrors impose much higher voltages to achieve improved confinement.

MCTrans++ solves simplified transport equations for the centrifugally confined plasma. The underlying model is discussed in the following sections, and it is primarily based on the assumptions that the plasma is well-confined, strongly magnetized and has low collisionality, all of which can be confirmed a posteriori.

2.1 Fundamentals

We first present the basic definitions needed to understand a centrifugal mirror. The flow velocity $\boldsymbol {u}$ is the same for all species and given by (Ellis et al. 2001)

(2.1)\begin{equation} \boldsymbol{u} = \omega R^2 \boldsymbol{\nabla}\phi, \end{equation}

where $\omega$ is the azimuthal angular velocity, $R$ is the radius and $\phi$ the azimuthal coordinate. The Mach number is defined as

(2.2)\begin{equation} M \equiv \frac{|\boldsymbol{u}|}{{c_{{s}}}} = \frac{\omega R_{\max}}{c_s} = \omega R_{\max}/ \sqrt{\frac{T_e}{m_i}}, \end{equation}

where we take $ {c_{{s}}} \equiv \sqrt {{T_e}/{m_i}}$, where $T_e$ and $m_i$ are the electron temperature and ion mass, respectively. We use the subscript ‘$\mathrm {max}$’ to denote the location of maximum radius in the square-well approximation (see figure 2). Note that, typically $T_e \sim T_i$ (Reid et al. 2014), but we take $ {c_{{s}}}$ as the cold ion limit of the sound speed, which turns out to be a convenient normalization factor that appears in the derivation of the confining potential (see § 2.2.1). Additionally, we assume that the magnetic field lies purely in the $R$–$z$ plane. In cylindrical coordinates $(R, \phi, z)$, the axisymmetric magnetic field can be expressed in terms of the poloidal flux $\psi$,

(2.3)\begin{equation} \boldsymbol{B} = \boldsymbol{\nabla}\psi\times\boldsymbol{\nabla} \phi, \end{equation}

where $\psi \equiv R A_\phi$, and $A_\phi$ is the azimuthal component of the vector potential. The electric field is dominantly perpendicular to the magnetic field, so that $\boldsymbol {E} \boldsymbol {\cdot } \boldsymbol {B} = 0$. To satisfy this condition, we introduce the electric field (also given by (12) in Abel et al. (2013)) as

(2.4)\begin{equation} \boldsymbol{E} ={-}\boldsymbol{\nabla} \psi \frac{{\rm d} \varPhi}{{\rm d}\psi} - \boldsymbol{\nabla} \varphi, \end{equation}

where $\varPhi$ is the part of the electrostatic potential associated with the plasma rotation (primarily the applied voltage) and $\varphi$ comprises all other pieces of the electrostatic potential.

Figure 2. Along a given field line, $R$ falls off along $z$. The blue $\boldsymbol {X}$ indicates the location of $R_\mathrm {max}$ (equivalently $B_\mathrm {min}$) and the red $\boldsymbol {X}$ the location of $R_\mathrm {min}$ (equivalently $B_\mathrm {max}$). This holds true for either the actual magnetic field or the square-well approximation. Thus, according to (2.9), the density along a flux surface will decrease exponentially as $|z|$ increases because $R$ decreases.

It is well known that (in the ideal, zero resistivity, limit) the field lines (or equivalently flux surfaces) in a plasma rotate as rigid bodies¹ (Ferraro 1937; Lehnert 1971). Mathematically, the angular frequency of any such surface $\omega$ obeys

(2.5)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\omega = 0, \end{equation}

and so $\omega = \omega (\psi )$ is a flux function. In strongly magnetized plasmas, any part of this flow perpendicular to $\boldsymbol {B}$ must be an $\boldsymbol {E}\times \boldsymbol {B}$ flow and so the angular frequency $\omega$ can be written in terms of the radial derivative of the electric potential $\varPhi$ (Lehnert 1971),

(2.6)\begin{equation} \omega ={-} \frac{{\rm d} \varPhi}{{\rm d}\psi}, \end{equation}

the radial electric field completely determining the flow profile and vice versa.

To estimate the size of $\varPhi$, we approximate the flow velocity as

(2.7)\begin{equation} |\boldsymbol{u}| \sim \frac{|E|}{|B|} \approx \frac{\varPhi}{a |B|} = \frac{e \rho_i \varPhi}{a m_i v_{th_i}}, \end{equation}

where $e$ is the electron charge, $\rho _i = {m_i v_{th_i}}/{e |B|}$ is the ion gyroradius and $v_{th_i} = \sqrt {{2 T_i}/{m_i}}$ is the ion thermal velocity. By assuming that the plasma is strongly magnetized (i.e. the ion gyroradius is small – ${\rho _i}/{a} \ll 1$) and our rotation is supersonic ($M \gg 1$), and substituting $M \equiv {|\boldsymbol {u}|}/{ {c_{{s}}}}$, we can compare the size of $\varPhi$ and $\varphi$, finding

(2.8)\begin{equation} \varPhi \sim M\frac{a}{\rho_i}\frac{T_e}{e} \gg \varphi \sim M^2 \frac{T_e}{e}, \end{equation}

where the Mach-number dependence of $\varphi$ comes from centrifugal effects (see (2.18)). In a centrifugal mirror, the dominant source of the potential $\varPhi$ is the voltage applied to the central conductor and $\varphi$ is determined by the fact that the plasma must remain quasineutral. For this reason we will often refer to $\varphi$ as the ambipolar potential.

If we assume that the plasma rotates at a large Mach number $M \gg 1$, and it has an ion temperature such that the plasma is fully ionized, then the plasma is well-confined and the confinement time is long compared with the collision time (this can be checked a posteriori). In such a situation, the plasma is locally Maxwellian (equivalently it is in local thermal equilibrium on a flux surface) and we can write the density $n$ of species $s$ as (see (96) in Catto, Bernstein & Tessarotto (1987), and Abel et al. (2013))

(2.9)\begin{equation} n_s = N_s (\psi) \exp\left( - \frac{\varXi_s}{T_s} \right) = N_s (\psi) \exp\left( - \frac{Z_s e \varphi}{T_s} + \frac{m_s \omega^2 R^2}{2 T_s}\right), \end{equation}

where $N_s$ is an arbitrary flux function with units of ${\rm m}^{-3}$, $\varXi _s$ is the potential energy of a particle (discussed further in § 2.2.1) and $Z_s$ is the charge number. We see that $n_s$ falls off exponentially along a field line because of its $R$-dependence (see figure 2). Because $\varXi _s$ is a potential, we can choose its zero to be anywhere. We therefore choose $\varXi _s = 0$ at the midplane ($z=0$, $R=R_{\mathrm {max}}$) so that $n_s = N_s$ at the midplane (see (2.9)). This choice simplifies solving for density because we can directly solve for the flux function $N_s$ at the midplane, and then easily compute $n_s$ along a field line.

MCTrans++ primarily solves the transport equations at the midplane, where $N_s = n_s$, (2.10)–(2.12). These have been derived from the low collisionality, strongly magnetized transport equations for an axisymmetric rotating plasma. Using the square-well approximation, there is no variation in $z$. The conservation of particles, energy and azimuthal angular momentum are thus,

(2.10)\begin{gather} {\frac{\partial N_s} {\partial t}} + \frac{1}{R} {\frac{\partial } {\partial R}} R{\varGamma}_s = S_{n,s}, \end{gather}

(2.11)\begin{gather}{\frac{\partial } {\partial t}} \left( \frac{3}{2} N_s T_s \right) + \frac{1}{R} {\frac{\partial } {\partial R}} R q_{s} = {{{\rm \pi}^{(R\phi)}_{s}}} {\frac{\partial \omega} {\partial R}} + Q_s + S_{E,s}, \end{gather}

(2.12)\begin{gather}{\frac{\partial } {\partial t}}(J\omega) + \frac{1}{R} {\frac{\partial } {\partial R}} R\sum_s {{{\rm \pi}^{(R\phi)}_{s}}} ={-} j_R R B + S_{\omega}. \end{gather}

In these equations, ${\varGamma }_s$ and ${q}_s$ are the radial particle and heat fluxes of each species, $J = \sum _s m_s N_s R^2$ is the moment of inertia of a flux surface at radius $R$, and ${ {{\rm \pi} ^{(R\phi )}_{s}}}$ is the radial flux of azimuthal angular momentum. MCTrans++ sets a constant density, but we leave $ {\frac {\partial N_s} {\partial t}}$ in this work to maintain generality. An external power source drives the radial current density $j_R$, which is how the plasma is spun up. The sources of heating in these equations come from viscous heating and the collisional equilibration between species (denoted by $Q_s$). We have also included arbitrary sources of particles $S_{n,s}$, energy $S_{E,s}$ and angular momentum $S_{\omega,s}$; these source terms will be used to account for additional energy sources such as alpha-particle heating. Currently, MCTrans++ models a particle source that is feedback-controlled to maintain the electron density at a fixed value. We use the SUNDIALS ARKODE library to solve these conservation equations via a fourth-order implicit Runge–Kutta scheme.

We often wish to consider the behaviour of a centrifugal plasma given a fixed input voltage (and thus fixed $\varPhi$). In this case, the momentum transport equation is used to determine the radial current drawn from the power supply, $I_R = 2{\rm \pi} R L j_R$.² We never need to solve (2.12) for $\omega$ under this fixed input. Instead, if we hold the input power or current fixed then (2.12) would be solved to find $\omega$.

In steady-state operation, MCTrans++ solves the time-dependent equations with internal time steps until the solution converges to balance heat generation and heat losses. Subsection 3.1 discusses the time-dependent capabilities. In the following subsections we will explain the approximations used to calculate those heat losses explicitly in terms of the system state.

2.2 Parallel transport

We begin with a discussion on parallel transport of particles, heat and momentum to feed into their respective equations, (2.10)–(2.12). The confining potential in the parallel direction primarily arises from the centrifugal potential, with a small (but not negligible) component from the electrostatic potential. In the following section, we derive said potential, $\varXi _s$, and then determine parallel losses in § 2.3.

2.2.1 Centrifugal potential

Contributions to the potential energy of a charged particle on a rotating field are two-fold: (i) the electrostatic potential $Z_s e \varphi$ and (ii) the centrifugal potential, which is given by $m_s \omega ^2 R^2 / 2$ for a particle of mass $m_s$ at radius $R$ rotating with velocity $\omega$. The potential energy of a particle of species $s$ is then

(2.13)\begin{equation} \varXi_s = Z_s e \varphi - \frac{m_s}{2} \omega^2 R^2. \end{equation}

The first term in this potential is confining (i.e. negative) for electrons, but deconfining for ions, whereas the second centrifugal term is confining for all species. We can reuse results from previous work (notably Pastukhov (1974)) which consider only electrostatic confinement, simply by making the substitution $Z_s e \varphi \rightarrow \varXi _s$. Physically, the first term serves to keep the electrons (which are light and barely affected by the centrifugal force) next to the ions, which are pushed to regions of large $R$ by the centrifugal force. Hence, the potential $\varphi$ can be found by insisting that the plasma is quasineutral along field lines and that the loss rate is ambipolar (Post 1961).

The ambipolar potential can be found by breaking it down into zeroth- and first-order components, $\varphi = \varphi _0 + \varphi _1$. The leading-order term is found by imposing quasineutrality assuming no losses, while the second term (found in § 2.3) is determined by insisting that the losses preserve quasineutrality.

If we now insist that the plasma is made up of ions (mass $m_i$ and charge $Z_i e$) and electrons (mass $m_e$ and charge $-e$), then the condition for quasineutrality is found by equating (2.9) for ions and electrons

(2.14)\begin{equation} Z_i N_i(\psi) \exp\left( - \frac{Z_i e\varphi_0}{T_i} + \frac{m_i}{2T_i} \omega^2 R^2 \right) = N_e(\psi) \exp\left( \frac{e \varphi_0}{T_e} + \frac{m_e}{2T_e} \omega^2 R^2 \right). \end{equation}

We eliminate $N_s$ by noting that $Z_i N_i = N_e$ on a given flux surface $\psi$. So, up to a possible constant offset, we have

(2.15)\begin{equation} \varphi_0 = \left( \frac{Z_i e}{T_i} + \frac{e}{T_e}\right)^{{-}1} \left(\frac{m_i}{2 T_i} - \frac{m_e}{2T_e}\right) \omega^2 R^2 \approx \left( \frac{Z_i e}{T_i} + \frac{e}{T_e}\right)^{{-}1} \frac{m_i}{2 T_i}\omega^2 R^2, \end{equation}

where we assume the temperatures are such that we can drop the term proportional to the electron mass in all further calculations. We have denoted this potential $\varphi _0$ because it is $O(M^2 T_e/e)$ and is, in fact, the leading-order term in an asymptotic series of $\varphi$ in $M^{-1}$.³ We end with a convenient expression for the potential drop from the centre of a flux surface (at $R = R_{\max }$) in terms of suitably normalized variables,

(2.16)\begin{equation} \frac{e \varphi_0}{T_e} = (Z_i + \tau)^{{-}1} \left(\frac{R^2}{R_{\mathrm{max}}^2} - 1\right) \frac{M^2}{2} + O(1), \end{equation}

where $\tau = T_i/T_e$ is the temperature ratio and we have changed the zero of $\varphi _0$ so that it vanishes on the midplane. As a consequence of flux conservation, one can approximate the ratio of radii in terms of the mirror ratio as follows:

(2.17)\begin{equation} \frac{1}{R_{\mathrm{mirror}}} \equiv \frac{B_{\mathrm{min}}}{B_{\mathrm{max}}} \approx \left(\frac{R_{\mathrm{min}}}{R_{\mathrm{max}}}\right)^2, \end{equation}

and so we can relate the ratio of the radius of the flux surface at throat, $R_\mathrm {min}$, and in the central cell, $R_\mathrm {max}$, to the mirror ratio, $R_{\mathrm {mirror}}$, given by the ratio of magnetic field strengths. This approximation is satisfactory for the vacuum field, but high speed rotation creates a self-consistent field that provides better confinement (Abel 2022).

In the simplest form of (2.16), where $Z_i = 1$, $T_i = T_e$ and $R = R_{\mathrm {min}}$, we end up with

(2.18)\begin{equation} \frac{e \varphi_0}{T_e} = \left(\frac{1}{R_{\mathrm{mirror}}} - 1 \right) \frac{M^2}{4} + O(1), \end{equation}

with the usual $M^2/4$ scaling for the potential drop (Ellis et al. 2001). We leave $Z_i$ as a free parameter elsewhere, but use $Z_i = 1$ simply to show agreement with Ellis et al. (2001). We will see in § 2.3 that this scaling is very important for minimizing parallel losses because it appears in an exponential term for the loss rate.

Equation (2.16) is only the leading-order term in the $M \gg 1$ expansion of $\varphi$, and so we collect all higher-order terms and denote them by $\varphi _1$. To compute the next-order terms in this series we need to know the particle loss rate and hence parallel transport (discussed in § 2.3). For electrons, the only term in the potential energy is the electrostatic potential,

(2.19)\begin{equation} \varXi_e = \frac{1}{Z_i + \tau} \left( 1 - \frac{R^2}{R_{\mathrm{max}}^2} \right) \frac{M^2}{2} T_e - e \varphi_1, \end{equation}

but for ions we need to include the centrifugal potential to obtain

(2.20)\begin{align} \varXi_i & = \frac{Z_i}{Z_i + \tau} \left( \frac{R^2}{R_{\mathrm{max}}^2} - 1 \right) \frac{M^2}{2} T_e + Z_i e \varphi_1 - \frac{m_i}{2} \omega^2 (R^2 - R_{\mathrm{max}}^2) \nonumber\\ & = \left[ \frac{Z_i}{Z_i + \tau} \left( \frac{R^2}{R_{\mathrm{max}}^2} - 1 \right) - \left( \frac{R^2}{R_{\mathrm{max}}^2} - 1 \right) \right] \frac{M^2}{2} T_e + Z_i e \varphi_1 \nonumber\\ & = \frac{\tau}{Z_i + \tau} \left( 1 - \frac{R^2}{R_{\mathrm{max}}^2} \right) \frac{M^2}{2} T_e + Z_i e \varphi_1 . \end{align}

Again, $\varphi _1$ is the higher-order part of $\varphi$ which comes from enforcing ambipolar losses, discussed in the following section.

2.3 Parallel losses

We assume that the time between particle collisions is long compared with the time a particle takes to travel along the mirror machine. This ‘low-collisionality’ assumption is true if the plasma is sufficiently hot with sufficiently low density. For reactor-grade plasmas this ratio, called $\nu ^* = \nu _{ii} L / { {v_{\mathrm {th}_s}}}$, can be as small as $1\times 10^{-5}$, but the assumption holds even in warm plasmas with temperatures of only 100 eV. The collisionality parameter $\nu ^*$ is an output of the code and validates the assumption a posteriori if it is significantly less than one.

The parallel losses are derived from Pastukov's work on low-collisionality plasmas taken in the case of a tandem mirror with an electron-confining electrostatic potential (Pastukhov 1974). Although originally found for electrons alone, we expand the results for a multispecies plasma. Catto (1981) found the parallel loss rates for a generic field shape ((40) of that work), and by taking the square-well limit, we obtain the particle loss rate in our notation⁴ ,

(2.21)\begin{equation} \left({\frac{\partial n_s} {\partial t}}\right)_{\text{End Losses}} ={-} \left( \frac{2 n_s \varSigma}{\sqrt{\rm \pi}} \right) \nu_s \frac{1}{\ln(R_{\mathrm{mirror}} \varSigma)} \frac{\exp(-\varXi_s/T_s)}{\varXi_s/T_s}, \end{equation}

where $\varSigma = Z_i + 1$ for electrons and $\varSigma = 1$ for ions, and $\nu _s$ is the appropriate collision frequency for species $s$ (see (2.5) in Braginskii (1965)).

The parallel heat losses can be calculated by multiplying the loss rate in (2.21) by the energy of a single particle. We approximate the total energy as the thermal and potential energy, $T_s + \varXi _s$. If $M \gg 1$, $\varXi _s \gg T_s$ following similar logic to (2.8).

Similarly, the parallel loss contribution to the total azimuthal angular momentum losses are due to the angular momentum of each particle $m_i \omega R^2$ when it is lost (with electrons carrying negligible momentum). As will be seen in § 4.2.2, this loss mechanism can be dominant. Because flux surfaces rotate rigidly, a particle will have less angular momentum when it is lost farther away from the midplane (i.e. at lower values of $R$, see figure 2). The choice of the exhaust radius, $R_\mathrm {exh}$, for MCTrans++ is therefore critical – the pessimistic assumption would be to assume that the ion is lost with angular momentum at the midplane ($R = R_\mathrm {max}$); whereas the optimistic assumption would be at the throat ($R = R_\mathrm {min}$). In reality, the momentum is lost when the density has dropped to the point where electrons can no longer shield parallel electric fields, i.e. $B \boldsymbol {\cdot }\boldsymbol {\nabla } \omega \neq 0$⁵. For this study, we have chosen $R_\mathrm {exh} = R_\mathrm {min}$, but § 4.2.2 discusses how varying $R_\mathrm {exh}$ changes results.

In (2.21) we see that $\varXi _s$ appears in the exponential. The leading-order part of $\varXi _s$ is $O(M^2 T_s)$, so this exponential is what strongly suppresses the collisional loss rate. Expanding this exponential, we have

(2.22)\begin{equation} \exp\left(-\frac{\varXi_s}{T_s}\right) \sim \exp[-(\dots) M^2] \exp\left(-\frac{Z_s e \varphi_1}{T_s}\right), \end{equation}

where $(\dots )$ represents the prefactor in (2.19) and (2.20) (which are independent of Mach number), and we see that even though $\varphi _1$ is small compared with the leading-order part of the potential it has an $O(1)$ effect on the loss rate and must be taken into account. MCTrans++ does not separate $\varphi _0$ and $\varphi _1$, and instead it solves the nonlinear equation for quasineutrality for $\varphi _1$ with the initial guess as $\varphi _0$, and $\varphi =\varphi _0+\varphi _1$.

We find $\varphi _1$ through a root-finding method by equating the electron and ion loss rates along the field line to enforce zero net charge loss. However, at low temperatures ($\lesssim$50 eV), (2.21) becomes very sensitive to changes in $T_s$, and produces poor confinement. MCTrans++ may be unable to find an equilibrium for these cases because quasineutrality cannot be satisfied⁶ . However, solutions may exist $\lesssim$50 eV and should be checked a posteriori for low collisionality ($\nu ^* = \nu _{ii} L_\parallel / {c_{{s}}} \ll 1$). Note that all plots in § 4 were checked and have low collisionality.

2.4 Perpendicular transport

We now proceed to determine particle, heat and momentum losses in the perpendicular direction. We make the assumption that turbulent transport will be fully suppressed by the flow shear (Huang & Hassam 2001), given that the velocity is everywhere perpendicular to the magnetic flux surfaces. A discussion of the literature on the suppression of turbulence by flow shear is given in appendix A. We thus assume that the only contributions to these fluxes are the classical collisional fluxes. These can be evaluated from the formulae in Braginskii (1965) or Helander & Sigmar (2002).

We begin by considering particle transport. Assuming that the plasma contains no impurities and has $Z_i = 1$, quasineutrality implies that the particle flux of ions must match that of electrons. Therefore (c.f. (5.6) of Helander & Sigmar (2002))

(2.23)\begin{equation} \varGamma_i = \varGamma_e ={-}D_e N_e \left[ \left( 1 + \frac{T_i}{Z T_e} \right) \frac{1}{N_e} {\frac{\partial N_e} {\partial R}} - \frac{1}{2T_e} {\frac{\partial T_e} {\partial R}} + \frac{1}{Z T_e} {\frac{\partial T_i} {\partial R}} \right], \end{equation}

where the classical electron diffusion rate $D_e$ is given by

(2.24)\begin{equation} D_e = \frac{T_e}{m_e \varOmega_{e}^2\tau_e}, \end{equation}

and $\varOmega _{e}$ is the electron cyclotron frequency, and $\tau _e$ is the electron–electron collision time.

We take the ion and electron collision times from Braginskii (1965) and convert them to SI units:

(2.25a,b)\begin{equation} \tau_i = \frac{6 \sqrt{2 m_i} {\rm \pi}^{3/2} T_i^{3/2} \epsilon_0^2}{n_i Z^4 e^4 \lambda_i};\quad \tau_e = \frac{6 \sqrt{2 m_e} {\rm \pi}^{3/2} T_e^{3/2} \epsilon_0^2}{n_i Z^2 e^4 \lambda_e}, \end{equation}

where $\lambda _i$ and $\lambda _e$ are the coulomb logarithms for ions and electrons, taken from Richardson (2019). These definitions differ by $\sqrt {2}$ from some other definitions of $\tau _s$.

To completely evaluate the transport equations we need expressions for $q_{s}$ and ${ {{\rm \pi} ^{(R\phi )}_{s}}}$. Similar to the particle flux, we only need to consider the classical collisional heat fluxes. The appropriate form of these fluxes is given by Braginskii (1965). For the ion heat flux, we use (2.14) and (2.16) from that work,

(2.26)\begin{equation} q_{i} ={-} 2 \frac{n_i T_i}{m_i \varOmega_{i}^2 \tau_i} \frac{{\rm d} T_i}{{\rm d} R}, \end{equation}

under the simplifying assumption that $\varOmega _{i} \tau _i \gg 1$ (as appropriate for our plasmas). As we do not evolve the radial temperature profile, we estimate this term by assuming that all the plasma profiles vary on a scale $a$, the half-width of the plasma (figure 1). Thus, we have that

(2.27)\begin{equation} \frac{1}{R} {\frac{\partial } {\partial R}} R q_{i} \approx 2\frac{n_i T_i^2}{m_i \varOmega_{i}^2 \tau_i a^2}. \end{equation}

For electrons, (2.13) from Braginskii (1965) gives

(2.28)\begin{equation} \frac{1}{R} {\frac{\partial } {\partial R}} R q_{e} \approx \frac{4.66 n_e T_e^2}{m_e \varOmega_{e}^2 \tau_e a^2}, \end{equation}

which is approximately $\sqrt {m_e/m_i}$ smaller than the ion heat loss, and usually small (though it is included in MCTrans++, nonetheless). The numerical coefficient in (2.28) is $Z_i$-dependent and takes the value given for $Z_i = 1$ (for the $Z_i$-dependence of this coefficient, see Braginskii (1965), Table 1). We leave $Z_i$ as a free parameter elsewhere. At the time of writing, $Z_i$-dependence of this coefficient is not implemented in MCTrans++. To these conductive heat fluxes, we add $(3/2)n_s T_s \varGamma _s$ to account for the non-zero particle flux.

To find the radial flux of angular momentum, ${ {{\rm \pi} ^{(R\phi )}_{s}}}$, we only compute the stress tensor for the ions because electron perpendicular viscosity is at least $m_e/m_i$ smaller than the ion viscosity and thus negligible. Similarly, convective transport of angular momentum by the ions is small compared with their viscous stress. Under these assumptions, from (2.23) of Braginskii (1965), we obtain

(2.29)\begin{equation} {{{\rm \pi}^{(R\phi)}_{i}}} = \frac{{\rm \pi}^{(\psi\phi)}_{i}}{BR} ={-} \frac{3}{10} \frac{n_i T_i}{\varOmega_{i}^2 \tau_i} R^2\frac{{\rm d}\omega}{{\rm d} R} \approx{-} \frac{3}{10} \frac{n_i T_i}{\varOmega_{i}^2 \tau_i} \frac{R^2\omega}{a}. \end{equation}

2.5 Assumptions and scope

Our model is valid when certain assumptions are met, some of which have been discussed already. The plasma must

1. (i) be strongly magnetized,

2. (ii) have low collisionality,

3. (iii) have a large applied bias in comparison with the ambipolar potential,

4. (iv) have a large mirror ratio,

5. (v) rotate supersonically.

Low-collisionality plasmas are required in MCTrans++ because in the collisional regime, it can no longer be assumed that the temperature is constant on a flux surface. Additionally, the calculations of parallel transport by Pastukhov (1974) and Catto (1981) are not valid for collisional plasmas. This applies to experiments like Ixion and MCX, or to the startup phase when the temperature is low and $\nu ^* \gtrsim 1$. Work is in progress to implement a collisional approach. The values to determine these operating conditions for a variety of devices are given in table 1, along with a other relevant plasma and experimental parameters. Details of prior experiments are given in § 4.1.1, and the results from MCTrans++ in § 4.2.

Table 1. Parameters (and their corresponding assumptions made in MCTrans++, if applicable) for prior experiments and projected conditions for CMFX and a reactor scenario. Also included are relevant plasma and experimental parameters for each experiment. Results for Ixion (Baker et al. 1961), PSP-2 (Volosov 2009) and MCX (Teodorescu et al. 2008; Reid et al. 2014) are from prior experiments, respectively, whereas those for CMFX and Reactor came from MCTrans++ predictive models.

3 Features

Our model contains several features, which can be turned ‘on’ or ‘off’ in MCTrans++. Typically MCTrans++ is operated in a steady-state mode, but (2.10)–(2.12) can be solved in a time-dependent mode. Additionally, models of neutral particles, continuum radiation sources and alpha heating are provided.

3.1 Time dependence

The system can be modelled as a circuit, where the plasma is charged by a capacitor bank and can be discharged through a crowbar into a dump resistor (figure 3). All of the passive circuit elements are static, but the plasma can be thought of as a variable resistor and capacitor in parallel.

Figure 3. Circuit model of CMFX. Here $C_{{\rm cap}}$ and $R_{{\rm cap}}$ are the capacitance and internal resistance of the capacitor bank, respectively; $R_l$ and $L_l$ are the line resistance and inductance, respectively. The plasma can be modelled as a variable resistor and capacitor in parallel. The dump resistor is in series to the plasma, and when the crowbar is switched, it is assumed all the stored energy from the plasma is transmitted through $R_{{\rm dump}}$ to ground.

The plasma has a voltage-dependent resistance because as it heats up, the current needed to rotate decreases, and thus the effective resistance increases. This is perhaps counter-intuitive as plasma resistance normally decreases with temperature in other devices like tokamaks; but in this case, the resistance increases with velocity shear because the potential gradient between neighbouring flux surfaces increases. Additionally, the plasma can be thought of as a capacitor that stores energy in its rotational momentum (Anderson et al. 1959).

Discharges in CMFX proceed as follows.

1. (i) The capacitor bank is charged to some nominal voltage.

2. (ii) Neutral gas is puffed into the chamber.

3. (iii) After some specified time, S1 closes and the capacitor bank discharges, applying high voltage to the central conductor.

4. (iv) A low-temperature plasma forms. The voltage on the capacitors drops because the current draw in this phase is relatively large.

5. (v) Rotational shear begins to heat up the plasma, and the voltage across the plasma reaches a quasi-steady-state with low current draw.

6. (vi) The circuit is crowbarred by S2, and stored energy in both the capacitor bank and plasma are discharged into a dump resistor.

Our model differs from discharges in CMFX in two ways. First, we assume that the electron density is constant, when really the density is time-dependent, especially as the plasma heats up. Second, a starting voltage must be specified, and it should be sufficiently large enough so that the plasma is not highly collisional (effectively skipping step (iv)). However, MCTrans++ does have the ability to model the crowbar sequence (step (vi)).

The high voltage circuit can be modelled with simple elements as in figure 3. The equations for voltage and current across the plasma are as follows:

(3.1)\begin{gather} \frac{\mathrm{d}V}{\mathrm{d}t} ={-} \frac{I_\mathrm{cap}}{C_\mathrm{cap}} - \frac{V_\mathrm{cap}}{R_\mathrm{cap} C_\mathrm{cap}} , \end{gather}

(3.2)\begin{gather}\frac{\mathrm{d}I}{\mathrm{d}t} = \frac{V_\mathrm{cap} - V - R_l I}{L_l}. \end{gather}

As a general rule of thumb, discharges with higher bank capacitance are able to sustain higher voltages across the plasma. As will be seen in § 4.2, higher plasma voltages almost always produce better performing plasmas.

3.2 Neutrals model

Ng & Hassam (2007) studied neutral penetration into centrifugal mirrors along the axis, finding that the neutral density drops exponentially along the field lines with good centrifugal confinement. However, even a neutral density that is orders of magnitude smaller than the plasma density can have a large effect on power losses (see § 4.2.2), so neutrals cannot be ignored. This section describes the basic model used to determine the neutral density.

To maintain a constant plasma density, neutrals must be supplied to the plasma at the same rate electrons are lost. We assume that a gas puff system provides an ambient source of cold neutrals. To calculate neutral density, we divide the electron loss rate (the sum of parallel and perpendicular losses) by the total ionization rate. Charge exchange is another important loss mechanism for ions. Therefore, we consider three neutral collisional processes: ion- and electron-impact ionization, and charge exchange.

Cross-sections for these collisions come from Janev & Smith (1993). Radiative recombination is negligible in the temperature range of interest. The plasma is assumed to be in coronal equilibrium (i.e. the atomic excitation frequency is much smaller than the de-excitation frequency), so excited states are not considered (Drawin & Emard 1976; Tallents 2018). We also do not consider wall recycling because, although it may decrease the necessary neutral source rate, it does not affect the steady-state neutral density in a 0-D model like MCTrans++.

We must first calculate the collision rates between neutrals $n$ and some species $s$ per unit volume, given by

(3.3)\begin{equation} R_{ns} (\boldsymbol{v_s}, \boldsymbol{v_n}) = n_s n_n v_r \sigma(v_r) f_s(\boldsymbol{v_s}) f_n(\boldsymbol{v_n})\,\mathrm{d}^3\,\boldsymbol{v_s}\,\mathrm{d}^3\,\boldsymbol{v_n}, \end{equation}

where $n_s$ and $n_n$ are the density of charged species and neutrals, respectively, and $v_r$ is relative velocity between species, $f_s$ is an arbitrary distribution function (normalized such that $\int f_s(\boldsymbol {v_s})\,\mathrm {d}^3\,\boldsymbol {v_s} = 1$, as is standard in atomic physics) and $\sigma$ is a collision cross-section. We choose a Maxwellian distribution for a rotating plasma with a bulk fluid velocity of $|\boldsymbol {u}| \equiv M {c_{{s}}} = \omega R_{\max }$ such that

(3.4)\begin{equation} f_s(\boldsymbol{v_s}) = \left( \frac{m_s}{2 {\rm \pi}T_s} \right)^{{3}/{2}} \exp \left( -\frac{m_s (\boldsymbol{v_s} - \boldsymbol{u})^2}{2 T_s} \right), \end{equation}

and the distribution function for the cold neutrals is

(3.5)\begin{equation} f_n(\boldsymbol{v_n}) = \delta( \boldsymbol{v_n} ), \end{equation}

where $\delta$ is the Dirac delta function. We then transform into spherical coordinates (for the sake of simplifying the integral, as is done in Appelbe & Chittenden (2011)) where the Jacobian is $\mathrm {d}^3 \boldsymbol {v_s} = v_s^2 \sin \theta _s \,\mathrm {d} v_s \,\mathrm {d} \theta _s\, \mathrm {d} \phi _s$. In order to perform the integration, we choose a coordinate system that is local to a single particle and align the $z$-axis of the transform with the fluid velocity. To be clear, this choice of coordinate system is not the global cylindrical coordinate system, where fluid flow is in the azimuthal direction,

(3.6)\begin{equation} \left.\begin{gathered} v_{sx} = v_s \sin \theta_s \cos \phi_s,\quad u_x = 0, \\ v_{sy} = v_s \sin \theta_s \sin \phi_s,\quad u_y = 0, \\ v_{sz} = v_s \cos \theta_s,\quad u_z = u, \end{gathered}\right\} \end{equation}

so that the integrand now becomes

(3.7)\begin{align} R_{ns} (v_s) & = n_s n_n \left( \frac{m_s}{2 {\rm \pi}T_s} \right)^{{3}/{2}} v_s^3 \sigma(v_s) \sin \theta_s \nonumber\\ & \quad \times \exp \left[ -\frac{m_s (v_s^2 + u^2 - 2 v_s u \cos \theta_s)}{2 T_s} \right] \delta( \boldsymbol{v_n} )\,\mathrm{d} v_s \,\mathrm{d} \theta_s \,\mathrm{d} \phi_s\, \mathrm{d}^3 \boldsymbol{v_n}. \end{align}

Only the delta distribution is a function of $\boldsymbol {v_n}$, so that term integrates out to 1. Completing the integral and taking the thermal velocity to be ${ {v_{\mathrm {th}_s}}} = \sqrt {2 T_s / m_s}$, we have

(3.8)\begin{align} R_{ns} & = \frac{2 n_s n_n}{u {{v_{\mathrm{th}_s}}} \sqrt{\rm \pi}} \exp \left( -\frac{u^2}{{{v_{\mathrm{th}_s}}}^2} \right) \nonumber\\ & \quad \times \int_0^\infty v_s^2 \sigma(v_s) \sinh \left( \frac{2 v_s u}{{{v_{\mathrm{th}_s}}}^2} \right) \exp \left( -\frac{v_s^2}{{{v_{\mathrm{th}_s}}}^2} \right) \mathrm{d} v_s. \end{align}

To simplify, we expand the $\sinh (\cdots )$ term and define a thermal Mach number $M_{th_s} = u / { {v_{\mathrm {th}_s}}} = M {c_{{s}}} / { {v_{\mathrm {th}_s}}}$ such that

(3.9)\begin{align} R_{ns} & = \frac{n_s n_n}{M_{th_s} {{v_{\mathrm{th}_s}}}^2 \sqrt{\rm \pi}} \int_0^\infty v_s^2 \sigma(v_s) \nonumber\\ & \quad \times \left\lbrace \exp \left[ -\left( M_{th_s} - \frac{v_s}{{{v_{\mathrm{th}_s}}}} \right)^2 \right] - \exp \left[ -\left( M_{th_s} + \frac{v_s}{{{v_{\mathrm{th}_s}}}} \right)^2 \right] \right\rbrace \mathrm{d} v_s. \end{align}

The computed rate coefficients for an arbitrary value of $M=4$ demonstrates the important role that rotation plays in decreasing charge exchange and increasing proton-impact ionization (figure 4).

Figure 4. Rate coefficients for a number of collisions involving neutrals. Solid lines are for a non-rotating plasma, while dashed lines are for a plasma with $M=4$. The electron-impact ionization rate is not affected by rotation, and in the limit $M \rightarrow 0$, the dashed lines equal the solid.

To give context to the collision rates with supersonic rotation, consider that prior rotating mirror experiments (like MCX Ellis et al. (2005)) produced plasmas with temperatures ${\sim }10^2$ eV, where the energy losses due to charge exchange actually increase. The target operating temperature of CMFX is ${\sim }10^3$ eV, where the charge exchange rate is roughly equal for both the rotating and non-rotating plasmas. Lastly, for reactor-scale rotating mirrors ${\sim }10^4$ eV, charge exchange losses are predicted to decrease in a rotating plasma by several orders of magnitude as evidenced in figure 4. Moreover, proton impact ionization increases by several orders of magnitude for all temperatures.

We pessimistically assume that when a charge exchange occurs, it produces a hot neutral at the rotational energy of the plasma that immediately exits the plasma. The mean free path of a hot neutral, $\lambda _{n^*}$, is given by

(3.10)\begin{equation} \lambda_{n^*} = \frac{v_{n^*}}{n_e \langle \sigma v \rangle_{n^*}} \approx \frac{|\boldsymbol{u}|}{n_e \langle \sigma v \rangle_{n^*}} = \frac{M \sqrt{\dfrac{T_e}{m_i}}}{n_e \langle \sigma v \rangle_{n^*}}, \end{equation}

where $v_{n^*}$ is the velocity of a hot neutral, which we assume comes from the kinetic energy of the ions (usually much greater than the thermal energy for rapidly rotating plasmas). For CMFX- and reactor-relevant plasmas $\lambda _{n^*} \gg a$, so the prompt loss assumption is typically valid. This model does also assume that the neutrals are supplied by an ambient gas puff source, whereas a method like neutral beam injection may be able to decrease charge exchange losses.

In fact, charge exchange can be the dominant mechanism for heat and momentum loss, especially at lower temperatures (see § 4.2.2). An increase in Mach number drastically decreases the charge exchange loss rate at a given temperature, so faster rotation is paramount to decreasing power draw.

3.3 Radiative losses

The MCTrans++ model includes continuum radiation from bremsstrahlung and cyclotron emission. Bremsstrahlung was modelled with the ‘lumped impurity’ assumption and is a function of $Z_{\mathrm {eff}}$, the effective charge of the plasma. However, we assume that the impurity ions do not dilute the main ion species, so quasineutrality is enforced by only considering the main ions and electrons.

Additionally, we assume that the plasma is in coronal equilibrium, meaning that the de-excitation frequency is much larger than the excitation frequency. Therefore, we only consider species in the ground state.

We have assumed that line radiation is negligible, which is true for a sufficiently hot and pure plasma. Future iterations of MCTrans++ may include this effect.

3.3.1 Bremsstrahlung and cyclotron radiation

Taking the formula for bremsstrahlung radiation from the NRL formulary ((62) in Richardson (2019)) and writing it in convenient units we have

(3.11)\begin{equation} \dot{Q}_{\mathrm{Brem}} = 5.34 \times 10^3 Z_{\mathrm{eff}} \left(\frac{n_e}{10^{20} \mathrm{m}^{{-}3}}\right)^2 \left(\frac{T_e}{1 \mathrm{keV}}\right)^{1/2}~\mathrm{W}\,\mathrm{m}^{{-}3}, \end{equation}

where we have used the definition

(3.12)\begin{equation} Z_{\mathrm{eff}} = \frac{1}{n_e} \sum_{s} Z_s^2 n_s \end{equation}

of the effective charge state of the plasma (summation taken over all positively charged species). Similarly, the cyclotron radiation by an electron in vacuum (from the formulary Richardson (2019)) is (in SI units)

(3.13)\begin{equation} \dot{Q}_{\mathrm{cyc,vac}} = 6.21\times 10^{{-}17} B^2 n_e T_e\,\mathrm{W}\,\mathrm{m}^{{-}3}. \end{equation}

The transparency of the plasma to cyclotron emission is calculated from the formulae in Tamor (1983). We have determined that all plasmas considered here are opaque to cyclotron emission and thus reabsorb most of the radiation. There are some amount of losses that occur at the surface of the plasma, and these are accounted for by assuming a reflection coefficient from the the vacuum vessel of $R= 95$% (modern vacuum vessels will typically exceed this reflectivity). The total power loss is $\dot {Q}_{\mathrm {cyc}} = k \dot {Q}_{\mathrm {cyc,vac}}$⁷, where

(3.14)\begin{gather} k = \frac{T_e^{3/2}}{200 l^{1/2}}, \end{gather}

(3.15)\begin{gather}l = \frac{2 a}{l_0 (1 - R)}, \end{gather}

(3.16)\begin{gather}l_0 = 1.66\times 10^{16} \frac{B}{n_e}\,\mathrm{m}; \end{gather}

and all quantities are in SI units.

3.4 Alpha particles

Our model for alpha particles is that they are all born with a delta-function distribution at the birth energy of $E_\alpha = 3.52$ MeV,

(3.17)\begin{equation} \left( {\frac{\partial f_\alpha} {\partial t}}\right)_{\mathrm{Source}} = \frac{S_\alpha \delta(v-v_*)}{4{\rm \pi} v_*^2}, \end{equation}

where $v_*$ is the birth velocity corresponding to $E_\alpha$ and $S_\alpha$ is the birth rate of alphas per unit time per unit volume, retrieved from the Maxwellian-averaged formula in Richardson (2019, p. 45).

Alphas are born at a high energy and lose energy via friction to the electrons. Eventually, they will reach a critical velocity, $v_c$, where they begin to scatter off ions. Helander & Sigmar (2002) gives an approximation of this critical velocity (or in this case, energy) as

(3.18)\begin{equation} \frac{m_\alpha v_c^2}{2} \sim 50 T_e. \end{equation}

This critical energy is still significantly greater than the centrifugal potential well (2.20), leading to the relation

(3.19)\begin{equation} \varXi_\alpha \sim \frac{M^2 T_e}{4} \ll \frac{m_\alpha v_c^2}{2} \sim 50 T_e, \end{equation}

as long as $M^2 / 4 \ll 50$ (which generally is true). So we can assume that the alphas are not centrifugally confined and generally do not deposit energy into the ions. However, energy is lost through drag to alpha–electron collisions. The assumption that all the alpha energy is deposited into the electrons is pessimistic because hotter electrons lead to worse ion confinement to maintain quasineutrality.

As alpha particles are born isotropically, a fraction of them are born directly into the unconfined region of velocity space. We thus model the loss region for energetic alphas as a cone, and only classical mirror confinement applies, thus an alpha is lost if

(3.20)\begin{equation} \mu_\alpha B_{\mathrm{max}} > E_\alpha, \end{equation}

where $\mu _\alpha$ is the magnetic moment. Integrating this over all velocity space, we see that the fraction of alphas that is lost is (Pastukhov 1987)

(3.21)\begin{equation} f_{\mathrm{lost}} = 1 - \sqrt{1 - \frac{1}{R_{\mathrm{mirror}}}}. \end{equation}

This fraction of alpha particles is lost promptly and is taken into account in assessing quasineutrality.

Currently, we do not have any explicit collisional losses of alpha particles. However, under the twin assumptions that of $R_{\mathrm {mirror}} \gg 1$ and that we can treat alphas like primary ions (see Ryutov 1988), the approximate lifetime of alpha particles in the machine is

(3.22)\begin{equation} \tau_\alpha \approx 0.4 \tau_{\alpha i} \ln R_{\mathrm{mirror}}, \end{equation}

where $\tau _{\alpha i}$ is the alpha–ion collision time, as alpha–electron collisions do not scatter the alpha particles (Helander & Sigmar 2002). Additionally, we do not account for dilution of the primary ion species due to the accumulation of helium ash, neither do we account for loss-cone-driven instabilities such as those discussed in Hanson & Ott (1984).

4 Results and discussion

MCTrans++ can be run in two steady-state modes: a single point in parameter space and batch mode to perform parameter scans. It also offers time-dependent options, including a capacitor bank discharge and ‘free-wheeling’ mode where a steady-state plasma spins down, i.e. is discharged through a dump resistor. The following sections discuss results from all these modes of operation.

4.1 Experimental comparison

Previously, there have been three centrifugal mirror experiments, including Ixion (Baker & Hammel 1961; Baker et al. 1961), PSP-2 (Volosov 2009) and MCX (Ellis et al. 2005). Results from MCTrans++ are compared with the available data of those experiments (table 2). The following paragraphs detail the assumptions made to model each experiment in MCTrans++.

Table 2. Benchmarking of MCTrans++ against previous experiments. Experimental results for Ixion, PSP-2 and MCX are taken from (Baker et al. 1961; Volosov 2009; Teodorescu et al. 2008), respectively. Some results did not report momentum confinement time, so *ion confinement time is reported instead. Collisionality is also calculated, and we find that both Ixion and MCX were highly collisional.

4.1.1 Prior experiments

In Ixion (Baker et al. 1961), preionized deuterium gas was injected into the chamber from the side and a negative bias was applied with short inner electrodes via capacitor bank discharges up to 20 kV (though usually 7.5 kV). When the gas was ionized, an axial column of plasma extended from one electrode to the other, thus creating a ‘plasma centre electrode’. The diameter of this plasma column at the midplane was $\sim$1.5 times larger than that of the electrode (6.4 cm), and we assume this is the innermost flux surface. The outer radius and axial extent of the plasma was inferred from voltage profile measurements (9.5 and 10 cm, respectively). The field at the midplane was 0.95 T, with a mirror ratio of 2.2. Baker et al. (1961) do note that poor vacuum quality could have meant that a typical discharge had as many impurity ions as deuterium, but without an actual measurement, we have assumed $Z_{\mathrm {eff}} = 3$. By modelling the plasma as a capacitor and discharging it into a known resistance, the characteristic ‘spin-down’ time can be related to the total charge stored, and therefore the ion density, which was estimated as $3\times 10^{21}\,\textrm {m}^{-3}$. Ixion produced highly collisional plasmas (see table 1), and because MCTrans++ assumes low-collisionality, the modelled electron temperature (and thus $M_A$) is significantly lower than the measured value. The momentum confinement time, $\tau _M$, is of similar magnitude.

Volosov (2009) wrote a review paper on the PSP-2 experiment which operated from 1975–1985, mostly detailing the work in Abdrashitov et al. (1991). The electric field was not generated via a centre electrode, but rather a series of matching ring electrodes which were charged up to 500 kV (positive bias), though typical discharges were 360 kV. The inner and outer plasma radii are given as 32 and 51 cm, respectively, and the length was measured by neutral detectors to be 40 cm. Hydrogen was simultaneously pumped through six valves spread azimuthally around the midplane. The typical midplane magnetic field was 0.99 T with a mirror ratio of 2.4. Usual discharge densities were of the order of $3\times 10^{17}\,\textrm {m}^{-3}$ and $1\times 10^{18}\,\textrm {m}^{-3}$ for ionized hydrogen and neutrals, respectively. The relative abundance of impurities is also given and $Z_{\mathrm {eff}}$ was calculated to be $\sim 2.4$. The energy confinement time was not calculated, so the ion confinement time ($\tau _D = 100\,\mathrm {\mu }$s in that paper) was used instead. Electron temperatures were indirectly measured based on assumptions relating the ambipolar potential and $T_e$, giving mean electron energies (non-Maxwellian) of 0.1–1 keV. Abdrashitov et al. (1991) reports the mean ion energy was up to 20 keV in the rotating frame; however, this was not a direct measurement of ion temperature, rather a measurement of fast neutrals from charge exchange. Ion drift velocities were found to be $2\times 10^6$ m s$^{-1}$. The primary difference between the experimental results and MCTrans++ is in the confinement time. This difference can be primarily attributed to the large neutral inventory (roughly an order-of-magnitude larger than the charged species density). By balancing loss rates and ionization, the neutral source in MCTrans++ is just large enough to keep a constant electron density such that, typically, $n_n \ll n_e$. However, if in reality $n_n \gg n_e$, the ion confinement time will be much smaller due to charge exchange losses.

MCX operated at the University of Maryland until 2012 (Ellis 2012), and Teodorescu et al. (2008, 2010) provided an overview of typical experiments. Hydrogen gas was pumped into the chamber until a base prefill pressure of $\sim$5 mTorr was achieved. Voltage to a central electrode was applied through a capacitor bank that was typically charged up to $-$10 kV. The magnetic field was such that $R_{\mathrm {mirror}} = 7.3$ and $B_{\mathrm {min}} = 0.23$ T. The inner conductor and chamber walls limit the inner and outer flux surfaces, at radii of 6 and 26 cm, respectively, with a plasma length of 1.3 m. We again assume $Z_{\mathrm {eff}}=3$. Interferometric methods measured densities in the range of $5\times 10^{20}$ m$^{-3}$, and thermal electron Bernstein emission provided peak electron temperatures of 100 eV (Reid et al. 2014). Teodorescu et al. (2008) briefly mentions that confinement times of $100\,\mathrm {\mu }$s were used for other calculations, as well. The performance of MCTrans++ shows good correspondence to MCX, despite the high collisionality (table 2).

Many of the experiments described above used a negative (rather than positive) bias to control plasma–surface interactions at the electrodes and aid in breakdown. This effect is only important when the plasma does not yet have a sufficient electron density and conductivity to completely screen imposed parallel electric fields. Once the plasma is fully ionized and weakly collisional, the parallel electric field is completely determined by the ambipolar potential from (2.16). At present, MCTrans++ does not model this early phase of the plasma and is thus agnostic to the direction of the electric field.

4.1.2 The CMFX comparison

At the time of writing, CMFX has demonstrated long-lived plasmas with discharges up to 40 kV (Schwartz et al. 2023). Unfortunately, detailed validation between the current model and the experiment is not yet possible because many crucial experimental measurements are still under way – crucially, $T_i$, $n_e$ and $L$. However, these discharges have provided some results for initial comparison. For this specific section, we will refer to the parameters as ‘CMFX (Expt.)’ to avoid confusion with the eventual goals of the experiment, which are referred to as simply ‘CMFX’ in the rest of the paper.

The following discussion compares MCTrans++ with preliminary results from CMFX. Capacitor discharges of nominally 25 kV resulted in steady-state deuterium plasmas at $8.1 \pm 0.1$ kV. The voltage and current traces provide estimates of some global variables like resistance, stored energy, capacitance and momentum confinement time. Table 3 displays the configuration parameters for MCTrans++ and a comparison with these preliminary experimental results. Due to unavoidable experimental uncertainties and unknowns, we find this result is consistent with our model.

Table 3. Nominal configuration parameters and outputs for a CMFX (Expt.) comparison. All the fields with ‘—’ are currently unavailable measurements or not applicable parameters. This comparison with recent experimental results should not be confused with the mention of CMFX in the rest of the paper, which considers the eventual operational goals of the device.

MCTrans++ provides optional variables to help fit experimental data when values like temperature or density are unavailable. Parallel loss factor and neutral density ($n_n$) were varied until the results approximately converged on the experimental data. Because MCTrans++ assumes a low-collisionality plasma, it over-predicts loss rates for collisional plasmas; thus, because the predicted value of $\nu ^* = 4.6 > 1$, we have set parallel loss factor to 0.1, i.e. the parallel loss rates have been artificially modified to 10 % of the value normally predicted by MCTrans++. Additionally, because gas is only puffed at the beginning of the experiment, not continuously throughout it, we expect the neutral inventory to be quite low. We therefore set $n_n$ to ${\sim }\times 10^{-6} n_e$.

4.2 The CMFX and reactor-scaling

To explore the parameter space of interest in CMFX and a reactor, the central field, electron density and applied voltage were varied. The fixed parameters provided to the input files for MCTrans++ are given in table 4. The ion density is assumed equal to the electron density to enforce quasineutrality.

Table 4. Nominal configuration parameters for CMFX and reactor scenarios considered here.

The Alfvén speed, $v_A = {B}/{\mu _0 n_i m_i}$, is the speed at which magnetic field perturbations propagate along the axial direction. The best performing plasmas occur when the Alfvén Mach number, $M_A \equiv {u}/{v_A}$ approaches unity. However, the results may appear sparse in some areas because configuration parameters with either $M_A > 1.25$ or $\rho ^* = {\rho _i}/{a} > 0.1$ were not plotted. As the Alfvén Mach number passes unity, the plasma is squeezed in the $z$-direction into a thin disk by the large centrifugal forces. Eventually, at $M_A \approx 1.25$, the magnetic field no longer has an equilibrium shape as the diamagnetic current and its associated field overcome the vacuum magnetic field, creating a magnetic null and causing the field lines close upon themselves (Abel 2022).

Surpassing the limit $\rho ^* \approx 0.1$ brings into question the strongly magnetized assumption that underlies (2.10)–(2.12) and the ordering $\varPhi \gg \varphi$ (2.8). However, as opposed to the $M_A$ limit, which results in a plasma with no equilibrium, the $\rho ^*$ limit is merely a constraint of the model. Future work should consider the regime where extended gyroradii exist – in particular the orbits of fusion-produced alpha particles.

Some lower temperature results are limited by the condition that $\nu ^* \leq 0.1$. Like the limit for $\rho ^*$ (and unlike the $M_A$ limit), this is not a physical limit, but merely a limitation of the model. The assumption that $T_s$ is constant on a flux surface is violated when the plasma becomes collisional, which may affect the results from § 2.2. Future work should consider plasmas in the collisional regime due to their importance during the start-up phase (step (iv) in § 3.1).

The results for the given CMFX and reactor configurations in table 5 provide some outputs of MCTrans++ for parameters of interest.

Table 5. Results predicted by MCTrans++ for CMFX and reactor configurations. See table 4 for the device parameters. Variables denoted by $P$ indicate angular momentum losses, and those by $\dot {Q}$ are heat losses. As described in § 4.2.2, viscous heating mediates between angular momentum and heat, in that it is a loss for the former and a gain for the latter – hence the crossover of two curly braces. Some values are not reported for CMFX because they are only relevant for devices with DD or DT fuel. Additionally, some parameters are listed in earlier tables, and are listed again for convenience.

4.2.1 Performance parameters

Some parameters of interest were recorded while central field strength, electron density and voltage were varied. Results for CMFX are in figure 5 and those for a reactor in figure 6. The results were only considered valid if $M_A \leq 1.25$, $\rho ^* \leq 0.1$ and $\nu ^* \leq 0.1$ as described in the prior section. Dashed lines indicating these limits are in each figure where appropriate.

Figure 5. Performance of a CMFX-like device for 3 T throat field with (a) a range of central fields ($B_\mathrm {min}$) and fixed electron density ($n_e = 1\times 10^{19}$ m$^{-3}$) and (b) a range of electron densities ($n_e$) and fixed central field ($B_{\mathrm {min}} = 0.3$ T). Results were cut off above the values of $M_A > 1.25$ and $\rho ^* > 0.1$, indicated by the magenta and green dashed lines, respectively.

Figure 6. Performance of a reactor-like device for 18 T throat field with (a) a range of central fields ($B_\mathrm {min}$) and fixed electron density ($n_e = 9\times 10^{19}$ m$^{-3}$) and (b) a range of electron densities ($n_e$) and fixed central field ($B_{\mathrm {min}} = 5$ T). Results were cut off above the values of $M_A > 1.25$ and $\rho ^* > 0.1$, indicated by the magenta and green dashed lines, respectively, and where appropriate.

Here $P_{\mathrm {in}}$ comes directly from the power supply and is used to drive the rotation. Section 4.2.2 provides a more in-depth explanation of power losses, but generally, because rotational velocity is inversely proportional to $B_{\mathrm {min}}$ ($u \propto E / B$), the power dissipated through viscous torque decreases with larger values of $B_{\mathrm {min}}$. We see that trend reflected in the results for CMFX in figure 5. There is a trade off in higher rotational velocity, however, between power lost to viscous torque and power saved by enhanced parallel confinement (2.21). For this regime, we see that higher rotational speeds lead to increased power draw.

For the same reason, we see higher ion temperatures for lower values of $B_{\mathrm {min}}$ (or higher values of $\varPhi$) because the viscous torque transfers more kinetic (and thus internal) energy to the plasma. We also find that supersonically rotating mirrors organically produce a hot-ion mode. The heating mechanism is viscous shear, which primarily heats the ions because they carry nearly all the rotational momentum in comparison with electrons. Direct heating of the ions is in fact one of the main benefits of supersonic rotation because it negates the need for auxiliary heating systems such as RF or neutral beams.

We also see that CMFX is limited by $\rho ^*$, not $M_A$, where the limit scales inversely with $B_{\mathrm {min}}$ and proportionally with $\sqrt {T_i}$ (recall that $\rho ^* = {\sqrt {2 m_i T_i}}/{ZeB}$). This limit does not mean that plasmas cannot exist for $\rho ^* > 0.1$, but rather the model is not applicable in that regime. We also see that our results are limited at lower temperatures, where $\nu ^* > 0.1$ and the plasma is collisional. Lastly, the overlapping lines in figure 5(b) show that ion and electron temperatures are approximately invariant with density in this temperature regime.

Similar plots for a reactor scenario (figure 6) demonstrate that the highest performance plasmas are at the limits of $M_A \leq 1.25$ and $\rho ^* \leq 0.1$. Moreover, there is competition between high scientific gain, $Q_{\mathrm {sci}}$, and high thermal power output, $P_{\mathrm {thermal}}$. Note this is dissimilar to tokamaks, where high fusion power correlates with high gain (Costley 2016).

Another general observation is that the best performing plasmas operate at ion temperatures that are around the peak location of $\sigma _{\mathrm {DT}} \sim 70$ keV. Additionally, we continue to achieve a hot-ion mode at lower voltages, but at higher voltages, a hot-electron mode develops. At these higher voltages, the alpha heating is significantly larger, and we assume that all the alpha energy is deposited into the electrons (see § 3.4). This is a pessimistic assumption because some amount of alpha energy will be transferred to the ions, with the exact proportion dependent on the precise conditions.

By observing the dashed lines in figure 6(a), we see that increasing the central field for a reactor scenario allows for much higher $Q_{\mathrm {sci}}$ (though it does require larger applied voltages). In contrast, the same dashed lines have negative curvature in the graph for $P_{\mathrm {thermal}}$, appearing to reach some asymptotic value. Thus, for a desired $P_{\mathrm {thermal}}$, achieving higher $Q_{\mathrm {sci}}$, requires both larger magnetic field and voltage. The $\nu ^*$ limit only affects the lowest temperature discharges in this regime.

The changing density in figure 6(b) is a more complicated story. The performance is primarily limited by $\rho ^* \leq 0.1$, except for very high values of density ${\gtrsim }1\times 10^{20}$ m$^{-3}$. Interestingly, there is a maximum value of $Q_{\mathrm {sci}}$ for $n_e \sim 5\times 10^{18}$ m$^{-3}$, but it has quite a low thermal power. While low densities are not limited by $M_A$ or $\rho ^*$, the output power is typically too low because fusion power is proportional to $n_s^2$ (Richardson 2019).

4.2.2 Power losses

We assume that a direct current (DC) power supply directly drives the rotation. The power balance in (2.12) is used to solve for the radial current density, $j_R$. The current draw from the power supply is balanced by angular momentum losses including particle losses and viscous heating. The viscous heating increases the thermal energy in (2.11), where the loss mechanisms are particle losses and radiation. Essentially, viscous heating acts like a mediator between momentum and heat – a loss term for the prior and a source term for the latter. A breakdown of loss mechanisms for both heat and momentum for a CMFX and reactor scenario shows the importance of each at various applied voltages (figure 7). As described in § 2.3, we optimistically choose $R_\mathrm {exh} = R_\mathrm {min}$, which affects the angular momentum lost by a particle in the parallel direction, $m_s \omega R_\mathrm {exh}^2$. For transparency, figure 7 also includes the pessimistic assumption that $R_\mathrm {exh} = R_\mathrm {max}$.

Figure 7. Power losses for (a) CMFX-like configuration ($n_e = 1\times 10^{19}$ m$^{-3}$ and $B_{\mathrm {min}} = 0.3$ T) and (b) reactor configuration ($n_e = 9\times 10^{19}$ m$^{-3}$ and $B_{\mathrm {min}} = 4$ T). Note the different scales on the horizontal axes for the applied voltage. The momentum lost due to parallel ion losses are a function of $R_\mathrm {exh}$, and the results for $R_\mathrm {exh} = R_\mathrm {min}$ are shown with solid lines and those for $R_\mathrm {exh} = R_\mathrm {max}$ with dashed lines.

In a CMFX scenario with lower voltages (corresponding to temperatures $\lesssim$1 keV), charge exchange is the dominant momentum loss mechanism, followed by viscous torque and then parallel ion losses. In the case where $R_\mathrm {exh} = R_\mathrm {max}$, the momentum lost due to ion losses is a factor of $R_\mathrm {mirror}$ larger. The rate coefficient for charge exchange is significantly higher than either electron- or ion-impact ionization (figure 4) for a given rotational speed, meaning that the ion loss rate due to charge exchange is large. However, there appears to be a local maximum in momentum losses from charge exchange, indicating that at sufficiently high temperatures, charge exchange ceases to be a limiting factor.

The heat loss mechanisms for CMFX monotonically increase, all being of equal magnitude except for radiation. The heat lost to charge exchange is not the dominant mechanism, as in momentum losses, because the rotational energy of lost ions is much greater than the thermal energy.

For a reactor scenario, viscous torque is the primary momentum loss mechanism, comparable to parallel ion losses at higher voltages. If we take $R_\mathrm {exh} = R_\mathrm {max}$, the momentum loss due to the parallel ion loss rate is dominant. Similarly, as the rotational speed (and therefore temperature) increases, we expect viscous torque to increase and charge exchange to decrease.

Heat loss in a reactor scenario is drastically different from CMFX in that charge exchange, and it is very small and radiation becomes important. Here, the primary losses are through parallel transport and radiation (roughly equal bremsstrahlung and cyclotron, see table 5), and charge exchange plays almost no role. The heat lost to electrons slightly dominates over parallel ion losses because the electron loss rate is larger and this device operates in a hot-electron mode.

As mentioned, the choice of $R_\mathrm {exh}$ plays a critical role in parallel ion losses for momentum balance. We can see exactly how that choice affects the required $P_\mathrm {in}$ in figure 8. Here $P_{\parallel i}$ increases by a factor of $R_\mathrm {mirror}$ from $R_\mathrm {min}$ to $R_\mathrm {max}$. For CMFX, this does not drastically alter the power draw because parallel ion losses are not the dominant loss mechanism of momentum (see figure 7, charge exchange is dominant). In contrast, $Q_\mathrm {sci}$ in a reactor decreases three-fold because parallel ion losses become dominant with choices of large $R_\mathrm {exh}$.

Figure 8. Variation of key parameters for different values of $R_\mathrm {exh}$ for (a) CMFX and (b) reactor scenarios.

4.2.3 Effect of some physical behaviours

MCTrans++ also has the option to turn off some physical behaviours, namely charge exchange and the ambipolar potential. While removing these options is entirely non-physical, the results are useful in demonstrating the importance of charge exchange and the ambipolar potential. The effect of these phenomena can be seen in figure 9 for both CMFX and reactor scenarios. In the case where we ‘turn off’ the ambipolar potential, we assume $\varphi = \varphi _0$, the centrifugal potential. There is then some unbalanced parallel current because we do not consider $\varphi _1$, and it is not compensated by a change in radial current. This option is present in the code mainly as a debugging tool. We use it here merely to demonstrate that we do indeed need to consider $\varphi _1$ and that its effects are not negligible.

Figure 9. Effect of turning off charge exchange and ambipolar potential for (a) CMFX-like configuration ($n_e = 1\times 10^{19}$ m$^{-3}$ and $B_{\mathrm {min}} = 0.3$ T) and (b) reactor configuration ($n_e = 9\times 10^{19}$ m$^{-3}$ and $B_{\mathrm {min}} = 4$ T). Note the different scales on the horizontal axes for the applied voltage.

For CMFX, a lack of charge exchange is beneficial because that is the major momentum loss mechanism at these voltages (figure 7a), which in turn increases the ion temperature. Conversely, when we turn off the ambipolar potential, there is a finite axial current (again, non-physical) that leads to higher electron loss rates and lower ion loss rates. Therefore, the viscous torque increases, requiring more power to sustain rotation. Despite the ion temperature decreasing slightly, the electron temperature decreases significantly more, establishing an even more pronounced hot-ion mode.

In contrast, for the reactor scenario, we see that charge exchange plays only a little role because it is the weakest loss mechanism (figure 7b). On the other hand, we see an exaggerated hot-ion mode appear with the absence of $\varphi$, leading to a $Q_{\mathrm {sci}}$ value that is quadrupled at some voltages. A lack of $\varphi$ allows for a finite parallel current, which produces higher ion confinement. This improved confinement is what sustains a hot-ion mode and therefore high $Q_{\mathrm {sci}}$ values.

4.3 Time dependence

As mentioned previously, there are two time-dependent modes of running: (i) capacitor bank discharge and (ii) free-wheeling into a load resistor. The $72\,\mathrm {\mu }$F capacitor bank discharge (figure 10a) starts with 100 kV across the plasma, and slowly draws current from the capacitors which decrease in voltage. The power draw is extremely nonlinear, with a large spike at the end due to an increase in current across the plasma at lower temperatures. The time scale of this behaviour is of the order of seconds because the plasma resistance is between 10 and 100 k$\Omega$ and therefore the resistor–capacitor (RC) time for a capacitor bank of this size is seconds long.

Figure 10. Time-dependent results for a CMFX configuration. (a) A $72\,\mathrm {\mu }$F capacitor bank, with 1 G$\Omega$ internal resistance, 10 $\Omega$ line resistance and $10\,\mathrm {\mu }$H line inductance is discharged into the plasma. (b) At 100 kV, the crowbar is closed and the plasma allowed to discharge into a 10 k$\Omega$ external resistor. Note the difference in time scales between the two scenarios.

The free-wheeling mode discharges a CMFX-like plasma at 100 kV into a 10 k$\Omega$ external dump resistor. The plasma stores energy like a hydromagnetic capacitor, meaning that its behaviour will be similar to a capacitor being discharged into the same load. The main difference is that the capacitance of the plasma is variable, and therefore the power draw does not fall exponentially as with a static capacitor. As the voltage, and therefore temperature, drops, charge exchange becomes the dominant loss mechanism. MCTrans++ predicts a plasma capacitance of 440 pF for a CMFX-like configuration at 100 kV and a density of $1\times 10^{19}\,\textrm {m}^{-3}$. With a dump resistance of 10 k$\Omega$, the RC time for this system is 4.4 ms. In reality, the voltage decay is slightly faster because the capacitance decreases as the plasma temperature falls.

5 Conclusions

MCTrans++ is a 0-D model for centrifugal mirrors. It is based on several assumptions, including a strongly magnetized and low-collisionality plasma, a large mirror ratio and supersonic rotation. The confining potential is due to both an ambipolar and centrifugal potential, the latter of which is confining for both ions and electrons. In the low-collisionality limit, perpendicular losses are classical, and continuum radiation losses like bremsstrahlung and cyclotron emission are included. A neutrals model has been implemented to determine the electron- and ion-impact ionization and charge exchange rates. The retained alpha particles are assumed to deposit all their energy into the electrons, and alphas are otherwise lost through a axisymmetric loss cone.

Comparison with prior results shows good agreement with MCX. MCTrans++ could not compute a solution for an Ixion-like configuration because the plasma was highly collisional. The comparison with PSP-2 showed differences likely due to the high neutral densities in that device.

Scaling of CMFX demonstrates that it is limited by $M_A$ at higher densities and increased performance as voltage increases. Additionally, supersonic rotation leads to a naturally occurring hot-ion mode, with direct ion heating resulting through momentum transfer to the ions via viscous shear. Scaling of the reactor scenario posits that a balance between low-field and high-density generates high $Q_{\mathrm {sci}}$ and $P_{\mathrm {thermal}}$. A hot-electron mode does appear, but this is due to the alpha particles depositing their energy into the electrons.

We find that at lower ion temperatures, the dominant power loss mechanism is through charge exchange, but at fusion-relevant temperatures, parallel losses become more important. Similarly, we see that radiation heat losses become more important at higher temperatures.

MCTrans++ can simulate time-resolved solutions, including a capacitor discharge and free-wheeling into an external dump resistor. The plasma behaves similarly to a capacitor, storing its energy in rotation. In fact, the time scales of parameter evolution are very close to the those of an RC circuit. In both cases, as the plasma temperature drops, we see a large spike in power draw due to charge exchange, while all the other loss mechanisms decrease.

Future work will implement line-radiation and collisional transport, as well as model the expander region near the insulators. Work is currently in progress to attain radial profiles and model magnetic equilibrium conditions.