2.1. Fokker–Planck transport and plasma heating using CQL3D-m and Freya

Mirror transport during high-performance operation where MHD and trapped particle flute-like modes as well as kinetic modes such as the drift-cyclotron loss-cone (DCLC) instability are not present is dominated by axial transport from collisions and neutral dynamics. To simulate transport under these conditions, a Fokker–Planck equation solver including plasma heating and fuelling models is employed. This work used the CQL3D-m Fokker–Planck code and its internal neutral beam deposition Monte Carlo solver Freya.

The CQL3D-m code (Harvey, Petrov & Forest Reference Harvey, Petrov and Forest2016) used here is a modified version of the well-known CQL3D Fokker–Planck code (Harvey & McCoy Reference Harvey and McCoy1992), and older versions of CQL3D-m were used for the design of WHAM and BEAM devices (Endrizzi et al. Reference Endrizzi, Anderson, Brown, Egedal, Geiger, Harvey, Ialovega, Kirch, Peterson and Petrov2023; Forest et al. Reference Forest, Anderson, Endrizzi, Egedal, Frank, Furlong, Ialovega, Kirch, Harvey and Lindley2024). Much like CQL3D, CQL3D-m solves the bounce-averaged Fokker–Planck equations for a coupled set of plasma species $\alpha , \beta , \ldots$ . The Fokker–Planck equation for plasma species $\alpha$ is

(2.1) \begin{equation} \frac {{\rm d}f_{\alpha 0}}{{\rm d}t} = C(f_{\alpha 0},f_{\beta 0},\ldots ) + Q(f_{\alpha 0},\boldsymbol{E}) + R(f_{\alpha 0}) + S_{+}(f_{\alpha 0}) - S_{-}(f_{\alpha 0}), \end{equation}

where $f_{\alpha 0}$ is the bounce-averaged distribution function of species $\alpha$ at the midplane, where $z=0$ (values taken at the midplane are denoted with subscript $0$ ), $f_{\beta 0}$ is the distribution function of species $\beta$ etc., $C(f_{\alpha 0},f_{\beta 0},\ldots )$ is the bounce-averaged nonlinear form of the Rosenbluth collision operator (Rosenbluth, MacDonald & Judd Reference Rosenbluth, MacDonald and Judd1957) considering collisions between all species in the simulation using their fully numerical distribution functions, $Q(f_{\alpha 0}, \boldsymbol{E}_{RF})$ is the bounce-averaged quasilinear diffusion coefficient associated with RF electric field $\boldsymbol{E}$ arising from either RF heating or an instability such as the DCLC, $R(f_{\alpha 0})$ is the radial diffusion (not considered in this work as losses from radial diffusion are ideally much smaller than axial losses in the end plug), and $S_+$ and $S_-$ are particle source and sink terms, respectively, arising from mirror losses, neutral fuelling and charge-exchange.

The CQL3D-m discretizes (2.1) in two velocity-space dimensions, $u = \gamma v$ the momentum-per-rest-mass (where $v$ is particle velocity and $\gamma$ is the relativistic factor) and $\vartheta$ the velocity-space pitch angle, as well as one spatial dimension $\psi$ the square-root poloidal flux (a radial-like coordinate in the magnetic mirror), and uses the implicit time advance solver described in Kerbel & McCoy (Reference Kerbel and McCoy1985). The coefficients $Q$ , $C$ , $R$ , $S_\pm$ in (2.1) have been bounce-averaged in magnetic mirror geometry assuming zero-orbit-width. The bounce-average of some coefficient $A$ is defined as

(2.2) \begin{equation} \langle \langle A \rangle \rangle = \frac {1}{\tau _b}\oint \frac {{\rm d}\ell _b}{\left |v_\parallel \right |} A \Bigg |_{\mu _m, \vartheta _0, \psi }. \end{equation}

and is performed over the equilibrium fields (vacuum fields are only used in the first iteration between CQL3D-m and Pleiades).

The mirror loss boundaries in CQL3D-m are calculated by the following method. Particles starting with a certain $\vartheta _0$ and $u_0$ at the midplane can be trapped by either the magnetic mirror force or the mirror’s ambipolar potential. These two trapping mechanisms have distinct loss boundaries that can be computed by the following methods. Potential trapping occurs when for species $\alpha$ , the quantity $v_p^2(z)$ becomes less than zero at some point in $z$ along the flux surface. The definition of this quantity is

(2.3) \begin{equation} v_p^2(z) = 2\frac {q_\alpha }{m_\alpha } (\phi _{p0} - \phi _p(z)), \end{equation}

where $\phi _p(z)$ is the ambipolar potential (Pastukhov Reference Pastukhov1974; Cohen et al. Reference Cohen, Rensink, Cutler and Mirin1978; Khudik Reference Khudik1997). This can be used to define a potential trapping boundary where

(2.4) \begin{equation} v_{\phi _p}(z) = \sqrt {-v_p(z)^2}. \end{equation}

For a species where $v_p^2 \gt 0$ along the entire flux surface in $z$ , no potential trapping takes place. Similarly, the mirror trapping boundary for species $\alpha$ with $v_0$ can be computed with equations

(2.5a) \begin{align} v_{\parallel , \textrm{mirr}} &= v_{0} \sin {\vartheta _{\textrm{mirr}}},\\[-10pt]\nonumber \end{align}

(2.5b) \begin{align} \vartheta _{\textrm{mirr}} &= \arcsin {\left (\sqrt {\frac {B(z)}{B_0} \frac {v_0^2}{v_0^2 + v_p(z)^2}}\right )_{\textrm{max},z}}, \end{align}

where the quantity in parentheses is evaluated at its maximum value in $z$ along a poloidal flux surface. If $\phi _p(z)$ is taken to be zero everywhere, this becomes the standard expression for the mirror loss cone. Components of the distribution function $f$ which have velocities less than the outermost loss boundary are removed by a sink term,

(2.6) \begin{equation} S_-^{{loss}}(f) = \frac {2v_\parallel }{L}f \end{equation}

where $L$ is the length of the mirror from mirror throat to mirror throat. The loss boundaries are then recomputed for the next time step based on an updated $\phi _p(z)$ profile calculated with the method below.

The CQL3D-m accounts for ambipolarity preservation by directly calculating the ambipolar potential response (Pastukhov Reference Pastukhov1974; Cohen et al. Reference Cohen, Rensink, Cutler and Mirin1978; Khudik Reference Khudik1997). This is done through an iterative relaxation solver which determines the potential required to meet the quasineutrality condition $n_e = \sum _{si} Z_{si} n_{si}$ . The ambipolar potential at a given $z$ location along a flux surface is computed by iterating two different calculation methods. The first method is an iterative solver which takes

(2.7) \begin{equation} \phi _p^{(n)}(z) = \phi _p^{(n-1)}(z) + f_{\textrm{rlx}} T_{e0} \frac {\sum _{si} Z_{si} n_{si}(z)-n_e(z)}{n_e(z) + \sum _{si} Z_{si} n_{si}(z)}, \end{equation}

where $\phi _p^{(n)}(z)$ is the updated ambipolar potential at the new time step, $\phi _p^{(n-1)}$ is the ambipolar potential at the previous time step, $f_{\textrm{rlx}} \approx 0.5{-}1.0$ is a relaxation factor added to improve numerical stability and $n_s$ refers to the local value of density along $z$ for species $s$ . If only (2.7) is used to calculate $\phi _p(z)$ , the code will eventually become unstable. To counteract this, the potential is occasionally relaxed with a Boltzmann-like response,

(2.8) \begin{equation} \phi _p(z) = T_{e0} \ln {\left (\frac {n_e(z) + \sum _{si} Z_{si} n_{si}(z)}{n_{e0} + \sum _{si} Z_{si} n_{si0}}\right )}, \end{equation}

which has been found empirically to ensure stability over many time steps. The generation of $\phi _p(z)$ leads to off-midplane trapping due to phenomena such as sloshing ions obtained by skew injected neutral beams (Kesner Reference Kesner1980) or sloshing electrons from directed electron cyclotron heating (ECH) (Baldwin & Logan Reference Baldwin and Logan1979). Such trapping is not included in a typical bounce-averaged formulation of the Fokker–Planck equation used by CQL3D-m. To account for these trapped populations, a Maxwellian component is added to the local off-midplane distribution function with a temperature corresponding to the local potential well depth if an insufficient density can be sourced from the remapped midplane distribution function.

The neutral beam solver Freya is unchanged from the typical version implemented in CQL3D and described in Harvey & McCoy (Reference Harvey and McCoy1992) and Rosenberg et al. (Reference Rosenberg2004). Neutral particles are injected using realistic neutral beam geometry and tracked until they are ionized in the plasma. The particle source list from Freya is used to calculate a bounce-averaged particle source of ions and electrons with the neutral beam velocity. A bounce-averaged particle sink corresponding to NBI charge exchange is also formulated. Similarly, the computation of the quasilinear diffusion coefficients from ray data in CQL3D-m works in effectively the same way as it did in CQL3D (Harvey & McCoy Reference Harvey and McCoy1992) with a few adjustments to the bounce-averaging procedure to account for mirror geometry.

An example of a CQL3D-m simulation of WHAM showing a sloshing ions distribution is shown in figure 2. This provides an example of many of the new features that are present in CQL3D-m. The plasma is fuelled by a $45^\circ$ NBI producing a sloshing-ions distribution function with strong density peaking about the off-midplane fast ion turning points at $B/B_0 \approx 2$ . The peaking causes a corresponding off-axis potential peak to enforce quasineutrality as the electrons are still effectively Maxwellian which in turn generates a midplane centred potential well which traps lower-energy warm-ions near the midplane.

Figure 2. An example of WHAM plasma profiles and distribution functions calculated with CQL3D-m. Panel (a) is a plot from CQL3D-m with logarithmic contours of the ion distribution in arbitrary units for the innermost $\sqrt {\psi _n} = 0.01$ normalized square root poloidal flux surface with the loss boundary shown in red. Panel (b) are the electron (solid red) and ion (dashed red) density profiles as well as the ambipolar potential (blue) versus axial distance along the mirror $z$ .

2.2. The MHD equilibrium using Pleiades

