2. Field configuration

The classic MM system comprises one central fusion cell and two MM sections (Post Reference Post1967; Budker et al. Reference Budker, Mirnov and Ryutov1971; Logan et al. Reference Logan, Brown, Lieberman and Lichtenberg1972a , Reference Logan, Lichtenberg, Lieberman and Makhijanib ; Mirnov & Ryutov Reference Mirnov and Ryutov1972; Makhijani et al. Reference Makhijani, Lichtenberg, Lieberman and Logan1974; Tuszewski et al. Reference Tuszewski, Lichtenberg and Eylon1977; Mirnov & Lichtenberg Reference Mirnov and Lichtenberg1996; Kotelnikov Reference Kotelnikov2007; Burdakov & Postupaev Reference Burdakov and Postupaev2016). Following Post (Reference Post1967), we model the MM section by a periodic magnetic field with an axial component of the form

(2.1) \begin{eqnarray} B_{z} = B_{0}\left [1+\left (R_{m}-1\right )\exp \left (-5.5\sin ^{2}\frac {\pi z}{l}\right )\right ]\! ,\end{eqnarray}

where $B_0$ is the minimal magnetic field, $R_m = B_{max} / B_0$ is the mirror ratio and $l$ is the length of each MM cell. Figure 1 illustrates the MM system and the axial component of the magnetic field, $B_z$ . To satisfy Gauss’s law $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{B}=0$ , we include $x,y$ components in the MM magnetic field such that the total (static) field near the mirror axis reads

(2.2) \begin{eqnarray} \boldsymbol{B}_{\text{mirror}}=-\frac {1}{2}\frac {\partial B_{z}}{\partial z}(x\hat {x}+y\hat {y})+B_{z}\hat {z} .\end{eqnarray}

Figure 1.

An illustration of one section of an MM system (top) and the amplitude of the axial magnetic field (2.1) of two MM cells with $R_m=3$ (bottom). The illustration shows the right half of the MM system, where the left half is assumed to lie to the left of the vertical dashed line.

As discussed in the introduction, we consider time-dependent TRMF and examine the effect of phase-space mixing and the resulting enhanced confinement with minimal plasma heating. First, a simple RMF in the mirror transverse plane ( $x-y$ ) but at a fixed axial $(z)$ location, can be written in phasor form as

(2.3) \begin{eqnarray} \boldsymbol{B}_{\mathrm{RMF}}&=B_{1}e^{-i\omega t} (\hat {x}-i\!\hat{y}), \end{eqnarray}

where $B_{1}$ is the field amplitude, $\omega$ is the RF angular frequency of the rotating field and the rotation direction is the same as that of the cyclotron motion of positively charged ions for $B_0\gt 0$ . The corresponding induced electric field, in this case, can be approximated to first order in $\omega r/c$ , with $r$ the distance from the mirror axis (in cgs convention) as

(2.4) \begin{eqnarray} \boldsymbol{E}_{\mathrm{RMF}}&=-B_{1}\omega e^{-i\omega t} (x-iy)\hat {z}. \end{eqnarray}

It is noted that the induced electrostatic field vanishes on the mirror axis and grows linearly with $r$ .

To extend the magnetic field to a travelling field along the axial $\hat {z}$ direction (i.e. a TRMF), we define the field’s phase as

(2.5) \begin{eqnarray} \varphi \left (z, t \right )=k z-\omega t. \end{eqnarray}

In this case, both the magnetic and the induced electric fields should include spatial corrections, which, to the first order, read

(2.6) \begin{align} \boldsymbol{B}_{\mathrm{TRMF}}& = B_{1}e^{i\varphi }\,[(\hat {x}-i\!\hat{y})+k(y+ix)\hat {z}] ,\\[-10pt]\nonumber\end{align}

(2.7) \begin{align} \boldsymbol{E}_{\mathrm{TRMF}}& = B_{1}\omega e^{i\varphi }\,[-(x-iy)\hat {z} + kxy(\hat {x}+i\!\hat{y})] ,\end{align}

