RF plugging of multi-mirror machines

One of the main challenges of fusion reactors based on magnetic mirrors is the axial particle loss through the loss cones. In multi-mirror (MM) systems, the particle loss is addressed by adding mirror cells on each end of the central fusion cell. Coulomb collisions in the MM sections serve as the retrapping mechanism for the escaping particles. Unfortunately, the conﬁnement time in this system only scales linearly with the number of cells in the MM sections and requires an unreasonably large number of cells to satisfy the Lawson criterion. Here, it is suggested to reduce the outﬂow by applying a traveling RF electric ﬁeld that mainly targets the particles in the outgoing loss cone. The Doppler shift compensates for the detuning of the RF frequency from the ion cyclotron resonance mainly for the escaping particles resulting in a selectivity effect. The transition rates between the different phase space populations are quantiﬁed via single-particle calculations and then incorporated into a semi-kinetic rate equations model for the MM system, including the RF effect. It is found that for optimized parameters, the conﬁnement time can scale exponentially with the number of MM cells, orders of magnitude better than a similar MM system of the same length but without the RF plugging, and can satisfy the Lawson criterion for a reasonable system size.


I. INTRODUCTION
Axial confinement is one of the main challenges for sustainable fusion in linear magnetically confined systems based on the mirroring effect.These open trap systems are advantageous for their simplicity, high-β (the ratio of plasma to magnetic pressures), and continuous mode of operations but suffer from radial instabilities and axial losses.Radial instabilities can be mitigated via passive, [1][2][3] radio frequency (RF), [3][4][5][6] or active 3,7,8 methods, while solutions for the outgoing flux through the loss-cone commonly involve geometrical modifications of the mirroring magnetic field.The variety of design concepts to reduce the loss-cone flux include tandem plugs with thermal barriers, [9][10][11][12][13][14][15][16][17][18] , diamagnetic confinement, 19,20 multi-mirrors (MM) systems, [21][22][23][24][25][26][27][28][29][30] moving multiple mirrors, 31,32 , helical mirror with rotating plasma, [33][34][35][36][37] , ponderomotive RF plugging, [38][39][40][41][42][43][44][45][46][47] and plugging using field reversed configuration (FRC) at the mirror throats. 48M systems are characterized by the two sections of multiple magnetic mirrors attached to each side of the central mirror cell, where fusion occurs.As a particle escapes the central cell into the right or left MM sections, it can be collisionally scattered out of the loss cone of one of the MM cells and later be scattered again back into the central cell.This process can be modeled as a one-dimensional diffusion dynamics, where we expect the axial flux to scale inversely with the number of cells in the MM section.The collisionality, which depends on the temperature and density profiles along the system, determines the number of cells needed for reducing the loss-cone flux to a desired level.Moreover, the collisionality in this scheme is required to be high enough i.e., Coulomb mean free path of the order of the mirror cell length, resulting in a requirement for lower temperatures and higher densities, which a) Electronic mail: ido.barth@mail.huji.ac.il are sub-optimal for fusion. 21,30,49M machines are commonly considered isothermal systems due to fast electron thermalization across the system, 24 resulting in a linear scaling of the confinement time with the system length.In a recent study, we developed a semi-kinetic rate equations model for MM systems that divides the ions in each mirror cell into three populations (trapped, escaping, and returning) and includes the processes of Coulomb scattering within each cell and the transmission between neighboring cells.The steady-state solution of the rate equations yields the density profile within the MM section and the outgoing axial flux.It was found that the scaling of the outgoing flux with the system's size depends on the thermodynamic scenario.The best confinement was obtained for isentropic systems, where the plasma adiabatically cools down with the density decreases within the MM section, and the mean free path drops along the MM section.][51] Therefore, new confinement enhancement methods should be developed in the quest for sustainable fusion in MM machines.
In this work, we suggest a new RF plugging method for MM systems.The idea is to apply an external RF electric field resonantly coupled with the ion cyclotron frequency in the moving frame of the outgoing particles.To this end, we employ a radially rotating electric field with a frequency slightly detuned from the exact ions' Larmor frequency and with a non-zero axial wave vector such that the Doppler shift compensates for the resonance mismatch only for outgoing particles.As a result, the RF field resonantly increases the perpendicular energy of the escaping particles only so they can be recaptured in the magnetic mirror cell.In contrast, it mildly affects the returning particles with the opposite axial velocities.This selection effect yields a significant confinement enhancement effect, which we study via single particle simulations and the semi-kinetic rate equations model.The structure of the paper is as follows.Sec.II introduces the considered static magnetic mirror field and the external RF field.Sec.III studies the effect of the RF fields on the time evolution of particles in phase space via single particle simulations.Sec.IV integrates the simulation results over the distribution of the fuel particles and evaluates the transition rates between different populations in phase space.Sec.V incorporates the effect of the RF fields into an extended rate equations model for MM systems and calculates the confinement enhancement resulting from the RF fields.Sec.VI summarizes the findings.