In order to calculate the equilibrium magnetic fields of the high beta plasmas under consideration in this work, a custom fork of the Green-function-based magnetic field solver code Pleiades introduced in Peterson et al. (Reference Peterson2019) was developed. Pleiades was used in a number of previous works on mirrors (Egedal et al. Reference Egedal, Endrizzi, Forest and Fowler2022; Endrizzi et al. Reference Endrizzi, Anderson, Brown, Egedal, Geiger, Harvey, Ialovega, Kirch, Peterson and Petrov2023; Forest et al. Reference Forest, Anderson, Endrizzi, Egedal, Frank, Furlong, Ialovega, Kirch, Harvey and Lindley2024), but the anisotropic equilibrium solver and self-consistent coupling found in this work represents a significant evolution of the code (these works included diamagnetic corrections to the magnetic equilibrium due to the plasma pressure but used a different method to solve for the corrections, and the corrections were not made self-consistently). This version of Pleiades solves for the MHD momentum balance in flux coordinates and outputs standard format eqdsk files containing the MHD equilibrium. The MHD momentum balance in equilibrium is defined as

(2.9) $$\boldsymbol{J} \times \boldsymbol{B} = \nabla \cdot\underline{\underline {\boldsymbol{P}}} ,$$

where $\boldsymbol{J}$ is the current density, $\boldsymbol{B}$ is the magnetic field and $\underline {\underline {\boldsymbol{P}}}$ is the anisotropic pressure tensor. To solve (2.9) using this technique, the divergence of the anisotropic pressure tensor is rewritten in component form,

(2.10) \begin{equation} \boldsymbol{J} \times \boldsymbol{B} = \boldsymbol \nabla p_\perp + (p_\parallel - p_\perp )\hat {\boldsymbol{b}} \boldsymbol\cdot \boldsymbol \nabla \hat {\boldsymbol{b}} + \boldsymbol{B} \boldsymbol \cdot \boldsymbol \nabla \left (\frac {p_\parallel - p_\perp }{B} \right ) \hat {\boldsymbol{b}}, \end{equation}

where $\hat {\boldsymbol{b}} = \boldsymbol{B}/|B|$ . Taking the cross product of $\boldsymbol{B}$ with both sides of (2.10) after some vector algebra yields

(2.11) \begin{equation} \left [1 + \frac {\mu _0 (p_\perp + p_\parallel )}{B^2} \right ] \boldsymbol{J}_\perp = \frac {\boldsymbol{B}}{B^2} \times \left [\boldsymbol \nabla p_\perp + \frac {\boldsymbol \nabla B}{B} (p_\parallel - p_\perp ) \right ]. \end{equation}

To solve the above equation in Pleiades, the correct choice of flux coordinates must be made such that derivatives of $p$ and $B$ are smooth and can be taken accurately. For an axisymmetric mirror, the most convenient choice of orthogonal coordinates for these purposes is the canonical Clebsch coordinate system ( $\psi$ , $\phi$ , $s$ ) (D’haeseleer et al. Reference D’haeseleer, Hitchon, Callen and Shohet1991). In this coordinate system, the poloidal flux $\psi$ is the radial-like coordinate, the poloidal angle $\phi$ is the symmetry direction and the distance along the field line $s$ is directed in the direction of the magnetic field $\hat {\boldsymbol{b}}$ . Applying these coordinates to the problem dramatically simplifies (2.11) providing an equation for the diamagnetic current response $J_{\textrm{dia}} \equiv J_\perp \equiv J_\phi$ given the spatial profiles of kinetic $p$ and $B$ ,

(2.12) \begin{equation} J_{\textrm{dia}} = \left [1 + \frac {\mu _0 (p_\perp + p_\parallel )}{B^2} \right ]^{-1} \left (\frac {1}{B} \frac {\partial p_\perp }{\partial \psi }\nabla _\perp \psi + \frac {p_\parallel - p_\perp }{B^2} \frac {\partial B}{\partial \psi } \nabla _\perp \psi \right ). \end{equation}

With (2.12), Pleiades can be used to calculate the magnetic equilibrium by the following method: first, Pleiades is run without the plasma included; the vacuum fields and poloidal flux, $\boldsymbol{B}_{\textrm{vac}}(R,Z)$ and $\psi _{\textrm{vac}}(R,Z)$ , produced by the current sources magnetic field coils are solved for using Green’s functions corresponding to the magnetic field coils’ currents. After the vacuum fields are determined, the anisotropic pressures versus $(\psi , s)$ are obtained from integrating the local distribution functions evaluated by CQL3D-m and mapped onto their corresponding $(\psi , s)$ in Pleiades. The $(\psi , s)$ grid is preferred here over a $(\psi , B)$ grid for the flux functions as it was found to be more numerically robust. The value of $J_{\textrm{dia}}$ is calculated using (2.12) at each location and the grid centred area $\Delta A (\psi , s)$ is multiplied by the local value of $J_{\textrm{dia}}(\psi , s)$ to provide a current source $I_{\textrm{dia}}(\psi , s)$ . Green’s functions associated with the $I_{\textrm{dia}}(\psi , s)$ current sources are calculated and used along with previously calculated Green’s functions associated with the magnetic coils to compute updated profiles of $\psi$ and $\boldsymbol{B}$ which account for diamagnetism. The anisotropic pressure and its derivatives are then mapped to the updated $\psi$ flux surface locations and the process is repeated until the change in $\psi$ over an iteration drops below some prescribed value $\epsilon = \Delta \psi / \psi$ . For $\epsilon = 10^{-6}$ , this typically takes ${\sim}10$ iterations. This flux advance and current redeposition scheme is very similar to the magnetic equilibrium solution scheme used in non-inductive current drive calculations with the tokamak code ACCOME (Devoto et al. Reference Devoto, Blackfield, Fenstermacher, Bonoli, Porkolab and Tinios1992). As the edge pressure profile is not solved for in the current version of the RealTwin, an exponential pressure drop at the plasma edge with an $e$ -folding length equal to the ion Lamor radius $\rho _i$ is used. If the value of $\epsilon$ is made small enough and a sufficiently refined ( $\psi$ , $s$ ) grid is used in the computation, the solution to the magnetic equilibrium obtained with this method will converge to the free-boundary pressure balance solution. While the Green’s function approach used by Pleiades is simpler to implement than a free-boundary Grad–Shafranov solver, it is more computationally expensive than Green’s functions corresponding to the diamagnetic current response must be recalculated each time $\psi$ is advanced within the code. However, as each Green’s function computation is independent, the diamagnetic currents and their associated Green’s functions can be computed in parallel, so a refined $(\psi , s)$ grid can be used without incurring long computational times.

The magnetic equilibrium solver in Pleiades has been quite successful in practice and is robust to high values of $\beta \lesssim 1$ . An example of a magnetic equilibrium in WHAM solved using Pleiades appears in figure 3 and shows the anticipated outward radial expansion of the flux surfaces consistent with a diamagnetic reduction in central magnetic field $B_0 \approx B_{\textrm{vac}}\sqrt {1-\beta }$ . The paraxial equilibrium condition from Ryutov et al. (Reference Ryutov, Berk, Cohen, Molvik and Simonen2011), also shown in this figure, is likewise well satisfied (though not perfectly as it is only approximate). Stability against the firehose and the mirror instabilities based on conditions (Grad Reference Grad1967)

(2.13) \begin{equation} \sigma \gt 0, \end{equation}

(2.14) \begin{equation} \frac {\partial }{\partial B}\left ( \sigma B\right ) \gt 0, \end{equation}

where $\sigma = 1/\mu _0 + (p_\perp - p_\parallel )/B^2$ , is checked at simulation runtime. Simulations which are unstable to these MHD modes or do not converge due to $\beta \gt 1$ are discarded. (An interesting note is that Pleiades simulations unstable to the mirror mode do not readily converge. This is due to the equilibrium bifurcation phenomenon described in Kotelnikov, Bagryansky & Prikhodko (Reference Kotelnikov, Bagryansky and Prikhodko2010) which causes the iterative equilibrium solver in Pleiades to fail.)

Figure 3. An example of a magnetic equilibrium for a high $\beta$ plasma in WHAM calculated with Pleiades. Panel (a) are the flux contours for the vacuum fields $\psi _{\textrm{vac}}$ in blue and the flux contours after diamagnetic evolution $\psi _{\textrm{new}}$ in red. Panel (b) are the kinetic pressure profiles and magnetic pressure profiles as well as the paraxial equilibrium condition in green which in paraxial equilibrium is equal to $B_{\textrm{vac}}^2/2\mu _0$ .

3.1. The POPCON for tandem mirrors

To generate tandem mirror POPCONs, the steady-state 0-D tandem mirror transport equations based on those found in Mirin et al. (Reference Mirin, Auerbach, Cohen, Gilmore, Pearlstein and Rensink1983) and Hua & Fowler (Reference Hua and Fowler2004) are taken to the steady state limit and all derivatives are taken to zero. However, consistent with the typical POPCON formulation, radial variation in plasma parameters (including radial variation in axial confinement) has been allowed. It is assumed that the tandem mirror is operating a classical tandem without thermal barriers (Fowler & Logan Reference Fowler and Logan1977). Under this assumption electrons in the central cell and end plug are in thermal contact and $T_e \equiv T_{ec} = T_{ep}$ . This assumption dramatically simplifies the system of equations which must be solved to make a POPCON. Under these assumptions, the primary energy and particle balance equations in the POPCON system are:

(3.1a) \begin{align} &\underline {\text{density}} \notag \\ &2 \pi \ell _c \int _0^{a_c} \frac {n_c}{\tau _c} \, r\mathrm{d}r = \Gamma _{\text{gas}}, \end{align}

(3.1b) \begin{align} &\underline {\text{ion energy}} \notag \\ &2 \pi \ell _c \int _0^{a_c} \left ( n_c\frac {\phi _i+T_{ic}}{\tau _c} + \mathcal{P}_{ie} -\mathcal{P}_{\text{fus},i} \right ) \, r\mathrm{d}r = P_{\text{RF},i}, \end{align}

(3.1c) \begin{align} &\underline {\text{electron energy}} \notag \\ &2 \pi \ell _c \int _0^{a_c} \left ( n_c\frac {\phi _e+T_{e}}{\tau _c} - \mathcal{P}_{ie} -\mathcal{P}_{\text{fus},e} + \mathcal{P}_{\textrm{rad}} \right ) \, r\mathrm{d}r = P_{\text{RF},e}, \end{align}

