Hostname: page-component-586b7cd67f-t8hqh Total loading time: 0 Render date: 2024-11-27T11:26:41.225Z Has data issue: false hasContentIssue false

Physics basis for the Wisconsin HTS Axisymmetric Mirror (WHAM)

Published online by Cambridge University Press:  20 September 2023

D. Endrizzi
Affiliation:
University of Wisconsin-Madison
J.K. Anderson
Affiliation:
University of Wisconsin-Madison
M. Brown
Affiliation:
Swarthmore College
J. Egedal
Affiliation:
University of Wisconsin-Madison
B. Geiger
Affiliation:
University of Wisconsin-Madison
R.W. Harvey
Affiliation:
CompX
M. Ialovega
Affiliation:
University of Wisconsin-Madison
J. Kirch
Affiliation:
University of Wisconsin-Madison
E. Peterson
Affiliation:
MIT
Yu.V. Petrov
Affiliation:
CompX
J. Pizzo
Affiliation:
University of Wisconsin-Madison
T. Qian
Affiliation:
University of Wisconsin-Madison Princeton University
K. Sanwalka
Affiliation:
University of Wisconsin-Madison
O. Schmitz
Affiliation:
University of Wisconsin-Madison
J. Wallace
Affiliation:
University of Wisconsin-Madison
D. Yakovlev
Affiliation:
University of Wisconsin-Madison
M. Yu
Affiliation:
University of Wisconsin-Madison
C.B. Forest*
Affiliation:
University of Wisconsin-Madison
*
Email address for correspondence: [email protected]
Rights & Permissions [Opens in a new window]

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.

Type
Research Article
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
Copyright © The Author(s), 2023. Published by Cambridge University Press

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 Reference Fowler1981; Post Reference Post1987), 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. Reference Beklemishev, Bagryansky, Chaschin and Soldatkina2010; Bagryansky et al. Reference Bagryansky, Anikeev, Beklemishev, Donin, Ivanov, Korzhavina, Kovalenko, Kruglyakov, Lizunov, Maximov, Murakhtin, Prikhodko, Pinzhenin, Pushkareva, Savkin and Zaytsev2011, Reference Bagryansky, Shalashov, Gospodchikov, Lizunov, Maximov, Prikhodko, Soldatkina, Solomakhin and Yakovlev2015; Yakovlev et al. Reference Yakovlev, Shalashov, Gospodchikov, Maximov, Prikhodko, Savkin, Soldatkina, Solomakhin and Bagryansky2018; Soldatkina et al. Reference Soldatkina, Maximov, Prikhodko, Savkin, Skovorodin, Yakovlev and Bagryansky2020). 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. Reference Whyte, Minervini, LaBombard, Marmar, Bromberg and Greenwald2016; Whyte Reference Whyte2019) 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. Reference Egedal, Endrizzi, Forest and Fowler2022). 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.13.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 Reference Bagryansky, Beklemishev and Postupaev2018), 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. Reference Yakovlev, Shalashov, Gospodchikov, Maximov, Prikhodko, Savkin, Soldatkina, Solomakhin and Bagryansky2018; Soldatkina et al. Reference Soldatkina, Maximov, Prikhodko, Savkin, Skovorodin, Yakovlev and Bagryansky2020) 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 Reference Boldyrev, Forest and Egedal2020). 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. Reference Egedal, Endrizzi, Forest and Fowler2022). 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. Reference Kremeyer, Flesch, Schmitz, Schlisio and Wenzel2020).

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. Reference Ialovega2023) 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. Reference Yakovlev, Shalashov, Gospodchikov, Solomakhin, Savkin and Bagryansky2017). 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 Reference Kesner1973; Simonen et al. Reference Simonen1983). In addition, angled injection decreases the pressure anisotropy compared with perpendicular injection, mitigating the Alfvén ion cyclotron (AIC) instability (Smith Reference Smith1984). 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 Reference Stix1992). 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. Reference Petty, Baity, de Grassie, Mau, Pinsker, Porkolab and Prater2001; Rosenberg et al. Reference Rosenberg, Menard, Wilson, Medley, Andre, Phillips, Darrow, LeBlanc, Redi, Fisch, Team, Harvey, Mau, Jaeger, Ryan, Swain, Sabbagh and Egedal2004) and also simulated by Fokker–Plank calculations (Harvey, Petrov & Forest Reference Harvey, Petrov and Forest2016).

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. Reference Lau, Berry, Jaeger and Bertelli2019). 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. Reference Green, Berry, Simpson and Younkin2018), while the AORSA code is also being adapted for the mirror geometry (Jaeger et al. Reference Jaeger, Berry, D'Azevedo, Batchelor and Carter2001).

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 Reference Barr, Moir and Hamilton1982). 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 Reference Bing and Roberts1961; Fowler & Rankin Reference Fowler and Rankin1962, Reference Fowler and Rankin1966). 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 Reference Killeen and Marx1970; Killeen, Mirin & Rensink Reference Killeen, Mirin and Rensink1976; Egedal et al. Reference Egedal, Endrizzi, Forest and Fowler2022) 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 Reference Ivanov and Prikhodko2013). 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. Reference Coensgen, Cummins, Logan, Molvik, Nexsen, Simonen, Stallard and Turner1975; Post Reference Post1987; Wurzel & Hsu Reference Wurzel and Hsu2022). 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. Reference Coensgen, Cummins, Logan, Molvik, Nexsen, Simonen, Stallard and Turner1975). 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. Reference Egedal, Endrizzi, Forest and Fowler2022), 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 Reference Beklemishev2016). 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 Reference Rosenbluth and Longmire1957). 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. Reference Ryutov, Berk, Cohen, Molvik and Simonen2011), including the use of good curvature (Mirnov & Ryutov Reference Mirnov and Ryutov1979; Bagryansky et al. Reference Bagryansky, Anikeev, Beklemishev, Donin, Ivanov, Korzhavina, Kovalenko, Kruglyakov, Lizunov, Maximov, Murakhtin, Prikhodko, Pinzhenin, Pushkareva, Savkin and Zaytsev2011), cusp boundaries (Prater Reference Prater1971; Ferron et al. Reference Ferron, Wong, Dimonte and Leikind1983; Kotelnikov et al. Reference Kotelnikov, Ivanov, Yakovlev, Chen and Zeng2019), ponderomotive effects (Ferron et al. Reference Ferron, Goulding, Nelson, Intrator, Wang, Severn, Hershkowitz, Brouchous, Pew, Breun and Majeski1987), rotating magnetic fields (Seemann, Be'ery & Fisher Reference Seemann, Be'ery and Fisher2018; Zhu et al. Reference Zhu, Shi, Yang, Zheng, Luo, Ying and Sun2019), close fitting conducting shells (Kesner Reference Kesner1985; Beklemishev Reference Beklemishev2016; Kotelnikov et al. Reference Kotelnikov, Zeng, Prikhodko, Yakovlev, Zhang, Chen and Yu2022) and feedback (Prater Reference Prater1971; Be'ery, Seemann & Fisher Reference Be'ery, Seemann and Fisher2014; Be'ery & Seemann Reference Be'ery and Seemann2015; Beklemishev Reference Beklemishev2017).

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. Reference Ryutov, Berk, Cohen, Molvik and Simonen2011)

