Centrifugal particle confinement in mirror geometry

Magnetic mirrors are often considered as configurations for fusion energy. A typical charged particle tightly gyro-orbits the magnetic field and traverses parallel to the field at speed $v∥$⁠. The parallel motion is impeded by the mirror force, $μ∇∥B$⁠. For large enough $v∥$⁠, the particle can escape the mirror; for smaller $v∥$⁠, the particle is trapped between the two high points of the B field, and mirrors back and forth. This condition describes a loss cone in velocity space, with $v∥/v⊥>Bmax/Bmin−1$ the condition for loss. Collisions, however, can transfer energy to $v∥$ at the expense of v_⊥, and all particles eventually escape. The mirror plasma thus has a particle confinement time of order the ion-ion 90° collision time, 1/ν_ii. To satisfy the Lawson criterion for net energy breakeven with a deuterium-tritium plasma at a temperature above 12 keV, one must have nτ > 1, where τ is the energy confinement time in seconds and n is in units of 10²⁰ m⁻³. If this is based on the ion loss time, nτ for a mirror trap is about 100 times below breakeven at an ion temperature of 10 keV. (Note that the electron heat losses are even faster, on account of the higher collision rate.)

A radial electric field, E, can enhance the confinement.^1–5 In the simplest description, a radial E field will induce both species to rotate azimuthally at the well-known drift velocity $u→E=E→×B→/B2$⁠. This plasma rotation creates a radially outward centrifugal force; the centrifugal force has a component parallel to the B field which always points towards the midplane, thus working to slow down $v∥$⁠. An estimate of the work done in the centrifugal potential is $muE2/2$⁠; thus, from parallel energy conservation, the condition $(1/2)mv∥2<(1/2)muE2$ is an estimate for the rotation speed needed to completely confine a particle. Thus, a mirror plasma rotating at supersonic speeds will effect fairly complete collisionless confinement, with exponentially small losses for a Maxwellian population. In addition, since there is no toroidal bending, the classical radial diffusion from the density gradient has no neoclassical enhancement although it could be turbulent.

This fusion concept has been investigated theoretically and experimentally at the University of Maryland^6–9 and shows promise for a steady state fusion reactor with good confinement, stable to magnetohydrodynamic (MHD) modes. Previous analysis has examined MHD stability, and although a rigidly rotating plasma is unstable to flute interchange modes, supersonic sheared flow stabilizes this mode.^10–13 In the experiment, the attainment of high Mach number rotation was achieved by biasing a metal rod running down the center of the plasma. In a reactor size device, an electric potential of several MeV is necessary, so the production and maintenance of this field is non-trivial. Alternatively, the rotation could be achieved using beam injection or possibly through the use of antennas. Kinetically, particle confinement has previously been studied using a guiding center approximation, in which a particle can be treated as being attached to a magnetic field line, behaving as a bead on a rotating wire.^5,6

This paper is concerned with quantifying better the dependence of the confinement time on the Mach number M, especially in the transition from supersonic to subsonic rotation, the dependence on particle collision rates, and the effect of the nonzero gyroradius, using a full cyclotron orbit code. The Mach number M refers to the rotation velocity with respect to the velocity of a particle at the thermal temperature of the Maxwellian distribution before the electric field producing rotation is applied. A full cyclotron orbit analysis reveals that the energy in the rotating frame is proportional to the square of the Mach number. The gyrofrequency is unchanged, but the gyroradius and velocity increase linearly with the Mach number. Thus, measured with respect to the particle velocity in the rotating frame, the Mach number never exceeds 1. Note that we are examining only ion transport. Electron losses, discussed later, are not fully addressed in this paper. In particular, the energetic ions would slow down on the electrons, with respect to the E × B frame, thus increasing the Mach number.

A simple analysis⁴ gives the asymptotic values of the parallel loss rate as ν_ii for M ≪ 1 and $νii/[M2 exp(M2/4)]$ for M ≫ 1, with the exponential behavior due to losses from the high energy Maxwellian tail. However, this analysis considers particles attached to a single field line rotating at Mach M, and the result must then be averaged over the rotation and density profiles. In addition, it ignores the effect of radial diffusion given by a nonzero gyroradius, and also, as we will see, the fact that the radial electric field increases v_⊥ in the rotating frame, essentially closing the loss cone. In this work, we examine the effects of the nonzero gyroradius and the particle collisions so as to obtain a more precise answer for the rotation rates required to produce effective classical confinement.