II. FIELDS CONFIGURATIONS
2][23][24][25][26][27][28][29][30] To study the confinement improvement of adding the MM sections, we consider only one side of the MM sections, say the right one, and model it by a periodic magnetic field with an axial component of the form 21 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.The other magnetic field components in cylindrical coordinates, (r, θ , z), are B θ = 0 and B r = −r ∂ ∂ z B z , to satisfy ∇ • B = 0. Fig. 1 illustrates the MM system and the axial component of the magnetic field, B z .The magnetic mirror trapping criterion with axial (v z ) and transverse (v ⊥ ) components that are measured in the minimum of B z , i.e., at the center of the magnetic mirror cell (z = l/2).Therefore, manipulating the velocities ratio, v ⊥ /v, for example, by RF fields, is an appealing path to enhance trapping, especially in MM systems where the effect on the pitch angle can be accumulated along the MM section.
The simplest approach is to apply an RF field near the ion cyclotron resonance on each MM.This approach acts symmetrically on both right-and left-going particles.However, in MM machines, an asymmetric drive would be preferable since one direction is inward (toward the fusion cell) while the other is outward (outside the system).Moreover, unlike in the central (fusion) cell, plasma heating in the MM sections does not contribute to the fusion rate.It may be even deleterious for confinement as it reduces the collision rate, which is the salient re-trapping mechanism in standard MM systems. 22,25,49Therefore, it is favorable to find ways to deliver transverse energy to outgoing particles while minimizing plasma heating.
An avenue that, as far as we know, has not yet been explored in MM systems is to apply traveling RF fields that are resonant only with the outgoing particles due to the combination of frequency detuning and Doppler shift.For each MM cell, we consider an electric field induced by external electrodes of the form where, E RF is the electric field amplitude, ω RF is the field's frequency, and k RF is the axial wave vector.Since the relevant used parameters obey ω RF r c 2 ≪ 1 (see Sec. III), we employ the electrostatic approximation and neglect all higher-order corrections to the electric and magnetic fields.Nonetheless, we have verified this assumption numerically by including the first-order correction to the time-dependent magnetic field for some of the calculated trajectories, which resulted in negligible differences.
The rotating electric field travels along the MM cell in velocity v RF = ω RF /k RF .Therefore, particles with axial velocity v z experience a Doppler-shifted RF frequency in their rest frame For k RF = 0 (spatially uniform RF field), particles of all velocities experience the same driving frequency.However, applying external fields with k RF ̸ = 0 allows differentiating between out-going and in-going particles with opposite v z since the resonance condition, ω rest = ω cyc , where ω cyc = qB/m is the ion cyclotron frequency, depends on the particle velocity, v z .For example, applying driving field with k RF < 0 and ω RF < ω cyc results in Doppler compensation for particles with positive axial velocity, v z > 0, such that they can meet the resonance condition, ω rest = ω cyc , for suitable values of v z .In contrast, the driving frequency in the rest frame of particles with v z < 0 is down-shifted and therefore gets farther away from resonance for the same driving parameters.Similarly, if we pick k RF > 0 and ω RF > ω cyc , particles with v z > 0 approach the resonance, while those with v z < 0 get farther away.
Naively thinking, both scenarios would exhibit a similar selectivity effect between ingoing and outgoing particles.Nonetheless, in the discussion above, we considered particles in a fixed location, e.g., in the mirror midplane where the magnetic field is minimal.But as a particle moves away from the center, ω cyc increases while |v z | decreases with the amplitude of the axial magnetic field that increases when the particle approaches the mirror throat.Namely, For k RF > 0 and ω RF > ω cyc , particles with v z < 0 also approach resonance, resulting in a reduction in the desired selection effect.On the other hand, when k RF < 0 and ω RF < ω cyc only particles with v z > 0 can experience Doppler compensation.Therefore, the latter scenario is more favorable for achieving selectivity.In Sec.III, we explore this selectivity effect analytically and numerically.Interestingly, similar compensation between density gradient and frequency detuning for maintaining a selective and directional resonance condition while suppressing unwanted noise amplification was employed in Raman amplifiers. 68inally, we note that the form of the rotating RF field defined in Eq. ( 2) is a simplified model for a real RF field configuration inside the plasma.In future work, we may study more complex yet more realistic electric and magnetic fields, for example, fields generated by the different types of Nagoya coils 41 or helicon antennas 69,70 , including the effects of plasma screening and field penetration.