(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. Reference Cho2005, Reference Cho2006) and vortex stabilization (Sakai, Yasaka & Itatani Reference Sakai, Yasaka and Itatani1993; Beklemishev et al. Reference Beklemishev, Bagryansky, Chaschin and Soldatkina2010; Bagryansky et al. Reference Bagryansky, Anikeev, Beklemishev, Donin, Ivanov, Korzhavina, Kovalenko, Kruglyakov, Lizunov, Maximov, Murakhtin, Prikhodko, Pinzhenin, Pushkareva, Savkin and Zaytsev2011; Yakovlev et al. Reference Yakovlev, Shalashov, Gospodchikov, Maximov, Prikhodko, Savkin, Soldatkina, Solomakhin and Bagryansky2018). 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 Reference Post1961; Kelley Reference Kelley1967) theoretically described (Fowler Reference Fowler1981) and experimentally measured (Post Reference Post1987). 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. Reference Severn, HershkowitE, Nelson and Pew1983). 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. Reference Sakai, Yasaka and Itatani1993; Bagryansky et al. Reference Bagryansky, Beklemishev and Postupaev2018). Applying the opposite polarity to reinforce the natural $E_r$ was experimentally challenging.

According to Beklemishev et al. (Reference Beklemishev, Bagryansky, Chaschin and Soldatkina2010), 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. Reference Yakovlev, Shalashov, Gospodchikov, Solomakhin, Savkin and Bagryansky2017) and GAMMA-10 (Shidara et al. Reference Shidara, Nakashima, Imai, Kariya, Minami, Ichimura and Yoshikawa2011) 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 Reference Smirnov, Harvey and Kupfer1994) 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. Reference Cho2006; Yakovlev et al. Reference Yakovlev, Shalashov, Gospodchikov, Maximov, Prikhodko, Savkin, Soldatkina, Solomakhin and Bagryansky2018). 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. (Reference Geiger, Stagner, Heidbrink, Dux, Fischer, Fujiwara, Garcia, Jacobsen, van Vuuren, Karpushov, Liu, Schneider, Sfiligoi, Poloskei and Weiland2020) 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. Reference Egedal, Endrizzi, Forest and Fowler2022) 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. Reference Egedal, Endrizzi, Forest and Fowler2022). (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.Footnote 1 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 Reference Baldwin1977; Ferron & Wong Reference Ferron and Wong1984). 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. Reference Francisquez, Rosen, Mandell, Hakim, Forest and Hammett2023; Petrov et al. Reference Petrov, Harvey, Forest and Anderson2023), 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 (Reference Kotelnikov, Chernoshtanov and Prikhodko2017), 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 Reference Ferron and Wong1984), 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 Reference Ferron and Wong1984), 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 Reference Yushmanov1966; Kesner Reference Kesner1973). 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 (Reference Smith1984)

(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. Reference Katsumata, Ichimura, Inutake, Hojo, Mase and Tamano1996). 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 Reference Ivanov and Prikhodko2013), and that our discussion here may be too pessimistic. Its likely that sloshing ion distributions may have negligible kinetic instabilities (Post & Ryutov Reference Post and Ryutov1995). 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. Reference Simonen2008). 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. Reference Horton, Fu, Ivanov and Beklemishev2010; Fowler, Moir & Simonen Reference Fowler, Moir and Simonen2017). 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. Reference Cho2005, Reference Cho2006).

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. Reference Schmitz2016). 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. Reference Zweben1997). 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 Reference Hastie, Hobbs and Taylor1969). As discussed by Cohen, Rowlands & Foote (Reference Cohen, Rowlands and Foote1978b), 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. Reference Cohen, Rowlands and Foote1978b). 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. Reference Cohen, Rensink, Cutler and Mirin1978a). 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. Reference Cohen, Rensink, Cutler and Mirin1978a). 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. Reference Harvey, Petrov and Forest2016), a special version of CQL3D (see Harvey & McCoy Reference Harvey and McCoy1992) for mirror geometry, is a bounce-averaged Fokker–Planck solver that, like FBIS, includes a self-consistent calculation of $E_\parallel$ (Petrov et al. Reference Petrov, Harvey, Forest and Anderson2023), 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.

Acknowledgements

For transparency of interests, the authors would like to state that J.K.A., C.B.F. and O.S. are cofounders of Realta Fusion. The authors would like to thank P. Catto for their help with the anisotropic equilibrium equations in Appendix A. The authors would also like to thank D. Ryutov for their helpful comments during the review process.

