## 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 Reference Post1987). 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 Reference Lehnert1971). Other concepts for improved confinement and stabilization in a magnetic mirror (or ‘open trap’) are outlined in Ryutov *et al.* (Reference Ryutov, Berk, Cohen, Molvik and Simonen2011).

Supersonic rotation has been demonstrated to improve axial confinement (Teodorescu *et al.* Reference Teodorescu, Young, Swan, Ellis, Hassam and Romero-Talamas2010) and stability (Huang & Hassam Reference Huang and Hassam2001), thereby reducing perpendicular losses too. This concept was first demonstrated with the Ixion device at Los Alamos National Lab (Baker & Hammel Reference Baker and Hammel1961; Baker, Hammel & Ribe Reference Baker, Hammel and Ribe1961), and variations have been constructed at the Institute of Nuclear Physics in Russia (PSP-2, Volosov (Reference Volosov2009)) and the University of Maryland (MCX, Ellis *et al.* (Reference Ellis, Case, Elton, Ghosh, Griem, Hassam, Lunsford, Messer and Teodorescu2005)), 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.* Reference Romero-Talamas, Beaudoin, Hassam, Koeth, Eschbach, Short, Schwartz, Kelly, Abel and Basu2022). 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 Reference Huang and Hassam2001; Romero-Talamás *et al.* Reference Romero-Talamás, Elton, Young, Reid and Ellis2012). 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.

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.* Reference Abdrashitov, Beloborodov, Volosov, Kubarev, Popov and Yudin1991), and current experiments, like GDT (Beklemishev *et al.* Reference Beklemishev, Bagryansky, Chaschin and Soldatkina2010) and TAE's C-2W (Binderbauer *et al.* Reference Binderbauer, Tajima, Steinhauer, Garate, Tuszewski, Schmitz, Guo, Smirnov, Gota and Barnes2015), use edge-biasing with concentric ring electrodes to establish a radial electric field profile in the plasma. Ryutov *et al.* (Reference Ryutov, Berk, Cohen, Molvik and Simonen2011) 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 Reference Huang and Hassam2001; Romero-Talamás *et al.* Reference Romero-Talamás, Elton, Young, Reid and Ellis2012). 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.* Reference Ellis, Hassam, Messer and Osborn2001)

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

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

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

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.

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 ^{1} (Ferraro Reference Ferraro1937; Lehnert Reference Lehnert1971). Mathematically, the angular frequency of any such surface $\omega$ obeys

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

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

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

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

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

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,

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 ^{2} 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

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

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

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 ^{3} 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,

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:

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

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,

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

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 ^{4} ,

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 ^{5}. 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

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 ^{6} . 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))

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

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:

*a,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,

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

For electrons, (2.13) from Braginskii (Reference Braginskii1965) gives

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.5 Assumptions and scope

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

(i) be strongly magnetized,

(ii) have low collisionality,

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

(iv) have a large mirror ratio,

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

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

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.

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

(ii) Neutral gas is puffed into the chamber.

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

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

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

(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:

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

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

and the distribution function for the cold neutrals is

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,

so that the integrand now becomes

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

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

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

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

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

where we have used the definition

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)

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 ^{7}, where

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,

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

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

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

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)

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

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

## 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 Reference Baker and Hammel1961; Baker *et al.* Reference Baker, Hammel and Ribe1961), PSP-2 (Volosov Reference Volosov2009) and MCX (Ellis *et al.* Reference Ellis, Case, Elton, Ghosh, Griem, Hassam, Lunsford, Messer and Teodorescu2005). 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++.

#### 4.1.1 Prior experiments