Test particle orbits are stepped forward in time in a system with a prescribed mirror magnetic field, $B→(x→)$⁠, and a model for the electric field, $E→(x→)$⁠. We use a full cyclotron motion code, stepping

We initiate the particle distribution as a Maxwellian deuterium plasma at a given temperature before applying the rotation. The temperature is taken to have the same radial profile as the rotation and the density, peaked at mid-radius.

The electric field arises as a collective plasma effect. As such, it will be the electric field in a magnetohydrodynamic (MHD), small Larmor radius plasma. This field is obtained as follows: Let the magnetic field be given by a flux function ψ(x) according to

We have, in steady state, $E→=−∇Φ, Φ(x)$ being the electric potential. Then, according to ideal MHD

Noting that Φ = Φ(ψ), we readily invert Eq. (3), using Eq. (2), to find

i.e., the rotation is that of a rigid rotor at angular frequency $Φ′(ψ)$ with respect to each ψ surface. Thus, a typical particle orbit is calculated according to

and

describes in simple terms the confining mirror force $−μ∇∥B$ and the centrifugal force $rΦ′2(ψ)∇∥r$⁠. We note that in the MHD ordering, $|E→⊥|≫E∥$⁠, the parallel electric field. Potential drops eΦ along $B→$ are of order T_e, while across the field they are much higher, of order (a/ρ_i)T_e. Here T_e is the electron temperature, ρ_i the ion Larmor radius, and e the electric charge. Thus corrections to eΦ should be allowed.

In Sec. II, we construct a model mirror field. Section III consists of a small gyroradius analysis of particle confinement in a rotating mirror plasma. Section IV discusses single particle orbits in a cylindrically symmetric magnetic field with a strong radial electric field, both analytically and by simulation using the model field. Section V is a numerical examination of orbit properties during reflection from the mirror. Section VI discusses collisionless loss as a function of mirror ratio and rotational Mach number. Section VII concerns analytic estimates of collisional particle and heat losses using the small gyroradius approximation. Section VIII gives simulation results for collisional loss. Section IX considers how these losses are modified by the electric field produced at the mirror ends by the ion density gradient. Section X considers alpha particle confinement in this system. Section XI discusses the use of a rotating mirror plasma to produce an ion thruster, and in Sec. XII are the conclusions.

For our analysis, we choose a simple cylindrical mirror field with cylindrical coordinates r, θ, z, given by

The magnetic flux inside radius r is

Choose the maximum radius at z = 0 to be r_w, giving the flux at this radius to be $ψw=rw2B0(1−A)/2$⁠. The bounding flux surface is then located at

shown in Fig. 1 for k = 6.25 × 10⁻¹⁰π, with the mirror at ±200 cm, A = 0.75, giving a mirror ratio of 7, and r_w = 100 cm. The parameters are chosen to be close to the reactor design given in Ref. 6. This analytic representation allows very accurate calculation of particle orbits and guarantees that $∇·B→=0$ and also $E→·B→=0$ to double precision. Note that additional harmonics can easily be added to change the shape of the bounding flux surface if desired. Also shown in Fig. 1 is the vertical centrifugal force vector and its component along $B→$⁠. The use of cos(kz⁴) in the form of B produces the extended flat domain about z = 0 seen in Fig. 1 and a strong component of the centrifugal force along $B→$ at z ≃ ±150 cm. This choice of coil design produces a large central confinement domain and strong confining centrifugal forces, and is an improvement over the more gentle change of B_z in the experiments carried out in Maryland.

Now introduce a radial electric field to model the plasma rotation. The experiments indicate a rotation profile approximately parabolic in minor radius, vanishing at the center and at the edge, where the plasma is in contact with confining surfaces. We model this by noting that ψ is proportional to r², so we take the potential to be

with V the potential drop from the plasma center to the edge. In fact, in the simulations the potential is taken negative at the center, the same choice as was made in the experiments. In the experiments, this is necessary because with a positive choice the electrons immediately escape out the ends and plasma breakdown is not achieved. The electric field as a function of radius is shown at z = 0 in Fig. 2. This is in close agreement with Fig. 2 from Ref. 7, giving an experimental determination of the plasma rotation profile. The plasma density and temperature are also observed to be approximately parabolic in radius and we model them simply at z = 0 through the parabola r(r_w – r). In general the density and temperature are functions of ψ, but we initiate all particles at z = 0 so this expression suffices.