Editor Per Helander thanks the referees for their advice in evaluating this article.

Declaration of interest

The authors report no conflict of interest.

Funding

The information, data or work presented herein was funded in part by the Advanced Research Projects Agency-Energy (ARPA-E), U.S. Department of Energy under Award Numbers DE-AR0001258, DE-AR0001261, and U.S. Department of Energy under Award Number DE-FG02-ER54744. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Appendix A. Equations for Anisotropic Equilibrium

For an anisotropic, low flow plasma in a magnetic field ${\boldsymbol B} = B {\boldsymbol b}$, MHD pressure balance ${\boldsymbol J}\times {\boldsymbol B} = \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol P}$ can be rewritten

(A1)\begin{align} {\boldsymbol J}\times {\boldsymbol B} & = \boldsymbol{\nabla} \boldsymbol{\cdot} [ p_{\|} {\boldsymbol b}{\boldsymbol b} + p_\perp ( {\mathbb 1} - {\boldsymbol b}{\boldsymbol b}) ] \nonumber\\ & = \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ p_\perp {\mathbb 1} + (p_{\|} - p_\perp)\frac{\boldsymbol B}{B}{\boldsymbol b} \right]\nonumber\\ & = \boldsymbol{\nabla} p_\perp+ {\boldsymbol B} \boldsymbol{\cdot} \boldsymbol{\nabla} \left[ \frac{p_{\|}-p_\perp}{B} {\boldsymbol b} \right] \nonumber\\ & = \boldsymbol{\nabla} p_\perp+ (p_{\|} - p_\perp) {\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol b} + {\boldsymbol b} {\boldsymbol B} \boldsymbol{\cdot} \boldsymbol{\nabla} \left(\frac{p_{\|}-p_\perp}{B}\right). \end{align}

The parallel component of the right-hand side is

(A2)\begin{equation} {\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} p_\perp+ {\boldsymbol B} \boldsymbol{\cdot} \boldsymbol{\nabla} \left(\frac{p_{\|}-p_\perp}{B}\right) = 0, \end{equation}

while the perpendicular current density ${\boldsymbol J}_\perp$ is

(A3)\begin{equation} {\boldsymbol J}_\perp= \frac{\boldsymbol B}{B^2} \times [ \boldsymbol{\nabla} p_\perp+ (p_{\|}-p_\perp) {\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol b}]. \end{equation}

Using the following property of vector calculus:

(A4)\begin{equation} ( {\mathbb 1} - {\boldsymbol b}{\boldsymbol b}) \boldsymbol{\cdot} \boldsymbol{\nabla} \times {\boldsymbol b} = {\boldsymbol b} \times [(\boldsymbol{\nabla} \times {\boldsymbol b}) \times {\boldsymbol b} ] = {\boldsymbol b} \times ({\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol b}), \end{equation}

and also this

(A5)\begin{align} ({\mathbb 1} - {\boldsymbol b}{\boldsymbol b}) \boldsymbol{\cdot} \boldsymbol{\nabla} \times {\boldsymbol b} & = \frac{{\mathbb 1} - {\boldsymbol b}{\boldsymbol b}}{B} \boldsymbol{\cdot} (\boldsymbol{\nabla} \times {\boldsymbol B} - \boldsymbol{\nabla} B \times {\boldsymbol b}) \nonumber\\ & = \frac{\mu_0 {\boldsymbol J}_\perp}{B} + \frac{{\boldsymbol b} \times \boldsymbol{\nabla} B}{B}, \end{align}

both can be combined into

(A6)\begin{equation} {\boldsymbol b} \times ({\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol b} ) = \frac{\mu_0 {\boldsymbol J}_\perp}{B} + \frac{{\boldsymbol b} \times \boldsymbol{\nabla} B}{B}. \end{equation}

As a result, ${\boldsymbol J}_\perp$ can be rewritten as

(A7)\begin{equation} \sigma {\boldsymbol J}_\perp{\equiv} \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_{\|}-p_\perp) \right]. \end{equation}

The perpendicular current density has $\boldsymbol{\nabla} p_\perp$ and $\boldsymbol{\nabla} B$ components. If you know $p_\perp = p_\perp (\psi,B)$, with the flux function $\psi$ satisfying ${\boldsymbol B}\boldsymbol{\cdot} \boldsymbol{\nabla} \psi = 0$, then ${\boldsymbol J}_\perp$ can be written as

(A8)\begin{equation} \sigma {\boldsymbol J}_\perp= \frac{\boldsymbol B}{B^2} \times \left[ \frac{\partial p_\perp}{\partial \psi} \boldsymbol{\nabla} \psi + \frac{\partial p_\perp}{\partial B} \boldsymbol{\nabla} B + \frac{\boldsymbol{\nabla} B}{B}(p_{\|}-p_\perp) \right], \end{equation}

which has separable $\boldsymbol{\nabla} \psi$ and $\boldsymbol{\nabla} B$ components

(A9)\begin{gather} \sigma {\boldsymbol J}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} \psi = \frac{1}{B} \left[ \frac{\partial p_\perp}{\partial B} + \frac{p_{\|} - p_\perp}{B} \right] {\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} B \times \boldsymbol{\nabla} \psi, \end{gather}
(A10)\begin{gather}\sigma {\boldsymbol J}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} B = \frac{1}{B} \frac{\partial p_\perp}{\partial \psi} {\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} B. \end{gather}

If you also know that $p_{\|} = p_{\|}(\psi,B)$, then parallel pressure balance gives

(A11)\begin{equation} \frac{\partial p_\perp}{\partial B} + B \frac{\partial }{\partial B} \left[\frac{(p_{\|}-p_\perp)}{B}\right] = 0 = \frac{\partial p_{\|}}{\partial B} - \frac{p_{\|}-p_\perp}{B}. \end{equation}

The current, of course, $\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol J} = 0$ with ${\boldsymbol J} = {\boldsymbol J}_\perp + J_{\|} {\boldsymbol b}$.