where subscripts $c$ and $p$ denote central cell and end plug quantities, respectively, subscripts $i$ and $e$ denote electrons and ions, respectively, $n_c \equiv n_{ec} = n_{ic}$ is the density assuming quasineutrality, $T_s$ is the temperature of species $s$ , $\ell _c$ is the central cell length, $a_c$ is the central cell radius, $\Gamma _{\text{gas}}$ is the gas fuelling rate, $\tau _c \approx \tau _p \approx \tau _E$ is the confinement time, $\mathcal{P}_{ie}$ is the collisional power density transfer from ions to electrons, $P_{\textrm{RF},s}$ is the applied RF heating power to the central cell, $\mathcal{P}_{\textrm{fus},s}$ is the fusion power density from alpha slowing down applied to species $s$ and $\mathcal{P}_{\textrm{rad}}$ is the power density radiated by bremsstrahlung and synchrotron radiation.

In order to create a set of POPCONs, (3.1) is solved for heating $P_{\textrm{RF}}$ (with some fraction directed at ions and electrons, respectively) and fuelling $\Gamma _{\textrm{gas}}$ at a given operating point with fixed $T_{ic}$ and $n_{c}$ . To solve this system, all variables must be written as a function of $T_{ic}$ and $n_c$ and fixed end plug and central cell parameters: $B_m$ , $a_m$ , $n_p$ , $\ell _c$ , $B_{0c}$ ( $T_e$ may be solved for as a dependent variable). This is done by applying the following definitions; the confinement time due to classical radial and axial transport (the dominant transport modes expected in the tandem mirror) $\tau _c$ is defined as

(3.2) \begin{equation} \tau _c = \left ( \frac {1}{\tau _{\text{Past}} + \tau _f} + \frac {1}{\tau _{\rho }}\right )^{-1}, \end{equation}

where the ambipolar confinement of ions in the central cell by the end plug potential is defined (Pastukhov Reference Pastukhov1974; Cohen et al. Reference Cohen, Rensink, Cutler and Mirin1978) as

(3.3) \begin{equation} \tau _{\text{Past}}=\frac {\sqrt {\pi }}{2}\tau _{ii}\frac {\phi _i}{T_{ic}}\exp \left (\frac {\phi _i}{T_{ic}}\right )\frac {G(R_{mc})}{1+T_i/2\phi _i-(T_{ic}/2\phi _i)^2}, \end{equation}

with $\tau _{ii}$ the central cell ion–ion pitch-angle scattering time and $\phi _i$ the ion confining potential between the central cell and the end plug resultant from their density difference,

(3.4) \begin{equation} \phi _i = T_{e} \ln \left (\frac {n_p}{n_c}\right ), \end{equation}

and function $G(R_{mc})$ of the central cell mirror ratio $R_{mc} = B_m / (B_{0c} \sqrt {1-\beta _c})$ defined as

(3.5) \begin{equation} G(x) = \sqrt {1+x^{-1}} \ln \left ( \frac {\sqrt {1+x^{-1}} + 1}{\sqrt {1+x^{-1}} - 1} \right ). \end{equation}

The confinement due to collisional trapping of ions, which can become important for large $\ell _c$ , is defined (Rognlien & Cutler Reference Rognlien and Cutler1980) as

(3.6) \begin{equation} \tau _f=\sqrt {\pi }R_{mc}\frac {\ell _c}{v_{\textrm{th}ic}}\exp \left (\frac {\phi _i}{T_{ic}}\right ), \end{equation}

where $v_{\textrm{th}ic} = \sqrt {T_{ic}/2 m_{ic}}$ is the ion thermal velocity in the central cell. The confinement time associated with classical radial transport is defined as

(3.7) \begin{equation} \tau _\rho = 0.25 \left (\frac {a_c}{\rho _{ic}}\right )^2 \tau _{ii}, \end{equation}

where $\rho _{ic} = v_{\textrm{th}ic}/\Omega _{ic}$ is the central cell ion Larmor radius. Returning to (3.1), the electron confining potential in the plug can be defined for classical tandem operation by solving the transcendental equation (Fowler & Logan Reference Fowler and Logan1977; Hua & Fowler Reference Hua and Fowler2004),

(3.8) \begin{equation} \frac {\phi _e}{T_e} \exp \left (\phi _e/T_e\right ) \approx \sqrt {\frac {m_i}{m_e}}\left (\frac {T_{ic}}{T_e} \right )^{3/2}\left (\frac {\phi _i}{T_{ic}}\right )\exp (\phi _i/T_{ic}). \end{equation}

The power density associated with ion to electron energy transfer is

(3.9) \begin{equation} \mathcal{P}_{ie} = -\mathcal{P}_{ei} = 3 \frac {m_e}{m_i} \frac {n_c}{\tau _{ei}}\left (T_{ic}-T_e\right ), \end{equation}

where $\tau _{ei}$ is the ion–electron collision time. The fusion power density is calculated for a deuterium–tritium plasma with alphas depositing their energy on species $s$ is calculated with

(3.10) \begin{equation} \mathcal{P}_{\text{fus},s} = f_{\alpha s}(x)\frac {E_\alpha }{4}n_c^2\langle \sigma v\rangle , \end{equation}

using the reactivity coefficients $\langle \sigma v\rangle$ from Bosch & Hale (Reference Bosch and Hale1992) and equation for the fraction of power deposited on each species $s$ by alpha particles from the analytic solution to the slowing down equation (Ott Reference Ott1997),

(3.11a) \begin{equation} f_{\alpha e}(x) = x^{-2}\left \{x^2+\frac {1}{3}\ln \left [\frac {(x+1)^2}{x^2-x+1}\right ]-\frac {2}{\sqrt {3}}\tan ^{-1}\left [\frac {1}{\sqrt {3}}(2x-1)\right ]-\frac {\pi }{3\sqrt {3}}\right \}, \end{equation}

(3.11b) \begin{equation} f_{\alpha i}(x) = 1 - f_{\alpha e}(x), \end{equation}

where $x=\sqrt {E_\alpha /E_{\text{crit}}}$ and the “critical” energy where the ion and electron energy deposition by the $\alpha$ particles is balanced $E_{\text{crit}}=[3\sqrt {\pi }/(4m_i\sqrt {m_e})]^{2/3}m_\alpha T_e=59.2 T_e /\mu ^{2/3}$ where $\mu = m_i/m_p$ . Finally, the radiated power density is the sum of the power densities from bremsstrahlung and synchrotron radiation. The radiated power densities, using the expression for synchrotron radiation assuming an infinite straight cylinder with finite wall reflectivity, are (Trubnikov Reference Trubnikov1979)

(3.12a) \begin{equation} \mathcal{P}_{\textrm{rad}} = \mathcal{P}_{\textrm{synch}} + \mathcal{P}_{\textrm{brem}}, \end{equation}

(3.12b) \begin{equation} \mathcal{P}_{\textrm{brem}} = 5.34\times 10^{-32} n_e^2 \sqrt {T_e}\,[\textrm{W}\, \textrm{m}^{-3}] ,\end{equation}

(3.12c) \begin{equation} \mathcal{P}_{\textrm{synch}} = 4.14\times 10^{-5} \sqrt {(1-R_w)\frac {n_e}{a_c}}\left (1+2.5\frac {T_e}{m_ec^2}\right )(T_e B_{0c}\sqrt {1-\beta _c})^{2.5}\,[\textrm{W}\,\textrm{m}^{-3}], \end{equation}

where $R_w$ is the wall reflectivity (taken to be 80 %, the same value used TGYRO tokamak transport simulations (Candy et al. Reference Candy, Holland, Waltz, Fahey and Belli2009)) and units of $T_e$ and $B_{0c}$ are taken to be in kiloelectronvolts and tesla, respectively.

Using this POPCON power-balance, it was possible to reproduce the tandem mirror operating points identified in Fowler et al. (Reference Fowler, Moir and Simonen2017) when the same assumptions were made, i.e. no classical transport (3.7), fusion heating (3.10) or radiation (3.12), and all central cell RF heating was taken to be electron directed. The $T_{ic} = 150\,\textrm{keV}$ , $\beta _c \sim 1$ operating points described in Fowler et al. (Reference Fowler, Moir and Simonen2017) were found to within a small factor of error in $T_i$ (<1.25) that was determined to be associated with a discrepancy in the fusion cross-sections used in their model and the Bosch–Hale cross-sections used here, they considered $T_{ic} \equiv T_{ec}$ when in reality $T_{ic} \sim T_{ec}$ at these parameters, as well as the sometimes poor convergence of root-finding algorithm used to evaluate the model here at very high values of $\beta _c$ . When the additional terms used here were included it was found that the operating points in Fowler et al. (Reference Fowler, Moir and Simonen2017) could not be attained. Radiation and classical transport significantly reduced performance (both terms become large at the very high $\beta \sim 1$ operating points considered there, though, it should be noted that the expression for classical transport used here is not valid when $\beta \sim 1$ and Fowler et al. (Reference Fowler, Moir and Simonen2017) assumed very high wall reflectivities $R_w$ ).

3.2. Identifying operating points with model end plugs

With the components of the POPCON system of (3.1) defined, it is possible to solve for the required heating and fuelling at each location in ( $n_c$ , $T_{ic}$ ) space for a given set of end plug radial density profiles and engineering parameters. The POPCONs with simplified end plugs were solved to establish target parameters for the end plug to be simulated with RealTwin. In these simplified end plugs, the density profile was taken to be flat, the ion directed fraction of $P_{\textrm{RF}}$ was taken to be $80\%$ . End plug parameters were taken to be fixed besides the plasma radius through the mirror throat $a_m$ . The central cell field was increased until either the number of alpha gyroradii across the plasma radius $N_\alpha \gt 5$ defined using

(3.13) \begin{equation} N_\alpha = \frac {a_{0c}}{\rho _\alpha } \approx 5.24 a_m \sqrt {B_m B_{0c}}, \end{equation}

was obtained or the central cell vacuum mirror ratio reached 8 (this limit was imposed in the interest of simplifying control of central cell interchange driven by trapped particle modes based on the results in Gerver et al. (Reference Gerver1989)). Average end plug density was varied until ignited operating points appeared. However, for illustrative purposes the POPCONs in this section all have the same density.

The results of the POPCON analysis described here are shown in figure 4. The POPCON analysis indicates that, for what will be shown later on are fairly reasonable end plug parameters, it is possible to create a classical tandem mirror with an ignited central cell. Ignition regions are present in all cases when the plasma radius at the mirror throat is larger than 0.15 m. These results indicate that a tandem mirror pilot plant designed with these principles should be able to satisfy the NASEM fusion pilot plant requirements at lengths under 100 m. Using the simplified equation for the net electric power:

(3.14) \begin{equation} P_{\textrm{net,elec}} = \left (P_{\textrm{fus},\alpha } + P_{\textrm{NBI}} + C_{\textrm{mult}} P_{\textrm{fus},n} \right )\eta _{\textrm{ele}} - \frac {P_{\textrm{NBI}}}{\eta _{\textrm{NBI}}}. \end{equation}

Figure 4. The POPCON plots for tandem mirrors with $\ell _c = 50\,\textrm{m}$ at four different plasma radii through the mirror throat, (a) 0.1 m, (b) 0.15 m, (c) 0.2 m and (d) 0.25 m, using $n_p = 1.5\times 10^{20}\,\textrm{m}^{-3}$ and $B_{0c} = 3.125\,\textrm{T}$ in all cases except for (a) where $B_{0c} = 5.0\,\textrm{T}$ . Operating points are shown as a function of central cell density $n_c$ divided by end plug density $n_p$ versus central cell ion temperature $T_i$ . Blue and red filled contours are of heating power to the central cell required at a given ( $n_c$ , $T_i$ ) operating point (blue regions indicate ignition). Regions which are off the scale for positive or negative and are inaccessible are denoted with red and blue hatching, respectively. Black contour lines are fusion power in the central cell. Green contour lines are central cell vacuum $\beta$ .

Assuming neutral beam efficiencies $\eta _{\textrm{NBI}} = 60 \%$ , an electrical conversion efficiency $\eta _{\textrm{ele}} = 50 \%$ assuming the use of a Brayton cycle (Schleicher, Raffray & Wong Reference Schleicher, Raffray and Wong2001), a blanket neutron energy multiplication factor $C_{\textrm{mult}} = 1.1$ and $P_{\textrm{NBI}} = 30\,\textrm{MW}$ into the end plugs; it was found a tandem mirror would need to produce $157.4\,\textrm{MW}$ of power to satisfy the NASEM $50\,\textrm{MWe}$ requirement. Such an operating point is achievable with $B_m = 25\,\textrm{T}$ , $a_m = 0.15\,\textrm{m}$ , $\langle n \rangle _p = 1.5\times 10^{20}\,\textrm{m}^{-3}$ , $\ell _c = 50\,\textrm{m}$ , $\beta _c \cong 0.6$ , $n_c/n_p \cong 0.55$ , $T_i \cong 45\,\textrm{keV}$ and $T_e \cong 125\,\textrm{keV}$ based on this POPCON analysis. The fusion power density of this case is $\sim 3.5\,\textrm{MW m}^{-3}$ which is competitive with high-field tokamaks (Frank et al. Reference Frank2022b ). Performance of such a system could be further enhanced by large $a_m$ or with technologies such as direct electrical conversion (Moir & Barr Reference Moir and Barr1973) or more efficient NBI utilizing advanced neutralization techniques and energy recovery (Fumelli, Jequier & Paméla Reference Fumelli, Jequier and Paméla1989; Grisham Reference Grisham2009; Kovari & Crowley Reference Kovari and Crowley2010). Values of $a_m \lt 15\,\textrm{cm}$ have difficulty obtaining good performance without very high plug densities $\langle n_p \rangle \gg 2.0\times 10^{20}\,\textrm {m}^{-3}$ due to radial central cell transport and enhanced synchotron radiation. It is noteworthy that good performance is achievable here despite using generally more conservative parameters than the DT classical tandem proposed in Fowler et al. (Reference Fowler, Moir and Simonen2017) which had $B_m = 24\,\textrm{T}$ , $a_m = 0.21\,\textrm{m}$ , $\langle n_p \rangle = 2.6\times 10^{20}\,\textrm{m}^{-3}$ , $\ell _c = 55\,\textrm{m}$ , $\beta _c \sim 1.0$ , $n_c/n_p \sim 0.35$ , $T_i = 150\,\textrm{keV}$ and $T_e = 150\,\textrm{keV}$ .

This POPCONs analysis exposed a number of surprising features of fusion relevant tandem mirror operation that were not previously well known. These discoveries are in part because, outside of the analysis in Fowler et al. (Reference Fowler, Moir and Simonen2017), axisymmetric tandem mirrors with high performance HTS magnets and high $R_m$ had never been considered. Older studies of tandem mirrors generally assumed end plug mirror ratios of order unity to account for minimum-B coil geometries and lower peak $B_m$ (Fowler & Logan Reference Fowler and Logan1977). This analysis builds upon that in Fowler et al. (Reference Fowler, Moir and Simonen2017), and very importantly, it includes the effect of fusion alpha particle heating on the electron temperature. The key result here is that, when high-field mirror coils are available, classical tandem central cells can and should be operated at much higher density fractions $n_c/n_p$ than previously assumed. Much of the previous classical tandem mirror literature used $n_c/n_p = 0.05{-}0.35$ (Dimov et al. Reference Dimov, Zadkaidakov and Kishinevsky1976; Fowler & Logan Reference Fowler and Logan1977; Hua & Fowler Reference Hua and Fowler2004). The analysis here suggests this value should be higher, $n_c/n_p = 0.35{-}0.6$ . Strong alpha heating of the electrons obtained by, increasing the central cell density, boosting confinement with large $R_m$ from axisymmetrization and directly heating the central cell ions with RF to enhance the fusion rate, allows the central cell to reach electron temperatures ${\sim}2 {-} 3\,T_i$ without supplemental electron heating as shown in figure 5. This operational mode takes advantage of tandem mirror confinement (3.3) improving more rapidly with increasing temperature than it degrades with increasing the density ratio.

Figure 5. Contours of $T_e$ for operating points at given $\langle n \rangle$ and $E_{\textrm{NBI}}$ satisfying the constraint equations (4.5). This plot uses the $a_m = 0.15$ m POPCONs case found in figure 4(b). Black contour lines are fusion power in the central cell in megawatts.

In conclusion, this scoping shows that for high-field axisymmetric mirrors, classical tandem operation can reach appreciable Q and is a valid alternative to operation with thermal-barriers (Baldwin & Logan Reference Baldwin and Logan1979). These findings are in agreement with similar analyses by Pratt & Horton (Reference Pratt and Horton2006) and Fowler et al. (Reference Fowler, Moir and Simonen2017). However, it has been found here that higher density operation than these previous analyses is possible thanks to alpha heating, somewhat relaxing the end plug performance criteria, but there are still some important caveats. The tandem mirror POPCON approach used here, much like the tokamak POPCON approach, is only able to identify operating points. It does not guarantee that an operating point it identifies is stable or achievable in actual operation. To determine this, a time dependent simulation which takes the mirror from an initial breakdown plasma to its final operating point must be performed. It is possible that, to actually reach the operational conditions proposed here, a transient heating source, such as a high-power ECH “sparkplug”, might be needed to boost $T_e$ initially. This would drive up the end plug potentials and improve confinement enough for the central cell to generate sufficient alpha heating power to reach ignition. Furthermore, the effects of spatial gradients or turbulent transport were not included in this analysis. Generally, the spatial density gradients in tandem mirrors are predicted to be small (Mirin et al. Reference Mirin, Auerbach, Cohen, Gilmore, Pearlstein and Rensink1983; Hua & Fowler Reference Hua and Fowler2004), and the effects of electron temperature gradient turbulence, the main form of turbulent energy transport predicted to be present in tandem mirrors, are anticipated to be more tolerable than the temperature gradient driven turbulent transport in tokamaks (Hua & Fowler Reference Hua and Fowler2004; Pratt & Horton Reference Pratt and Horton2006). These topics will be treated in future work where gyrokinetic turbulence simulations at mirror central cell conditions will be performed to evaluate the growth rates of gradient driven instabilities. A time dependent solver based on those found in Mirin et al. (Reference Mirin, Auerbach, Cohen, Gilmore, Pearlstein and Rensink1983) and Hua & Fowler (Reference Hua and Fowler2004) in order to evaluate the accessibility of these operating points under more realistic assumptions will also be implemented. In addition to these confinement issues, macrostability of the tandem system, specifically its stability against trapped particle modes and MHD, have not been addressed here and will be the topic of a future publication.

One potentially interesting feature of the relatively high $n_c/n_p$ ratio used in this study is the question of whether outflow from the central cell may have a stabilizing effect on kinetic modes in the end plugs. The flux from the central cell through the end plug midplane can be determined by particle balance of the outflow along the flux tube:

(3.15) \begin{equation} \Gamma _{c,\textrm{out}} = \frac {n_c \ell _c}{2 \tau _c} \left (\frac {R_{mc}}{R_{mp}} \right ). \end{equation}

If this outflow flux becomes large enough, it has been shown to be stabilizing to DCLC (Drake et al. Reference Drake, Casper, Clauser, Coensgen, Correll, Cummins, Davis, Foote, Futch, Goodman, Grubb, Hornady, Nexsen, Simonen and Stallard1981), but large amounts of central cell outflow could also reduce an end plug’s classical performance by increasing the ion–ion collisional scattering rate. Taking the estimate for the outflow required to stabilize the plug from Correll et al. (Reference Correll1980), we have

(3.16) \begin{equation} \Gamma _{p, \textrm{stab}} \gt 1.2 \times 10^5 \frac {n_p \phi _e^{3/2}}{(R_{mp}\sin ^2\theta _{\textrm{NBI}}-1) \sqrt {\mu } E_{ip}} \left (\frac {a_m\sqrt {R_{mp}}}{\rho _{i0}}\right )^{-4/3}\,[\textrm{m}^{-2}\,\textrm{s}^{-1}], \end{equation}

where $\theta _{\textrm{NBI}}$ is the angle of the NBI in the plug, $\mu$ is the average ion mass in the plug in AMU, $\phi _e$ is the plug potential in kilovolts, $E_{\textrm{NBI}}$ is the plug NBI energy in kiloelectronvolts (all other quantities are assumed to be in MKS) and $\rho _{i0}$ is ion Larmor radius at the midplane of the plug. Using these two equations, the minimum $n_c/n_p$ density ratio at which DCLC stabilization in the plug can be obtained:

(3.17) \begin{equation} \frac {n_c}{n_p} \gt 2.4 \times 10^5 \frac {\tau _c \phi _e^{3/2}}{\ell _c (R_{mp}\sin ^2\theta _{\textrm{NBI}}-1) \sqrt {\mu } E_{ip}} \left (\frac {R_{mp}}{R_{mc}} \right ) \left (\frac {a_m\sqrt {R_{mp}}}{\rho _{i0}}\right )^{-4/3}. \end{equation}