For small gyroradius, in the collisionless limit, the particles can be considered attached to the field line.⁶ Thus the particles radial position is r = r(z), and the Hamiltonian is

where w is the azimuthal rotation frequency. Let z_l be the value at the mirror throat, and z₀ the center. Also $B(zl)=Bmax, r(zl)=rmin, B(z0)=Bmin, r(z0)=rmax$⁠.

Equating H at the turning point and at the midplane yields an equation for the loss cone, now modified by the presence of the plasma rotation

Assuming that the rotation frequency from the E × B drift (w_E) is independent of z, yields

Further using $ż0=v∥$⁠, we have

For B_max/B_min ≫ 1 and r_max/r_min ≫ 1, the centrifugal force associated with a V_E(r) traps all the particles with parallel velocities less than V_E(r). For the effect of the rotation to exceed the confining effect of the mirror then requires $VE(r)/v⊥>Bmax/Bmin$⁠. Note that in this guiding center analysis the rotation does not narrow the loss cone, and it gives a restriction on the value of $v∥2$⁠. Thus it limits loss to the high energy tail of the Maxwellian, independent of the ratio $v∥/v⊥$⁠. In fact, we will see that the rotation increases v_⊥ in the rotating frame so that the loss cone is closed.

In one collision time, the fraction of lost particles predicted by the guiding center analysis is

with E₀/T = M²,

Furthermore, the collision frequency scales as 1/v³ ∼ 1/M³, so the loss in a given interval of time for large M is

We will see in Sec. IV that the guiding center analysis is importantly modified at significant Mach numbers, in two ways. First, the electric field increases the gyroradius and thus finite Larmor radius effects on classical transport need consideration. Second, radial excursions of the ions from a flux surface, via collisions, produce an essential modification of the particle energy and confinement, and parallel loss of particles in the Maxwellian tail is not the primary loss channel. The perpendicular loss is, of course, present in any magnetized plasma confinement system, and must be accounted for as a concomitant loss channel. In this sense, we regard Eq. (17) as a baseline loss which accounts for the parallel transport only.

Particle orbits using full cyclotron analysis show some very important differences from what a guiding center treatment gives. The profile of the electric field means that one cannot transform away its effects by moving to a locally rotating frame. Particle orbits, due to the cyclotron radius, sample different parts of the field, and $E→·v→$ is not zero.

In the central confinement domain, near z = 0 the radial B field is negligible and the field can be approximated by a constant B_z. In this case, the particle motion can be found analytically.

We study the motion of a charged particle in a constant magnetic field $B=B ẑ$ in the presence of a radial electric field $E=E(r) r̂$⁠. The equation of motion Eq. (1) is decomposed in cylindrical components as

In what follows, we ignore the parallel motion along the z-axis and focus our attention on the motion in the plane perpendicular to the magnetic field.

We begin our investigation by, first, considering the case of a constant radial electric field. Using cylindrical coordinates in straight magnetic-field lines $B=B ẑ$⁠, the equations describing gyromotion in the plane perpendicular to the magnetic field are

where E₀ denotes the constant electric field and the terms $(r θ̇,−ṙ)$ represent the effects of the magnetic force, with Ω = eB/m.

If we multiply Eq. (22) with r, we obtain the conservation law

Using the initial conditions

we find

from which we obtain the angular velocity

After substituting this solution into the radial equation (21), we obtain

where V₀(r) denotes the effective radial potential in the presence of a constant electric field.

The effective radial potential V₀(r) has a single minimum at a radial position r = R₀, which is defined by the relation $V0′(R0)=0$⁠, which yields the identity

where ν ≡ 2ω/Ω and E₀ = E(R₀) is the (constant) electric field at r = R₀. The dimensionless parameter ϵ ≡ (E₀/B)/(R₀Ω) is defined as the ratio of the E × B velocity (which is positive in the clockwise angular direction) to the cyclotron tangential velocity R₀Ω (which is in the clockwise angular direction for positively charged particles). Hence, the sign of ϵ is determined by the sign of E₀ (e.g., ϵ > 0 when E₀ > 0). Another interpretation of the dimensionless parameter ϵ = (E₀/BΩ)/R₀ is given as the ratio of the radial polarization shift (E₀/BΩ) to the radial position R₀.

We now investigate the radial motion in the vicinity of R₀ by writing r = R₀ + ρ, with $|ρ|≪R0$⁠. Hence, Eq. (27) becomes the simple-harmonic radial equation $ρ¨=−ω02 ρ$⁠, where the simple-harmonic frequency is

