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. Reference Ellis, Hassam, Messer and Osborn2001)

(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. Reference Reid, Romero-Talamás, Young, Ellis and Hassam2014), 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. (Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013)) 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 bodiesFootnote ¹ (Ferraro Reference Ferraro1937; Lehnert Reference Lehnert1971). 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 Reference Lehnert1971),

(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 (Reference Catto, Bernstein and Tessarotto1987), and Abel et al. (Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013))

(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$.Footnote ² 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 (Reference Pastukhov1974)) 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 Reference Post1961).

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}$.Footnote ³ 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 Reference Abel2022).

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. Reference Ellis, Hassam, Messer and Osborn2001). We leave $Z_i$ as a free parameter elsewhere, but use $Z_i = 1$ simply to show agreement with Ellis et al. (Reference Ellis, Hassam, Messer and Osborn2001). 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 Reference Pastukhov1974). Although originally found for electrons alone, we expand the results for a multispecies plasma. Catto (Reference Catto1981) 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 notationFootnote ⁴ ,

(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 (Reference Braginskii1965)).

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$Footnote ⁵. 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 satisfiedFootnote ⁶ . 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 Reference Huang and Hassam2001), 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 (Reference Braginskii1965) or Helander & Sigmar (Reference Helander and Sigmar2002).

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 (Reference Helander and Sigmar2002))

(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 (Reference Braginskii1965) 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 (Reference Richardson2019). 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 (Reference Braginskii1965). 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 (Reference Braginskii1965) 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 (Reference Braginskii1965), 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 (Reference Braginskii1965), 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 (Reference Pastukhov1974) and Catto (Reference Catto1981) 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. Reference Baker, Hammel and Ribe1961), PSP-2 (Volosov Reference Volosov2009) and MCX (Teodorescu et al. Reference Teodorescu, Clary, Ellis, Hassam, Lunsford, Uzun-Kaymak and Young2008; Reid et al. Reference Reid, Romero-Talamás, Young, Ellis and Hassam2014) are from prior experiments, respectively, whereas those for CMFX and Reactor came from MCTrans++ predictive models.

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. Reference Anderson, Baker, Bratenahl, Furth and Kunkel1959).

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 (Reference Ng and Hassam2007) 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 (Reference Janev and Smith1993). 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 Reference Drawin and Emard1976; Tallents Reference Tallents2018). 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 (Reference Appelbe and Chittenden2011)) 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. (Reference Ellis, Case, Elton, Ghosh, Griem, Hassam, Lunsford, Messer and Teodorescu2005)) 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 (Reference Richardson2019)) 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 (Reference Richardson2019)) 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 (Reference Tamor1983). 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}}$Footnote ⁷, 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 (Reference Richardson2019, 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 (Reference Helander and Sigmar2002) 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 Reference Pastukhov1987)

(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 Reference Ryutov1988), 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 Reference Helander and Sigmar2002). 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 (Reference Hanson and Ott1984).