where the two small parameters in the approximation are $\omega r/c \ll 1$ and $kr \ll 1$ (see the Appendix). This expression for the induced electric field assumes vacuum conditions, whereas in plasmas, the electric field may be suppressed or penetrate less effectively than the magnetic field. Thus, we also analyse the case in which it is absent, realised by setting it to zero in part of the single-particle calculations in § 3. Notably, in this scenario, the TRMF does not transfer energy to the particle and therefore does not contribute to plasma heating. Hereafter, we refer to this scheme as TRMF–noE.

For completeness and comparison, we will also consider and simulate the TREF plugging method that was first proposed in Miller et al. (Reference Miller, Be’ery, Gudinetsky and Barth2023). The expressions (in cgs convention) for the electric and the (induced) magnetic fields in the TREF method, which are analogous to those of the TRMF method, but with a role exchange of the magnetic and electric fields, read

(2.8) \begin{align} \boldsymbol{E}_{\mathrm{TREF}}&=E_1e^{i\varphi }\,[(\hat {x}-i\!\hat{y})+k\,(y+ix)\hat {z}] ,\\[-10pt]\nonumber\end{align}

(2.9) \begin{align} \boldsymbol{B}_{\mathrm{TREF}}&=E_1\frac {\omega }{c^2} e^{i\varphi }\,[(x-iy)\hat {z} - kxy\,(\hat {x}+i\!\hat{y})], \end{align}

where $E_1$ is the TREF electric-field amplitude. Here, we also include the first-order temporal ( $\omega r/c$ ) and spatial ( $kr$ ) correction terms for both electric and induced magnetic fields, which were neglected in Miller et al. (Reference Miller, Be’ery, Gudinetsky and Barth2023). Nonetheless, the forthcoming simulation results show that this assumption is justified for TREF.

In the next section, we consider three scenarios: TRMF with and without an induced electric field, and TREF with an induced magnetic field, although for the parameters considered, its effect on the dynamics is negligible. The two TRMF scenarios, i.e. with and without an induced electric field, can serve as models for the limiting behaviour of dilute laboratory plasmas (where the electric field can penetrate readily) and for dense fusion plasmas (where the electric field is strongly suppressed). A crucial caveat is that, in the high-density plasma limit, the magnetic field itself might be suppressed or modified by the plasma, thereby affecting the efficiency of TRMF plugging. A detailed study of this collective effect is left for future work.

3. Single-particle simulations

To quantify the RF effect on particles in the mirror system, we employ the single-particle approximation and perform a Monte Carlo analysis, where the initial velocities were sampled from a thermal distribution with an ion temperature of $k_{B} T_i=10\,\mathrm{keV}$ and a uniformly distributed random direction. The ions were initialised at the mirror midplane ( $z = l/2$ ), with radial positions, $r_0$ , sampled uniformly from the range $[0, 10]\,\mathrm{cm}$ to represent a realistic plasma profile. The azimuthal angle (in the $x$ – $y$ plane) was also sampled uniformly over the interval $[0, 2\pi ]$ to simulate a random phase relative to the phase of the RF field. The static mirroring magnetic field modelled one MM cell with axial magnetic field as defined in (2.1), where $l=1\,\mathrm{m}$ , $B_0=1\,\mathrm{T}$ and $R_m=5$ . For the time-dependent RF fields, we considered 3 cases of RF fields, as described in § 2: (a) TRMF ((2.6) and (2.7)), (b) TRMF–noE ((2.6) solely) and (c) TREF ((2.8) and (2.9)).

For TRMF, i.e. cases (a) and (b), the magnetic field amplitude was $B_{RF}=0.05\,\text{T}$ while for TREF of case (c), the electric-field amplitude was $E_{RF}=50\,\text{kV m}^{-1}$ . In case (c), we found that the effect of the induced magnetic field (2.9) is negligible, so we do not show a case of TREF without the induced magnetic field. Neglecting collisions and collective effects, we calculate the trajectories of deuterium and tritium (D-T) ions under the influence of the Lorentz force $\boldsymbol{F}=q(\boldsymbol{E}+\boldsymbol{v}\times \boldsymbol{B})$ . We employed a symplectic scheme, which preserves phase-space volume (He et al. Reference He, Sun, Liu and Qin2015), to solve the particles’ trajectories under the influence of the static and RF fields. The simulation time for all particles of the same isotope was identical, chosen to be the thermal passage time of a particle along one MM cell length, $\tau _{\textit{th}} \equiv l/v_{\textit{th}}$ , where $v_{\textit{th}}$ is the ion’s thermal velocity (different for D-T ions).