where we used identity (28). Equation (26) for the azimuthal motion, on the other hand, becomes

where Eq. (28) is used once again, with the assumption $|ρ|≪R0$⁠. We now return to the initial condition $θ̇0=ω$ to find

where we defined ρ₀ ≡ r₀ – R₀. Hence, we find the dimensionless parameter $ν(ϵ,ρ0)≡2 ω(ϵ,ρ0)/Ω$ expressed as

and Eq. (30) becomes

We note that, if $|ϵ|≪1$⁠, then Eq. (32) yields $ν(ϵ,ρ0)≃−2 (ϵ+ρ0/R0)$⁠, which vanishes if $ϵ=−ρ0/R0≡ϵ0$⁠, so that $ω0(ϵ0)≃Ω$ and ω(ϵ₀) ≃ 0.

The solution for ρ(t) and θ(t) is now simply expressed in terms of trigonometric functions as the polar-cycloid solution

where ω₀ and ω are defined in Eqs. (29) and (31), respectively, with ν defined in Eq. (32).

The physical motion of the charged particle in the (x, y)-plane perpendicular to the magnetic field is therefore expressed as

which is shown in Fig. 3 for R₀ = 49.1 cm and ρ₀ = 0.9 cm (so that r₀ = 50 cm), and ϵ = –0.019 is chosen so that ν ≡ 0 in Eq. (32).

We note that the solution Eq. (36) exhibits a net radial variation Δr = 2 ρ₀ ≃ 2 cm. In addition, the choice of ν = 2ω/Ω ≡ 0 implies that the counterclockwise cycloidal motion in Fig. 3 exhibits no retrograde motion and the radial excursion Δr is directly proportional to the constant radial electric field E₀.

To find the orbit and kinetic energy in the rotating frame, simply drop the secular terms in the expression for θ in Eq. (35). In Fig. 4 is shown the orbit and kinetic energy of a particle in the rotating frame for the parameters given above. The effect of the radial electric field is to increase the gyroradius, approximately linearly in the Mach number, with very little variation of the particle energy around the orbit.

We now consider the case of a radial electric field with constant gradient

where $E0′$ denotes the radial electric-field gradient at r = R₀. The effective radial potential defined in Eq. (27) is now modified to include the effect of the radial electric-field gradient

where R₀ is still defined as the minimum of the new potential V(r), i.e., Eqs. (28) and (33) are still valid.

The addition of a constant radial electric-field gradient modifies the simple-harmonic frequency (29)

while the form of the polar-cycloidal solution (34)–(35) remains the same. Here, we note that the dimensionless ratio $E0′/(BΩ)$ is defined as the gradient of the radial polarization shift E(r)/(BΩ).

Simulations are made using a fourth order Runge-Kutta code for advancing Eq. (1).

At Mach 6, the gyroradius for the parameters chosen for our simulations with nominal gyroradius of 0.17 cm is 1 cm, shown in Fig. 5, in agreement with Figs. 3 and 4. The gyroradius increases linearly with M. This is only one of several confirmations of the numerical simulation carried out by comparison with analytic theory. For the parameters chosen, Mach 1 requires a peak field at center of 250 kV and Mach 6 a field of 1.5 MV. The increase in the gyroradius gives an increase in the radial diffusion rate as well as an increase in the particle energy in the rotating frame to over 12 keV. For this reason, we have chosen a design with a large radius and a sufficiently strong field to make the radial diffusion loss less important. For fusion breakeven, both the parallel and perpendicular confinements must be satisfactory. The fact that the radial electric field causes a value of v_⊥ linear in the Mach number also means that there are no particles in the loss cone. Shown in Fig. 6 are the particle distributions in the local rotating frame for Mach 3 and Mach 6, with a mirror ratio of 7. The only particles near the loss cone are those at large r where the rotation is weak.

It is impossible for a particle to have small v_⊥ in the presence of the rotation. This increase in the particle energy to a velocity given by the rotation is similar to the phenomenon of “pick up ions” studied, for example, in planetary atmospheres^14,15 and in the solar corona.¹⁶ 