In Ixion (Baker *et al.* Reference Baker, Hammel and Ribe1961), 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.* (Reference Baker, Hammel and Ribe1961) 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 (Reference Volosov2009) wrote a review paper on the PSP-2 experiment which operated from 1975–1985, mostly detailing the work in Abdrashitov *et al.* (Reference Abdrashitov, Beloborodov, Volosov, Kubarev, Popov and Yudin1991). 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.* (Reference Abdrashitov, Beloborodov, Volosov, Kubarev, Popov and Yudin1991) 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 Reference Ellis2012), and Teodorescu *et al.* (Reference Teodorescu, Clary, Ellis, Hassam, Lunsford, Uzun-Kaymak and Young2008, Reference Teodorescu, Young, Swan, Ellis, Hassam and Romero-Talamas2010) 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.* Reference Reid, Romero-Talamás, Young, Ellis and Hassam2014). Teodorescu *et al.* (Reference Teodorescu, Clary, Ellis, Hassam, Lunsford, Uzun-Kaymak and Young2008) 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.* Reference Schwartz, Romero-Talamas, Hassam, Beaudoin, Abel, Koeth, Eschbach, Short and Kelly2023). 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.

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.

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 Reference Abel2022).

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.

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

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 Reference Costley2016).

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 Reference Richardson2019).

#### 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}$.

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}$.

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

For CMFX, a lack of charge exchange is beneficial because that is the major momentum loss mechanism at these voltages (figure 7*a*), 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 7*b*). 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 10*a*) 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.

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.

## Acknowledgements

*Editor C. Forest thanks the referees for their advice in evaluating this article.*

## Funding

This work was supported by ARPA-E award no. DE-AR0001270.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Review of shear flow stabilization

Shear flow stabilization of turbulence is key to the success of centrifugally confined mirrors. The assumption that perpendicular transport is classical hinges on this behaviour. While this phenomenon has been demonstrated for rotating mirrors both theoretically and experimentally, the concept is also pervasive in the tokamak literature.

There are two primary instabilities which must be suppressed in order to eliminate turbulent transport. The interchange (also known as ‘flute’) instability occurs when the magnetic field curvature and pressure gradients point in the same direction, and is similar to the fluid Rayleigh–Taylor instability, where the effective outward gravity comes from centrifugal rotation (Goldston Reference Goldston2020). The interchange instability occurs on the same scale as the pressure gradient ($\boldsymbol {\nabla } p / p$). This instability was the primary one observed on older mirror machines (Post Reference Post1987), which led to the creation of the ‘minimum-B’ configuration that included complicated magnetic coil shapesFootnote ^{8} to combat the poor magnetic curvature.

The Kelvin–Helmholtz (K-H) instability can occur when there is shear flow and thus relative motion between fluid surfaces. For a strongly magnetized plasma, the scale of this instability is much larger than the ion gyroradius (Horton, Tajima & Kamimura Reference Horton, Tajima and Kamimura1987). The concern in a centrifugal mirror is that the shear flow between radial layers might cause K-H vortices to form.

#### A.1 Theoretical

Burrell (Reference Burrell1997) provides a review of the theoretical underpinnings of shear flow stabilization. A highly idealized approximation demonstrates nonlinear decorrelation of turbulence, whereby radial turbulent flux can be reduced, even if the plasma starts in a highly turbulent configuration. A more complex analysis of linear stabilization is mode dependent, but the general theme is that unstable modes are coupled to stable modes through velocity shear, which in turn dampen the unstable modes. The key takeaway is that, if the turbulent decorrelation time is long compared with the shearing rate, turbulence is suppressed.

Some work has been specifically applied to stabilization of centrifugally confined plasmas. For example, Huang & Hassam (Reference Huang and Hassam2001) used a three-dimensional ideal magnetohydrodynamic (MHD) model to demonstrate that supersonic shear flows suppress growth of the interchange instability. Huang & Hassam (Reference Huang and Hassam2004) considered four primary instabilities – interchange and K-H in the low $\beta$ limit, and magnetorotational and Parker in the high $\beta$ limit. The conclusions for each are as follows: (i) the interchange mode can be suppressed with sufficiently high $M$, the limit of which can be lowered for elongated plasmas, (ii) the K-H instability is marginally MHD stable provided there are no inflection points in the density and rotation profile. Particle sources play a key role in attaining stability for these two, and must be considered further. (iii) Because $M_A \lesssim 1$, the magnetorotational instability does not occur, and (iv) the Parker instability may simply create new equilibria, but further works needs to be done.

