Physics basis for the Wisconsin HTS Axisymmetric Mirror (WHAM)

Abstract

The Wisconsin high-temperature superconductor axisymmetric mirror experiment (WHAM) will be a high-field platform for prototyping technologies, validating interchange stabilization techniques and benchmarking numerical code performance, enabling the next step up to reactor parameters. A detailed overview of the experimental apparatus and its various subsystems is presented. WHAM will use electron cyclotron heating to ionize and build a dense target plasma for neutral beam injection of fast ions, stabilized by edge-biased sheared flow. At 25 keV injection energies, charge exchange dominates over impact ionization and limits the effectiveness of neutral beam injection fuelling. This paper outlines an iterative technique for self-consistently predicting the neutral beam driven anisotropic ion distribution and its role in the finite beta equilibrium. Beginning with recent work by Egedal et al. (Nucl. Fusion, vol. 62, no. 12, 2022, p. 126053) on the WHAM geometry, we detail how the FIDASIM code is used to model the charge exchange sources and sinks in the distribution function, and both are combined with an anisotropic magnetohydrodynamic equilibrium solver method to self-consistently reach an equilibrium. We compare this with recent results using the CQL3D code adapted for the mirror geometry, which includes the high-harmonic fast wave heating of fast ions.

1. Introduction

The Wisconsin high-temperature superconductor axisymmetric mirror (WHAM) is intended to operate in the weakly collisional or classical mirror (CM) regime in which particle confinement is governed by ion–ion pitch angle scattering rather than electron drag on fast ions. While this was typical of the earlier non-axisymmetric experiments (see figure 1), WHAM will extend recent advances from the gas dynamic regime, characterized by strong sonic losses, to the CM regime, where ion pitch angle scattering controls the particle confinement time.

Figure 1. The historical record of fast ion confinement time for beam driven simple mirrors compared with (3.4) and for thermal ions in the tandem configuration.

Axisymmetric linear devices are unstable to magnetohydrodynamic (MHD) interchange modes, with fast growth rates $\gamma _{\rm MHD}$ of the order of an inverse ion bounce time, $\tau_{\rm bounce} = L_p/v_{\rm thi}$ (Fowler 1981; Post 1987), the mirror half-length $L_p$ over the ion thermal velocity. However, experiments and simulations have shown that interchange modes with $m \geq 2$ can be stabilized by the ion's finite Larmor radii (FLR) when the experiment is long enough, and also that the $m=1$ mode can be shear-flow stabilized at low amplitude when the shearing rate exceeds the growth rate by roughly a factor of 2 (Beklemishev et al. 2010; Bagryansky et al. 2011, 2015; Yakovlev et al. 2018; Soldatkina et al. 2020). A combination of a high-mirror-ratio axisymmetric magnet system and electron cyclotron heating of a sufficiently dense plasma has led to electron temperatures of 1 keV and dispelled the myth that electrons are always ‘dead cold’ in mirrors. Separately, the development of high-temperature superconductor (HTS) magnets appropriate for fusion applications (Whyte et al. 2016; Whyte 2019) now allows for larger mirror ratios and better classical confinement, while axisymmetry precludes neoclassical effects. Together, these significant developments motivate this re-examination of fusion mirror devices.

WHAM will demonstrate several key technological milestones that show the viability of fusion mirrors: retiring the engineering risks of constructing circa 20 T magnets; use of electron cyclotron heating for breakdown, localized $T_e(r)$ control and temperatures of $T_e\geq 1$ keV; extending MHD stabilization of axisymmetric plasmas to collisionless environments; and demonstrating high-harmonic fast wave (HHFW) heating on injected neutral beam ions. Perhaps most important, the proper modelling of the physics of these features and benchmarking of computational methods against experimental results is essential to its mission as a development platform for future linear devices. This paper lays out performance predictions for the WHAM experiment, and the methods in this paper will be used in a concurrent paper designing a larger break-even axisymmetric mirror device.

We start with § 2 and give an overview of the engineering design of WHAM, with summaries of the magnets and field geometry (§ 2.1), the vacuum vessel and fuelling system (§ 2.2), the electron cyclotron heating (ECH) system (§ 2.3), the neutral beam injection (NBI) system (§ 2.4), the radio-frequency (RH) HHFW heating system (§ 2.5) and the biased limiter and end-ring system (§ 2.6). Then, in § 3, we discuss the performance modelling of WHAM, beginning first with a general discussion of the mirror physics. Section 3.1 addresses the most pressing issue for all axisymmetric mirror experiments, MHD stability, and briefly discusses FLR effects and vortex stabilization. Section 3.2 explains the several roles of the ECH system: first, in ionizing the target plasma and in driving the transition to the collisionless mirror regime, and in localized heating to produce sheared flow profiles. Section 3.3 looks at the feasibility of fuelling via NBI on WHAM given its limited opacity, focusing on the outsize role of charge exchange collisions in determining the ion distribution function. Section 3.4 discusses the anisotropic equilibrium calculation method using the semi-analytic results from the fast beam-ion solver code (FBIS) (see Egedal et al. 2022). Section 3.5 discusses the anticipated role of kinetic instabilities, as well as enhanced radial transport from turbulence. Section 3.6 looks at the role of fast-ion adiabaticity in limiting device performance at low fields. Finally, § 3.7 concludes by integrating §§ 3.1–3.6 into a holistic upper-limit prediction for WHAM performance, presenting results from both the FBIS simulations, but also recent results from the newly developed code CQL3D-mirror.

With the exception of §§ 3.5 and 3.6, all of these methods are purely classical in nature. These results then represent an absolute best case scenario for a steady-state plasma, one unperturbed by MHD and kinetic instabilities, uninhibited by gradient driven turbulence and unbothered by impurities. Despite the many fewer mirror experiments as compared with tokamaks, these issues have each individually been addressed theoretically and managed experimentally. However, solving them simultaneously remains to be demonstrated at reactor relevant scales. Modelling of each of these of substantial topics will be presented in future papers.

2. Overview of WHAM and subsystems

Basic operation of the WHAM device works as follows: deuterium gas is puffed into the central cell close to the ECH injection point. A total of 1 MW of 110 GHz crosses the 4 T resonant surface with X-mode polarization to ionize a target plasma of radius $a=0.1$ m and density $1-3\times 10^{19}$ m$^{-3}$. This target is next struck with a 10 ms pulse of 25 keV, 40 A, 45$^\circ$-inclined deuterium neutral beam producing a sloshing ion distribution. Biasing of the limiters and rings in the expander region will impose a transverse electric field $\boldsymbol{E}$ which interacts with the magnetic field $\boldsymbol{B}$ to produce an $\bf \textit{E}\times \textit{B}$ drift in the edge, providing the vortex confinement to saturate MHD instabilities at small amplitude. After a year of operation, a 1 MW RF antenna will use second harmonic heating of the sloshing ions, accelerating these particles above their injection energy. Figure 2 shows an engineering schematic of the central region of WHAM for reference.

Figure 2. Overhead section of WHAM with contours of field strength. The colour contours represent the 4 T ECH resonance (yellow), the fast ion turning point at 2$B_0 = 1.72$ T (green) and the expansion ratio for electron thermal confinement $B_{\rm mir}\sqrt {m_e/m_i}$ (blue). Numbered: (1) ECH injection port, (2) NBI beam path, (3) HHFW antenna, (4) limiter.

2.1. Magnets and Field Geometry

The 17 T, 2 kA steady state, 5.5 cm warm bore HTS mirror magnets, custom built by Commonwealth Fusion Systems, are centred at $z=\pm 98$ cm. In addition, two pulsed copper coils, repurposed from the Wendelstein-7A experiment, will be located near the midplane at $z=\pm 20$ cm. They will be driven with Transrex power supplies with a 1 s rise time to increase the central field from 0.32 T to a maximum of 0.86 T. At the highest fields, the compressive force between the HTS magnets will reach 60 tons, with three 304 stainless steel struts on the central cell bearing this load.

As seen in the Gas Dynamic Trap (GDT) experiment (see Bagryansky, Beklemishev & Postupaev 2018), a sufficiently large expansion ratio in the field outside of the mirrors, $R_{\exp }\equiv B_{\rm mir}/B \geq \sqrt {m_i/m_e}$, leads to increased electron thermal confinement (Yakovlev et al. 2018; Soldatkina et al. 2020) when compared with prior mirror experiments. This is due to the parallel ambipolar potential that preserves quasineutrality, trapping the thermal electron population (Boldyrev, Forest & Egedal 2020). In the less collisional WHAM case, the parallel potential drop is predicted to happen in the central cell rather than the expander region (Egedal et al. 2022). Still, the transition from gas dynamic to a CM will benefit from the large expansion ratio by reducing thermal losses from the electrons to the end rings on startup.