In Fig. 7 is shown the orbit and kinetic energy of a particle initially launched at 500 eV, rotating at Mach 6. Above Mach 1, there is no retrograde motion in the lab frame. The gyroradius is about 1 cm. Note the gyroperiod is unchanged. This is in agreement with Eq. (39) since the orbit is at r = 50 cm, where $E′=0$⁠. At other radii, the gyrofrequency is slightly modified according to this equation. The radial extent of the orbit and the frequency agrees with the analytic results in Sec. IV A. By Mach 6, we mean that the potential is six times the potential required to produce Mach 1, where the gyroradius is not significantly altered. Since the gyroradius, and hence the gyroorbital velocity, increases linearly with the Mach number, we find that, in fact, if the Mach number is defined with respect to the actual velocity in the rotating frame, one can never significantly exceed Mach 1.

The analytic sections above give two key results, confirmed in the simulations. First the radial electric field, besides resulting in plasma rotation, stretches the gyroradius linearly in the applied field, while the gyrofrequency is relatively unchanged. Hence the particle kinetic energy around the gyro-orbit in the local rotating frame increases as the square of the Mach number. These results have significant consequences for the analysis of confinement and fusion rates. The increased radius means an increase in the classical radial diffusion rate, and the energy in the local rotating frame is the relevant energy for fusion reactions. Thus one can initiate a 500 eV plasma, which upon rotation at Mach 6 acquires a thermal distribution above 12 keV.

In this section, we examine the trajectory of a particle during reflection from the mirror. The particles are initiated with $v∥/v=1$ so the gyroradius and value of v_⊥ are due solely to the rotation. One example is shown in Fig. 8, which was for a particle of 25 keV with Mach 6 rotation. The last plot shows that the deviation from the flux surface, smaller than the flux variation over a gyro-orbit, changes sign with the change in the sign of $E→$⁠. Several different energies and Mach numbers were used. Detailed plots at different times were used to numerically determine the gyroradius and frequency at different values of B_z. For all cases, within numerical error, ρ ∼ 1/B and ω ∼ B, so in the rotating frame v_⊥ was found to be constant during reflection.

To study collisionless loss launch particles with a full range of r and pitch $λ=v∥/v$⁠. Particles are deposited at z = 0 with the radial density profile given in Fig. 2, uniform in pitch, with a local Maxwellian energy distribution. The temperature is taken to have the same radial profile as the rotation and the density. An initial distribution, peaked at midradius at 500 eV before rotation, is deposited in the local rotating frame, which upon rotation produces a Maxwellian of over 12 keV in this frame. Particles are then advanced in time for a large number of mirror bounces since losses may depend on particle gyro-phase during reflection. The Mach number refers to the maximum value, near the midpoint in radius, r = 50 cm, with the rotation profile given by Eq. (2). The collsionless loss as a function of mirror ratio is shown in Fig. 9 for plasma rotation at Mach 0, 2, 4, and 6. At Mach 0, the losses are through the conventional loss cone. Mirror ratios above 2 are very effective in improving confinement. Note that the effectiveness of the rotation depends on the slope of the magnetic field near the mirror, at the point where the centrifugal force is most effective, see Fig. 1. Different results for the dependence of loss on mirror ratio can be expected for different mirror field geometries.

We examine here analytic estimates of parallel particle and heat losses to the insulators obtained using guiding center analysis. For high Mach numbers, the ions are very tightly confined, as has been discussed. Neglecting cross field transport, energetic ions on the very tail of the distribution may escape the centrifugal potential well. This leads to particle as well as momentum loss. Once these energetic ions are lost, further losses are incurred only if the confined ions are scattered into this exponentially small loss cone, via ion-ion collisions. This loss occurs by ion-ion collisions, but at a rate based on the energetic ions; in addition, the number of energetic ions is exponentially small for a Maxwellian distribution. The resulting loss rate can be estimated to be of the order^4,17 given by Eq. (17)

where the M² factor is the well-known Pastukhov factor,¹⁷ which results from the deep centrifugal potential well as well as an electrostatic potential well (discussed below). However, these estimates neglect what we find to be an important competing loss channel, consisting of radial diffusion to the edge where the rotation is weak.

The lighter electrons, on the other hand, are not centrifugally confined. Nonetheless, their tendency to escape is arrested by the build-up of a space charge electrostatic potential along the field line. By quasineutrality, the electron particle loss rate must equal the ion loss rate. The electron heat loss rate, however, is larger and the more important concern here is that energetic electrons can reach the insulators, transfer heat, and return as cold electrons. The parallel electron heat loss is thus based on ν_ee and scales as⁴ 