Theoretical work on PSP-2 posits that the combination of shear flow and line-tying (achieved by concentric conducting rings to set the electric field profile) are necessary for MHD stability (Volosov Reference Volosov2009). However, that hypothesis is only true for particular radial electric field profiles. The PSP-2 profile is relatively flat (see figure 14 in that work), and at any point where the gradient is near zero, there will be no shear and therefore no shear stabilization. In contrast, a singly peaked profile (as has been achieved with a centre conductor (Huang & Hassam Reference Huang and Hassam2001; Romero-Talamás *et al.* Reference Romero-Talamás, Elton, Young, Reid and Ellis2012)) and sufficiently large rotational shear, can stabilize centrifugally confined plasmas in the absence of line-tying. Ryutov *et al.* (Reference Ryutov, Berk, Cohen, Molvik and Simonen2011) also states that the centrifugal force alone is sufficient to achieve stability.

Recent theoretical work on shear flow stabilization in tokamaks has demonstrated suppression of all linear turbulence in the condition of zero magnetic shear (Highcock *et al.* Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010). Large enough shear flow can suppress nonlinear turbulence due to ion temperature gradients (Highcock *et al.* Reference Highcock, Barnes, Schekochihin, Parra, Roach and Cowley2010, Reference Highcock, Barnes, Parra, Schekochihin, Roach and Cowley2011). For subsonic shear flows, heat is transported neoclassically and momentum by turbulence. However, if the flow is supersonic, the parallel velocity gradient instability can cause a transition back to turbulence (Newton, Cowley & Loureiro Reference Newton, Cowley and Loureiro2010; Highcock *et al.* Reference Highcock, Barnes, Parra, Schekochihin, Roach and Cowley2011). Parra *et al.* (Reference Parra, Barnes, Highcock, Schekochihin and Cowley2011) concludes that, for a given tokamak configuration, there is an optimal amount of perpendicular momentum injection to suppress turbulence. While tokamak plasmas may become turbulent when shear flow parallel to the magnetic field is too large, the centrifugal mirror does not have this problem because the flow is strictly perpendicular.

#### A.2 Experimental

Burrell (Reference Burrell1997) also includes an extensive review of experimental results that demonstrate shear flow suppression of turbulence. This work only considers tokamaks, but lists experiments like JET and TFTR that demonstrate neoclassical diffusion in the core with shear flow, and like DIII-D and JT-60 U that show similar behaviour across the entire minor radius. Experiments in JET have demonstrated that higher ion temperature gradients can be achieved with larger shear flows, and show good correlation to several gyrokinetic codes (Mantica *et al.* Reference Mantica, Strintzi, Tala, Giroud, Johnson, Leggate, Lerche, Loarer, Peeters and Salmi2009). More recent work on DIII-D has demonstrated access to a stable ‘super-H’ mode with large stored energy, primarily due to sheared flow (Knolker *et al.* Reference Knolker, Evans, Snyder, Grierson, Hanson, Jaervinen, Jian, McClenaghan, Osborne and Paz-Soldan2021).

The review from Volosov (Reference Volosov2009) also provides a discussion of experimental results, which show good agreement with the MHD theory underpinning shear flow stabilization in centrifugal mirrors. Results from MCX demonstrated suppression of the interchange instability with rotational speeds an order of magnitude faster than interchange growth times (Ghosh *et al.* Reference Ghosh, Elton, Griem, Case, DeSilva, Ellis, Hassam, Lunsford and Teodorescu2006).

## Appendix B. Symbols, subscripts and superscripts

The following table presents a list of symbols, subscripts, and superscripts which have been used in this work. The table includes the symbol, it's definition, and the units (if applicable).