A.1 Rotation effects

To add the effects of ${\boldsymbol E}\times {\boldsymbol B}$ rotation into the MHD equilibrium, we assume axisymmetry to allow the flow to be sonic. Sonic rotation adds a term to the right side of the MHD pressure balance equation of the form $\boldsymbol{\nabla} \boldsymbol{\cdot} (m_i n_i {\boldsymbol vv}) = \rho _i {\boldsymbol v}\boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol v}$, with $n_i$ the ion density, $m_i$ the ion mass, $\rho _i$ the mass density and ${\boldsymbol v}$ the ion flow velocity. Using cylindrical $r$, $\theta$, $z$, coordinates and writing the magnetic field as

(A12)\begin{equation} {\boldsymbol B} = \boldsymbol{\nabla} \psi(r,\theta) \times \boldsymbol{\nabla} \theta = \frac{1}{r} \frac{\partial \psi}{\partial r} \boldsymbol{\nabla} z - \frac{1}{r}\frac{\partial \psi}{\partial z} \boldsymbol{\nabla} r, \end{equation}

then the ion flow velocity is

(A13)\begin{equation} {\boldsymbol v} = \omega(\psi,\ell) r^2 \boldsymbol{\nabla} \theta. \end{equation}

We now assume the rotation is a result of an electrostatic field

(A14)\begin{equation} {\boldsymbol E}={-}\boldsymbol{\nabla} \varPhi(\psi,\ell) ={-} \frac{\partial \varPhi(\psi,\ell)}{\partial \psi}\boldsymbol{\nabla} \psi - \frac{\partial \varPhi(\psi,\ell)}{\partial \ell} \boldsymbol{b}, \end{equation}

where $\boldsymbol {b} \equiv \boldsymbol {B}/|B|$, thus

(A15)\begin{equation} {\boldsymbol E}_\perp={-} {{\boldsymbol v} \times {\boldsymbol B}} ={-}\omega r^2 \boldsymbol{\nabla} \theta \times (\boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} \theta) ={-}\omega \boldsymbol{\nabla} \psi, \end{equation}

and the angular velocity becomes

(A16)\begin{equation} \omega = \frac{\partial \varPhi(\psi,\ell)}{\partial \psi}. \end{equation}

We note here that a ideal and resistive MHD do not hold since $E_\parallel = -({T_e}/{e n}) \boldsymbol {b} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol {n}$

Let us expand the ${\boldsymbol v} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol v}$ term

(A17)\begin{align} {\boldsymbol v} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol v} & = \omega^2 r^2 \boldsymbol{\nabla} \theta \boldsymbol{\cdot} \boldsymbol{\nabla} (r^2 \boldsymbol{\nabla} \theta ) \nonumber\\ & = \omega^2 r^3 \boldsymbol{\nabla} \theta \boldsymbol{\cdot} [ (\boldsymbol{\nabla} r)(\boldsymbol{\nabla} \theta) - (\boldsymbol{\nabla} \theta)(\boldsymbol{\nabla} r) ] ={-}\omega^2 r \boldsymbol{\nabla} r. \end{align}

Now separate the rotation term into its parallel component

(A18)\begin{equation} \rho_i {\boldsymbol v} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol v} \boldsymbol{\cdot} {\boldsymbol b} ={-} \frac{\rho_i \omega^2 r}{B} {\boldsymbol B}\boldsymbol{\cdot} \boldsymbol{\nabla} r = \left( \frac{\rho_i \omega^2}{B} \right) \left( \frac{\partial \psi}{\partial z} \right), \end{equation}

and its perpendicular component

(A19)\begin{equation} {\boldsymbol b} \times [ (\rho_i {\boldsymbol v} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol v}) \times {\boldsymbol b} ] ={-} \rho_i \omega^2 r ({\mathbb 1} - {\boldsymbol b}{\boldsymbol b}) \boldsymbol{\cdot} \boldsymbol{\nabla} r. \end{equation}

Using the new term, parallel force balance becomes

(A20)\begin{equation} {\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} p_\perp+ {\boldsymbol B} \boldsymbol{\cdot} \boldsymbol{\nabla} \left[ \frac{p_{\|} - p_\perp}{B}\right] = \rho_i \omega^2 r {\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} r = \rho_i {\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} \left( \frac{\omega^2 r^2}{2}\right). \end{equation}

Substituting the new term in the expression for perpendicular current

(A21)\begin{align} \frac{\boldsymbol B}{B^2} \times (\rho_i {\boldsymbol v} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol v}) & ={-} \frac{\rho_i \omega^2 r}{B^2} {\boldsymbol B}\times \boldsymbol{\nabla} r \nonumber\\ & ={-} \frac{\rho_i \omega^2 r}{B^2} \left( \frac{\partial \psi}{\partial r} \right) \boldsymbol{\nabla} \theta, \end{align}

so that the perpendicular current becomes

(A22)\begin{equation} {\boldsymbol J}_\perp= \frac{\boldsymbol B}{B^2} \times [ \boldsymbol{\nabla} p_\perp+ (p_{\|}-p_\perp) {\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol b} - \rho_i \omega^2 r \boldsymbol{\nabla} r ]. \end{equation}

Again using (A6)

(A23)\begin{equation} {\boldsymbol b} \times ({\boldsymbol b} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol b}) = \frac{\mu_0 {\boldsymbol J}_\perp}{B} + \frac{{\boldsymbol b} \times \boldsymbol{\nabla} B}{B}, \end{equation}

and substituting

(A24)\begin{align} \sigma {\boldsymbol J}_\perp & \equiv \left[ 1 + \frac{\mu_0}{B^2} (p_\perp{-} p_{\|}) \right] {\boldsymbol J}_\perp \nonumber\\ & = \frac{\boldsymbol B}{B^2} \times \left[ \boldsymbol{\nabla} p_\perp+ \frac{p_{\|} - p_\perp}{B} \boldsymbol{\nabla} B - \rho_i \omega^2 r \boldsymbol{\nabla} r \right]. \end{align}