The ion cyclotron frequencies of D-T at the mirror midplane are $\omega _{0,D} = eB_0/m_D = 2\pi \times 7.6$ MHz and $\omega _{0,T} = eB_0/m_T = 2\pi \times 5$ MHz, respectively. Note that in the figures, the RF frequency $\omega$ is usually normalised by $\omega _{0,T}$ . The time step in the simulations was set to be $\Delta t=\tau _{cyc}/50$ where $\tau _{cyc}=2\pi / \omega _{0D,0T}$ is the cyclotron period for the relevant isotope of each calculation (tested for numerical convergence). Also note that although the calculations were performed separately for D-T, both present in a D-T fusion reactor but exhibiting distinct dynamics due to their different masses, we show only the tritium results in some of the figures to avoid clutter. The key results related to the plugging enhancement, however, are presented for both isotopes.

For each RF scheme, we survey a range of RF parameters, $k$ and $\omega$ , and for each configuration, we perform $10\,000$ Monte Carlo realisations of the initial conditions to obtain reliable statistics on the kinetic transition rates. The range selected for the survey in the $k$ – $\omega$ parameter space was chosen to remain consistent with the validity limits of the first-order approximation in (2.6)–(2.7), i.e. $\omega r / c \ll 1$ and $k r \ll 1$ , and with the system’s physical dimensions. Consequently, the RF frequency $\omega$ was chosen in the range $\omega \in [0.5, 2]\,\omega _{0,T}$ , roughly centred around the tritium cyclotron frequency, while the wave-vector range was $k \in 2\pi \times [-1,\, 1]\,\mathrm{m}^{-1}$ , corresponding to wavelengths bounded by the MM cell length. The resulting restrictions on the particle position are $r\ll c/\omega _{max}=5$ m and $r \ll 1/k_{max}=16$ cm. These conditions are satisfied by the particles’ initial radii in our simulations, which were set to be at most $10$ cm.

In the semi-kinetic rate-equation model (see § 4 below), we classify the particles into three populations according to their initial conditions: (a) confined particles, (b) right-going particles and (c) left-going particles. In our convention (see figure 1), right-going particles escape outward from the MM system, while left-going particles propagate inward toward the central fusion cell (the roles would be reversed for the left half of the MM system). This classification is performed at the beginning of the simulation, by checking the loss-cone condition at the midplane $( v_{\perp }/v )^2 \lt B_{\textit{min}} / B_{\textit{max}}$ , where $v_{\perp }$ is the perpendicular velocity and $v$ is the total velocity.

We denote the number of particles in each population by $N_{c}$ for confined particles, $N_{r}$ for right-going (escaping), and $N_{l}$ for left-going (incoming) particles. During the simulations, we take $50$ equidistant time snapshots and check, for each particle, the local loss-cone condition $( v_{\perp }/v )^2 \lt B / B_{\textit{max}}$ , where $B$ is the local mirror magnetic field at the particle’s location at a given time. We then track the time evolution of each particle’s classification, which is updated whenever the particle crosses one of the loss-cone boundaries. The critical parameter required to estimate the overall efficiency of RF plugging in MM machines via the rate-equation model (see § 4) is the number of ‘converted’ particles, say, that originated at the right loss cone but ended up confined $\Delta N_{r \rightarrow c}$ . Let us define the time-dependent population-conversion metric

(3.1) \begin{align} \Delta \bar {N}_{\textit{ij}} \left ( t \right ) \equiv \frac {\Delta N_{i \rightarrow j}(t)}{N_i (t=0)}, \qquad ij \in \{rc, cr, lc, cl, rl, lr\}, \end{align}

for the six possible transitions. These quantities are normalised by the original number of particles in the source population, so by definition $\Delta \bar {N}_{\textit{ij}} ( t )\in [0,\, 1]$ .