For tandem mirror parameters similar to the ignited operating points in figure 4(b): $\tau _c = 5\,\textrm{s}$ , $\phi _e = 400\,\textrm{keV}$ , $\ell _c = 50\,\textrm{m}$ , $R_{mp} = 8$ , $R_{mc} = 13.3$ (accounting for finite $\beta \sim 0.6$ in both mirror ratios), $\theta _{\textrm{NBI}} = 45^\circ$ , $E_i = 500\,\textrm{keV}$ , $\mu = 2.5\,\textrm{AMU}$ and $a_m \sqrt {R_{mp}}/\rho _{i0} = 25$ , it is calculated that $n_c/n_p \gt 1000$ indicating DCLC likely cannot be stabilized with central cell outflow in the tandem mirrors under consideration here. As $\Gamma _{c,\textrm{out}}$ is insufficient for DCLC stabilization, DCLC stability in § 4 will be ensured by the suppression of the diamagnetic drift wave.

Finally, there is the matter of alpha particle removal from the central cell in steady-state operation. Alpha particles are well confined in the central cell by the end plug’s ambipolar potential. Without a method to deconfine them, they will accumulate. This will need to be addressed with a technique such as drift pumping, which removes thermalized alphas by applying a RF perturbation to the background magnetic field that causes them to drift out of the plasma after enough central cell transits (Dimov et al. Reference Dimov, Zadkaidakov and Kishinevsky1976; Logan Reference Logan1983).

4.1. Assumptions and 0-D constraints

When performing simulations with the RealTwin model, several assumptions about end plug operation were made due to the limitations of present mirror modelling tools (a discussion of ways in which these capability gaps may be closed in the future appears in § 5). In addition, a number of constraints related to engineering feasibility and plasma stability were used to do significant optimization work prior to even running the RealTwin.

Beginning with the modelling assumptions and fixed parameters; two large physics assumptions were made in this study to overcome the limitations of present modelling tools. First, it was assumed that the end plug, for the purposes of equilibrium and transport modelling, could be treated as a standalone simple mirror with certain correction factors. The justification of the first assumption lies in the fact that the ion populations in the plug and the central cell are to lowest order “detached”. Ions confined in the central cell do not enter the plugs and ions confined in the plugs do not enter the central cell. Similarly ions that are not confined in the plug are also not confined in the central cell and vice versa. This property of tandem mirrors ends up being useful for reasons that will be discussed momentarily. In fusion-relevant conditions, confinement is good in both the plugs and the central cell and the unconfined population of ions flowing through either region is small, acting as a (relatively) low-density stream which passes outwards until it reaches the expander region (Cohen et al. Reference Cohen, Bernstein, Dorning and Rowlands1980). However, in a classical tandem mirror the electron population is able to transit between the plug and central cell, and alpha particle heating from the central cell is sufficient to enable us to reach high electron temperatures. In fact, the predicted electron temperatures in a tandem are higher than they otherwise would be in a standalone end plug due to alpha heating of electrons in the central cell. The initial simulations in § 4.2 of the standalone simple mirror end plug utilize a dynamic electron temperature. However, later simulations in § 4.3 use a fixed values of $T_e$ and ambipolar potential $\phi _e$ corresponding to the values computed in the tandem POPCON.

The second important physics assumption is that MHD (Rosenbluth & Longmire Reference Rosenbluth and Longmire1957) and trapped particle mode (Berk & Lane Reference Berk and Lane1986 a,b) were not considered in detail. While kinetic stability to the DCLC instability and MHD stability to the mirror and firehose instabilities are considered here, stability against flute modes and rigid $m=1$ displacements driven by MHD and trapped particle modes have not been assessed. While it is acknowledged that stability to these modes is vital to developing an operating mirror fusion device, in order to effectively assess stability and develop working stabilization actuators, self-consistent mirror equilibrium parameters, like those calculated here, must be obtained first. This work encompasses that important first step. The IPS-driven RealTwin simulation framework built here was constructed for the express purpose of coupling transport, heating and equilibrium simulations self-consistently to stability simulations. An early demonstration of this capability, in which CQL3D-m/Pleiades simulations are coupled to hybrid vector particle-in-cell (hybrid-VPIC) simulations in order to assess extended MHD and kinetic stability. This is presented in this article’s companion article, Tran et al. (Reference Tran2025). Realta Fusion plans to continue to investigate kinetic and MHD stability utilizing a full complement of simulation tools coupled to IPS driven simulations including: linear MHD (FLORA) (Cohen, Freis & Newcomb Reference Cohen, Freis and Newcomb1986), nonlinear MHD (NIMROD) (Sovinec et al. Reference Sovinec, Glasser, Gianakon, Barnes, Nebel, Kruger, Schnack, Plimpton, Tarditi and Chu2004), particle-in-cell methods (VPIC, WARP-X) (Bird et al. Reference Bird, Tan, Luedtke, Harrell, Taufer and Albright2021; Fedeli et al. Reference Fedeli2022; Le et al. Reference Le, Stanier, Yin, Wetherton, Keenan and Albright2023) as well as gyrokinetic and multifluid codes (Gkeyll) (Hakim Reference Hakim2008; Francisquez et al. Reference Francisquez, Rosen, Mandell, Hakim, Forest and Hammett2023). These simulation tools can also account for extended effects: such as charge separation’s impact on trapped particle modes (Kesner & Lane Reference Kesner and Lane1983), finite Larmor radius effects’ impact on MHD stability (Ryutov et al. Reference Ryutov, Berk, Cohen, Molvik and Simonen2011) and sloshing ions as well as neutral particle fuelling’s impact on DCLC (Baldwin Reference Baldwin1977; Tran et al. Reference Tran2025). As stability actuators and these effects can interact in ways which are nonlinear and difficult to predict analytically, high-fidelity simulation tools and self-consistent simulations are vital to developing a complete picture of stability.

This simulation model will eventually allow a systematic evaluation of axisymmetric mirror stability to trapped particle modes, MHD and kinetic modes, using realistic equilibrium parameters to be performed under the influence of stability actuators such as end-ring induced shear flows (Sakai, Yasaka & Itatani Reference Sakai, Yasaka and Itatani1993; Beklemishev et al. Reference Beklemishev, Bagryansky, Chaschin and Soldatkina2010; Bagryansky et al. Reference Bagryansky2011), localized ECH-induced shear flows (Cho et al. 2005, Reference Cho2006), active feedback control (Lieberman & Wong Reference Lieberman and Wong1977; Kang, Lieberman & Sen Reference Kang, Lieberman and Sen1988), high- $\beta$ (Berk, Wong & Tsang Reference Berk, Wong and Tsang1987) field ripple in the central cell in combination with conformal conducting walls (Li, Kesner & LoDestro Reference Li, Kesner and LoDestro1987), kinetic stabilizers (Post Reference Post2001) and line-tying using cold plasma mantels (Fornaca, Kiwamoto & Rynn Reference Fornaca, Kiwamoto and Rynn1979; Molvik et al. Reference Molvik, Breun, Golovato, Hershkowitz, McVey, Post, Smatlak and Yujiri1984). In addition to MHD modes, there is interest in the impact of these actuators on trapped particle modes in the tandem mirror. Kinetic stabilizers (Post Reference Post2001) have been shown to be potentially ineffective at stabilizing trapped particle modes (Berk & Pratt Reference Berk and Pratt2011; Fowler et al. Reference Fowler, Moir and Simonen2017). However, as trapped particle modes in tandem mirrors cause fluid-like interchange which has similar growth rates and is phenomenologically similar to MHD, it may be possible to stabilize them with the same feedback control, line tying and shear flow actuators as MHD (Kesner Reference Kesner2018). Another potential avenue for trapped particle mode stabilization is the difference in the electron and ion bounce points in the central cell which can have a stabilizing effect on trapped particle modes in classical tandem mirrors like those here (Kesner & Lane Reference Kesner and Lane1983).

As mentioned earlier, the detached nature of the end plug and central cell ion distributions can be used to improve performance. In a tandem, the end plug ion composition need not be the same as the central cell as confined ion populations in the end plug and central cell are not exchanged (Fowler & Logan Reference Fowler and Logan1977; Cohen et al. Reference Cohen, Bernstein, Dorning and Rowlands1980). This presents a unique opportunity to enhance end plug physics performance and simplify end plug engineering. To this end, end plugs fuelled with only tritium are used in conjunction with a deuterium and tritium fuelled central cell. This combination maximizes end plug performance and reduces the nuclear engineering challenges presented by the end plug. The physical argument for this operational mode is as follows: ion confinement in the end plug is dominated by two processes: ion–ion scattering and electron drag. The particle confinement times associated with these two processes are (Logan, Mirin & Rensink Reference Logan, Mirin and Rensink1980)

(4.1a) \begin{align} \tau _s &= 2.76\times 10^{18}\mu ^{1/2} E_p^{3/2} \frac {\log _{10}(R_m \sin ^2\theta _{\text{NBI}})}{Z^4 n_i\ln \Lambda _{ii}} ,\\[-12pt]\nonumber \end{align}

(4.1b) \begin{align} \tau _d &= 1\times 10^{14}\mu T_{e,p}^{3/2}\frac {\ln (E_{\textrm{NBI}} / \langle E_L \rangle )}{Z^2 n_e\ln \Lambda _{ei}}, \\[-12pt]\nonumber \end{align}

(4.1c) \begin{align} \tau _p &\approx \left (1/\tau _s + 1/\tau _d \right )^{-1}, \end{align}

where the average ion loss energy $E_L$ can be defined as

(4.2) \begin{equation} E_L \approx \frac {1}{1+(\tau _s/\tau _d)}\left [E_{\textrm{NBI}}+E_h\left (\frac {\tau _s}{\tau _d}\right )\right ], \end{equation}

with ambipolar hole loss energy,

(4.3) \begin{equation} E_h =\frac {\phi _e}{R_m\sin ^2\theta _{\text{NBI}}-1}. \end{equation}