Using axisymmetry this becomes

(A25)\begin{equation} \sigma {\boldsymbol J}_\perp= \frac{\boldsymbol{\nabla} \theta \boldsymbol{\nabla} \psi}{B^2} \boldsymbol{\cdot} \left[ \boldsymbol{\nabla} p_\perp+ \frac{p_{\|}-p_\perp}{B} \boldsymbol{\nabla} B - \rho_i \omega^2 r \boldsymbol{\nabla} r \right]. \end{equation}

Notice that because of axisymmetry, if $p_\perp = p_\perp (\psi,B)$, then $\sigma {\boldsymbol J}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} \psi = 0 = \sigma {\boldsymbol J}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla} B$.

Appendix B. Charge Exchange in Fast Beam-Ion Solver

At the NBI energy of $E_{\rm nbi} = 25$ keV, the CX cross-section is substantially larger than the sum of the ion and electron impact ionization collisions. We parameterize this as $\chi \equiv \sigma _{\rm cx} / (\sigma _{\rm cx} + \sigma _{i/e})$, where at 25 keV $\chi \simeq 0.8$. For each ion deposited by electron/ion impact ionization, $\chi / (1-\chi )$ ions are removed from elsewhere in the distribution function and redeposited by charge exchange at the injection energy and pitch angle. (59) of Egedal et al. (Reference Egedal, Endrizzi, Forest and Fowler2022) can be rewritten as

(B1)\begin{align} & {-}\beta \lambda_j \frac{v_c^3 \tilde{g}_1}{v^3} f_j + \frac{1}{v^2} \frac{\partial}{\partial v}\left[ (v^3 + v_c^3 \tilde{h}) f_j + \frac{v_c^3 \tilde{g}_2}{2} \frac{\partial f_j}{\partial v}\right] \nonumber\\ & \quad {-} C w_j f_j ={-}\tau_s S_j \frac{\delta(v-v_0)}{v^2}. \end{align}

The final term on the left-hand side reflects CX losses and is proportional to the distribution function weighted by the normalized axial dependence of each eigenfunctions’ ‘residence time’ within the NBI beam path. To solve for the constant $C$, we reason that the sum over $j$ of the integral of this term must add up to $\chi$ times the total number of injected particles

(B2)\begin{equation} C \sum_j w_j \int 4\pi v^2 \, {\rm d} v \int \, {\rm d} \eta f_j = \chi \sum_j \tau_s S_j \frac{\delta(v-v_0)}{v^2}. \end{equation}

The right-hand side can be rewritten

(B3)\begin{equation} \chi \tau_s \sum_j S_j \frac{\delta(v-v_0)}{v^2} = \chi \tau_s S_0 =\chi \tau_s \frac{{\rm d}N}{{\rm d}{\boldsymbol x} \,{\rm d}t} = \chi \tau_s \frac{I/e}{V}. \end{equation}

The left-hand side can be integrated up into a weighted density $n^*$ and then using

(B4)\begin{equation} \tau_{\rm cx} \equiv \frac{n^* V}{I/e}, \end{equation}

our constant $C$ becomes

(B5)\begin{equation} C = \chi \frac{\tau_s}{\tau_{\rm cx}}. \end{equation}

This is calculated numerically in two steps: first solving the equation with the case where $\tau _{\rm cx} = \tau _s$, and then solving again once $\tau _{\rm cx}$ is known. We use

(B6)\begin{equation} N = V \int 4 \pi v^2 \, {\rm d} v \int_0^1 f(v,\eta) \, {\rm d} \eta. \end{equation}

It is worth noting that the effect here is precisely the concept conceived of to pump out ions in thermal barrier plasmas (Carlson et al. Reference Carlson1981) that rely upon a sloshing ion distribution and prevent buildup of plasma in the Yushmanov potential (Yushmanov Reference Yushmanov1966; Kesner Reference Kesner1973) by pumping out the weakly confined particles.

References