In figures 2, 3 and 4, we plot the time evolution of the population-conversion metric for tritium for different values of the RF-parameters $k$ and $\omega$ , and for the different RF schemes, TRMF, TRMF–noE and TREF, respectively. Each subplot in the figures corresponds to a specific pair of $k$ and $\omega$ values, where the horizontal axis shows the time over the interval $[0,\, \tau _{\textit{th}}]$ , and the vertical axis displays the population-conversion metric over the range $[0,\, 1]$ . The corresponding plots for deuterium are similar to those for tritium, so we omit them here. The fluctuations observed in some of the trajectories arise from particles switching identities multiple times, oscillating between different populations. To estimate the statistical error associated with the finite number of Monte Carlo particle samples, we apply a bootstrap procedure in which $50$ resampled subsets (with repetitions) are drawn from the original set of $10,000$ particles, and the population-conversion metric is recomputed for each subset. In the population-conversion plots, the solid lines denote the median, while the shaded regions indicate the $\pm 2 \sigma$ bounds at each time point.

Figure 2.

Population-conversion $\Delta \bar {N}_{\textit{ij}}$ plots for tritium in TRMF for different values of $k,\omega$ . Colours indicate different transitions between the three populations (see legend). In each subplot, the horizontal axis represents time over the interval $[0,\, \tau _{\textit{th}}]$ , while the vertical axis shows the population-conversion metric over the range $[0,\, 1]$ .

Figure 3.

Population-conversion $\Delta \bar {N}_{\textit{ij}}$ plots for tritium in TRMF–noE for different values of $k,\omega$ . Colours indicate different transitions between the three populations (see legend). In each subplot, the horizontal axis represents time over the interval $[0,\, \tau _{\textit{th}}]$ , while the vertical axis shows the population-conversion metric over the range $[0,\, 1]$ .

Figure 3. Long description

The image contains multiple line graphs arranged in a grid. Each subplot represents population-conversion metrics for tritium in TRMFnoE for different values. The horizontal axis in each subplot represents time over the interval, while the vertical axis shows the population-conversion metric over the range. The graphs are color-coded to indicate different transitions between the three populations, as shown in the legend. The colors used are blue, green, red, and yellow, each representing different transitions. The trends in the graphs show how the population-conversion metric varies over time for different transitions and values.

Figure 4.

Population-conversion $\Delta \bar {N}_{\textit{ij}}$ plots for tritium in TREF for different values of $k,\omega$ . Colours indicate different transitions between the three populations (see legend). In each subplot, the horizontal axis represents time over the interval $[0,\, \tau _{\textit{th}}]$ , while the vertical axis shows the population-conversion metric over the range $[0,\, 1]$ .

For $k=0$ , the RF fields are symmetric in the $\pm \hat {z}$ directions, so particles travelling in opposing axial directions are expected to behave similarly. This expectation is consistent, within the numerical accuracy indicated by the shaded regions, with the results presented in figures 2, 3 and 4). For example, this symmetry is evident in the comparison of $\Delta \bar {N}_{rc}$ (blue) and $\Delta \bar {N}_{lc}$ (green) in the middle column of the figures. Due to symmetry, for a fixed $\omega$ , the solutions corresponding to opposing wave vectors $\pm k$ are expected to be antisymmetric under the transformation $\hat {z} \rightarrow -\hat {z}$ . This behaviour is clearly visible in the curves in the panels with positive $k$ when compared with the corresponding panels with negative $k$ .

The next step is to evaluate from the single-particle simulations the transition rates between the different populations, which will serve as the transition coefficients in the rate-equation model to be introduced in § 4. To this end, we extract the population-conversion values at the end of each simulation $\bar {N}_{\textit{ij}} \equiv \Delta \bar {N}_{\textit{ij}} ( \tau _{\textit{th}} )$ , defined as the average over the final five time frames. To get RF-induced transition rates in units of inverse time, we divide by $\tau _{\textit{th}}$ , the characteristic time scale for an average non-trapped particle to travel a single mirror cell

(3.2) \begin{equation} \nu _{\textit{RF},\textit{ij}} \equiv \frac {\bar {N}_{\textit{ij}}}{\tau _{\textit{th}}} = \frac { \Delta \bar {N}_{\textit{ij}} \left ( \tau _{\textit{th}} \right ) }{\tau _{\textit{th}}}. \end{equation}