At typical tandem operational parameters, confinement is assumed to be electron drag dominated, $\tau _d \ll \tau _s$ (Logan et al. Reference Logan, Mirin and Rensink1980; Mirin et al. Reference Mirin, Auerbach, Cohen, Gilmore, Pearlstein and Rensink1983; Hua & Fowler Reference Hua and Fowler2004; Fowler et al. Reference Fowler, Moir and Simonen2017). As the end plug ion population composition need not be the same as the central cell, performance can be enhanced by choosing ions with large $\mu /Z^2$ . In addition, choosing an ion mix to suppress fusion neutron rates in the end plug can substantially simplify engineering by reducing radiation shielding requirements for components such as magnets and NBI. Suppressing fusion rates in the end plugs is acceptable as the total fusion power produced by the end plug is a small component of the overall fusion power in a tandem, since the end plug volume is much less than central cell volume, $V_p \ll V_c$ . If taken to its logical conclusion, this argument suggests that a tritium-only end plug is optimal with respect to classical transport and neutron damage. Tritium has the highest $\mu /Z^2$ ratio of any standard fusion fuel ion, and the tritium–tritium fusion cross-section is nearly an order of magnitude lower than the deuterium–deuterium fusion cross-section, and tritium–tritium fusion produces neutrons with energies of $\sim 3\,\textrm{MeV}$ (Brown et al. Reference Brown2018). Thus, a tritium-only end plug provides a factor of $\sim 1.2$ increase in confinement performance and a many-order of magnitude reduction in the neutron rate as well as a significant reduction in neutron energy relative to a deuterium–tritium end plug (tritium fuelled plugs will have inferior microstability to deuterium or hydrogen plugs as the tritium Larmor radius is larger, making it easier to excite diamagnetic drift waves. However, it will be shown that it should be possible to avoid diamagnetic drift waves with large $a/\rho _i$ even with tritium plugs).

Another assumption is the use of only negative-ion neutral beams for end plug heating in steady-state (a discussion about the limitations associated with high-harmonic RF heating coupled with positive-ion NBI, similar to what is planned in WHAM, and the reason it is not used in this work is discussed in Appendix A). The use of ECH has also been limited to achieving breakdown (assuming a scheme like that found in Yakovlev et al. (Reference Yakovlev, Bagryansky, Gospodchikov, Shalashov and Solomakhin2016)). Presently, the use of any significant quantity of ECH in CQL3D-m causes the potential profile to invert. This behaviour is actually promising as it indicates that it may be possible to generate thermal barriers (Baldwin & Logan Reference Baldwin and Logan1979), but at the moment it causes significant problems with the bounce-averaging and ambipolar potential calculations in CQL3D-m. Thus, further investigation of ECH is left for future work. Neutral transport into the tandem mirror plasma has also been neglected. This effect is anticipated to be of some importance in small relatively low-density mirrors such as WHAM, where neutrals recycled from the edge can enter the plasma and charge exchange with well-confined beam ions causing them to be lost and replacing them with a promptly lost cold ion. However, fusion-grade plasmas like those under consideration here are opaque to edge neutrals, and the effect of these neutrals on transport can be ignored. The final assumption made here is that transport in the end plug will remain classical during operation. To this end, it has been insured that the DCLC stability condition is well satisfied in the plasmas optimized here.

In addition to physics assumptions, some engineering constraints were imposed on the mirror design here. First, the field of the two large-bore high-field mirror coils placed at either end of the end plug (denoted HF) was fixed at $25\,\textrm{T}$ . Feasibility studies of HTS coils at these fields have found that such designs should be viable at bores sizes appropriate for fusion applications, $\sim 1$ m warm bore diameter. A single low-field coil (denoted LF) was placed at the centre of the end plug to give control of the central field. The magnetic configuration for an example case is shown in figure 6. The minimum length $\ell$ (where subscript $p$ on end plug quantities in this section has been dropped for simplicity), defined as the mirror throat to mirror throat distance. This limit was chosen not by the finite Larmor radius stabilization condition for $m \geqslant 2$ modes, but rather by engineering constraints. The substantial size of the neutral curvature central cell in addition to the good curvature expanders is expected to provide significant stabilization against high $m$ modes (Roberts & Taylor Reference Roberts and Taylor1962; Rosenbluth, Krall & Rostoker Reference Rosenbluth, Krall and Rostoker1962; Rosenbluth Reference Rosenbluth1965; Pearlstein Reference Pearlstein1966; Anikeev et al. Reference Anikeev, Bagryansky, Ivanov, Kuzmin and Salikova1992; Ryutov et al. Reference Ryutov, Berk, Cohen, Molvik and Simonen2011; Fowler et al. Reference Fowler, Moir and Simonen2017). Active stabilization schemes are anticipated to provide for remaining MHD stabilization needs. Here, $\ell \geqslant 4.5\,\textrm{m}$ has been chosen as with values of $\ell \lt 4.5\,\textrm{m}$ require large reverse polarity currents in the LF coil to cancel out the HF field. Additionally, very short end plugs will almost certainly impede NBI port access and require the use of inefficient NBI launch angles that limit beam path length, causing reduced NBI absorption. The optimal NBI configuration for maximizing beam path length, and maintaining good confinement (noting there is an effective mirror ratio reduction associated with NBI launch angles smaller than $90^\circ$ relative to the magnetic axis) occurs at, or very close to, a neutral beam with $\theta _{\textrm{NBI}} = 45^\circ$ targeted at the midplane. To maintain an effective $45^\circ$ ion velocity-space angle for different beam launch angles, the beam’s target location must be moved to different values of $B/B_0$ . The reduction factor in beam path length incurred by this move is in a simple mirror is $\sqrt {(1+\cot ^2 \theta _{\textrm{NBI}})/4}/\sin \theta _{\textrm{NBI}}$ . To simplify this study, only the near optimal case of $\theta _{\textrm{NBI}} = 45^\circ$ with the beam directed at the midplane was examined in detail (the NBI here used a $0.4\,\textrm{m}$ by $0.8\,\textrm{m}$ source with no beam divergence as a placeholder until more realistic NBI specifications are known). In later work, with better defined, beam, magnet and radiation shielding parameters, this angle may be slightly adjusted to prevent beam-line interference with components.

Figure 6. An example of the magnetic configuration used in the RealTwin end plug simulations. The high-field (HF) mirror coils are shown in red and the low-field (LF) central coil is shown in black. Contours of poloidal flux are shown in blue and the limiting surface is shown with a black dashed line.

After fixing some engineering parameters, a basic physics constraint-based optimization was used to identify a parameter space which seemed likely to produce attractive results. To do this, the primary parameters of an negative-ion NBI driven plug, $a_m$ , $\ell$ , $B_m$ , $B_0$ , $E_{\textrm{NBI}}$ , $\langle n \rangle$ , $B_0$ and $P_{\textrm{NBI}}$ were defined. From these parameters, most other major parameters in the plasma can be defined. Fixed parameters were $B_m = 25\,\textrm{T}$ and $\ell = 4.5\,\textrm{m}$ . From the POPCON analysis in § 3, the acceptable values of $\langle n \rangle$ were known to range from $1 \times 10^{20}\,\textrm{m}^{-3}$ to $2 \times 10^{20}\,\textrm{m}^{-3}$ . By fixing $\langle n \rangle$ and $E_{\textrm{NBI}} \sim E_i$ to some desired values, it is possible to presume that the electron energy is $\sim 0.1 E_{\textrm{NBI}}$ which Egedal et al. (Reference Egedal, Endrizzi, Forest and Fowler2022) and Endrizzi et al. (Reference Endrizzi, Anderson, Brown, Egedal, Geiger, Harvey, Ialovega, Kirch, Peterson and Petrov2023) have shown is roughly accurate for an NBI heated simple mirror plasma. The beam path was estimated to be $\sim 2\sqrt {2} a_m$ and the plug volume was taken to be $V \sim 2/3 \pi R_m a_m^2 \ell$ where the factor of 2/3 is added based on simulation results to account for the field curvature reducing the volume from a simple cylinder. With these assumptions in place, it was possible to constrain the rest of the mirror’s parameters by simultaneously solving the system,

(4.4a) \begin{align} &\underline {\text{NBI absorption}} \notag \\ &\exp {\left (-\int n \sigma _s \,\text{d}l\right )} \lt 0.2, \\[-10pt]\nonumber\end{align}

(4.4b) \begin{align} &\underline {\text{DCLC stability}} \notag \\ &\frac {a_0}{\rho _i} \gtrsim 25 ,\end{align}

for $a_m$ and $B_0$ . The NBI absorption condition (4.4a ) uses the analytic NBI stopping cross-sections $\sigma _s$ from Janev, Boley & Post (Reference Janev, Boley and Post1989) and (4.4b) is a slightly modified version of the high- $\beta$ drift-cyclotron wave stability criterion which appears in Tang, Pearlstein & Berk (Reference Tang, Pearlstein and Berk1972) and Baldwin (Reference Baldwin1977). According to these conditions, cool plasma trapped by sloshing ions or high- $\beta$ and large plasma size can be used to stabilize a plasma against DCLC. In the end plug here, stabilization of DCLC will come not from sloshing ions, but primarily from large plasma size and high- $\beta$ as the POPCONs in § 3 predict that $T_e$ will be too high for effective sloshing ion trapping of cool plasma due to $\tau _{ii} \ll \tau _{d}$ . (This stabilization method has been demonstrated experimentally in the LAMEX mirror (Ferron & Wong Reference Ferron and Wong1984).) Using the assumptions outlined above, (4.4) may be rewritten (taking (4.4b) to be equal to 25) as

(4.5a) \begin{align} &\underline {\text{DCLC stability}} \notag \\ &B_0^2 - 0.16\left (\frac {\mu E_{i,100\,\textrm{keV}}}{B_m a_m^2} \right ) - 3.0 E_{i,100\, \textrm{keV}} \langle n_{20} \rangle = 0, \\[-5pt]\nonumber \end{align}

(4.5b) \begin{align} &\underline {\text{NBI absorption}} \notag \\ &6.7\times 10^{-21} = \langle n_{20} \rangle \frac {\sqrt {1+\cot ^2(\theta _{\textrm{NBI}}/4)}}{\sin \theta _{\textrm{NBI}}} a_m \sqrt {\frac {B_m}{B_0 (1-\beta )^{1/2}}} \sigma _s, \end{align}