where eϕ/T_e is the self-consistent electrostatic potential drop along $B→$⁠. The Pastukhov factor¹⁷ again applies; in addition, the potential well is deep, and it can be shown that it scales as the square of the Mach number. Note that this well is deeper than for classical mirrors and, to an extent, controllable. The estimate for the electron heat loss rate is thus of the same order as Eq. (40) with the ion-ion collision rate replaced by the electron-electron collision rate. Based on this, the Lawson criterion can be attained against electron heat losses, at T_e = 10 K for a Mach number M > 6. Experimentally it has been shown that high electron temperature can be reached, and 900 eV values were obtained in mirror geometry.¹⁸ However, the form of the ion density in the mirror throat and the nature of the electron response are essential to understand electron heat transport, and is thus beyond the scope of the present simulations.

To study collisional loss, we load all particles with a parabolic density profile, with a Maxwellian at 500 eV, which upon rotation assumes an energy of over 12 keV in the local rotating frame. The temperature is taken to have the same radial profile as the rotation and the density. Particles are then advanced in time for a large number of mirror bounces, typically for a full millisecond in time or for a full collision time.

The collisions take place among particles in the rotating frame. Numerically each particle is transformed to the local rotating frame to carry out energy conserving pitch angle scattering, and then transformed back to the lab frame.

In introducing collisions, we also introduce viscous drag impeding the particle rotation. Classical viscous drag is much smaller than the collision rate, given by ν_iiρ²/a² with ν_ii the ion-ion collision rate. The direct heating of the ions by the rotation means that the conventional heating of ions through the electrons, occurring in tokamak fusion scenarios, may not be necessary. In the experiments, the observed heating was thought to be due to viscous drag, but it is apparently a natural effect of the rotation.

Particles with different values of pitch and radius reflect at different z values, but the predominant reflection takes place in the region of strong centrifugal force component parallel to $B→$⁠. In Fig. 10 is shown the particle distribution in z after one collision time. It is seen that the confinement is greatly improved by the rotation, with no particles arriving beyond z ≃ ±100, where the centrifugal force begins to have a large component along $B→$⁠. See Fig. 1. In the Maryland experiments, approximately exponential decay of the density in z was observed. The exact nature of the density near the mirror is important because the ions produce an electrostatic confining field for the electrons, which control plasma heat loss. It appears that the coil design we are using gives much smaller ion density in the mirror throat region than what was observed experimentally, but the nature of electron heat loss is beyond the scope of the present investigation.

In Fig. 11, we see that it indeed requires small collisionality to obtain good confinement, but the necessary small value is not unreasonable for typical reactor densities and temperatures. For low collisionality and Mach numbers above 2, the confinement time, limited by end losses, is more than 1 s. Radial diffusion loss is also seen as end loss since upon diffusing radially to a point near the edge particles encounter low rotation and escape through the end.

In Fig. 12 is shown the fractional loss in 1 ms for an initial 500 eV Maxwellian distribution with Mach 6 rotation versus the mirror ratio. The mirror ratio increases the confining component of the centrifugal force, as is clear from Fig. 1. The exponential improvement of the confinement is simply understood and can be seen through an examination of the energy distribution of the lost particles. As the mirror ratio increases, the centrifugal confining force causes higher and higher energy particles to be confined.

Collisional losses are still parallel, but they involve the intermediate path of radial diffusion to large r. This is shown in Fig. 13, giving the distribution in r of mirror end losses after 1 ms with ν_ii = 70/s. All end losses are at the maximum radius at the mirror, r = 35 cm, showing that the loss path consists of radial diffusion to a point near the edge, where the rotation is weak. This is also a confirmation that at radial locations where the rotation is strong, the loss cone is unoccupied.

The Lawson criterion for a 12 keV D-T plasma is approximately nτ = 10²⁰ s/m³. A density of 10²⁰ then requires a confinement time of 1 s, obtainable with Mach numbers of 4 or more and low collisionality, provided the energy confinement time is approximately that given by the ions.

An ion density gradient gives rise to an electric field of the form

with p the plasma pressure and n the density. To examine the effect of this field on the density, the space of r, z is divided into domains, and a running time average of the local density is calculated. This density is then periodically used to calculate the electric field due to the local gradient, and this field is added to the equilibrium field producing rotation. The calculation is continued until there is a convergence of the local density. The code used 10 steps per gyro-orbit, and 10⁵ particles were stepped. The density in each of the 50 × 50 bins in r, z was averaged for 10⁴ steps, giving a value of the local density with approximately 4 × 10⁵ points per bin. This density was then used to calculate a new value of $E→$⁠, and the process repeated until it converged well enough to reliably give the particle density in the presence of a self-consistent field. The convergence of the quantity