For each RF scheme, we plot the RF-induced transition rates (six panels per figure) for tritium, obtained from the single-particle Monte Carlo simulations. The schemes presented are TRMF (figure 5), TRMF–noE (figure 6) and TREF (figure 7). For visualisation purposes, we plot the dimensionless RF rates $ \bar {N}_{\textit{ij}}= \nu _{RF,ij} \tau _{\textit{th}}$ after applying a Gaussian filter to smooth the RF rate maps and reduce numerical noise. For comparison, we also conducted a full set of Monte Carlo simulations to obtain the transition rates for deuterium, the second component of the D-T fuel. The results are similar to those of tritium, and for illustration, we present them in figure 8 for the case of TRMF–noE.

Figure 5.

Smoothed and dimensionless RF rates, $\bar {N}_{\textit{ij}}$ , as a function of $k,\omega$ for tritium in TRMF (with induced electric field). The overlaid dashed black lines indicate the theoretical resonance condition for right- and left-going particles.

Figure 5. Long description

Six heat maps depict the dimensionless RF rates for tritium in TRMF with induced electric field. Each heat map shows the rates as a function of k/(2πm−¹) on the horizontal axis and ω/ω_b,τ on the vertical axis. The color scale on the right indicates the magnitude of the rates, with red representing higher values and blue representing lower values. The overlaid dashed black lines indicate the theoretical resonance condition for right- and left-going particles. Panel A: N_rc (T) shows a range from 0.2 to 0.8. Panel B: N_cr (T) shows a range from 0.00 to 0.08. Panel C: N_rl (T) shows a range from 0.00 to 0.12. Panel D: N_lc (T) shows a range from 0.2 to 0.8. Panel E: N_cl (T) shows a range from 0.00 to 0.08. Panel F: N_lr (T) shows a range from 0.00 to 0.12. Each panel displays distinct patterns and trends in the dimensionless RF rates.

Figure 6.

Smoothed and dimensionless RF rates, $\bar {N}_{\textit{ij}}$ , as a function of $k,\omega$ for tritium in TRMF–noE. The overlaid dashed black lines indicate the theoretical resonance condition for right- and left-going particles.

Figure 7.

Smoothed and dimensionless RF rates, $\bar {N}_{\textit{ij}}$ , as a function of $k,\omega$ for tritium in TREF. The overlaid dashed black lines indicate the theoretical resonance condition for right- and left-going particles.

Figure 7. Long description

Six heat maps depict the dimensionless RF rates for tritium in TREF as a function of wave number (k) and frequency ratio (ω/ωb,T). Each heat map has a color scale indicating the magnitude of the RF rates. The x-axis represents the wave number (k) in units of 2πnm^-1, ranging from -1.0 to 1.0. The y-axis represents the frequency ratio (ω/ωb,T), ranging from 0.6 to 1.8. Overlaid dashed black lines indicate the theoretical resonance condition for right- and left-going particles. Panel A (Nrc(T)) shows a heat map with RF rates ranging from 0.2 to 0.8. Panel B (Ncr(T)) displays RF rates from 0.00 to 0.08. Panel C (Nrl(T)) shows RF rates from 0.00 to 0.12. Panel D (Nlc(T)) depicts RF rates from 0.2 to 0.8. Panel E (Ncl(T)) shows RF rates from 0.00 to 0.08. Panel F (Nlr(T)) displays RF rates from 0.00 to 0.12.

Figure 8.

Smoothed and dimensionless RF rates, $\bar {N}_{\textit{ij}}$ , as a function of $k,\omega$ for deuterium in TRMF–noE. The overlaid dashed black lines indicate the theoretical resonance condition for right- and left-going particles.

The RF-induced rates for the four dominant transitions, $\bar {N}_{rc},\bar {N}_{cr}, \bar {N}_{lc}, \bar {N}_{cl}$ , are found to be concentrated along straight lines in the ( $k,\omega$ ) parameter space. These lines are associated with the cyclotron resonance condition, which results from the compensation between the Doppler frequency detuning (in the moving frame of each population) and the RF wave phase velocity, as discussed in § 2. The overlaid dashed black lines in the figures indicate the resonance condition, with $v_z$ replaced by the mean axial velocity within the loss cones

(3.3) \begin{eqnarray} \bar {v}_{z,LC}=\frac {\intop _{LC}\,f_{\textit{MB}}(\boldsymbol{v})\,v_{z} \,\text{d}^{3}v}{\intop _{LC}f_{\textit{MB}}(\boldsymbol{v})\,\text{d}^{3}v}. \end{eqnarray}