2.2. Vacuum Vessel and Fuelling System

The central cell vessel (manufactured by Nor-Cal) is bounded on either end by the HTS mirror magnets, which are in turn bookended by large expansion tank end cells (repurposed from the LANL CTX experiment). The strong stray magnetic fields require turbopumps to be several metres from the HTS magnets: two turbos will be placed on long pumping ducts leading to the central cell, and two on the end cells. Including the cryopumps, the total centre cell pumping is 4800 l s$^{-1}$ for 200 l, and the total end cell pumping is 8400 l s$^{-1}$ each for 3900 and 4400 l. This does not include the additional pumping in the NBI source and titanium gettering in the beam dump, described later.

WHAM will be fuelled by piezoelectric gas valves placed close to the plasma inside the central cell, with an actuation time of under 1 ms. A maximum of four injectors will be available and each can be supplied with an arbitrary mixture of hydrogen, deuterium, helium and argon gas for fuelling and diagnostic purposes. Two injectors on pivoting mounts will supply gas fuelling at both ends of the central cell near the 4 T ECH resonance surfaces, one injector on a telescopic mount will be available at the mid-plane for diagnostics such as gas puff imaging and the remaining injector will be fixed on one of the end cells. The Wisconsin in situ Penning gauge will provide fast measurements of neutral gas pressure near the plasma edge (Kremeyer et al. 2020).

Successful plasma performance relies on the removal of neutral hydrogen atoms in the plasma edge to reduce energy losses from charge exchange events. In WHAM, a cold-sprayed tantalum coating on the inner surface of the central cell vacuum vessel will be tested as a neutral pump. The high melting temperature, excellent room-temperature ductility, outstanding corrosion resistance, low sputtering yield by hydrogen isotopes and exothermic hydrogen absorption are beneficial properties of tantalum, alongside the enhanced retention capacity of the engineered cold spray coatings. Our recent experimental work (Ialovega et al. 2023) has shown that tantalum coatings sprayed on SS 316L substrates maintain their mechanical integrity and physical properties upon high dose ($3\times 10^{25}$ D m$^{-2}$) deuterium irradiation and subsequent thermal annealing up to 1100 K. The absorbing tantalum first wall interface can therefore be regenerated by heating during wall conditioning to free up the trapping sites for hydrogen. This same technology may later be used for coating the plasma-facing components in the expansion tank rings.

2.3. The ECH System

The ECH system is built around a Gycom gyrotron and various ancillary equipment decommissioned from both the DIII-D and KSTAR tokamak experiments. The remainder of the ECH system includes a modern field programmable gate array (FPGA) control system, a General Atomics designed two-mirror aiming system that utilizes an electrically biased split waveguide to cross through the cyclotron resonance to the high-field side where it is launched and a full transmission line including two power monitor miter bends, two polarization shift miter bends, a 5 kV DC voltage break, pumpout tee, waveguide switch (to dummy load or machine), waveguide vacuum valve and all necessary straight lengths and couplers. Figure 4 contains a rendering of the in-vessel section, including the split waveguide.

Figure 3. Vacuum fields of WHAM at low field (a) and high field (b), not to scale. While the ECH resonant surface stays nearly constant axially, the fast ion turning point at $R_m \equiv B/B_0 = 2$ moves 20 cm and the radius decreases from 15 to 9 cm.

Figure 4. Transmission line leading from low stray field area of the gyrotron to WHAM device. (a) Shows the first section between the matching optics unit and the dummy load, (b) shows the long waveguide run from the gyrotron cage to the vacuum vessel and (c) shows a CAD rendering of the in-vessel transmission line.

The WHAM experiment plans to use fundamental high-field side X-mode ECH for two purposes. The first is to breakdown plasma and build up target density prior to NBI heating, similar to how the GDT device initiates discharges using ECH (Yakovlev et al. 2017). The second is to provide additional heating to the electrons, beyond the collisional heating from fast neutral beam ions, to transition from the gas dynamic regime to the CM regime and improve overall confinement. Without ECH, the initially cold electrons will act as a drag on the fast ions, preventing them from moving into the high ion temperature and good confinement regime of the CM. The very large ECH power density (1 MW/40 l) should be enough to allow this transition. This configuration is optimal as high-field side injection is strongly absorbed at low plasma density (see § 3.2) and is thus good for building up density. And with the higher frequency (110 GHz, compared with 54.5 GHz on the GDT), the waves can propagate without diffraction at four times higher density.

Accurately predicting the electron temperature achievable with our ECH system involves a significant multi-factor analysis, which must include effects of collisions with background neutrals, slowing down of fast ions on electrons, ion loss rate and bias ring voltage profile among other factors. Until such an analysis is complete, we will compare with the results of the GDT device to get a rough estimate of our target electron temperature. GDT was able to achieve a stable electron temperature of 600 eV using ECH with a power density of 0.4 MW/300 l. The ECH power density in WHAM is 1 MW/40 l, therefore we predict we will be able to achieve an electron temperature of at least 1 keV.

2.4. The NBI System

The WHAM experiment will use a 25 keV, 40 A neutral beam originally built for TAE Technologies by the Budker Institute of Nuclear Physics. The beam uses an inductively driven RF plasma source, with characteristic primary-, half- and third-energy components of 70 %, 20 % and 10 %, respectively. The neutralizer for the NBI is 80 cm long and will operate at a pressure of 5.2 mTorr during a shot. In order to prevent this neutral gas from entering the central cell, the NBI is equipped with 2 cryopumps, each capable of a hydrogen pumping speed of $7 \times 10^4$ l s$^{-1}$. NBI injection is through the origin at a 45$^{\circ }$ angle with respect to the machine axis (see figure 2). The $45^\circ$ injection angle sets the turning point for injected fast ions at the position where the magnetic field strength is twice the central field $B = 2B_0$ (see figure 3). The first phase of the experiment will begin with power supplies capable of 10 ms of injection.

One NBI system acts as an unbalanced $m=1$ perturbation on the plasma column. In general, when the shear-flow stabilization is operational then the $E\times B$ drifts will rapidly make symmetric the fast-ion content and this perturbation should be negligible. Although one NBI system will be sufficient for full plasma performance in WHAM (see § 3.3), the vessel has space for a second NBI system in case this assumption is not correct.

If injected beam particles slow without pitch angle scattering, they form a ‘sloshing ion’ distribution with several important features. Such a distribution has an axial density peak at the turning points and consequently forms a potential well that can trap a warm thermal ion population. The trapped thermal ions provide kinetic stability against the drift cyclotron loss cone instability (DCLC) by filling in the ambipolar hole in the ion distribution function (Kesner 1973; Simonen et al. 1983). In addition, angled injection decreases the pressure anisotropy compared with perpendicular injection, mitigating the Alfvén ion cyclotron (AIC) instability (Smith 1984). The sloshing ion population will be resolved with a proton detector array to measure the axial and radial intensity of 3.02 MeV D-D fusion protons. The degree of peaking is sensitive to the collisionality, with cold $T_e$ plasmas more strongly peaked as fast ions rapidly slow without diffusing in pitch angle. These sloshing ions are also critical to the HHFW RF heating described below.

2.5. The RF Heating System

A primary physics mission of WHAM is to accelerate high-energy beam-born sloshing ions by ion cyclotron absorption at the second to fourth harmonics of deuterium. With the 0.86 T central field, the 2$\varOmega _{ci}$ deuterium ion cyclotron frequency at the turning point is 26 MHz. A well-suited 1 MW Thales transmitter has been transported from the Archimedes device and is currently being commissioned pending substantial laboratory infrastructure upgrades. A complete RF transmission line utilizing components from both the LDX and C-Mod experiments at MIT is on site and ready for assembly. This includes two long stub tuners, bi-directional couplers with detectors, suitable pre-amplifiers and downstream electronics, several crates and cabinets full of miscellaneous elbows, straight lengths, tee adaptors and vacuum windows.

For the first phase of the experiment, RF heating will be done via an $m=0$ single strap antenna collocated with the sloshing ions’ turning point, where their density and residence times are maximized. As higher harmonic cyclotron damping is an FLR effect, fast ions will be preferentially heated, further improving their confinement (Stix 1992). The synergy between HHFW heating and neutral beam heating for enhancing the number of high-energy ions has been robustly observed and demonstrated on tokamaks (Petty et al. 2001; Rosenberg et al. 2004) and also simulated by Fokker–Plank calculations (Harvey, Petrov & Forest 2016).