where ion energies are in hundreds of kiloelectronvolts and densities are in $10^{20}\,\textrm{m}^{-3}$ (other units are standard MKS). At high values of $\langle n \rangle$ and $E_{\textrm{NBI}}$ , the constraint on $B_0$ set by $\beta$ rather than DCLC stability. In this case $B_0$ was constrained with condition $\beta \leqslant 0.5$ . Once the constraint equations were solved, from the simple mirror confinement time scaling law $n_{20}\tau _p = 0.25 E_{\textrm{NBI},100\, \textrm{keV}}^{3/2} \log _{10}(R_m)$ from Egedal et al. (Reference Egedal, Endrizzi, Forest and Fowler2022) (this scaling law is for deuterium and does not include all the details of the scaling laws presented in (4.1), but is sufficient for the constraint exercise here), one can calculate the upper bound on the absorbed $P_{\textrm{NBI}}$ required to reach the target value of $\langle n \rangle$ using fuelling rate balance:

(4.6) \begin{equation} P_{\textrm{NBI}} = \frac {E_{\textrm{NBI}} \langle n \rangle V}{\tau _p}. \end{equation}

From the POPCON analysis performed in § 3, it was identified that plasma radius $a_m \gtrsim 0.15\,\textrm{m}$ and density $\langle n \rangle \sim$ 1.0–1.5 $\times 10^{20}\,\textrm{m}^{-3}$ . Plotting contours of $P_{\textrm{NBI}}$ and $a_m$ versus $E_{\textrm{NBI}}$ and $\langle n \rangle$ produces a plot like that shown in figure 7. This plot indicates an operating region accessible with reasonable $P_{\textrm{NBI}}$ values $\lt 20\,\textrm{MW}$ exists between $E_{\textrm{NBI}}$ of $200\,\textrm{keV}$ and $300\,\textrm{keV}$ . However, the method used to constrain NBI absorption with (4.4a) is suboptimal. It has not accounted for that, despite not being absorbed as readily at high energies, the beam may still be more efficient at fuelling, and simplified scaling laws were used here rather than an actual Fokker–Planck calculation. Furthermore, the electron temperature in this analysis is implicitly taken to be $\sim 0.1 E_{\textrm{NBI}}$ in accordance with the scaling law from Egedal et al. (Reference Egedal, Endrizzi, Forest and Fowler2022). This is lower than the values predicted by the POPCONs in § 3. These discrepancies will be corrected in the next parts of this section with more detailed simulations using the RealTwin.

Figure 7. Contours of $a_m$ (red lines) and $P_{\textrm{NBI}}$ (filled contours) for operating points at given $\langle n \rangle$ and $E_{\textrm{NBI}}$ satisfying the constraint equations (4.5). The grey region represents operation regimes requiring greater than 20 MW of NBI power to access.

4.2. Parametric simulations of standalone simple mirror performance

With the constraints in the previous section established, a parametric scan in $B_0$ , $P_{\textrm{NBI}}$ , $E_{\textrm{NBI}}$ and $a_m$ was carried out with the RealTwin model described in § 2 was specified. This scan used the end plug parameters shown in table 1 to simulate a standalone simple mirror. For this work, points in the scan were selected over a linear range for each of the parameters. This scoping method is coarse, but substantial constraints were already placed on the problem through the use of the POPCON analysis in § 3 and the analytic arguments used to constrain the physics and engineering earlier in this section. This made the coarse high-fidelity scan sufficient to deliver an optimized simple mirror that exceeded the end plug requirements identified by the POPCON analysis. A scan in $\ell$ was also performed to determine if scaling power density with $\ell$ versus $P_{\textrm{NBI}}$ was substantially different (changing the aspect ratio of the plug can slightly change NBI path length and modify NBI absorption). Varying $\ell$ , however, had effectively the same impact on power density and plug performance as varying $P_{\textrm{NBI}}$ . In the rest of the discussion here will focus on the highest performing $\ell = 4.5\,\textrm{m}$ case.

Table 1. Mirror parameter ranges used in the parametric scan in § 4.2 and in the machine learning aided optimization in § 4.3.

The CQL3D-m simulations performed in this section used 24 flux surfaces evenly spaced from 0.01 to 0.9 in normalized square-root poloidal flux $\sqrt {\psi _n}$ , and the axial $z$ grid used to formulate the ambipolar response along the field line had 64 evenly spaced points extending from the midplane to the mirror throat. The distribution functions were discretized over 500 points in momentum-per-rest-mass $u_0$ with a grid which used increased resolution at low energies as to accurate capture the ion distributions and 400 evenly spaced points in the pitch angle $\vartheta$ . The maximum energy in the momentum-per-rest-mass grid was set to correspond to a $250\,\textrm{keV}$ electron energy. The CQL3D-m simulations used 800 time steps. The time advance was initially set to $1\,\textrm{ms}$ for the first 200 time steps, then increased to $2\,\textrm{ms}$ for the next 200 time steps, and then further increased to $4\,\textrm{ms}$ for the remaining 400 time steps providing an overall simulation runtime of $2.2\,\textrm{s}$ , several $\tau _p$ . Ramped time stepping was used to maintain numerical stability, as despite being an implicit time-advance code, the iterations of CQL3D-m with Freya and Pleiades and the internal ambipolar potential solve are explicit. The Freya NBI deposition calculation was run at each CQL3D-m time step using $2.5\times 10^6$ particles, and the Pleiades equilibrium was updated after every 100 time steps.

The results of the simulations using the RealTwin model are shown in figure 8. Simulations which failed to satisfy the DCLC constraint (4.4b ) or in which the equilibrium failed to converge due to $\beta \gt 1$ or the onset of the mirror instability were removed from this analysis (this is why there is an empty region in the figure 8 plots; the onset of the mirror instability as the source of a $\beta$ limit in mirror plasmas will be discussed further in a future publication). A potentially surprising result of these scans is that there is little difference in peak density versus $E_{\textrm{NBI}}$ , however, remembering the analysis in § 4.1, this makes sense. Neutral beam absorption drops with higher $E_{\textrm{NBI}}$ , but the NBI fuelling efficiency increases $\propto E_{\textrm{NBI}}^{1/2}$ . Another observation is that $\langle n \rangle$ improves superlinearly with increasing $P_{\textrm{NBI}}$ . This is due to a twofold nonlinearity in the NBI absorption equation (4.4a) with NBI power. As $P_{\textrm{NBI}}$ increases, the $\langle n \rangle$ increases and therefore the neutral beam absorption fraction increases. In addition, the plasma $\beta$ also increases as $P_{\textrm{NBI}}$ increases. This causes the plasma to expand outwards from diagmagnetism, increasing the NBI absorption path length, further increasing NBI absorption.

Figure 8. Contour plots of $\langle n \rangle$ for different values of $a_m$ and $B_0$ at three different $E_{\textrm{NBI}}$ for an end plug with fixed length $\ell = 4.5\,\textrm{m}$ with three different NBI heating powers $10\,\textrm{MW}$ (a), $15\,\textrm{MW}$ (b) and $20\,\textrm{MW}$ (c).

As anticipated in § 4.1, the target end plug $\langle n\rangle$ values identified in § 3 were achieved at the lowest $P_{\textrm{NBI}}$ for $E_{\textrm{NBI}}$ of $240\,\textrm{keV}$ and $360\,\textrm{keV}$ . Sufficient end plug conditions for fusion relevant central cell operation were achievable with $B_0$ of $4 - 6\,\textrm{T}$ , $P_{\textrm{NBI}} \sim 15\,\textrm{MW}$ , $E_{\textrm{NBI}} \sim 240-360\,\textrm{keV}$ and $a_m$ of $0.15 - 0.2\,\textrm{m}$ . Under these conditions effectively full NBI absorption (>95 %) was obtained with a $240\,\textrm{keV}$ beam, so lowering NBI energy further would only serve to reduce performance. The plasma profiles of the most optimal mirror identified in this scan are shown in figure 9. The mirror had parameters $a_m = 0.15\,\textrm{m}$ , $B_0 = 4\,\textrm{T}$ , $E_{\textrm{NBI}} = 240\,\textrm{keV}$ and $P_{\textrm{NBI}} = 15\,\textrm{MW}$ and achieved $\langle n \rangle \gt 1.5\times 10^{20}\,\textrm{m}^{-3}$ , more than sufficient to drive a $50\,\textrm{m}$ long central cell to $Q \gt 5$ conditions based on the analysis in § 3.

Figure 9. Density $n_p$ and average ion energy $\langle E_i \rangle$ , electron energy $\langle E_e \rangle \equiv 1.5 T_e$ and ambipolar potential $e\phi$ , versus the square root of the normalized poloidal flux $\sqrt {\psi _n}$ . The profiles come from a simulation with parameters $a_m = 0.15\,\textrm{m}$ , $B_m = 25\,\textrm{T}$ , $B_0 = 4\,\textrm{T}$ , $\ell = 4.5\,\textrm{m}$ , $E_{\textrm{NBI}} = 240\,\textrm{keV}$ and $P_{\textrm{NBI}} = 15\,\textrm{MW}$ .

4.3. Optimization for a self-consistent end plug operating condition

While the end plug-like simple mirror identified in § 4.2 appears to be sufficient to deliver a viable central cell, the scan those conditions were derived from was coarse. It also did not include the effects of increased electron temperature in the end plug as a result of alpha particle heating from the central cell. To address this, a secondary set of RealTwin simulations were created. These simulations used a modified CQL3D-m model in an attempt to better capture the tandem mirror physics. The electron distribution was assumed to be fixed by the classical tandem mirror central cell dynamics and the CQL3D-m simulations in this section used a fixed $T_e$ Maxwellian electron background where $n_e=n_i$ and a fixed ambipolar potential $\phi _e$ corresponding to those computed in the tandem mirror POPCON using (3.1) and (3.8). The operating point for the tandem mirror was chosen to be at a thermally stable ignited central cell (i.e. operating on the right-hand side of the ignited region) operating point with $\beta _c \sim 0.6$ with the same plug density as the RealTwin simulations. The POPCON and the RealTwin simulations were iterated until a self-consistent condition where $T_e$ , $\phi _e$ and $\langle n_p \rangle$ in the POPCON and RealTwin models matched. These RealTwin simulations consistent with the POPCON conditions provide a more detailed picture of steady-state end plug operation in a tandem mirror. Most notably, the larger ambipolar hole in the ion distribution increases the optimal values of $E_{\textrm{NBI}}$ and $B_0$ .