Here, $f_{\textit{MB}}(\boldsymbol{v})=\pi ^{-3/2}\, v_{th}^{-3} \exp (-\boldsymbol{v}^2/v_{th}^2)$ is the Maxwell–Boltzmann distribution function, where the integral is taken only over the loss-cone region of the velocity space. One finds that

(3.4) \begin{eqnarray} \bar {v}_{z,LC}=\pi ^{-1/2} \Big( 1 +\sqrt {1-R_m^{-1}} \Big) v_{th}, \end{eqnarray}

where for the considered mirror ratio $R_m=5$ the solution is $\bar {v}_{z,LC}=1.07 \,v_{th}$ . Notably, the right and left loss cones have mean axial velocities in opposite directions, resulting in resonance lines with opposite slopes.

It can also be seen that the theoretical lines appear at slightly lower frequencies compared with the peak rate values. This can be explained by the increase in magnetic field experienced by ions as they move away from the mirror midplane, leading to higher cyclotron frequencies. For deuterium, the transition-rate maps in figure 8 are shifted to higher frequencies relative to tritium, reflecting the inverse dependence of the cyclotron frequency on ion mass.

4. Rate-equation model

In previous work, we generalised the semi-kinetic rate-equation model for the MM system (Miller et al. Reference Miller, Be’ery and Barth2021), which included only particle transmission between neighbouring cells through the loss cones and Coulomb scattering within each cell, and added additional terms that represent driven transport induced by external RF fields (Miller et al. Reference Miller, Be’ery, Gudinetsky and Barth2023). Here, we adapt the rate-equation model and further generalise it to include the transitions between left- and right-going populations, which were neglected in (Miller et al. Reference Miller, Be’ery, Gudinetsky and Barth2023). The resulting rate model, therefore, includes all six RF-induced transitions $\nu _{RF,ij}$ as defined in (3.2) and estimated from the single-particle Monte Carlo simulations.

The generalised model then reads

(4.1) \begin{align} \dot n_{c}^i &= \nu _{s} \big[(1-2 \alpha ) \big(n_{l}^i + n_{r}^i\big) - 2 \alpha n_{c}^i\big]\nonumber \\ & \quad - ( \nu _{RF,cl} + \nu _{RF,cr} ) n_{c}^i + \nu _{RF,lc} n_{l}^i + \nu _{RF,rc} n_{r}^i , \end{align}

(4.2) \begin{align} \dot n_{l}^i & = \nu _{s}\left [\alpha \big(n_{r}^i+n_{c}^i\big) - (1-\alpha ) n_{l}^i\right ] - \nu _{t} n_{l}^i + \nu _{t} n_{l}^{i+1} \nonumber\\ & \quad - ( \nu _{RF,lc} + \nu _{RF,lr} ) n_{l}^i + \nu _{RF,cl} n_{c}^i + \nu _{RF,rl} n_{r}^i , \end{align}

(4.3) \begin{align} \dot n_{r}^i & = \nu _{s} \left [\alpha \big(n_{l}^i +n_{c}^i\big) - (1-\alpha ) n_{r}^i\right ]- \nu _{t} n_{r}^i + \nu _{t} n_{r}^{i-1}\nonumber \\& \quad - ( \nu _{RF,rc} + \nu _{RF,rl} ) n_{r}^i + \nu _{RF,cr} n_{c}^i + \nu _{RF,lr} n_{l}^i , \; \end{align}

where $\nu _s$ is the ion–ion Coulomb scattering rate that roughly scales with density and temperature as $\propto n/T^{3/2}$ ; $\nu _t = f_t v_{th}/l$ the inter-cell transmission rate, where $v_{th}$ is the thermal velocity of the ions and $l$ is the length of the mirror cell and the ambipolar factor is $f_t=( T_i + T_e )/T_i=2$ ; and finally $\alpha$ is the normalised loss-cone solid angle. For the mirror ratio considered $R_m=5$ the loss-cone angle is $\theta _{LC} = \arcsin ( R_{m}^{-1/2} ) \approx 26.5 ^{\circ }$ giving the normalised loss-cone solid angle $\alpha = \sin ^{2} (\theta _{LC}/2) \approx 0.053$ ; Therefore, the relative size of each loss cone is $\alpha$ and the confined section $1-2\alpha$ . Among the different thermodynamic scenarios studied in Miller et al. (Reference Miller, Be’ery and Barth2021), we considered here the simple isothermal scenario, in which $\nu _s,\nu _t$ are constant across all cells, since the external RF fields dominate the scattering rather than collisions.