A COMSOL Multiphysics finite-element model has been constructed using the RF module in the frequency domain to simulate antenna performance. To model wave propagation in the plasma, a fully anisotropic cold plasma dielectric tensor including artificial collisional damping was used (Lau et al. 2019). Both two-dimensional (2-D) axisymmetric and full 3-D simulations were run to inform the preliminary design of the first antenna. The 3-D RF electric field magnitudes with moderate, uniform collisional damping are shown in figure 5 together with a rendering of the antenna geometry. A single strap antenna fed by parallel plate transmission lines has been shown to couple to a cavity eigenmode in the limit of low single pass damping over a large range of core densities. When collisional damping is increased, the loading impedance was found to be strongly dependent on the edge density profile, but was of order 0.1–$1$ Ohm.

Figure 5. A three-dimensional COMSOL simulation of RF electric field strength (V m$^{-1}$) with 100 kW incident power at 28 MHz, with artificial collision rate at $0.01\omega _{\rm ci}$.

There is an ongoing effort to incorporate hot plasma effects such as higher harmonic ion cyclotron damping in self-consistent, full-wave simulations. We are exploring the possibility of coupling the KineticJ code to the cold plasma finite-element COMSOL model via kinetic correction currents (Green et al. 2018), while the AORSA code is also being adapted for the mirror geometry (Jaeger et al. 2001).

2.6. Biased limiter and end-ring system for controlling sheared flow

Magnetic mirrors have both radial and axial boundary conditions. In WHAM, the radial boundary condition is set by a tungsten limiter positioned in the central cell, and the axial boundary condition is set by an ensemble of concentric molybdenum end rings installed in both end cells. Tungsten and molybdenum were chosen for their high melting temperatures, resistance to sputtering by hydrogen isotopes and for manufacturing considerations. Figure 6 shows a photograph of the end rings, and also the central cell mounted limiter, which defines the potential on the outermost confined flux surface of the plasma. Both the limiter and end rings are mounted on linear axial translation stages to account for changes in the shape of the magnetic flux surface (see figure 2).

The limiter will be positioned to screen the HTS mirror coil cryostat from particle fluxes. Each end-ring ensemble consists of 19 rings, 10 of which can be independently biased at voltages up to 2.5 keV to create a radial electric field and hence a sheared rotation at the radius of choice. The biased rings are spaced with 9 isolated (floating) rings to prevent arcing. A set of programmable power supply modules has been modified to power the biased rings. The expected bias potential is comparable to the electron temperature, $T_e \sim 1$ keV, while the expected current draw of several hundred amps is set by ion saturation current and is far below the 2 kA power supply rating.

Figure 6. Panel (a) (photo) shows one end-ring ensemble mounted in the end cells for radial field control and panel (b) (rendering) a tungsten limiter to bias the outermost confined flux surface. The limiter is electrically insulated from the rest of the assembly and is linearly translatable.

A sub-goal of the WHAM experiment is to combine the expander bias end rings for vortex confinement with a direct energy converter (DEC). This will have substantial plasma–material interaction and engineering challenges. Single and multi-stage converters have long since demonstrated more than 50 % efficiency on tandem mirror devices (Barr, Moir & Hamilton 1982). The ambipolar accelerated distribution makes the use of a DEC promising, as even a single stage DEC should be able to recapture $5 T_e/ E_0$ of the injected neutral beam energy. figure 7 shows the exhaust ion energy distribution for the same dataset used later in § 3.7.

Figure 7. An example FBIS histogram of the on-axis lost ion distribution at the segmented end rings, after acceleration by the $5 T_e/e$ potential.

3. Performance Modelling of WHAM

Classical confinement in the simple mirror is governed by pitch angle scattering of ions into a loss cone modified by the ambipolar field that limits the loss of electrons (Bing & Roberts 1961; Fowler & Rankin 1962, 1966). The collisional relaxation of the ion distribution function is governed by two time scales

(3.1)\begin{gather} n_{20} \tau_s = 5 \frac{T_{e,{\rm keV}}^{3/2} \mu }{ Z^2}\ \mbox{ms}, \end{gather}

(3.2)\begin{gather}n_{20} \tau_{ii} = \frac{\langle E_{i,{\rm keV}}\rangle^{3/2} \mu}{8 Z^2}\ \mbox{ms}, \end{gather}

where $\tau _s$ represents slowing and energy loss to electrons, $\tau _{ii}$ is the ion–ion pitch angle scattering time, $\mu$ is the ion mass in units of proton masses and $Z$ is the ion charge. It is easy to see that these time scales are similar when $T_e=\langle E_i \rangle / 40^{2/3}\sim \langle E_i \rangle / 10$.

Consider the idealized experiment that fuels plasma with a high-energy neutral beam: when the electron temperature is low, $T_e \leq \langle E_i \rangle / 10$, the ions slow down before scattering into the loss cone, and the confinement, which is still determined by the ion temperature and pitch angle scattering rate, is essentially set by the electron slowing time. At higher electron temperatures, the situation reverses and the confinement becomes independent of the electrons, dependent only on the ion energy. This confinement gets very long at high ion energies, as $\tau _{ii} \propto E_i^{3/2}$. The fusion gain $Q$, however, does not continue improving at higher temperatures due to increased radiative losses, as can be seen in figure 1.

The much lighter electrons pitch angle scatter more rapidly, on a time scale associated with electron–electron collisions

(3.3)\begin{equation} n_{20} \tau_{\rm ee} = 5.8 T_{e,{\rm keV}}^{3/2} \ \mathrm{\mu} {\rm s}. \end{equation}

These electrons would be lost very quickly (compared with the ions) were it not for the development of a self-consistent ambipolar potential that retains the electrons and increases ion losses. The result is that electrons essentially thermalize and only the very high-energy electrons with energies above $5 T_e$ leave the system. The role of this ambipolar potential is profound: it ensures the electron confinement and also leads to plasma flows (both an outflow but also, as will be shown, rotation). However, this potential also creates an ambipolar ‘hole’ in the ion distribution function that can lead to kinetic instabilities.

For the simple case of NBI-only plasmas, a number of authors (see Killeen & Marx 1970; Killeen, Mirin & Rensink 1976; Egedal et al. 2022) have considered how the classical confinement scales with beam injection energy by solving the Fokker–Planck equation for the self-consistent ambipolar potential, electron temperature and ion distribution function. The results are all generally in agreement with

(3.4)\begin{equation} n_{20} \tau_p =250 E_{b,100\ \text{keV}}^{3/2} \log_{10} R_m\ \mbox{ms}, \end{equation}

although the exact coefficient depends on details, such as the shape the magnetic field, the injection angle and the degree to which electrons are heated by alpha particles in a deuterium-tritium (DT) device. Contrast this with the gas dynamic (GD) or collisional confinement time of

(3.5)\begin{equation} \tau_{\rm GDT} = \frac{ R_m L_p}{c_s} = 5.2 R_m L_p T_{e,{\rm keV}}^{{-}1/2} \ \mathrm{\mu}{\rm s}, \end{equation}

where $R_m$ is the mirror ratio and $L_p$ is again the plasma half-length (Ivanov & Prikhodko 2013). Mirrors might be expected to transition from the GD to CM regime as the ions become less collisional. For this reason (among others), care must be taken to self-consistently model the exact magnetic coil geometry, the ion distribution function and also the resulting effect of the plasma currents on magnetic equilibrium. Equation (3.4) is well supported by the historical record of fast ion confinement shown in figure 1 in minimum-B mirrors (Coensgen et al. 1975; Post 1987; Wurzel & Hsu 2022). Historically, the data collection ended when the US mirror program was cancelled in 1986. Perhaps most relevant to this paper are the 2XIIB results at mirror ratio 2 that achieved the first 10 keV ion plasma (Coensgen et al. 1975). The GDT is not included here as its pulse length is shorter than the expected fast-ion confinement time, and because it is a two component plasma with both warm and hot ions.

In the modelling that follows, our starting point is provided by the FBIS code (Egedal et al. 2022), with the FBIS quasi-analytic solution for the fast ion distribution function and electron temperature, the self-consistent ambipolar potential and density along a single flux tube. In the FBIS solver, electron drag on the fast-ion population heats the electrons, and the electrons do work sustaining the ambipolar potential which accelerates the exhaust ions. An additional input heating term, $p_{\rm aux.}$, is added to describe ECH heating, radiative cooling or some other source of power to the electron species. To extend these results, we combine a set of 1-D solutions from the FBIS code of $T_e, p_{\|}(B/B_0),\ p_\perp (B/B_0),$ $n(B/B_0)$ and $\phi (B/B_0)$ to construct the 2-D equilibrium for a given magnetic geometry, where again $B_0$ is the midplane field. The solution includes both magnetic and electric equilibria including the plasma flows.