III. 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 = 10keV and a random direction.We consider one MM cell with axial magnetic field given in Eq. ( 1), where, l = 1m, B 0 = 1T, and R m = 3. Neglecting collisions and collective effects, we calculate the trajectories of deuterium (D) and tritium (T) ions under the influence of the Lorentz force F = q (E + v × B).The static magnetic field and the time-dependent RF electric field are given in Eqs. ( 1) and ( 2), respectively, while we used a symplectic scheme 71 , which preserves phase space volume.We calculate for both fusion species because a specific RF field might affect the two hydrogen isotopes differently due to the mass difference, yet both would exist in a D-T fusion reactor.
First, let us demonstrate the dynamics of a charged ion in velocity space under the influence of a traveling rotating electric field and the resulting selectivity effect for four sets of parameters, detailed in Table I.For all the parameter sets, we picked the RF field amplitude E RF = 50kV/m, which is large but realistic and results in an appreciable effect, as will be shown below.The values of ω RF and k RF in the different sets were chosen to represent typical scenarios of the diverse dynamical behavior of the system.In Table I, the frequencies, ω RF , are normalized by the ion cyclotron frequencies of deuterium and tritium at the mirror midplane, ω 0,D = eB 0 /m D = 4.75•10 7 s −1 and ω 0,T = eB 0 /m T = 3.17•10 7 s −1 , respectively.The wave numbers, k RF , are given in the units of 2πm −1 .In the simulations, all particles were initialized at the center of a magnetic mirror cell (minimum axial magnetic field).The simulation time was twice the characteristic passage time of the slower isotope, tritium, across one mirror cell, 2l/v th,T , where v th,T = 2k B T /m T = 7.97 • 10 5 m/s is the thermal velocity of tritium.
Since in the absence of the external RF field, the magnetic moment of both bouncing and passing particles is adiabatically conserved, the transverse and axial degrees of freedom exchange energy such that the total energy is conserved.Therefore when the particle experiences a field B at a position z, its perpendicular and axial velocity components, v ⊥ and v z , respectively, can be mapped into those at the center of the mirror, ṽ⊥ and ṽz , where the magnetic field is B 0 via the transformation, Here, we defined s = sign (v z ) to indicate the initial axial direction.Although in the presence of an external RF field, this transformation is not precise, we use this approximated transformation in Fig. 2 to visualize the particle dynamics in velocity phase space under the influence of the RF.In each panel of the figure, we plot the dynamics in the projected mid-plane, ( ṽz , ṽ⊥ ), for many different initial velocities.The different panels are associated with the different sets of RF parameters of Table I for both deuterium (D) and tritium (T).The dashed black lines indicate the approximated loss cones boundaries, ṽ⊥ / ṽz = (R m − 1) −1/2 .The large dots represent the initial velocities v ⊥,0 , v z,0 and the faint lines connected to them depict the time-evolution in the projected mid-plane according to the transformation in Eq. ( 4).The color indicates how much the RF fields affect each particle during the simulation time, according to the metric defined as the maximal displacement in the projected velocity plane ( ṽz , ṽ⊥ ) where the maximum is taken over the simulation time.In the figure, one can see that the maximal effect (red dots) in all figures is focused mainly around vertical lines around a specific axial velocity v z , associated with the Doppler shifted ioncyclotron resonance.
Next, we analytically analyze the Doppler compensation mechanism and develop a simple model for the resonant region of initial conditions in the (v z,0 , v ⊥,0 ) plane.The idea is that particles will be most affected by the RF if they happen to be at a location z with an axial velocity, v z such that in their rest frame, the Doppler shifted RF field, ω rest , which is defined in Eq. ( 3), meets the resonance condition, where, ω c (B) = ω c,0 B/B 0 depends on z through Eq. ( 1).Here, we neglect the radial components of the magnetic fields that are small for particles near the mirror's axis.As a particle travels down the mirror, its Larmor frequency increases with the axial mirror field within the range 1 ≤ B/B 0 ≤ R m , but the Doppler shifted RF frequency also changes since it depends on the particle velocity.Similarly to the transformation (4), from the conservation of energy and magnetic moment, one finds, where, s 0 = sign (v z,0 ).Therefore, the approximated condition for resonance to take place at some field B reads Substituting v z (B) from Eq. ( 8) and squaring the resonance condition yields a quadratic equation for B. An acceptable solution exists if the determinant is non-negative and at least one of the roots falls in the range 1 ≤ B/B 0 ≤ R m .It is noted that the interaction time with the RF is larger near the mirror center, where the magnetic field gradient is smaller.Therefore, we further constrain the effective resonant region to the range 1 ≤ B/B 0 ≤ 1.25 even though in our simulations R m = 3.Note that for a particle with initial velocities in the loss cone, we also require that the resonant axial velocity is in the same direction as the initial axial velocity because such an escaping particle does not have the opportunity to reverse its direction.In other words, the frequency compensation depends on the direction of the escaping particles (inward or outward in the MM system), giving rise to selectivity.
In Fig. 2, we colored (in grey) the acceptable resonant regions in mirror mid-plane velocity space.Notably, the theoretical resonant regions fit well with the single-particle simulation results, meaning all the most affected particles (colored red) are inside the theoretical (gray) regions.However, the 0 2vth, T max(ΔvΔ FIG. 2. single particle simulations for 1000 different initial conditions sampled from the Maxwell-Boltzmann distribution, according to Eq. ( 4).The upper-left text in each panel indicates the set of RF parameters used (see Table I) and particle type used, deuterium (D) or tritium (T).The color (rainbow colormap going from blue to green to red) of each initial condition point represents the magnitude of the RF effect as defined in Eq. ( 5).The grey areas are the resonant regions found in Sec.III.magnitude of the RF effect in velocity space, and even more importantly, the RF trapping probability, are not predictable by this simple model and thus will be addressed numerically in the next section.

IV. RF TRAPPING PROBABILITY
At any given time, each particle can be characterized as being inside or outside the loss cones by checking the local loss cone condition (v ⊥ /v) 2 < B max /B z , where B z is the local axial mirror magnetic field at the particle's location.Here, we assume that the particles are near the vicinity of the mirror axis, so the radial field components are small.Of course, this condition is valid only when there is no external RF, so energy and magnetic moment are conserved.In other words, it correctly determines the particle's final state if the RF were to be instantly turned off.We divide the particles into three populations by their initial conditions being either trapped particles or, right or left, non-trapped particles, where right means escaping outside the MM system, and left means the opposite, i.e., towards the central fusion cell (see Fig. 1).We define the number of particles in each population as N c for captured, N r for right-going, and N l for left-going.Then we track each particle's identity as a function of time i.e., whether or not it crosses one of the loss cone lines, as depicted in Fig. 2 (dashed lines).
The critical parameter required to estimate the overall efficiency of RF plugging in MM machines is the number of "converted" particles, say, that originated at the right loss cone but ended up trapped.We define four transitions: N r→c , N l→c , N c→r , and N c→l .The other two possible transitions, left to right and right to left, were negligible under the influence of the external RF field, and therefore are omitted from the generalized rate equation model (see Sec. V).
To get relative transition quantities, we normalize each converted particle quantity by the initial number of particles in that population.These normalized transition quantities are plotted in Fig. 3 for the different RF parameter sets of Table I.For example, the blue line represents the number of particles that originated in the right lose cone that at time t were in the capture region in velocity space, i.e., N r→c /N r,0 , where N r,0 is the initial number of particles in the right lose cone in the simulation.It can be seen from the figure that, typically, the number of converted particles grows with time at early times but may fluctuate as particles can change their identity multiple times.In most cases, the number of converted particles saturates at some value after the transient increase.Therefore, for each case, we define and calculate the saturated value by averaging over the second half of the simulated time.It is noted that in some of the cases, the saturation has not been reached yet, so the approximation is not as good as in other cases.We will address this caveat in the next section when studying the overall RF plugging effect.We denote the saturation values by Nrc , Nlc , Ncr , Ncl and plot them by dashed lines in Fig. 3, where the values are recorded in the legend of each panel.
The next step is to study the RF effect for a wide k RF , ω RF parameter space.We calculate the saturation values of con- verted particles for the four relevant transitions and plot the results in Fig. 4 for both deuterium and tritium.It can be seen that the maximal transitions from the two loss-cones populations to the captured population, i.e., Nrc and Nlc (panels (a), (b), (e), and (f)) are concentrated along a straight line in k RF , ω RF space.These lines correspond to the cyclotron resonance condition that depends on the Doppler compensation between the frequency detuning and the RF wave velocity as described in Sec.II.In dashed black lines we overlaid the resonance condition of Eq.( 3), where v z is replaced by the mean axial velocity in the loss cones Here, is the Maxwell- Boltzmann distribution function and the integral is taken only over the loss cone section of the velocity space.One finds vz,LC = where for the considered mirror ratio, R m = 3, the result is vz,LC = 1.025 v th .In panels (a,b) of Fig. 4, the thermal velocity, v th = 2k B T /m, is calculated with the deuterium mass, and in panels (e,f) with the tritium mass.Notably, the right and left loss cones have the opposite mean axial velocity directions and therefore the opposite slope of the resonance line in the figure.We can see that the theoretical lines match quite well with the peak transition rates found in simulations.Also, one can see that the results for deuterium (a-d) and tritium (eh) look very similar, except the latter are downshifted toward lower values of ω RF .This is because the resonance condition is associated with the cyclotron frequency that is inversely related to the ion's mass, and tritium is heavier than deuterium.
Since RF plugging is based on the asymmetric effect of leftand right-going particles, it is elucidating to focus on the ratio, Nrc / Nlc .We call this ratio the right-left selectivity, and we draw it for both deuterium and tritium in Fig. 5 for the same parameter space as Fig. 4. One can again see that, due to the mass difference, the results for tritium (Fig. 5b) look similar but downshifted in frequency compared to those of deuterium (Fig. 5a).In the figure, we also indicated by numbered blue circles the locations (in k RF , ω RF plane) of the parameter sets used in the examples of Table I and Figs. 2 and 3.
Finally, we repeated the single particle simulations for a half-amplitude RF electric field, i.e., 25kV/m.We found that although the transition rates in the optimal parameters regions decrease by a factor of two, the selectivity, Nrc / Nlc , changed by a few percent only.We conclude that the selectivity effect is not sensitive to the RF amplitude because of its resonant nature.However, the overall plugging efficiency might be more sensitive and should be studied and optimized for a specific system in future research.

V. A GENERALIZED RATE EQUATIONS MODEL
In previous work, we developed a semi-kinetic rate equations model the MM system, which includes Coulomb scattering within each cell and transmissions between neighboring cells through the loss 49 Next, the model's definitions and assumptions and introduce new terms for a driven transport induced by the external RF The ion-ion Coulomb scattering is by ν s and roughly scales with density and as ∝ n/T 3/2 .The inter-cell transmission rate is approximated by ν t = v th /l, where v th is the ion thermal velocity l is the mirror cell length.The particles are divided into three populations per cell, captured particles due to the mirroring effect, and and left-going particles through the right and left loss cones, respectively.The densities of the three populations in the i'th cell are denoted n i c , n i r , and n i l .The steady-state solution is defined by ṅ = 0 for all cell populations.The outgoing flux between two neighboring cells (e.g., i . In steady state, definition, inter-cell are the and denoted φ ss .overall system time inversely lated to steady-state flux, namely, ∝ 1/φ , which aim to maximize for satisfying the criterion in fusion system. 50,51In the absence RF fields, confinement time expected to scale linearly with the of MM cells N, in with simple one-dimensional models. 22,25,49even for optimized system parameters optimistic thermodynamical regimes, fusion scales an impractical number of cells. 49Therefore, as discussed we propose apply an RF that asymmetrically affects the particle transport in MM systems, recaptures more escaping (right-going) than returning (left-going) particles.
Here, we study the accumulative plugging by including the RF-induced transition terms in an rate equation model.In Sec.IV, we quantified the amounts of particles that transport between the different populations of one MM cell as Nrc Nlc , Ncr , Ncl under the influence of traveling rotating electric field for an extensive parameter range (see Fig. 4).The characteristic time for a non-trapped particle to pass across one mirror cell depends on its thermal velocity via τ th = l/v th .Therefore, estimate the RF rates by the ratio between the relative conversion amounts and the acteristic time scale, i.e., ν RF,rc = Nrc /τ th , ν RF,lc = Nlc /τ th , ν RF,cr = Ncr /τ th , and ν RF,cl = Ncl /τ th .We have generalized the rate equations model, 49 to include these four induced transition rates such that the generalized model reads The normalized loss-cone solid angle is α = sin 2 (θ LC /2), where the loss-cone angle satisfies sin . Our convention is that the left side of the MM section is connected to the main cell, and the right side is the exit.The boundary conditions that close the rate equations are a constant density in the left cell (n 1 c + n 1 l + n 1 r = n 0 =const., which practically reads n 1 r = n 0 − n 1 c − n 1 l since n 1 l and n 1 c are determined by the rate equations themselves), where n 0 is the central fusion cell density, and a free flow boundary condition on the right cell (ν N+1 t n N+1 l = 0), i.e., no left going flux from outside the system.The Coulomb scattering rate, ν s , depends on the thermodynamic scenarios along the MM section. 49Here, for simplicity, we focus on the common isothermal scenario, which can be justified by electron thermalization since electrons are much faster than ions.
The classic MM operates best when the Coulomb scattering mean free path (MFP), λ = v th /ν s , is of the order of the mirror cell length, l, where the assumption of one-dimensional diffusion between cells is valid.However, for plasma parameters commonly considered for D-T fusion machines, n = 10 21 m −3 and k B T = 10keV, the MFP is of the order of kilometers, making MM systems impractical.It is possible to tweak the plasma parameters (increasing the density and reducing the temperature) to get a lower MFP for the same magnetic pressure.0][51] It is notable that for symmetric RF transition rates, i.e., when ν RF,rc = ν RF,lc , the role of the RF terms in Eqs. ( 12)−( 14) is equivalent to that of the Coulomb scattering.Therefore, by using RF, one can circumvent the difficulty of choosing suitable plasma parameters (i.e., temperature and density), as the RF transition rates depend on the externally applied field rather than on the plasma parameters.Furthermore, when the RF trapping rates are not symmetrical, we can get a right-left selectivity such that the trapping becomes more favorable for escaping particles than for returning particles.Exploiting this effect, the plugging of the MM section can be based on the RF effect rather than on collisions.As a result, one can choose plasma parameters that are preferable for fusion in the central cell without worrying about the collision rates in the MM sections.
We demonstrate the RF plugging effect by solving Eqs.(12-14) for a steady-state flux.In Fig. 6, we present the density profiles and the steady state flux obtained in simulations for different numbers of cells N and different sets of RF parameters.Each set of RF parameters corresponds to four different values of the RF rates that also vary for deuterium and tritium (see Table I).We note that in the rate equations model, we can solve only for a single particle type, which does not describe well a D-T system that composes both deuterium and tritium.However, in the cases where the Coulomb scattering rates are negligible compared to the RF, we can assume both species do not interact.We, therefore, solve the rate equations for deuterium and tritium separately.where V is the volume of the central fusion cell and is the minimal required confinement time due to Lawson's criterion for ignition.0][51] The factor of 1/2 in Eq. ( 15) stands for the two sides of the systems.Further higher-order adjustments to the Lawson criterion for RF-plugged MM systems, such as the RF energy deposited in the escaping particles, are left to a future study.For a practical example, if the main cell length is 100m and its diameter is D = 0.5m (plasma volume in the central cell of V ≈ 20m 3 ), then φ Lawson ≈ 3 • 10 22 s −1 .One can see that the flux drops exponentially with N in the cases with high rightleft selectivity.This exponential drop allows the system to reach the fusion relevant regime φ ss ≪ φ Lawson with a reasonable MM section size, i.e., tens of cells and less than 100m total system length.This should be contrasted with Fig. 5 in Ref. 49, which shows that without the RF terms, the steady-state flux drops inversely linear with N, leading to an impractical number of cells for reducing the flux to a required value.
Finally, the mass difference between deuterium and tritium results in a frequency shift in the particles' response to the RF fields, as shown in Figs. 4, 5. Consequently, the optimal parameters for RF plugging are also different, as illustrated in the examples of Table I and in the rate equations results presented in Fig. 6.Importantly, there is an overlap region in parameter space, so a single RF field can plug both isotopes.For example, in sets 3 and 4, the selectivities for D and T are of the same order (less than a factor of two).The differences in the flux suppression between the two species (solid and dashed lines for D and T, respectively) are minor in these examples (green and red sets in Fig. 6).On the other hand, for set 1 (blue), the flux suppression effect is more dominant in tritium than in deuterium and vice versa for set 2 (black).Interestingly, the scenario of set 1 can be advantageous for fusion reactors fueled and heated by injecting a high-energy neutral deuterium beam.This is most of the fusion power comes from the reactions between the high-energy, injected deuterium and the thermal tritium.Therefore, it might be desirable in this case to remove the slowed-down thermal deuterium particles while plugging the tritium solely. 15,18,72.CONCLUSIONS Axial confinement in MM fusion machines can be significantly enhanced by applying external RF fields that are resonant mostly with the outgoing particles.The considered plugging field is a rotating electric field with a frequency slightly detuned from the ion cyclotron frequency and with a non-zero axial wave vector, such that in the particles' rest frame, the Doppler shift compensates the frequency detuning.The resonant interaction would then increase the transverse energy of the escaping fuel particles and thus move them outside of the loss cone, i.e., recapture them.On the other hand, the incoming particles move in the opposite direction, from the MM section toward the central fusion cell.Therefore, they experience an opposite Doppler shift, i.e., away from resonance with the external field, resulting in an insignificant change in their transverse energy.
We studied this asymmetric effect by sampling initial conditions from a thermal distribution and employing single particle simulations to find the trajectories under the mutual influence of the mirror magnetic field and the RF electric field.We developed theoretical criteria for the resonant region in velocity space, which agreed well with the single particle calculations.From the Monte Carlo simulations results, we extracted the net transport rate of particles between the three populations of the rate equations model, trapped, escaping, and returning particles.Because the re-trapping effect depends on the particle's mass and electric charge, we calculate the transition rates for each set of RF parameters for both species of D-T plasma.We incorporated the RF transition rates into an extended semi-kinetic rate equations model for MM and solved it for steady state.
The main result is that adding RF plugging significantly changes the scaling of the axial confinement time with the number of MM cells from linear (without RF) to exponential (with RF).Remarkably, this flux reduction could help satisfy the Lawson criterion for a reasonable system size.We also demonstrated that although the plugging effect is different for deuterium and tritium, there are RF parameters for which the results are similar for both species.However, it is noted that a confinement enhancement of tritium alone could be advantageous for fusion machines that include high-energy, neutral deuterium beams.Despite the promising results, a few serious caveats remain for us to address in future studies.First, the rotating electric field, which increases the kinetic energy of the particles, also heats the plasma to temperatures that might be deleterious for confinement and stability in some scenarios.Second, the measure of the external electric field penetration into the plasma, the excitation of various collective modes and waves, and the issue of wave backscattering from the plasma should be carefully examined theoretically and experimentally.Third, the problem of radial stability, which is crucial for linear machines, is out of the scope of this paper.Finally, the amount of plugging that can be induced by RF fields generated by realistic antennas like Nagoya coils 41 or helicon antennas 69,70 and the influence of interactions between the particles on the RF plugging are also left for future study along with other types of RF fields, such as a rotating magnetic field.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

FIG. 1 .
FIG.1.An illustration of one side of a MM system (top) and the axial magnetic field of two MM cells used in this study (bottom).

FIG. 4 .
FIG.4.The saturation values (as in Fig.3) for a wide range of k RF , ω RF values for both deuterium and tritium.Nrc in panels (a, e), Nlc in panels (b, f), Ncr in panels (c, g), and Ncl in panels (d, h).Note the color scale is different for the saturation metrics pairs Nrc , Nlc and Ncr , Ncl .The black dashed lines in panels (a,b,e,f) indicate the theoretical resonant lines described in Sec.IV.

7 FIG. 5 .
FIG. 5. Right-left selectivity, Nrc / Nlc , as a function of the RF parameters k RF and ω RF for deuterium (a) and tritium (b).The locations of the parameter sets of table I are indicated by numbered (from one to four) blue dots.

Fig. 6 (
a) shows the steady state density profiles for MM systems with N = 30 cells.Fig. 6(b) shows the steady state flux as a function of the total number of cells, N. The flux in the figure is normalized by the Lawson flux, φ Lawson = nV 2τ Lawson ,

FIG. 6 .
FIG. 6. Steady-state simulation results of the MM rate equations models, including (a) the ion density profiles as a function of cell number for an MM section with N = 30 cells, and (b) flux normalized by the maximal Lawson flux φ ss /φ Lawson as a function of system size.Presented are different sets of RF parameters and different isotopes (D and T), as described in the legend.

TABLE I .
The four sets of RF parameters (frequency and k-vector) in this work.The rows denoted by D and T refer to the isotopedependent saturation values of the transition probabilities, Nrc,lc,cr,cl , and to the selectivity parameter, s, calculated in Sec.IV.
N r→c /N r,0 (blue), N l→c /N l,0 (green), N c→r /N c,0 (red) and N c→l /N c,0 (orange).The dashed lines indicate the mean of each curve in the second half of the simulation time.These saturation values, Nrc , Nlc , Ncr , and Ncl are recorded in the legend and summarized in TableI.The upper-right text in each panel indicates the RF parameters set number and the particle type (see TableI).