We close the rate-equation model by imposing the following boundary conditions: a constant total density in the throat of the central fusion cell (the leftmost cell)

(4.4) \begin{align} n_c^1+n_{l}^1+n_{r}^1 =n_0 =\mathrm{const}, \end{align}

where $n_0$ is the density in the central cell, and a free-flow boundary condition at the exit of the outermost cell

(4.5) \begin{align} \nu _t n_l^{N+1}=0, \end{align}

corresponding to the absence of left-going flux from outside the system.

We solve the rate equations ((4.1)–(4.3)) for the steady-state solution, i.e. $\dot {n}=0$ for all cell populations in all cells. The outgoing flux between two neighbouring cells (e.g. $i\rightarrow i+1$ ) is proportional to $\phi _{i,i+1} \propto v_{th} ( n_{r}^i - n_{l}^{i+1} )$ . In a steady state, by definition, all inter-cell fluxes are the same and denoted by $\phi _{ss}$ . The system confinement time scales inversely with the steady-state flux, $\tau \propto 1/\phi _{ss}$ , which we aim to maximise to satisfy the Lawson criterion (Lawson Reference Lawson1957; Wurzel & Hsu Reference Wurzel and Hsu2022). It is noted that in the absence of RF fields, the confinement time scales linearly with the number of mirror cells. However, even under optimised and optimistic conditions, fusion-relevant confinement would require an impractically large number of cells (Logan et al. Reference Logan, Lichtenberg, Lieberman and Makhijani1972b ; Makhijani et al. Reference Makhijani, Lichtenberg, Lieberman and Logan1974; Miller et al. Reference Miller, Be’ery and Barth2021).

Figure 9.

Steady-state flux in an MM system with $N=50$ cells as a function of $k$ and $\omega$ for deuterium (left) and tritium (right). The RF schemes are TRMF (top), TRMF–noE (centre) and TREF (bottom). The overlaid dashed black lines indicate the theoretical resonance condition for right- and left-going particles.

Figure 9. Long description

Panel A: A heat map showing the steady-state flux in an MM system with cells for deuterium under the TRMF scheme. The x-axis represents k/(2πm−¹) ranging from -1.00 to 1.00, and the y-axis represents ω/ω_{0 r} ^T ranging from 0.6 to 1.8. The color scale indicates the log10(φ_{s s}/φ₀) values, with red representing higher values and blue representing lower values. The dashed black lines indicate the theoretical resonance condition for right- and left-going particles. Panel B: A heat map showing the steady-state flux in an MM system with cells for tritium under the TRMF scheme. The x-axis represents k/(2πm−¹) ranging from -1.00 to 1.00, and the y-axis represents ω/ω_{0 r} ^T ranging from 0.6 to 1.8. The color scale indicates the log10(φ_{s s}/φ₀) values, with red representing higher values and blue representing lower values. The dashed black lines indicate the theoretical resonance condition for right- and left-going particles. Panel C: A heat map showing the steady-state flux in an MM system with cells for deuterium under the TRMFnoE scheme. The x-axis represents k/(2πm−¹) ranging from -1.00 to 1.00, and the y-axis represents ω/ω_{0 r} ^T ranging from 0.6 to 1.8. The color scale indicates the log10(φ_{s s}/φ₀) values, with red representing higher values and blue representing lower values. The dashed black lines indicate the theoretical resonance condition for right- and left-going particles. Panel D: A heat map showing the steady-state flux in an MM system with cells for tritium under the TRMFnoE scheme. The x-axis represents k/(2πm−¹) ranging from -1.00 to 1.00, and the y-axis represents ω/ω_{0 r} ^T ranging from 0.6 to 1.8. The color scale indicates the log10(φ_{s s}/φ₀) values, with red representing higher values and blue representing lower values. The dashed black lines indicate the theoretical resonance condition for right- and left-going particles. Panel E: A heat map showing the steady-state flux in an MM system with cells for deuterium under the TREF scheme. The x-axis represents k/(2πm−¹) ranging from -1.00 to 1.00, and the y-axis represents ω/ω_{0 r} ^T ranging from 0.6 to 1.8. The color scale indicates the log10(φ_{s s}/φ₀) values, with red representing higher values and blue representing lower values. The dashed black lines indicate the theoretical resonance condition for right- and left-going particles. Panel F: A heat map showing the steady-state flux in an MM system with cells for tritium under the TREF scheme. The x-axis represents k/(2πm−¹) ranging from -1.00 to 1.00, and the y-axis represents ω/ω_{0 r} ^T ranging from 0.6 to 1.8. The color scale indicates the log10(φ_{s s}/φ₀) values, with red representing higher values and blue representing lower values. The dashed black lines indicate the theoretical resonance condition for right- and left-going particles.