As pressure and therefore $\beta = 2 \mu _0 p/ B_0^2$ increases, plasma diamagnetism excavates the magnetic field, increasing the mirror ratio and improving confinement

(3.6)\begin{equation} R_m = R_m^{\rm vac}(1 - \beta)^{-{1}/{2}}.\end{equation}

As has been noted by a number of authors, this improvement in confinement at high pressures could lead to very high $Q$ values, even in the simple mirror (Beklemishev 2016). Hence, getting the equilibrium right (§ 3.4) is critical for assessing performance. However, very high $\beta$ values such as those observed in 2XIIB, would seem to be inconsistent with mirror confinement due to adiabaticity issues, and so equilibria must be evaluated against this constraint as well. Furthermore, (3.4) and (3.6) make clear that in going from $\beta =0$ to $\beta =0.9$, the gain in $n\tau$ is only 50 %.

With this background, we now address individual subtopics before integrating them in § 3.7.

3.1. MHD Stability

The principal challenge in developing the axisymmetric mirror is, and has always been, stabilizing the flute or interchange mode that is universally unstable in non-cusped axisymmetric systems (Rosenbluth & Longmire 1957). The strategy for WHAM will be to use the known physics of FLR and vortex stabilization first, and later to test other promising techniques that may extrapolate better to future fusion reactors. Many ideas have been suggested and demonstrated to work over the past 50 years, (summarized well by Ryutov et al. 2011), including the use of good curvature (Mirnov & Ryutov 1979; Bagryansky et al. 2011), cusp boundaries (Prater 1971; Ferron et al. 1983; Kotelnikov et al. 2019), ponderomotive effects (Ferron et al. 1987), rotating magnetic fields (Seemann, Be'ery & Fisher 2018; Zhu et al. 2019), close fitting conducting shells (Kesner 1985; Beklemishev 2016; Kotelnikov et al. 2022) and feedback (Prater 1971; Be'ery, Seemann & Fisher 2014; Be'ery & Seemann 2015; Beklemishev 2017).

FLR stabilization of $m>1$: hot temperatures and energetic ions in mirrors are beneficial to MHD stability because FLR effects stabilize interchange for azimuthal mode numbers $m>1$. For this to occur, a magnetic mirror must be long enough to satisfy (Ryutov et al. 2011)

(3.7)\begin{equation} m > 2 \frac{a^2}{L_p \rho_i} = 2 \frac{a/\rho_i}{L_p/a}, \end{equation}

which sets a critical plasma half-length $L_p$, given ion gyroradius $\rho _i$ and plasma radius $a$. Ryutov's scaling is derived from a dispersion relation of the form $\omega ^2 + m\omega (v_{\rm Ti} \rho _i/a^2) + \varGamma _0^2$, where $\varGamma _0 \sim v_{\rm Ti} / L_p$. For WHAM-0.86 with $a/\rho _i=4$ and $L_p/a=10$, all modes with $m\geq 2$ should be FLR stabilized.

Vortex flow stabilization of $m=1$: the most successful technique for countering $m=1$ interchange thus far has been to use sheared flows (Cho et al. 2005, 2006) and vortex stabilization (Sakai, Yasaka & Itatani 1993; Beklemishev et al. 2010; Bagryansky et al. 2011; Yakovlev et al. 2018). Vortex stabilization relies on an edge region rotating around a stationary central plasma, while sheared flow refers to strong azimuthal flow shear not necessarily localized to the edge. In both cases, growing instabilities are saturated at small amplitude by a strong shearing rate. We focus on vortex stabilization here and sheared-flow stabilization will be mentioned in § 3.2.

The ambipolar potential that maintains quasineutrality between the electron and ion species in mirrors has long been recognized (Post 1961; Kelley 1967) theoretically described (Fowler 1981) and experimentally measured (Post 1987). When paired with a radial temperature gradient, this ambipolar potential naturally leads to positive $E_r$ and an azimuthal $E\times B$ flow in the same direction as the ion diamagnetic drift (see § 3.4). This is, in general, destabilizing for interchange (Severn et al. 1983). With external biasing, both the HIEI and GDT experiments applied a positive bias to the limiter and edge, and a negative or grounded potential to the low radius expander rings, producing a negative $E_r$ in the edge and reversing the direction of the $E\times B$ drift (Sakai et al. 1993; Bagryansky et al. 2018). Applying the opposite polarity to reinforce the natural $E_r$ was experimentally challenging.

According to Beklemishev et al. (2010), this experimental technique works if, firstly, the plasma has sufficient line-tying to the end plates to neutralize the polarization electric field of the flute mode and ensure closed flow lines, and secondly, if the width of the flow layer is smaller than the plasma radius. These two conditions ((20) and (22) in their text) are

(3.8)\begin{gather} \frac{e \phi}{T_e} \geq 4 R_m^2 \left(\frac{\rho_*^3 L_p^2}{a^2 L_\kappa^3} \right) \left(\frac{p_i}{p_e} + 1\right)^{{3}/{2}}, \end{gather}

(3.9)\begin{gather}R_m \leq 0.7 \left(\frac{a}{\rho_*}\right)^2 \frac{L_\kappa}{L_p} \left(\frac{p_i}{p_e} +1 \right)^{-{1}/{2}}, \end{gather}

where $\phi$ is the applied bias, $\rho _*\equiv c_s / \varOmega _i = \sqrt {m_i T_e} / q_e B$ is a length scale given by the sound speed and ion cyclotron frequency, $L_p$ the plasma half-length and $L_\kappa$ is the field curvature scale length. While having a large mirror ratio $R_m$ is in general good for confinement, it makes satisfying both of these conditions more difficult. An independent assessment by the Budker team has estimated that vortex stabilization will likely be sufficient for WHAM at high fields where $T_e,B \sim 1\ \text {keV}, 0.86\ \text {T}$, while it will likely not be for WHAM at low fields where $T_e,B \sim 0.5\ \text {keV}, 0.32\ {\rm T}$. However, their derivation uses gas dynamic confinement to model the outflowing ion current and may not apply in the CM regime.

3.2. High-Density Target Formation through ECH Breakdown and Electron Temperature Control

To ionize the plasma and achieve the high electron temperatures necessary for CM performance, the WHAM ECH system will use fundamental X-mode heating launched from the high-field side. Due to its high absorption at low densities, both the GDT (Yakovlev et al. 2017) and GAMMA-10 (Shidara et al. 2011) experiments used similar high-field launched X-mode ECH. Each observed that breakdown leads to an initial low-density, high-energy electron population, which then causes rapid and near complete ionization of neutral background gas on 1 ms time scales. Rapid ionization is important as WHAM currently does not have power resources for longer than a 20 ms shot duration.

The GENRAY ray-tracing code from CompX (Smirnov, Harvey & Kupfer 1994) predicts similar behaviour for ECH in WHAM. After an initial buildup of a low-density plasma through neutral ionization (the physics of which is not properly captured in ray tracing), ECH is calculated to be 20 % absorbed even at $2.5\times 10^{16}\ {\rm m}^{-3}$ and over 99 % absorbed above $3\times 10^{17}\ {\rm m}^{-3}$, as shown in figure 8. The electron temperature is also shown to have a very minor effect on absorption. Thus, there should be no issues with the ECH coupling to the low-temperature electrons produced though ionization of neutrals.

Figure 8. Absorbed ECH power vs plasma density in WHAM-0.32, using a centrally launched ECH configuration at different electron temperatures, calculated using the GENRAY ray-tracing code with ion temperature fixed at 10 eV.

As ECH breakdown leaves relatively few background neutral particles compared with a plasma gun discharge, it is less of a drain on the fast ion population through charge exchange. However, with plasma radius of 0.1 m, WHAM will be a modest experiment without a clear distinction between the edge and core plasma, and being so small it will have substantial NBI shine-through. The question of whether a steady-state beam-only operation can be achieved on WHAM is discussed in § 3.3.

After forming a seed plasma, the second role of the ECH system is to control the radial electron temperature profile for both the transition to the weakly collisional CM regime and for $E_r$ control. If a gradient in electron temperature is established though ECH, a radial electric potential gradient will form and interact with the potential defined by the limiter and bias rings. Both GAMMA-10 and the GDT have successfully used localized ECH to create sheared flows and internal transport barriers (Cho et al. 2006; Yakovlev et al. 2018). The WHAM experiment will control the beam direction with a rotatable mirror, which allows for the ECH beam to be aimed at any radial position, see figure 9. By targeting the beam at the plasma edge, the electron temperature will increase with radius and give the desired ambipolar potential profile, assuming that the absorbed power is successfully transferred to thermal electrons. While the beam diameter in waveguide is roughly 2 cm, it will spread after exiting the launcher and the extent of control over the deposition profile is not yet known. The effect of off-axis ECH injection on the radial density profile is also yet to be seen, but data from GDT suggest that, despite the localized deposition profile at the ECR surface, the resulting seed plasma will have broad and uniform radial profile. This means that, in most cases, the launching optics can be tuned for optimal ECH at later stages of the discharge, without much regard for the initial breakdown stage.

Figure 9. (a) Location and absorption percentage of ECH power vs mirror angle in WHAM-0.32. (b) The proposed in-vessel ECH beam path including a fixed mirror mounted on the end of the waveguide run and a second, rotatable mirror mounted on the vessel.

3.3. NBI fuelling at CX dominated energies

An ideal neutral beam heated mirror experiment would begin with a target plasma dense enough to fully attenuate the neutral beam, and the neutral beam alone would fuel and heat the plasma. Unfortunately, WHAM, GDT, TAE and all previous mirror experiments operate in a regime in which charge exchange (CX) collisions (which do not fuel the plasma) dominate over impact ionization collisions (which do). The challenge is that at 25 keV for deuterium, the CX cross-section $\sigma _{\rm cx}$ is roughly four times the combined cross-sections for electron and ion impact ionization, $\sigma _{e+i}$. We represent this with the variable $\chi \equiv \sigma _{\rm cx} / (\sigma _{\rm cx} + \sigma _{e+i})$, shown in figure 10. The consequence is that a 25 keV neutral deuterium atom is more likely to CX, contributing nothing to the overall density, but increasing the temperature while pumping out warm plasma.

Figure 10. Deuterium cross-sections for charge exchange and cumulative impact ionization, including the resultant $1-\chi$. Note that at beam energies of 25, 70 and 120 keV, the fraction of impact ionization collisions $1-\chi$ is roughly 20, 50 and 80 %.

To better understand these CX and fuelling challenges, NBI fuelling on WHAM has been studied using FIDASIM. Specifically, we use the python-based pyFIDASIM routines written by Geiger et al. (2020) to model neutral beam deposition. FIDASIM, like TRANSP, calculates probabilities of all the atomic processes and uses a Monte Carlo method to track the fuelling rates of both the incident neutral beam and the reionization of halo neutrals generated by the CX process. In concert with the FBIS solution particle confinement time, we use FIDASIM in an iterative scheme to calculate the evolution of the density profile achieved by NBI fuelling. This assumes some pre-existing ECH produced target plasma, but does not include additional fuelling that would come from cold neutrals being ionized at the plasma edge.

This iterative procedure begins with a simple differential equation

(3.10)\begin{equation} \frac{{\rm d} n}{{\rm d} t} = S_b - \frac{n}{\tau_p}, \end{equation}

where the density $n$ is the average density along a field line, or equivalently the total number of particles divided by the flux surface volume $N/V(\psi )$. pyFIDASIM provides the fuelling source rate $S_b \equiv (1/\Delta V) {\rm d}N/{\rm d}t$ for each flux surface. The density is then iterated as

(3.11)\begin{equation} n^{t+1} = n^t + \Delta t \left[ S^t - \frac{n^t}{\tau_p^t}\right].\end{equation}

The FBIS solution gives the steady-state particle confinement time $\tau _p^*$, as well as the analytic result that $n \tau _p$ is constant, so that $\tau ^t_p \equiv \tau _p^* n^{*}/n^{t}$ and the above equation can be solved. The results from this technique are plotted in figure 11, which predicts after 10 ms a density of around $3\times 10^{19}/m^3$ with 65 % shine-through. Thus, despite the short path length and unfavourable $\chi$, WHAM should be able to achieve its target density solely via neutral beam fuelling on a tenuous background plasma. It will not reach a steady state, however, as the particle confinement time is much longer than the experiment duration.

Figure 11. Density profile evolution including CX losses predicted from pyFIDASIM and FBIS in an iterative method, for deuterium NBI with primary/secondary fractions of 0.6/0.4.

Role of CX on fast-ion distribution: since CX processes will dominate over impact ionization collisions, the ion distribution function is likely to be strongly modified due to ‘pump out’ across the entire distribution. To reiterate, while CX still increases the total plasma energy, it does not affect the fuelling. One conclusion from figure 10 is that, for high efficiency fuelling, a beam energy above 100 keV is required.

Thus far, the FBIS Fokker–Planck solver presumes only losses through the trapped– passing boundary. In principle, the CX process can be represented by a source term with the energy and pitch angle of the NBI and an equal loss term that is distributed over the entire distribution. Accurately including the effects of CX into FBIS requires weighting the losses by the relative time $\tilde {t}_{\rm nbi}$ that particles of a given pitch angle spend in the beam path, with deeply trapped ions being most likely to be lost.

For WHAM, the NBI beam path is defined by a region which extends to $\pm z_{\rm nbi}$ around the midplane. Particles with turning points inside of $\pm z_{\rm nbi}$ spend all of their time in the beam path. For all other particles, the fraction of time spent in the path $\tilde {t}_{\rm nbi}$ can be calculated as

(3.12)\begin{equation} \tilde{t}_{\rm nbi}(L_{\rm bounce}) = \frac{ \displaystyle \int_0^{z_{\rm nbi}} \, {\rm d} z / v_{\|}(z) }{\displaystyle \int_0^{L_{\rm bounce}} \, {\rm d} z / v_{\|}(z) }. \end{equation}

The result for WHAM 0.86 is plotted in figure 12.

Figure 12. Relative time $\tilde {t}$ spent in the NBI beam path for particles as a function of $\eta$, the bounce-averaged pitch angle variable. Beyond the trapped–passing boundary (red dashed line), $\tilde {t}$ is treated as 0.

Properly incorporating this weighting into the analytic FBIS solver for CX losses requires a modification to the pitch angle eigenfunctions (and loss rate eigenvalues), which is beyond the scope of this work. As an approximation, we assume the eigenfunctions are unchanged and allow FBIS to solve for the ion distribution function normally. We then incrementally remove density proportional to $\Delta t_{\rm nbi}$ weightings until the losses through the trapped–passing boundary equal $(1-\chi )$ of the injected current. In other words, we treat all ionized particles as part of the injected current, but remove CX losses until the current into the loss cone equals the impact ionization current.

There are several problems with this method. The first is that, if the CX replacement time constant $\tau _{\rm CX}^{\rm rep.} \equiv N_{\rm tot.} / I_{\rm nbi}^{\rm CX}$, or the time to replace all the confined particles $N_{\rm tot.}$ via CX fuelling $I_{\rm nbi}^{\rm CX}$ is substantially less than the particle confinement time, then the distribution will be replaced by beam ions faster than it can slow or scatter. This ratio $\tau _{p}/\tau _{\rm CX}$ is simply $\chi / (1-\chi )$, which for $\chi =0.8$ is 4, so unfortunately larger than one. The second problem is that FBIS will calculate the scattering as a result of both the CX and impact ionized populations, when it should only include the impact ionized part. Ignoring both concerns, as a simple check we can use the $n \propto \sqrt {I}$ dependence (Egedal et al. 2022) and predict that the total particle inventory should be reduced roughly as $\sqrt {1-\chi } = 0.45$ of its initial value. This is very close to the 0.47 value observed in figure 13, and this approximation improves at smaller values of $\chi$. The initial and final distribution functions using this method are plotted in figure 13.

Figure 13. (a) The unmodified ion distribution function in $(\eta, v)$ space, where $\eta$ is a bounce-averaged pitch angle variable (see Egedal et al. 2022). (b) The same distribution after including pitch angle weighted CX losses for $\chi =0.80$. (c) The axial density profiles for the two cases, including the ratio of total particle inventory. (d) The final velocity space distribution function, showing the reduction in particles with large $v_\perp /v_{\|}$. The red dots represent the NBI injection location, and the red dashed lines the trapped–passing boundaries.

Aside from the lower fuelling efficiency, two features illustrate the importance of including this effect when evaluating the overall performance of WHAM. First, deeply trapped particles are most effectively pumped, resulting in a ‘loss wedge’ feature at ${v_\parallel ^2}/{v^2} < 1 - {B_0}/{B(z_{\rm nbi})}$. Second, the distribution function is considerably more sloshing, with a stronger density enhancement at the turning point. While we expect these results will be an important feature of the experiment, we emphasize the approximate nature of these estimates. Furthermore, the FBIS solver is still not configured to model the half and third energy components of a real neutral beam.

3.4. Anisotropic MHD Equilibrium and Ambipolar Potential

The modelling of how the vacuum magnetic field is modified for high-$\beta$, anisotropic plasmas with flow is done with an iterative Grad–Shafranov solver. First, the vacuum flux function $\psi _{\rm vac}=\psi (r,z)$ and field ${\boldsymbol B}={\boldsymbol B}(r,z)$ is calculated from the coil and current parameters using an internally developed, open-source Green's function solver named Pleiades.¹ Next, the ion distribution function solution is calculated from the FBIS code, giving a parallel and perpendicular pressure as a function of the normalized magnetic field strength $\tilde {B} = B/B|_{z=0}$. Extending the 1-D FBIS results to the 2-D domain requires a radial pressure profile $p(\tilde {\psi })$, where $\tilde {\psi } \equiv \psi /\psi _{\lim }$ is the flux function normalized to the flux value at the limiter, and $p(\tilde {\psi }>1) \equiv 0$. This radial profile can either be chosen a priori or a realistic $p_{\rm aux}$ profile representing increased cooling of the electrons at the edge can be specified in FBIS .

The combined 2-D ion pressure profiles of $p_\perp (\tilde {B},\tilde {\psi })$ and $p_\parallel (\tilde {B},\tilde {\psi })$ are then interpolated onto the $RZ$-grid and used to calculate the plasma currents from the diamagnetic (gradient and curvature) drifts as derived in Appendix A

(3.13)\begin{equation} {\boldsymbol J}_{{\perp}} ={-}\frac{{\boldsymbol b} \times \boldsymbol{\nabla} p_\perp}{B} - (p_{\|}-p_\perp) \frac{{\boldsymbol b} \times {\mathbb \kappa } }{B} - n m_i \omega^2 r^2 \frac{{\boldsymbol b} \times \hat{\boldsymbol r} }{B}, \end{equation}

where ${\boldsymbol b} = {\boldsymbol B}/|B|$ is the magnetic field unit vector, $p$ is the pressure, $B$ the magnetic field strength and ${\boldsymbol \kappa } = {\boldsymbol b} \boldsymbol{\cdot} {\boldsymbol{\nabla} }{\boldsymbol b}$ the magnetic curvature. The last term, which is associated with plasma rotation, is included here for completeness but is in general small for subsonic rotation. Although not yet implemented in this code, biased limiters and end rings can be used to modify both the rotation and centrifugal confinement from the boundary

(3.14)\begin{equation} \frac{e \phi(\ell,\psi)- e V_{\rm bias}(\psi) }{T_e} = \ln{\frac{n(\ell)}{n(0)}}, \end{equation}

where $V_{\rm bias}(\psi )$ is an applied potential.

The flux $\psi _{\rm plas}$ from the current in (3.13) is computed using the Pleiades Green's function method and is added to the vacuum magnetic flux solution $\psi _0$. From the sum $\psi _0 + \psi _{\rm plas}$, the magnetic field ${\boldsymbol B}$ is recalculated. With the updated field geometry, the pressure profiles are again mapped to the $RZ$-grid and the process is repeated. A solution is reached after the difference between iterations, summed over the entire domain, falls below specified small threshold. A result for WHAM 0.86 at $\beta =0.2$ with an imposed $\cos (\tilde {\psi })^2$ pressure profile is shown in figure 14. It is worth reiterating here that this axisymmetric solution remains MHD unstable to perturbations as, unlike the minimum-B configurations of earlier experiments, axisymmetric mirrors are always MHD unstable, though they may be stabilized by FLR or sheared-flow effects.

Figure 14. Anisotropic equilibrium solution for WHAM at $\beta =0.2$ corresponding to a density of $n=0.3\times 10^{20}$ m$^{-3}$ and an average ion energy of 10 keV at the midplane. The parallel and perpendicular pressure profiles at left and the azimuthal plasma current at right. Colour contours are identical to figure 3.

Parallel force balance is controlled by

(3.15)\begin{equation} {\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} p_\perp+ {\boldsymbol B}\boldsymbol{\cdot} \boldsymbol{\nabla} \left(\frac{ p_\parallel{-} p_\perp}{B}\right) - n m_i {\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} \left( \frac{\omega^2 r^2}{2}\right) = 0, \end{equation}

where we recognize the second term as the mirror force and the third term as a centrifugal confining force. The FBIS solution provides parallel force balance self-consistently as a solution of the drift-kinetic solution that results in local values of $p_\perp$ and $p_\parallel$. The resulting density variation along the magnetic field is then used to compute the electric potential, since the electrons are electrostatically confined for many $\tau _{\rm ee}$ with a nearly Boltzmann distribution. The centrifugal confining potential, shown on the right of figure 15, is predicted to be small for NBI-only plasmas with shallow radial $T_e$ gradients.

Figure 15. At left, the azimuthal $E\times B$ and diamagnetic flow speeds out to the limiter, showing nearly solid body rotation ($v_\phi \propto R$). Centre, the potential distribution after the FBIS+Pleiades+pyFIDASIM iterative procedure described in § 3.7. At right, the axial ambipolar, centrifugal and total potentials along the dashed flux surface.

3.5. Non-Classical Confinement

After § 3.1, we have been ignoring other non-ideal effects. In reality, once the problem of MHD interchange has been solved, other limits may degrade the mirror performance. The following are softer critical thresholds that generally increase velocity space diffusion, but are not disruptive like interchange. The two principal kinetic instabilities are the DCLC and the AIC (Baldwin 1977; Ferron & Wong 1984). A rigorous analysis of these instabilities would evaluate stability from Particle-in-cell (PIC) or gyrokinetic simulations. While substantial progress has been made addressing the computational challenges of high-field mirrors (see Francisquez et al. 2023; Petrov et al. 2023), in the qualitative analysis that follows, the fast-ion distributions’ outputs from the FBIS solver are used to assess the severity of each kinetic instability. Drift waves and convective cells can ruin the assumption that parallel transport losses dominate over radial losses, and these are addressed below. Fast-ion adiabaticity matters at high energies or low magnetic field strengths and will be addressed in § 3.6.

DCLC stabilization: the DCLC instability is the result of an ion-Bernstein-like mode driven by the free energy associated with the non-thermal ambipolar hole coupling to an electron drift wave associated with the density gradient. Following the work done by Kotelnikov, Chernoshtanov & Prikhodko (2017), DCLC is stabilized if the radial density gradient scale length $L_r \equiv n/\boldsymbol{\nabla} n$ is much larger than the ion Larmor radius $\rho _i$. When $L_r$ is approximately the plasma radius $a$, then $N_\rho \equiv a/\rho _i$ must exceed approximately 100 for stability (Ferron & Wong 1984), while in WHAM $N_\rho$ will be 3 or 4. An astute reader may observe the apparent contradiction between FLR stabilization (which benefits from small $N_\rho$), and DCLC (which is stabilized by large $N_\rho$). As a preventative strategy, one could produce a pressure profile with a flat core and a steep edge: DCLC will be stable in the core where the pressure gradient drive is zero, and in the edge the strong vortex flow shear will saturate DCLC at small amplitudes. As noted in (Ferron & Wong 1984), this localization of the DCLC will limit its effects to regions of strong gradients and also saturate at manageable levels.

A different strategy involves the addition of low energy ions into the ambipolar hole, which inhibits DCLC by reducing the positive gradient in the distribution function. Kotelnikov provides a quantitative threshold

(3.16)\begin{equation} \alpha\equiv \frac{n_c}{n_h} \gtrapprox 4.22 \left( \frac{\mu}{N_\rho} \right)^{3/2}, \end{equation}

where $\mu$ is the ion atomic mass number, and $\alpha$ represents the fraction of cold to hot ions. In past experiments, $n_c$ has been introduced by cold plasma streams or by trapping in the Yushmanov potential associated with sloshing ions (Yushmanov 1966; Kesner 1973). The GDT, for example, is likely stabilized by the cold target plasma at first and then later by the sloshing ions. When evaluated at values expected for WHAM, this constraint unfortunately demands a cold ion population comparable with the fast-ion population.

A new mechanism for maintaining DCLC stability uses the centrifugal potential associated with plasma rotation to trap low-energy ions. As described earlier, the radial electron temperature profile and the axial ambipolar potential together produce an $E \times B$ flow and consequently a centrifugal confining potential $\phi _{\rm cent}$ shown in figure 15. To enhance this potential, we assume a strong but realistic radial electron temperature gradient of $\Delta T_e / \rho _i \sim 300$ eV cm$^{-1}$ has been made in the edge with off-axis ECH. This potential has been included self-consistently into the FBIS solution of the distribution function, and the result is shown in figure 16. With this substantial radial potential gradient, the combination of the trapping potential between fast-ion turning points and the centrifugal potential increases the population of ‘cold’ ions by 4 %. This may not be sufficient in the small $N_\rho$ of WHAM but could likely play a role in a larger next step device.

Figure 16. The change in the midplane ion distribution functions in the region of steep edge gradients at $\varPsi = 0.5 \varPsi _{{\rm lim.}}$ without (left) and with (right) the centrifugal effects. Note the difference in the dashed pink line representing the loss cone at the midplane.

AIC stabilization: we can similarly make predictions of the degree of AIC instability from FBIS distribution functions. The AIC is driven by the ion temperature anisotropy and depends on the plasma $\beta$ in a way similar to mirror and firehose instabilities. An absolute instability threshold was identified by Smith (1984)

(3.17)\begin{equation} \beta_{{\perp}}\left(\frac{T_{{\perp}}}{T_{\|}} \right)^2 > 3.5. \end{equation}

However, the threshold for convective instability can be much lower: GAMMA-10 observed AIC when $\beta _{\perp }({T_{\perp }}{T_{\|}})^2$ exceeded only 0.3 (Katsumata et al. 1996). For comparison, in the high $\beta =0.2$ case of figure 14 this parameter approaches 0.4 at the turning points. At present, we do not have a good estimate for the precise value of AIC onset in WHAM, although this comparison indicates it will not be a problem until WHAM exceeds the performance prediction of § 3.7.

We conclude our discussion of kinetic instabilities by noting that the GDT device has seen no evidence of either DCLC and minimal evidence of AIC (Ivanov & Prikhodko 2013), and that our discussion here may be too pessimistic. Its likely that sloshing ion distributions may have negligible kinetic instabilities (Post & Ryutov 1995). Regardless, without a fully kinetic simulation test, understanding the extent of kinetic instabilities will be left to the experiment.

Cross-field transport: drift waves, trapped particle modes, convective cells: the basic assumption of the CM scaling is that parallel transport dominates over perpendicular transport and that energy is lost via particle losses, i.e. the Lawson parameter $n\tau$ is set by the particle confinement time $\tau _p$. In addition to the kinetic instabilities described above, the presence of traditional drift waves could mix plasma radially via turbulence, and thereby allow energy to diffuse outward, just as in toroidal systems (including trapped electron modes, electron temperature gradient modes and ion temperature gradient modes). The lack of toroidal curvature removes an important driver of these modes from mirrors (Simonen et al. 2008). Much less work has been carried out to evaluate which instabilities (if any) will grow and be most problematic, but numerous authors have speculated that trapped electron modes may be present in the long central cells of tandem mirrors (Horton et al. 2010; Fowler, Moir & Simonen 2017). Radial transport associated with neoclassical transport was inferred in TMX-U, trapped particle modes were observed in Tara and convective cells measured directly in GAMMA-10. In a remarkable set of follow-on papers, GAMMA-10 went on to document how sheared flow created by the natural ambipolar potential fully stabilized these convective modes and restored the tandem to having purely axially losses (Cho et al. 2005, 2006).

Beneficially, fast-ion FLR effects may mitigate turbulence as well as interchange modes by averaging over small fluctuations ($k_\perp \rho \gg 1$). An absence of turbulence was observed by TAE that was attributed to fast ions(Schmitz et al. 2016). In tokamaks as well, experimental measurements on the Tokamak Fusion Test Reactor (TFTR) show that $\alpha$-particles exhibit classical diffusion due to the FLR effect (Zweben et al. 1997). Mirrors tend to optimize with $T_i/T_e \gg 1$ and this is known to stabilize ion temperature gradient (ITG) driven turbulence in toroidal devices, e.g. TFTR supershots. On the other hand, $T_i/T_e > 1$ is known to destabilize electron temperature gradient (ETG) modes. Sheared flow associated with differential rotation long wavelength streamer-like ion-scale turbulence such as ITG and TEMs. Electron-scale ETG may not benefit from this effect depending on the radial correlation lengths. Further evaluation of these instabilities is warranted.

3.6. Fast-Ion Adiabaticity

Non-adiabatic effects are always present as finite jumps $\Delta \mu$ in the magnetic moment as trapped particles pass through the field minima (Hastie, Hobbs & Taylor 1969). As discussed by Cohen, Rowlands & Foote (1978b), fast-ion adiabaticity is often a more restrictive condition on $\beta$ than the onset of the mirror mode. Conservation of the first adiabatic invariant $\mu = m v^2_\perp / 2B$ fails if a particle experiences a large change in magnetic field strength in one gyro-period. This is represented by the adiabaticity parameter $\alpha \equiv L_B / \rho _{\|}$, where $\rho _{\|} \equiv v_{\|}/\omega _{\rm ci}$ and $L_B \equiv |B|/|\boldsymbol{\nabla} B|$, where large $\alpha$ means the magnetic field gradient scale length is larger than the parallel distance travelled in one gyroperiod. The size of the jumps scales as $\Delta \mu \propto \, {\rm e}^{-\alpha }$ (Cohen et al. 1978b). In general, particles are well confined only if $\alpha \gg 1$, and we find non-adiabatic effects become significant when $\alpha \leq 10$.

In WHAM, this adiabaticity condition is not satisfied everywhere. Figure 17 shows $\alpha$ for the two magnetic configurations. In WHAM at low field (0.3 T), NBI deuterons are born with $\alpha \simeq 10$ and experience appreciable changes in $\mu$ each bounce. In contrast, for the 0.8 T high-field case, $\alpha$ is much larger than 10 all the way to the loss cone. In short, non-adiabatic effects may be present in WHAM at low field but are not expected at high field.

Figure 17. Numerical calculation of the adiabaticity parameter $\alpha \equiv L_B/\rho _{\|}$ in velocity space for the two WHAM configurations. Speeds are normalized to the 25 keV deuterium injection energy, with the red dot locating the expected injected velocity.

To more carefully explore when adiabaticity fails in our specific geometries, we implemented a standard Buneman–Boris algorithm for computing single particle trajectories. As this algorithm is known to not conserve energy, simulations are run with enough temporal resolution to keep the cumulative error in energy below 0.1 %. To model fast-ion slowing behaviour without pitch angle scattering (the cold $T_e$ case), an electron drag term is added so that at each timestep a scalar factor $\exp (-{\rm d}t/\tau _s)$ is multiplied to the particle velocity vector. This technique is appropriate in the limit that one bounce time is less than the ion slowing down time, which is less than the ion pitch angle scattering time, i.e. $\tau _b \ll \tau _s \ll \tau _i$.

To understand some of the behaviour observed, we use conventions from early mirror work to note that fast ions in WHAM will be ‘superadiabatic’, meaning that the gyrophase angle at the midplane is well correlated between bounces (Cohen et al. 1978a). This is largely due to the small ratio of bounce period to bounce-averaged gyroperiod, $\tau _b/\tau _{\rm ci} \leq 10$ (Cohen et al. 1978a). In contrast, for longer or more collisional experiments, the gyrophase is essentially random from bounce to bounce. Because of this, WHAM may encounter gyro-bounce resonance effects. At low field in WHAM, the bounce period for 25 keV D ions injected at 45$^\circ$ is six times the mean gyroperiod, allowing for a 6 : 1 resonance.

Figure 18 shows numerical calculations of the gyro-bounce ratio of injected NBI deuterons into WHAM at 0.32 T and 0.86 T. For the high-field case, the particle slows along the red line while preserving its initial pitch angle, as expected. In contrast, the low field case shows the particle tracking the 6 : 1 gyro-bounce resonant surface. We call these ‘resonant superadiabatic’ particles, as this resonance between bounce and gyrofrequency can shift the pitch angle of the ion. The resonant condition breaks when the bounce-to-bounce fluctuations in magnetic moment $\Delta {\mu }$ (which are proportional to $\alpha$, shown in figure 17) can no longer compensate for the decrease in $\mu$ due to the slowing. This effect might be measured in a less peaked fusion proton axial intensity profile as fast-ions are deflected towards larger $v_\perp /v_{\|}$.

Figure 18. Numerical calculation of the ratio of bounce period to bounce-averaged gyroperiod across velocity space. The pink line tracks the particles at each bounce as they slow over 1 ms with $\tau _s = 2$ ms.

3.7. Performance Predictions

The FBIS code can be combined with the methods in §§ 3.3 and 3.4 to produce an upper-limit classical prediction for WHAM performance, with results displayed in figure 19. The inputs to the model are shown in the top panel, with the mirror-ratio profile calculated from the vacuum magnetic field, and an arbitrary auxiliary power profile representing strong electron cooling at the edge. With that magnetic geometry and cooling profile, results from the FBIS code are interpolated to each flux surface. pyFIDASIM next calculates the deposition rate given an initial density profile, and uses the FBIS particle confinement time to calculate the steady-state density. After this, the magnetic equilibrium is recalculated using this new pressure profile and the process is iterated. Finally, the bottom panel of figure 19 predicts the upper-limit Deuterium-Deuterium (DD) fusion neutron rate for this configuration of $6\times 10^{13}\ {\rm s}^{-1}$.

Figure 19. The FBIS, Pleides and pyFidasim performance predictions. Assumptions and results for WHAM with $B_0=0.86$ T, with 1 MW of 25 keV NBI onto an ECH target plasma.

As FBIS calculates the equilibrium solution, it does not provide information on the time to reach this equilibrium state other than the particle confinement time. We know that, for large $\chi >0.5$, CX will dominate over the fuelling and this method will not be accurate. It will be more useful in calculating a next step experiment with high-energy neutral beams.

3.7.1. CQL3D-mirror Results

In addition to not modelling temporal behaviour, the FBIS code also does not include a model for the high-harmonic ion cyclotron heating. CQL3D-m (Harvey et al. 2016), a special version of CQL3D (see Harvey & McCoy 1992) for mirror geometry, is a bounce-averaged Fokker–Planck solver that, like FBIS, includes a self-consistent calculation of $E_\parallel$ (Petrov et al. 2023), but does not yet include the potential drop in the expanders outside the mirror. Rather, a simple potential sheath drop to maintain ambipolar losses is used. It has been used here to model a combined neutral beam and second harmonic RF heating scenario, including subtleties like the half and third energy components of the neutral beams. While the details of the implementation will be given in a subsequent publication, initial benchmarking has shown good agreement between the semi-analytic FBIS and CQL3D-m for NBI only, predicting similar ion confinement times and parallel electrostatic potential drops. CQL3D-m also has implemented an advanced energy-dependent cross-section for CX such that the low-energy ions are pumped out of plasma by CX with injected neutrals, while very fast ions (that are accelerated by RF above injection energy) are less likely to CX with the neutrals. This is separate but similar to the spatially dependent effect detailed in § 3.3. Important features of the simulation, including the NBI deposition from NFREYA and ray tracing from GENRAY, are shown in figure 20

Figure 20. CQL3D-m uses NFREYA to model NBI fuelling and GENRAY-C to track ray trajectories. For clarity, only one ray is shown here, with reduced number of bounces. Labelled are the first through fourth ion cyclotron harmonics.

Preliminary simulations were initialized with a flat density profile at $6\times 10^{19}\ {\rm m}^{-3}$, and a nominal 100 kW each of NBI and RF heating. The runs evolve with NBI fuelling balancing the particles lost through the mirror loss cone and due to CX, resulting in a final density of $2\times 10^{19}$ m$^{-3}$. The temporal evolution of the absorbed power, plasma $\beta$ and DD fusion neutron rate is shown in figure 21. Only a small fraction of the injected neutral beam is captured, and the fast wave is mostly absorbed by Landau damping and transit time magnetic pumping on the electrons, with only 10 %–20 % deposited into the fast-ions (mostly at second harmonic resonance), although the RF power partition is gradually shifting towards ions. As the accelerated fast-ions have long confinement times, vastly reduced NBI and RF powers compared with experiment values help prevent the simulation from rapidly becoming unrealistic. Indeed, by 200 ms, the plasma has already reached $\beta = 0.5$.

Figure 21. Time-dependent simulation of absorbed RF and NBI power, and DD neutron yield. Central $\beta$ corresponds to the right-hand axis. Note the magnetic equilibrium uses only the vacuum field and is not self-consistently calculated, and so only results for small beta are reasonable.

As the second harmonic resonance condition is collocated with the turning points of the 45$^\circ$ injected NBI ions, the RF preferentially increases the energy of the sloshing ion distribution. The sloshing ions provide a confining potential in the central region, visible in the axial plot of figure 22. The resultant midplane ion distribution functions are shown in figure 23, demonstrating substantial diffusion to high energies at the 45$^\circ$ pitch angle of the injected neutral beams, and the notable lack of an ambipolar hole. These electrostatically trapped ions reduce the severity of the DCLC instability, although to what extent will be addressed in future work.

Figure 22. A snapshot of the spatial profiles from the CQL3D-m simulation at 200 ms. Radial profiles (left) are taken at the midplane. Axial profiles (right) are shown for flux surface with $\rho /a=0.17$. In the right panel plot $z/z_{\max } = 0$ corresponds to the midplane, while $z/z_{\max } = 1.0$ is at the mirror coil position.

Figure 23. Ion (a,b) and electron (c,d) distribution functions from a CQL3D-m simulation after 20 ms (a,c) and 200 ms (b,d). Distributions are shown at $\rho =0, z=0$ midplane point. Dashed black lines correspond to trapped–passing boundaries in the absence of the axial electric potential, for reference. The solid black lines include the effect of $\phi (z)$. Note that the confining region in the ion distribution formed by an increasing population of sloshing ions grows in time. The fast ion birth energy 25 keV is at the top of the green region in $f_D(u_{\|},u_\perp )$ plots, at $u/c \sim 0.005$.

These very promising initial results show that $\beta =0.5$ can be reached in approximately 200 ms, but that further evolution would continue to even higher densities and pressures. It should be noted that in these simulations the gyro-radius cutoff was implemented to reduce the growth of high-energy tails, however, the effective ion temperature is still increasing because of the tails getting wider (more populated). This suggests in an experiment that the power might be feedback controlled and lowered to not exceed an effective $\beta$ limit. Unlike the combined results at the beginning of § 3.7, these results do not yet include the diamagnetic modification to the equilibrium that will increase the mirror ratio and lead to even better confinement. Furthermore, the ray-tracing approximations may not be quite valid as the rays fill up the plasma volume with no interference effects included. However, by using the initial azimuthal spectrum, we were able to obtain reasonably wide profiles of power deposition. A next step will be to use a full-wave code like AORSA3D (special version for cylinder geometry) to couple to the Fokker–Planck solver and investigate if this approximation is reasonable. Other next steps include modelling the electron cyclotron resonance heating, including the role of CX halo particles and edge effects.

3.8. Concluding Remarks

MHD stability is paramount: forming a stable plasma is critical for assessing confinement and establishing a platform for testing the viability of alternative stabilization concepts like feedback or mantel line-tying. Stable plasmas will be a prerequisite for a number of additional goals, such as demonstrating the HHFW acceleration of fast ions, and the direct energy conversion of real exhaust distributions. Unfortunately, using biased limiters for $m=1$ stability in WHAM (duplicating the GDT method in the CM regime) does not fully retire the MHD stability risk for a next step device when plasma confinement improves even further. While the physics principle of sheared-flow stabilization has already been shown to work experimentally, implementing it in the next step may require a different set of techniques or actuators appropriate to larger, better confined plasmas. Because of its increased size, such a device will have access to additional stabilization strategies appropriate for a reactor grade plasma including: using a low-temperature mantel plasma with good electrical connectivity to the ends to assist line-tying; the use of close fitting shells theoretically expected to stabilize/ larger plasmas at high $\beta$; saddle feedback stabilization or tail wagging; more precise control of $T_e$ and $E_r$ profiles, using steerable ECH deposition as an actuator.

Scientific validation of the fast-ion confinement against models is not trivial: simply measuring the total ion content and dividing by the neutral beam injection current is not particularly informative because of the strong role that CX plays in the neutral beam deposition and losses. Measurements of the ion distribution function with nuclear product diagnostics, fast-ion D-$\alpha$ emission and CX neutrals will focus on observing these effects. For these reasons, benchmarking WHAM results against hybrid-PIC calculations and gyrofluid calculations that can include CX effects is essential. Validating a CQL3D-mirror, FBIS, FIDASIM model will also require very careful measurements of beam shine-through in order to accurately assess confinement.

In summary, the WHAM device is a prototype for a future $Q\sim 1$ class simple axisymmetric mirror. It will operate in the low collisionality ‘classical mirror’ limit in which ion pitch angle scattering dominates over electron slowing. We have presented initial modelling predictions for the WHAM experiment that couples fast-ion modelling from the FBIS package with anisotropic MHD equilibria from Pleiades and beam-ion fuelling from pyFIDASIM to predict the upper limits of instability-free performance. We have also included results from CQL3D-m showing the efficacy of second harmonic ICRF heating.