where the brackets indicate an average over the bins, is shown in Table I, with n_k,t the time averaged number of particles in bin k at iteration t. With the introduction of the electric field, there is initially a strong modification of the ion density, but after seven iterations the density remains fixed except for a small residual statistical variation.

The result is shown in Fig. 14. The only significant modification of the local density consists of a flaring out of density in the regions of strong gradient, and there is no increase in the number of lost particles. The sharp peaks in density occur at the rotation induced reflection points, where particles reverse direction.

Figure 15 shows the density distribution in z, including the electric field due to density gradient, for Mach 6 rotation. This is compared with Fig. 10, the same simulation without the electric field. An exponential flaring of the density outwards is produced by the electric field, but no significant particle loss. This is understandable because the regions of significant electric field, at the large density gradients, are also the regions where the centrifugal force is very strong, and so it is difficult for the density to extend itself very much in these regions.

At birth, alpha particles within the loss cone are immediately lost because their velocity is much too large to be affected by the rotation. But particles not in the loss cone slow down on the electrons without significant pitch angle scattering, so they remain confined until the pitch diffuses them into the loss cone. The collision operators are, for slowing down dv_α/dt = –ν_sv_α, and for transverse diffusion $d(vα)⊥2/dt=ν⊥vα2$⁠, with $νs=(1+mα/mβ)ψ(x)ν0$⁠, and $ν⊥=2[(1−1/2x)ψ(x)+ ψ′(x)]ν0$⁠, where α refers to the test particle and β the background particle with which it is colliding, $ν0=4πeα2eβ2Lnβ/mα2vα3, x=mβvα2/(2kTβ)$⁠, and L = ln(Λ) is the Coulomb logarithm. So, the pitch diffusion is

The Rosenbluth function

has an approximation asymptotically valid in the limits of both small and large x and with a relative error of dψ/ψ < 10⁻² over the whole range of x given by¹⁹ 

with $p=x1.5[a0+a1x+a2x2+a3x3+a4x4]$⁠, and $a0=0.75225, a1=−0.238, a2=0.311, a3=−0.0956, a4=0.0156$⁠.

Shown in Fig. 16 is the time evolution of the alpha particle energy, normalized to the initial energy of 3.5 MeV, and the evolution of the pitch λ at Mach 6 with respect to the 12 keV ion distribution with n = 10²⁰/m³. Initially the collision operator is dominated by slowing down, and only at lower energies the perpendicular scattering begins to take effect. An alpha particle loses significant energy before scattering into the loss cone, with the velocity arriving below Mach 1 at 6 ms.

An actual simulation following an initial uniform pitch distribution of 3.5 MeV alpha particles reveals that about half of the alphas are lost in the process of slowing down before being confined by the rotation, due to pitch angle diffusion of that part of the population near the loss cone and diffusion to larger r values where the rotation is weak, the same primary loss channel for ions. Of the total alpha particle energy initially present, 34% is lost and 64% is absorbed by the electrons, with only 2% remaining as kinetic energy of the remaining confined alpha particles.

As is well known,²⁰ the loss cone distribution is also subject to an electrostatic instability that can repopulate the loss cone and hence lead to complete particle loss. The dispersion relation indicates wave absorption at the mirror ends with no reflection. This implies that standing waves are not possible, and only the convective instabilities need to be considered. Neglecting mode reflection at the ends leads to a critical device length L, with rapid particle loss only for $L>104m/Mρi$ with m and M the electron and ion masses, and ρ_i the ion gyroradius. For our system L ≃ 10³ cm, much larger than the device, but a more detailed examination of the mode structure with the actual particle distribution and the field configuration is necessary to determine whether this instability is actually operative.

To produce a thruster geometry with a rotating mirror, it is sufficient to introduce asymmetry in z. We take as model

with the parameters chosen to produce asymmetric mirrors at z = z₁, z₂. In Fig. 17 is shown the configuration of a rotating mirror as a thruster along with points showing the test particles in the distribution. We examine a small thruster, with z₁ = –20 cm, z₂ = 20 cm, w₁ = 0.05, w₂ = 0.1, a = 0.525, b = 0.55. The bounding lines give the location of the limiting flux surface. Also shown is the magnitude of B_z(z), with an asymmetry of the mirrors producing one directional loss for thrust, with mirror ratios of approximately 10 and 6. The rotation rate is Mach 2 with a 100 eV proton Maxwellian distribution. The required potential is 900 eV.