It is noted that because our rate-equation model describes a single ion species, it cannot directly capture the dynamics of a deuterium–tritium plasma. However, in regimes where Coulomb scattering rates are negligible compared with the RF-induced transition rates, the two ion species may be treated as non-interacting. Under this assumption, the rate equations are solved separately for D-T, using RF transition rates obtained independently for each species from the single-particle simulations presented in § 3. Figure 9 shows the results for D-T under different RF schemes and for various values of $k$ and $\omega$ at a fixed number of cells, $N=50$ . The theoretical resonance lines, as defined in § 3, are indicated by dashed lines. The fluxes shown in the figures are normalised by the single-mirror flux, $\phi _0 = n_0 v_{th}$ , which provides a convenient reference for illustrating the flux reduction achieved by the various plugging schemes considered in the simulations. This normalisation differs from that used in Miller et al. (Reference Miller, Be’ery and Barth2021, Reference Miller, Be’ery, Gudinetsky and Barth2023), where the flux was normalised by the maximum value compatible with satisfying the Lawson criterion (Lawson Reference Lawson1957; Wurzel & Hsu Reference Wurzel and Hsu2022). Here, however, we use a simpler normalisation to emphasise the primary objective of flux reduction relative to the single-mirror value.

From the figure, it can be seen that the confinement improvement, quantified by the flux reduction, in all three RF scenarios over the parameter range explored in the simulations lies between one and approximately $2.5$ orders of magnitude. Along the resonant line corresponding to the right-going (outgoing) particles, both the TRMF and TREF configurations exhibit the strongest confinement enhancement, approaching a reduction of about $2.5$ orders of magnitude. In contrast, for the TRMF–noE configuration, the optimal parameter region lies along the opposite resonant line, corresponding to resonance with right-going (incoming) particles. This non-intuitive result is discussed in detail in the following section.

Figure 10.

Steady-state flux for tritium as a function of the number of MM cells for three different parameter sets (see legend). Line styles denote TREF (solid), TRMF (dashed) and TRMF–noE (dotted).

In addition, outside the resonant lines, both the TRMF and TREF schemes provide stronger confinement enhancement than the TRMF–noE configuration, resulting in greater suppression of the steady-state flux, as discussed in the following section. As a rough estimate, the axial flux must be reduced by at least four orders of magnitude relative to $\phi _0$ for an MM system to become viable as a fusion device (Miller et al. Reference Miller, Be’ery, Gudinetsky and Barth2023). Consequently, while the present technique provides substantial confinement enhancement, additional plugging mechanisms will be required for fusion-relevant applications, or alternatively, operation in a different parameter regime or the adoption of a modified system design.

Figure 10 examines the scaling of the steady-state flux with the number of MM cells, $N$ , for selected RF configurations, as indicated in the legend. In the resonant cases, the flux decays approximately exponentially with $N$ . However, the decay rate associated with resonance of right-going (outgoing) particles, corresponding to the TREF and TRMF schemes (red curves), is significantly larger than that associated with resonance of left-going particles in the TRMF–noE configuration (green curve). This indicates that, even under optimal resonant conditions for each scheme, the TRMF–noE configuration is less effective at suppressing the axial flux than the other two RF schemes. Nonetheless, despite its weaker confinement performance, the TRMF–noE configuration offers a critical practical advantage, i.e. dramatically lower power consumption, as discussed in the next section.