In addition to fixed $T_e$ and $\phi _e$ , these simulations employed several other modifications and refinements on the model versus § 4.2. Rather than the coarse parametric scan used in the previous section, the BO technique described in § 2 was used to efficiently find values of $a_m,$ $B_0,$ $E_{\textrm{NBI}}$ and $P_{\textrm{NBI}}$ which optimize the normalized Lawson criterion $\langle n \tau \rangle (P_{\textrm{abs}}/P_{\textrm{NBI}})$ , where $P_{\textrm{abs}}/P_{\textrm{NBI}}$ is the absorbed fraction of the total NBI power. This metric can be thought of as a measure of NBI fuelling efficiency. The normalized Lawson criteria is the preferred objective function quantity as it simultaneously optimizes density per unit of absorbed NBI power and ensures good beam absorption. Rather than the simple $a_0/\rho _{i0} \gt 25$ DCLC constraint used previously, the $\beta$ adjusted DCLC constraint from Tang et al. (Reference Tang, Pearlstein and Berk1972),

(4.7) \begin{equation} N_\rho ^{\textrm{Tang}} \equiv a_0/\rho _{i0} \approx \frac {1}{2}\frac {1}{\sqrt {\beta }\sqrt {m_e/m_i + \Omega ^2_i /\omega _{pi}^2}}, \end{equation}

was employed in order to provide a more accurate estimate of the DCLC stability. This constraint was implemented automatically in the objective function through the use of penalty function $f(a_0/\rho _i) = 1/2 \tanh {\{[a_0/\rho _i - (N_\rho ^{\textrm{Tang}} - 5)]/5\}}$ which multiplies $\langle n \tau \rangle (P_{\textrm{abs}}/P_{\textrm{NBI}})$ .

The optimization in this section was initially performed using tritium to improve classical confinement. However, the higher $E_{\textrm{NBI}}$ required to overcome the ambipolar hole at high $T_e$ caused end plugs using tritium as the main ion species to not satisfy the DCLC constraint. Because of this, a switch to deuterium was made to reduce $\rho _i$ . The confinement reduction factor resultant from this switch was $\sim \sqrt {3/2}$ rather than $\sim 3/2$ as the elevated $T_e$ in these simulations caused the plasma to be ion–ion scattering dominated rather than electron drag dominated as it was in § 4.2.

In order to recover some classical confinement performance in deuterium operation, the NBI angle $\theta _{\textrm{NBI}}$ was shifted from the value of $45^\circ$ used in § 4.2 to $50^\circ$ , increasing the effective mirror ratio. This small shift was acceptable as $\theta _{\textrm{NBI}}$ should not be perpendicular enough to excite the Alfven ion-cyclotron instability (Smith Reference Smith1984) or the mirror instability (2.14), sloshing ions are not required for DCLC stability, and the NBI path length reduction is small enough that the shift in $\theta _{\textrm{NBI}}$ does not significantly reduce ${P_{abs}}/{P_{\textrm{NBI}}}$ .

End plug simulations, especially those with large $E_{\textrm{NBI}}$ , were found to exhibit an instability tied to $\beta \rightarrow 1$ near the magnetic axis. The higher NBI power density due to smaller flux surface volumes as $\sqrt {\psi _n}\rightarrow 0$ (low $\sqrt {\psi _n}$ surfaces absorb a similar amount of NBI power to large $\sqrt {\psi _n}$ surfaces for an NBI width substantially less than $a_0$ ) caused $\beta$ on the inner flux surfaces to grow quickly. Higher $\beta$ increased $R_m$ and $\tau _p$ . This in turn increased $n$ leading to a further increase to the NBI power density on the inner flux surfaces. This feedback loop rapidly caused $\beta$ on the inner flux surfaces to exceed unity resulting in a loss of magnetic equilibrium. To combat the instability, a radially varying classical diffusion coefficient was added to CQL3D-m’s Fokker–Planck equation in order to flatten the pressure profiles in the end plug. A typical value for this diffusion coefficient is $D_{\rho} \sim 1\times 10^{-4}\,\textrm{m}^{2}\,\textrm{s}^{-1}$ . This addition caused a reduction in performance due to the introduction of radial losses, but classical diffusion successfully resolved the $\beta \rightarrow 1$ instability in most cases and flattened the radial end plug density profiles. An example of the profile flattening from radial diffusion is shown in figure 10.

Figure 10. Plots of density $n$ (blue) and ion temperature $E_i$ (red) versus normalized square root poloidal flux $\sqrt {\psi _n}$ obtained from CQL3D-m simulations of a fixed $\beta =0.6$ end plug with $a_m=0.25\,\textrm{m}$ , $B_m = 25\,\textrm{T}$ , $B_0 = 6\,\textrm{T}$ , $E_{\textrm{NBI}} = 400\,\textrm{keV}$ , $P_{\textrm{NBI}} = 20\,\textrm{MW}$ and $T_e = 120\,\textrm{keV}$ both with (solid) and without (dashed) classical radial diffusion. Radial diffusion has a significant effect on confinement in these cases. Diffusion notably reduces the total density in addition to flattening the density. The $E_i\gt E_{\textrm{NBI}}$ seen in the plot is the result of a combination of up-scattering in ion energy and the preferential loss of low energy ions.

However, when the magnetic equilibrium was evolved self-consistently with Pleiades using the pressure profiles from the CQL3D-m simulations with diffusion, numerical noise developed which could corrupt the solution. The noise appeared to be related to the radial boundary conditions used in CQL3D-m, and there will be attempts to resolve this noise in future work so that fully self-consistent CQL3D-m/Pleiades optimizations of end plug performance may be performed. To combat the noise here, fixed $\beta _p = 0.6$ Pleiades equilibria were used during the optimization to prevent failures which can impair the GPBO’s convergence. To maintain a degree of self-consistency, the optimization imposed a Gaussian envelope penalty function (centred at $\beta _p = 0.6$ with $\sigma = 0.1$ ) multiplying $\langle n \tau \rangle (P_{\textrm{abs}}/P_{\textrm{NBI}})$ to ensure that optimized end plugs had $\beta _p\sim 0.6$ . After optimal end plug conditions were identified, the optimal cases were rerun in standalone simulations which self-consistently included field evolution to confirm that the same results could be obtained.

The CQL3D-m simulations in this section could be evolved for a longer period of time at lower resolution as the explicit potential solver was disabled and electrons were considered as a background species rather than undergoing nonlinear evolution with the ions. These CQL3D-m simulations used 24 flux surfaces in the radial grid even spaced from $\sqrt {\psi _n} = 0.01{-}0.95$ and 64 axial points evenly spaced in $\hat z$ between the midplane and the mirror throat. The distribution function was discretized over 300 points in momentum-per-rest-mass $u_0$ and 200 evenly spaced points in the pitch angle $\vartheta$ . The maximum energy in the momentum-per-rest-mass grid was set to correspond to $1000\,\textrm{keV}$ deuterons. The CQL3D-m used 100 ms time steps and Freya NBI deposition calculations were run at each time step using $2.5\times 10^6$ particles. Simulations were evolved for $10\,\textrm{s}$ , several $\tau _p$ , to time converge the end plug transport. Rather than self-consistently optimize both the POPCON and the end plug in a single simulation loop, optimal end plug conditions were calculated at a number of different values of $T_e$ independently using the BO procedure described in § 2. At each value of $T_e$ , five batches of 10 simulations were run with parameters selected using the CL-min acquisition function. The initial three batches used a more exploratory acquisition function with $\kappa = 10$ and final two batches used a lower $\kappa = 5$ which allowed the optimization to quickly find a global optimum.

As anticipated in § 4.2, when self-consistent electron temperatures were accounted for, $E_{\textrm{NBI}}$ , $B_0$ and $a_m$ had to increase to ensure that the plasma remained below the $\beta$ limit, had good confinement, had adequate NBI absorption and remained stable to DCLC. The optimization procedure selected the largest end plug radii $a_m=25\,\textrm {cm}$ set by engineering feasibility as it both improved $\langle n\tau _p\rangle (P_{\textrm{abs}}/P_{\textrm{NBI}})$ and allowed the plasma to more easily satisfy the DCLC stability condition (4.7). With end plug $T_e = 80{-}120\,\textrm{keV}$ , $E_{\textrm{NBI}} \approx 400\,\textrm{keV}$ and $P_{\textrm{NBI}} \approx 15{-}20\,\textrm{MW}$ , desirable $1.5\times 10^{20}\,\textrm{m}^{-3}$ densities were achievable. These densities were achievable despite the $a_m = 25\,\textrm{cm}$ plugs’ higher volume because of the increased fuelling efficiency obtained with high-energy NBI and the reduction in electron drag from using self-consistent values for $T_e$ . The best end plug performance was obtained at $T_e=100\,\textrm{keV}$ as it provided a good balance of increased ambipolar losses and decreased electron drag. However, for the best overall system performance occurred at $T_e=120\,\textrm{keV}$ , and it was chosen for the optimized operating point as it increased the achievable central cell $P_{fus}$ by $\sim 100\,\textrm{MW}$ . The optimized tandem mirror POPCON operating point based on this optimization process is shown in figure 11 and the parameters of the optimized tandem mirror are given in the second column labelled “Optimum consistent $T_e$ ” in table 2. It is noted that, in future evolutions of this design, slightly higher $\beta$ may be desirable in order to achieve high $\beta$ stabilization using conducting walls (Berk et al. Reference Berk, Wong and Tsang1987; Li et al. Reference Li, Kesner and LoDestro1987) and higher $P_{\textrm{fus}}$ . This may be achieved under the same parameters with elevated $P_{\textrm{NBI}}$ , or can be optimized for in a new optimization. The optimized tandem mirror operating point in table 2 has self-consistent $T_e$ , uses technologically feasible components, obeys the DCLC stability constraint (4.7) (the Tang constraint specified $a_0/\rho _i \gt 27$ for diamagnetic drift wave stability) and satisfies the NASEM requirements for a fusion power core (Hawryluk et al. Reference Hawryluk2021).

Figure 11. A POPCON of the tandem mirror operating point described the second column labelled “Optimum consistent $T_e$ ” (table 2). Heating power $P_{\textrm{RF}}$ is shown in the filled red (positive) and blue (negative, ignited) contours, fusion power $P_{fus}$ in black contours, $\beta _c$ in green contours, and the operating point in the table is marked with a star.

Table 2. Optimized tandem mirror parameters based on the results of parametric scans and POPCON analysis using the end plug conditions obtained using the standalone simulations with CQL3D-m and Pleiades performed in § 4.2 in column “Optimum Standalone” and the CQL3D-m and Pleiades simulations with electron temperature fixed to the central cell value performed in § 4.3 in column “Optimum Consistent $T_e$ ”.