In the nozzle, with z > 20 cm, the flux surfaces fan out radially as a function of z. This causes a large increase in the particle gyroradius and also the particle energy as it exits the device. In Fig. 18 is shown the z component of the velocity distribution in the exhaust, about 6 × 10⁵ m/s. The thermal velocity in the distribution is only 1.4 × 10⁵ m/s. This increase is simply understood as being due to the combination of the mirror and the rotation. The mirror restricts the particles to have high pitch and the rotation, in this case only at Mach 2, restricts the exhaust to consist of only a part of the Maxwellian tail. Also shown is the particle distribution across the exit. In addition to the selection of high velocity particles imposed by the mirror and the rotation, the gradient of B in the nozzle acts to increase the particle velocity during exit. In Fig. 19 is shown a spiraling loss orbit, gaining energy as it moves outward in the exhaust region.

The confinement time is 2 ms. Assuming a plasma density of 10²⁰/m³, the mass loss rate is

and the exit velocity is 6 × 10⁵ m/s, giving a thrust F = vdm/dt of about 0.1 N. This would accelerate a one ton vehicle at 10⁻⁵g. This is about the same value produced by a NSTAR thruster.²¹ However the present design would be much heavier and would consume much more power.

Because of the large weight and power requirements, the best solution for a thruster appears to be a large combination fusion energy and thrust producing device, or even two separate devices. Taking a reactor with one cubic meter of volume and a plasma density of 10²⁰/m³ and a confinement time of 3 ms gives

The velocity of the deuterium at 12 keV is 10⁶ m/s, but there is a factor of about 5 coming from the mirror, the rotation, and the nozzle, so the force is about 500 N. For a mass of 1000 kg, this will give an acceleration of 0.5 m/s². Accelerating for a month gives a velocity of 10⁶ m/s, with a fuel consumption of 100 kg of deuterium. Jupiter is currently about 500 × 10⁶ km from Earth. So the flight time to Jupiter would be about one month.

The simplicity of a rotating mirror geometry makes it a very attractive alternative to tokamak or stellarator designs. Since there is no current profile, there are no disruptions and since there is no toroidal bending, the cross field classical radial diffusion is not enhanced by neoclassical geometrical effects, and there are no runaway electrons. The coil design is trivial as compared to either tokamaks or stellarators. The limited experiments conducted at Maryland and the theoretical work done thus far have revealed no significant large scale MHD instabilities, with the sheared rotation profile thought to be responsible for the stability. Rotation at Mach 6 appears to severely limit the mirror end losses to the contribution of those particles which manage to diffuse radially to a point where the rotation is small.

Confinement times with the design examined are greater than 1 s for mirror end loss and are negligible for radial diffusive loss, but it is in fact the radial diffusion to large r values where the rotation is small which dominates the loss. The plasma radius and the field strength were chosen large enough to make the radial diffusive loss manageable, and for reasonable collision rates, the Lawson criterion for ions can be achieved, using fusion rates associated with a 12 keV plasma. Note that these simulations include only classical radial transport, and thus the existence of any additional turbulent radial transport would modify these results. A full nonlinear two fluid MHD simulation is necessary to investigate this. A high resolution nonlinear MHD simulation is also necessary to test velocity shear stabilization at high Reynolds number. Estimates indicate that the loss cone instability is not important, but also this needs to be further investigated.

The rotation rate is far too slow to affect alpha particle confinement, so they are initially lost through the usual mirror loss cone. After this, assuming that there is no significant excitation of the electrostatic loss cone instability, slowing down on the electrons without significant pitch angle scattering insures that over half of the alpha energy is transferred to the electrons. A detailed analysis of the loss cone instability for the given particle distribution and field configuration is necessary to establish this result. Further work including also electron dynamics is necessary to confidently predict heat confinement times. In addition, the nature of MHD stability at high plasma beta must be addressed.

Introducing small asymmetry of the design produces very asymmetric loss, leading to a possible design of a thruster, but large weight and power requirements limit its use to a large space vehicle.

This work was partially supported by the U.S. Department of Energy Grant DE-AC02-09CH11466. We gratefully acknowledge encouragement from and multiple discussions with Curt Bolton, DoE, and conversations with Greg Hammett.