Bagryansky, P.A., Anikeev, A.V., Beklemishev, A.D., Donin, A.S., Ivanov, A.A., Korzhavina, M.S., Kovalenko, Y.V., Kruglyakov, E.P., Lizunov, A.A., Maximov, V.V., Murakhtin, S.V., Prikhodko, V.V., Pinzhenin, E.I., Pushkareva, A.N., Savkin, V.Y. & Zaytsev, K.V. 2011 Confinement of hot ion plasma with $\beta \le 0.6$ in the gas dynamic trap. Fusion Sci. Technol. 59 (1T), 3135.CrossRefGoogle Scholar
Bagryansky, P.A., Beklemishev, A.D. & Postupaev, V.V. 2018 Encouraging results and new ideas for fusion in linear traps. J. Fusion Energy 38 (1), 162181.CrossRefGoogle Scholar
Bagryansky, P.A., Shalashov, A.G., Gospodchikov, E.D., Lizunov, A.A., Maximov, V.V., Prikhodko, V.V., Soldatkina, E.I., Solomakhin, A.L. & Yakovlev, D.V. 2015 Threefold increase of the bulk electron temperature of plasma discharges in a magnetic mirror device. Phys. Rev. Lett. 114 (20), 205001.CrossRefGoogle Scholar
Baldwin, D.E. 1977 End-loss processes from mirror machines. Rev. Mod. Phys. 49 (2), 317339.CrossRefGoogle Scholar
Barr, W.L., Moir, R.W. & Hamilton, G.W. 1982 Experimental results from a beam direct converter at 100 kV. J. Fusion Energy 2 (2), 131143.CrossRefGoogle Scholar
Be'ery, I. & Seemann, O. 2015 Feedback effect on flute dynamics in a mirror machine. Plasma Phys. Control. Fusion 57 (8), 085005.CrossRefGoogle Scholar
Be'ery, I., Seemann, O. & Fisher, A. 2014 Active feedback stabilization of multimode flute instability in a mirror trap. Plasma Phys. Control. Fusion 56 (7), 075002.CrossRefGoogle Scholar
Beklemishev, A.D. 2016 Diamagnetic ‘bubble’ equilibria in linear traps. Phys. Plasmas 23 (8), 082506.CrossRefGoogle Scholar
Beklemishev, A.D. 2017 Tail-waving system for active feedback stabilization of flute modes in open traps. Fusion Sci. Technol. 59 (1T), 9093.CrossRefGoogle Scholar
Beklemishev, A.D., Bagryansky, P.A., Chaschin, M.S. & Soldatkina, E.I. 2010 Vortex confinement of plasmas in symmetric mirror traps. Fusion Sci. Technol. 57 (4), 351360.CrossRefGoogle Scholar
Bing, G.F. & Roberts, J.E. 1961 End-losses from mirror machines. Phys. Fluids 4 (8), 1039.CrossRefGoogle Scholar
Boldyrev, S., Forest, C. & Egedal, J. 2020 Electron temperature of the solar wind. Proc. Natl Acad. Sci. USA 117 (17).CrossRefGoogle ScholarPubMed
Carlson, G.A., et al. 1981 Tandem mirror reactor with thermal barriers. Nucl. Engng Des. 63 (2), 233250.CrossRefGoogle Scholar
Cho, T., et al. 2005 Recent progress in the GAMMA 10 tandem mirror. Fusion Sci. Technol. 47 (1T), 916.CrossRefGoogle Scholar
Cho, T., et al. 2006 Observation and control of transverse energy-transport barrier due to the formation of an energetic-electron layer with sheared $e \times b$ flow. Phys. Rev. Lett. 97, 055001.CrossRefGoogle Scholar
Coensgen, F.H., Cummins, W.F., Logan, B.G., Molvik, A.W., Nexsen, W.E., Simonen, T.C., Stallard, B.W. & Turner, W.C. 1975 Stabilization of a neutral-beam – sustained, mirror-confined plasma. Phys. Rev. Lett. 35 (22), 15011503.CrossRefGoogle Scholar
Cohen, R.H., Rensink, M., Cutler, T.A. & Mirin, A.A. 1978 a Collisional loss of electrostatically confined species in a magnetic mirror. Nucl. Fusion 18, 1229.CrossRefGoogle Scholar
Cohen, R.H., Rowlands, G. & Foote, J.H. 1978 b Nonadiabaticity in mirror machines. Phys. Fluids 21 (4), 627.CrossRefGoogle Scholar
Egedal, J., Endrizzi, D., Forest, C.B. & Fowler, T.K. 2022 Fusion by beam ions in a low collisionality, high mirror ratio magnetic mirror. Nucl. Fusion 62 (12), 126053.CrossRefGoogle Scholar
Ferron, J.R., Goulding, R., Nelson, B.A., Intrator, T., Wang, E.Y., Severn, G., Hershkowitz, N., Brouchous, D., Pew, J., Breun, R.A. & Majeski, R. 1987 Electrostatic end plugging accompanied by a central-cell density increase in an axisymmetric tandem mirror. Phys. Fluids 30 (9), 28552869.CrossRefGoogle Scholar
Ferron, J.R. & Wong, A.Y. 1984 The dependence of the drift cyclotron loss cone instability on the radial density gradient. Phys. Fluids 27 (5), 1287.CrossRefGoogle Scholar
Ferron, J.R., Wong, A.Y., Dimonte, G. & Leikind, B.J. 1983 Interchange stability of an axisymmetric, average minimum-b magnetic mirror. Phys. Fluids 26 (8), 22272233.CrossRefGoogle Scholar
Fowler, T.K. 1981 Mirror theory. In Fusion, vol. 1, chap. 6, pp. 291–356. Academic Press.CrossRefGoogle Scholar
Fowler, T.K., Moir, R.W. & Simonen, T.C. 2017 A new simpler way to obtain high fusion power gain in tandem mirrors. Nucl. Fusion 57 (5), 056014.CrossRefGoogle Scholar
Fowler, T.K. & Rankin, M. 1962 Plasma potential and energy distributions in high-energy injection machines. J. Nucl. Energy. C, Plasma Phys. Accel. Thermonucl. Res. 4 (5), 311320.CrossRefGoogle Scholar
Fowler, T.K. & Rankin, M. 1966 Fusion energy balance in mirror machines. J. Nucl. Energy. C, Plasma Phys. Accel. Thermonucl. Res. 8 (2), 121128.CrossRefGoogle Scholar
Francisquez, M., Rosen, M.H., Mandell, N.R., Hakim, A., Forest, C.B. & Hammett, G.W. 2023 Towards continuum gyrokinetic study of high-field mirrors. arXiv:2305.06372.Google Scholar
Geiger, B., Stagner, L., Heidbrink, W.W., Dux, R., Fischer, R., Fujiwara, Y., Garcia, A.V., Jacobsen, A.S., van Vuuren, A.J., Karpushov, A.N., Liu, D., Schneider, P.A., Sfiligoi, I., Poloskei, P.Z. & Weiland, M. 2020 Progress in modelling fast-ion d-alpha spectra and neutral particle analyzer fluxes using fidasim. Plasma Phys. Control. Fusion 62 (10), 105008.CrossRefGoogle Scholar
Green, D.L., Berry, L.A., Simpson, A.B. & Younkin, T.R. 2018 Kinetic-j: a computational kernel for solving the linearized vlasov equation applied to calculations of the kinetic, configuration space plasma current for time harmonic wave electric fields. Comput. Phys. Commun. 225, 122127.CrossRefGoogle Scholar
Harvey, R.W. & McCoy, M.G. 1992 The CQL3D code. In Proceedings of the IAEA TCM on Advances in Simulation and Modeling of Thermonuclear Plasmas pp. 489–526 (see also CQL3D manual, with corrections https://www.compxco.com/cql3d_manual.pdf).Google Scholar
Harvey, R.W., Petrov, Y.V. & Forest, C.B. 2016 3D distributions resulting from neutral beam, ICRF and EC heating in an axisymmetric mirror. In AIP Conference Proceedings 1771, p. 040002. https://doi.org/10.1063/1.4964187CrossRefGoogle Scholar
Hastie, R.J., Hobbs, G.D. & Taylor, J.B. 1969 Non-adiabatic behavior of particles in inhomogeneous magnetic fields. In Plasma Physics and Controlled Nuclear Fusion Research, vol. 1. IAEA.Google Scholar
Horton, W., Fu, X.R., Ivanov, A. & Beklemishev, A. 2010 Parameter optimization studies for a tandem mirror neutron source. J. Fusion Energy 29 (6), 521526.CrossRefGoogle Scholar
Ialovega, M., et al. 2023 Initial study on thermal stability of cold spray tantalum coating irradiated with deuterium. Phys. Scr. (in review).Google Scholar
Ivanov, A.A. & Prikhodko, V.V. 2013 Gas-dynamic trap: an overview of the concept and experimental results. Plasma Phys. Control. Fusion 55 (6), 063001–32.CrossRefGoogle Scholar
Jaeger, E.F., Berry, L.A., D'Azevedo, E., Batchelor, D.B. & Carter, M.D. 2001 All-orders spectral calculation of radio-frequency heating in two-dimensional toroidal plasmas. Phys. Plasmas 8 (5), 15731583.CrossRefGoogle Scholar
Katsumata, R., Ichimura, M., Inutake, M., Hojo, H., Mase, A. & Tamano, T. 1996 Eigenmode excitation of Alfvén ion cyclotron instability. Phys. Plasmas 3, 44894495.CrossRefGoogle Scholar
Kelley, G.G. 1967 Elimination of ambipolar potential-enhanced loss in a magnetic trap. Plasma Phys. 9 (4), 503.CrossRefGoogle Scholar
Kesner, J. 1973 Inverse ambipolar potential in a magnetic mirror configuration. Plasma Phys. 15 (6), 577.CrossRefGoogle Scholar
Kesner, J. 1985 Axisymmetric, wall-stabilized tandem mirrors. Nucl. Fusion 25 (3), 275282.CrossRefGoogle Scholar
Killeen, J. & Marx, K.D. 1970 Solution of the Fokker-Planck equation for a mirror-confined plasma. Meth. Comput. Phys. 9, 421.Google Scholar
Killeen, J., Mirin, A.A. & Rensink, M.E. 1976 The solution of the kinetic equations for a multispecies plasma. Meth. Comput. Phys.: Adv. Res. Appl. 16, 389431.Google Scholar
Kotelnikov, I.A., Chernoshtanov, I.S. & Prikhodko, V.V. 2017 Electrostatic instabilities in a mirror trap revisited. Phys. Plasmas 24 (12), 122512.CrossRefGoogle Scholar
Kotelnikov, I.A., Ivanov, A.A., Yakovlev, D.V., Chen, Z. & Zeng, Q. 2019 Divertor for a steady-state gas-dynamic trap. Nucl. Fusion 60 (1), 016008.CrossRefGoogle Scholar
Kotelnikov, I., Zeng, Q., Prikhodko, V., Yakovlev, D., Zhang, K., Chen, Z. & Yu, J. 2022 Wall stabilization of the rigid ballooning $m = 1$ mode in a long-thin mirror trap. Nucl. Fusion 62 (9), 096025.CrossRefGoogle Scholar
Kremeyer, T., Flesch, K., Schmitz, O., Schlisio, G. & Wenzel, U. 2020 Wisconsin in situ penning (wisp) gauge: a versatile neutral pressure gauge to measure partial pressures in strong magnetic fields. Rev. Sci. Instrum. 91 (4), 043504.CrossRefGoogle ScholarPubMed
Lau, C., Berry, L.A., Jaeger, E.F. & Bertelli, N. 2019 Cold plasma finite element wave model for helicon waves. Plasma Phys. Control. Fusion 61 (4), 045008.CrossRefGoogle Scholar
Mirnov, V.V. & Ryutov, D.D. 1979 Linear gasdynamic system for plasma confinement. Sov. Tech. Phys. Lett. 5 (11), 678682.Google Scholar
Petrov, Y.V., Harvey, R.W., Forest, C.B. & Anderson, J.K. 2023 Calculations of WHAM2 mirror neutron rates and FI transport using the GENRAY/ CQL3D-M and MCGO-M codes. In Proceedings of 49th European Conference on Plasma Physics, July 3–7 2023, Bordeaux, France.Google Scholar
Petty, C.C., Baity, F.W., de Grassie, J.S., Mau, T.K., Pinsker, R.I., Porkolab, M. & Prater, R. 2001 Fast wave current drive at high ion cyclotron harmonics on DIII-D. Plasma Phys. Control. Fusion 43 (12), 17471758.CrossRefGoogle Scholar
Post, R.F. 1961 Equilibrium ambipolar potentials in a mirror machine. Phys. Fluids 4 (7), 902.CrossRefGoogle Scholar
Post, R.F. 1987 The magnetic mirror approach to fusion. Nucl. Fusion 27 (10), 15791739.CrossRefGoogle Scholar
Post, R.F. & Ryutov, D.D. 1995 Mirror fusion research: update and prospects. Comm. Plasma Phys. Control. Fusion 16 (6), 375399.Google Scholar
Prater, R. 1971 Feedback suppression of a large-growth-rate flute mode. Phys. Rev. Lett. 27 (3), 132135.CrossRefGoogle Scholar
Rosenberg, A.L., Menard, J.E., Wilson, J.R., Medley, S.S., Andre, R., Phillips, C.K., Darrow, D.S., LeBlanc, B.P., Redi, M.H., Fisch, N.J., Team, N.S.T.X., Harvey, R.W., Mau, T.K., Jaeger, E.F., Ryan, P.M., Swain, D.W., Sabbagh, S.A. & Egedal, J. 2004 Fast ion absorption of the high harmonic fast wave in the national spherical torus experiment. Phys. Plasmas 11 (5), 24412452.CrossRefGoogle Scholar
Rosenbluth, M.N. & Longmire, C.L. 1957 Stability of plasmas confined by magnetic fields. Ann. Phys. (N.Y.) 1 (2), 120140.CrossRefGoogle Scholar
Ryutov, D.D., Berk, H.L., Cohen, B.I., Molvik, A.W. & Simonen, T.C. 2011 Magneto- hydrodynamically stable axisymmetric mirrors. Phys. Plasmas 18 (9), 092301.CrossRefGoogle Scholar
Sakai, O., Yasaka, Y. & Itatani, R. 1993 High radial confinement mode induced by dc limiter biasing in the HIEI tandem mirror. Phys. Rev. Lett. 70, 40714074.CrossRefGoogle ScholarPubMed
Schmitz, L., et al. 2016 Suppressed ion-scale turbulence in a hot high-$\beta$ plasma. Nat. Commun. 7, 13860.CrossRefGoogle Scholar
Seemann, O., Be'ery, I. & Fisher, A. 2018 Stabilization of magnetic mirror machine using rotating magnetic field. J. Plasma Phys. 84 (5), 905840502.CrossRefGoogle Scholar
Severn, G., HershkowitE, N., Nelson, B. & Pew, J. 1983 Radial Potential Modification Experiments in the Phaedrus Tandem Mirror Machine. Bulletin of the American Physical Society.Google Scholar
Shidara, H., Nakashima, Y., Imai, T., Kariya, T., Minami, R., Ichimura, M. & Yoshikawa, M. 2011 Ec wave plasma breakdown on gamma 10 device. Fusion Engng Des. 86 (6–8), 913915.CrossRefGoogle Scholar
Simonen, T., et al. 2008 The axisymmetric tandem mirror: a magnetic mirror concept game changer magnet mirror status study group. LLNL Tech. Rep..CrossRefGoogle Scholar
Simonen, T.C., et al. 1983 Operation of the tandem-mirror plasma experiment with skew neutral-beam injection. Phys. Rev. Lett. 50, 16681671.CrossRefGoogle Scholar
Smirnov, A.P., Harvey, R.W. & Kupfer, K. 1994 A general ray tracing code genray. Bull. Am. Phys. Soc. 39 (7), 1626.Google Scholar
Smith, G.R. 1984 Alfvén ion-cyclotron instability in tandem-mirror plasmas. I. Phys. Fluids 27 (6), 14991513.CrossRefGoogle Scholar
Soldatkina, E.I., Maximov, V.V., Prikhodko, V.V., Savkin, V.Y.A., Skovorodin, D.I., Yakovlev, D.V. & Bagryansky, P.A. 2020 Measurements of axial energy loss from magnetic mirror trap. Nucl. Fusion 60 (8), 086009.CrossRefGoogle Scholar
Stix, T.H. 1992 Waves in Plasmas. American Institute of Physics.Google Scholar
Whyte, D. 2019 Small, modular and economically attractive fusion enabled by high temperature superconductors. Phil. Trans. A. Math. Phys. Engng Sci. 377 (2141), 20180354.Google ScholarPubMed
Whyte, D.G., Minervini, J., LaBombard, B., Marmar, E., Bromberg, L. & Greenwald, M. 2016 Smaller and sooner: exploiting high magnetic fields from new superconductors for a more attractive fusion energy development path. J. Fusion Energy 35, 4153.CrossRefGoogle Scholar
Wurzel, S.E. & Hsu, S.C. 2022 Progress toward fusion energy breakeven and gain as measured against the Lawson criterion. Phys. Plasmas 29 (6), 062103.CrossRefGoogle Scholar
Yakovlev, D.V., Shalashov, A.G., Gospodchikov, E.D., Maximov, V.V., Prikhodko, V.V., Savkin, V.Y.A., Soldatkina, E.I., Solomakhin, A.L. & Bagryansky, P.A. 2018 Stable confinement of high-electron-temperature plasmas in the GDT experiment. Nucl. Fusion 58 (9), 094001.CrossRefGoogle Scholar
Yakovlev, D.V., Shalashov, A.G., Gospodchikov, E.D., Solomakhin, A.L., Savkin, V.Y.A. & Bagryansky, P.A. 2017 Electron cyclotron plasma startup in the GDT experiment. Nucl. Fusion 57 (1), 016033.CrossRefGoogle Scholar
Yushmanov, P.P. 1966 Confinement of slow ions of plasma with positive potential in a mirror trap. Sov. Phys. JETP 22, 409.Google Scholar
Zhu, G., Shi, P., Yang, Z., Zheng, J., Luo, M., Ying, J. & Sun, X. 2019 A new method to suppress the Rayleigh-Taylor instability in a linear device. Phys. Plasmas 26 (4), 042107.CrossRefGoogle Scholar
Zweben, S.J., et al. 1997 Alpha-particle physics in the tokamak fusion test reactor DT experiment. Plasma Phys. Control. Fusion 39 (5A), A275.CrossRefGoogle Scholar
Figure 0

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.

Figure 1

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.

Figure 2

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 3

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.

Figure 4

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

Figure 5

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.

Figure 6

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.

Figure 7

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.

Figure 8

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.

Figure 9

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

Figure 10

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.

Figure 11

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.

Figure 12

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.

Figure 13

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.

Figure 14

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.

Figure 15

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.

Figure 16

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.

Figure 17

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.

Figure 18

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.

Figure 19

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.

Figure 20

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.

Figure 21

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 22

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