Formation of cylindrical plasma equilibria with β > 1

The possibility to exceed the MHD limit $β=1$ seems very attractive for the controlled fusion problem, since it increases the efficiency of magnetic plasma confinement and makes it possible to achieve plasma fusion parameters at weaker magnetic fields. The purpose of this work is to answer the question of how strongly the $β$ parameter can exceed the MHD limit if the currents that flow along the plasma boundary and are responsible for displacing the magnetic field from its volume are concentrated in a transition layer of sub-ion scale. For simplicity, we will limit ourselves to the same idealized problem of confining an infinitely extended cylindrical plasma column in a uniform vacuum magnetic field that was studied in the paper.¹² A significant decrease in the width of the transition layer compared to Kotelnikov's predictions [ $λ≈(6−8)ρi$] can occur if a jump in the electric potential is formed in the layer.¹³ This jump is the result of charge separation caused by deeper penetration of ions into the strong magnetic field. In this case, the electric field not only reflects the ions from the transition layer, preventing them from full Larmor rotation but also forces the electrons to drift in crossed $E$ and $B$ fields, thereby creating the azimuthal current necessary for equilibrium. This kind of equilibrium (with electrostatic confinement of ions and electron $E×B$–drift current) has been recently observed in magnetic holes of sub-ion scale (⁠ $<ρi$⁠).⁶ In addition, recent particle-in-cell (PIC) simulations of the formation of high- $β$ equilibrium in a multicusp trap¹⁴ have shown that the appearance of a jump in the electric potential localizes the electron $E×B$-drift current at a size comparable to the electron gyroradius (⁠ $λ≈4ρe$⁠). It seems obvious that decreasing the width of the transition layer should increase the non-gyrotropic effect and the associated maximum value of $β$ because of the stronger influence of inhomogeneity on the Larmor rotation of ions. However, to check this, we need to answer the following questions.

First, it is necessary to find out whether the size of the current sheet in the considered problem can really be reduced to several Larmor radii of electrons. The fact is that currently there are different opinions on this issue. For example, in PIC simulations of plasma inflow into the Earth's magnetosphere¹⁵ it was concluded that the width of the transition layer is first set at the level of $λ≈ρeρi$ and then increases to $λ≈ρi$ because of the development of drift ion-cyclotron instability. The previously mentioned PIC simulations of plasma injection into a “picket fence” trap¹⁴ demonstrated the possibility of forming a transition layer with the electron gyroradius scale (⁠ $λ≈4ρe$⁠); however, due to the axial symmetry of the problem, these simulations could not resolve instabilities with non-zero azimuthal numbers m. In our recent work,¹⁶ the same type of equilibrium with a strong jump in the electric potential and a dominant electron current, which was observed inRef. 14, has been shown to form even if one starts from the Kotelnikov equilibrium¹² with the zero electric field and the ion current localized on a large scale of $8ρi$⁠, but this restructuring of the equilibrium occurred without changing the width of the transition layer. It should be noted, however, that in Ref. 16 we used a greatly reduced mass ratio $mi/me=16$⁠, at which the scale $10ρe$ already became indistinguishable from $ρi$⁠.

The second question is related to the stability of the equilibrium in which the current sheet is localized on the scale $λ<ρi$⁠. Obviously, this question strongly depends on the curvature of the magnetic field lines, which has a different sign in cusp and mirror traps. However, in this work, we will be interested in those specific instabilities that develop in high- $β$ plasma even in a straight magnetic field lines. Here, it is more important for us not to reproduce the entire spectrum of instabilities that occur in real magnetic field configurations, but to find out how their development affects the maximum achievable value of $β$⁠. It has long been known that in a low- $β$ plasma the difference in the velocities of ion and electron azimuthal drifts under the influence of $∇n$ leads to the buildup of flute instabilities both at low frequencies $ω<Ωi$¹⁷ and on the ion-cyclotron frequency harmonics $nΩi$⁠.¹⁸ In a plasma with an electric field directed along the density gradient, similar instabilities develop in the frequency range around the lower hybrid resonance.¹⁹ The reason for the unstable plasma behavior in this case is the difference in the electron and ion $E×B$-drift velocities. In equilibria with $β=1$⁠, similar instabilities of drift waves were observed in a highly inhomogeneous current layer,¹⁶ where the radial gradient is present not only in the plasma density but also in the magnetic field, due to which ion-cyclotron harmonics in this layer overlap with the lower hybrid frequency. The difference in the velocities of ion and electron drifts in crossed fields is also the cause of the instability growth both in weakly ionized plasma, where electrons and ions have different collisional rates (Simon–Hoh instability^20,21) and in the collisionless case, when the ion $E×B$ drift velocity is greatly reduced, since the ion Larmor radius exceeds the characteristic size of the electric field localization (modified Simon–Hoh instability^22–24). In low-pressure plasma, instabilities of this kind were observed during PIC modeling of a Penning discharge.^25,26 In this work we will show how strongly the development of instabilities (both drift and MHD) affects the width of the transition layer at limiting values of $β$⁠.

To study all of these questions, we will simulate the problem of plasma injection into a selected cylindrical region of space with an initially uniform vacuum magnetic field, using a semi-implicit electromagnetic 3D PIC code.²⁷ We will study the axisymmetric plasma injection model because it is the simplest case when the non-gyrotropy of plasma pressure can play a significant role in establishing equilibrium. This knowledge will help us to differ non-gyrotropic effects in straight and curved magnetic field lines in future full-fledged 3D simulations of real confinement schemes (mirror and cusp). Note that taking into account the non-gyrotropy of plasma pressure is of fundamental importance not only for fusion but also for problems of magnetic reconnection and turbulence.^28,29

To simulate the formation of a cylindrical plasma column, infinitely extended along the magnetic field, we use our own semi-implicit 3D PIC code with exact conservation of energy and charge.^27,30 The problem is modeled in Cartesian coordinates. To eliminate as much as possible the influence of longitudinal inhomogeneities, periodic boundary conditions are used in the z-direction, and the length of the region has the minimal size (five cells). The layout of the computational domain in the transverse $(x,y)$ plane is shown in Fig. 2. At the initial moment of time, we set only a uniform vacuum magnetic field $B(0)=(0,0,Bv)$ with $Bv=0.2$ in the whole domain. Then plasma particles are injected into a central region with a radius of $R=15$ at a constant rate of $dn/dt=1/τ$⁠. This means that at each time step a certain number of ion–electron pairs are thrown into this region with random coordinates, so that by the time $τ=6000$ the density of particles of each type should reach the value $n=1$⁠. This injection rate is chosen so that the ions have time to react adiabatically to changes in the magnetic field. In our case, the time interval $τ$ approximately corresponds to two periods of their cyclotron rotation. The injected particles have a Maxwellian velocity distribution with different ion and electron temperatures (⁠ $Ti=10$⁠, $Te=1$ keV). This situation is typical for experiments in mirror traps, where hotter ions appear during neutral injection. If we choose the value $n0=1013$ cm⁻³ as the unit density, then the vacuum magnetic field will correspond to the conditions of the CAT experiment (Compact Axisymmetric Toroid).³¹ To reduce the complexity of calculations, a reduced mass ratio $mi/me=100$ is used. In the peripheral region $r>120$ there is an absorbing layer, which eliminates particles that reach it and absorbs all perturbations of electromagnetic fields.

The time step is chosen to resolve the Larmor rotation of electrons $Δt=1.5 ωpe−1=0.3 Ωe−1$⁠, where $Ωe=eBv/(mec)$ is the electron cyclotron frequency in a vacuum magnetic field. The spatial steps are equal to $dx=dy=dz=0.5 c/ωpe$⁠, which is about 2 electron gyroradii. An acceptable noise level is achieved when the unit density $n=1$ is represented by 2000 computational particles per cell.

Let us present the results of modeling the continuous injection of hot ion plasma into a vacuum. Furthermore, we will say that the magnetic field is displaced from the plasma if its dimensionless magnitude $b=B/Bv$ becomes lower than unity. To completely displace the magnetic field from the injection region, it is necessary to create a plasma with a pressure of $Bv2/2$⁠. At a constant temperature and without plasma expansion, such pressure would be achieved by the time $t≈τ$ at $n=1$⁠. In reality, due to the conservation of the magnetic flux, field displacement is accompanied by the expansion of the flux tube and a corresponding decrease in plasma density. In addition, the temperature also does not remain constant during such expansion, since the plasma spends its energy by doing work to displace the magnetic field. As a result, the complete displacement (⁠ $b→0$⁠) occurs at a higher density [Fig. 3(g)] and over a longer time [in Fig. 3(h) pressure $1−b2$ reaches a close vicinity of 1 at $t=4τ$].

Before the detailed study of the dynamics of this process, as well as the possibility for plasma pressure to exceed the pressure of the vacuum magnetic field, we will first answer the question about the width of the current sheet. As already mentioned in the Introduction, if the radial electric field of charge separation is for some reason absent in the transition layer (e.g., short-circuited on the conducting plates of plasma receivers), then Kotelnikov's kinetic theory¹² predicts that the current necessary for equilibrium is created by a diamagnetic ion drift, and the width of the transition layer should be determined by the ion gyroradius (⁠ $λ≈8ρi$⁠). PIC simulations of $β=1$ equilibrium in a multicusp trap with electron Larmor radius resolution¹⁴ showed that the appearance of a jump in the electric potential $Δφ∼T/e$ in the transition layer contributes to the reflection of ions. As a result, the current is carried by electrons drifting in crossed fields, and its typical size of localization is limited by the electron gyroradius scale (⁠ $λ≈4ρe$⁠). These results contradict other studies,^15,16 in which equilibrium with a strong electric field was formed with a much more extended transition layer (tens of $ρe$⁠). Since the mass ratio in Refs. 15 and 16 was significantly underestimated compared to the real value (⁠ $mi/me≲100$⁠), the observed layer size was interpreted by the authors as $λ∼ρi$⁠. A feature of the simulations carried by Park et al.¹⁴ was the axially symmetric geometry incapable of reproducing drift microinstabilities at harmonics of the ion-cyclotron frequency, which were considered in Refs. 15 and 16 as the main reason for the broadening of the layer to the scale of the ion gyroradius.

To clarify this issue, let us follow the evolution of the current sheet size in our simulations. Figure 3(i) shows that if the profiles of the electric current at different moments of time (up to the moment of complete displacement of the field) are normalized to their maximum values and the positions of these maxima are combined, then they all fit well onto the Gaussian-like curve in which one wing is elongated toward the plasma periphery. From the two-dimensional maps of the current density shown in Figs. 3(a)–3(f), it is clear that the current is indeed mainly created by electrons, and the growth of an asymmetric wing occurs in the region where drift oscillations are excited. In this case, the presence of drift microturbulence apparently does not affect the width of the profile, the upper part of which remains symmetrical and rather smooth [Fig. 3(i)]. If the radial coordinate is expressed in units of the electron Larmor radius, calculated from the vacuum magnetic field and injection temperature [ $ρe=Te/(mec2)/Bv$], then the width of the current layer reaches the value $λFWHM≈32ρe$⁠. Even if the width of this current profile is measured in electron gyroradii calculated using the mean magnetic field in the middle of the layer (⁠ $Bz≈Bv/2$⁠), as is done in the paper,¹⁴ then it will still turn out to be 4 times larger than that observed by Park et al. The observed thickness $λ≈32ρe$ coincides with the gyroradius of ions with the mass $mi=100 me$ [ $ρi=miTi/(meTe)ρe≈31.6ρe$]; however, we will further show that the conclusion made in Refs. 15 and 16, that the layer width is related to the ion gyroradius in the presence of a jump in the electric potential $Δφ∼Ti/e$⁠, is erroneous.

To do this, we will carry out an additional simulation in which the ion Larmor radius of the injected plasma will coincide with the electron gyroradius (⁠ $ρi=ρe$⁠). In this case, isothermal (⁠ $Ti=Te=1$ keV) electron–positron plasma (⁠ $mi=me$⁠) is injected into a uniform vacuum magnetic field of the same magnitude $Bv$⁠. This set of parameters is not dictated by interest in any particular physical problem, but is used here only as an example of a plasma with equal gyroradii. A comparison of the simulation results in the cases of $ρi=32ρe$ and $ρi=ρe$ at times corresponding to the same value of injected energy is shown in Fig. 4. Despite the 30-fold difference in ion gyroradii in these simulations, no significant changes in the size of the transition layer are observed [Fig. 4(g)]. Indeed, if we combine the normalized electric current density profiles [Fig. 4(i)], it becomes clear that the entire effect of a more extended ion gyroradius is reduced to the growth of the asymmetric wing, which does not affect the width of this profile at half maximum. Figures 4(h) and 4(a) show that when the mass and temperature of different types of particles are equal, a jump in the electric potential at the plasma boundary does not form. Under these conditions, the current is equally created by the diamagnetic drift of both ions and electrons, and the current sheet remains stable with respect to drift perturbations [Figs. 4(b) and 4(c)]. In a plasma with heavy ions, charge separation occurs. The electric potential is established at the level $Δφ≈8Ti/e$ [Fig. 4(d) and 4(h)], sufficient to reflect almost all ions, and the current necessary for equilibrium is created exclusively by electrons due to their $E×B$ drift [Fig. 4(f)]. The strong difference in the drift motions of ions and electrons leads to buildup of the same microinstabilities in the range of the ion-cyclotron frequency harmonics [Fig. 4(e)], which have already been studied in detail in our recent work.¹⁶ 

Thus, we can conclude that in $β=1$ equilibria with a strong jump in the electric potential $Δφ∼Ti/e$⁠, the size of the transition layer, in accordance with the work of Park et al.,¹⁴ is determined by electrons, but is limited by tens rather than units of their Larmor radii. The microinstabilities observed in Refs. 15 and 16 at the harmonics of the ion-cyclotron frequency actually expand the region through which the equilibrium current flows, but have almost no effect on the width of its profile at half maximum. Note that in hydrogen plasma the observed width of the transition layer $λ=32ρe$ should be considered as the sub-ion scale (⁠ $λ<ρi$⁠).

Let us see whether Eq. (8) is satisfied in our simulations and what sign $Δβ$ has in the problem of continuous plasma injection. The radial profiles of $βr$⁠, $βϕ$⁠, $Δβ$⁠, and $βb$ averaged over the azimuthal angle at different times are presented in Fig. 5 (bottom row). It can be seen that the pressure of the displaced magnetic field $βb(r)$ (violet curve) coincides well with the plasma pressure $βr(r)−Δβ(r)$ (red curve) accounting for the non-gyrotropic addition, therefore, at the selected injection rate, the equilibrium (8) can be considered established at each moment of time. At the initial stage of injection (⁠ $t<3τ$⁠), when $βr(0)$ has not yet approached unity, $Πrr<Πϕϕ$ over most of the transition layer, therefore the non-gyrotropic addition to pressure is negative (⁠ $Δβ<0$⁠). Approaching the MHD limit, $Δβ$ changes sign, so the pressure in the center of the plasma column $βr(0)=βϕ(0)$ continues to grow above the MHD limit due to the growth of the non-gyrotropic additive. The dynamics of relative pressure growth in the center of the plasma is shown in Fig. 6. It can be seen that, at the moment $t≈5τ$⁠, the plasma pressure on the axis reaches a maximum value of $β=1.15$ and does not change further. From the two-dimensional maps of $Πrr(x,y)$⁠, $Jϕi(x,y)$⁠, and $Jϕ(x,y)$ presented in Fig. 5, one can see that the stabilization of pressure growth coincides in time with the rise of the flute-like MHD instability with an azimuthal number $m=3$⁠. The most likely reason for pressure saturation is the radial removal of the azimuthal momentum during the growth of this unstable mode, which leads to the need to take into account the term $(Πϕr/r)ϕ′$ in the equilibrium equation (1).

Now we will determine the main characteristics of the observed MHD instability and find out the reason for its buildup. From the perturbation of plasma shape in Fig. 5, it is already clear that the dominant role in this instability is played by the mode with the azimuthal number $m=3$⁠, the wavenumber of which has the form $kϕ=m/r$⁠. Figure 7(a) shows the time evolution of the azimuthal electric field $Eϕ(ϕ,t)$ at a fixed radius $r=25$⁠. Using Fourier filtering, we select in this field the contribution of the most unstable mode $m=3$ [Fig. 7(b)]. By the displacement of the wave crest position, we determine the real part of the frequency, which turns out to be significantly lower than the ion-cyclotron frequency in a vacuum field $ω∼0.06 Ωvi$ (the inscribed dashed line has a slope $ω/kϕ$⁠). The growth rate of this instability $Γ∼0.05 Ωvi$ has approximately the same value, which is determined from the stage of the exponential growth of the perturbation amplitude [Fig. 7(c)].

Identification of the observed instability is complicated by the fact that the considered equilibrium is not stationary and changes noticeably during its development. For example, at the moment $t=6τ$ the magnetic field reverses, so that even the electrons become unmagnetized inside the injection region. Until this moment, instability develops under conditions when electrons are magnetized, but ions are not (they only reflect from the plasma boundary), and an azimuthal current flows inside the injection region (⁠ $r<10$⁠), to which both electrons and ions make comparable contributions [Figs. 7(d) and 7(e)]. At the stage of nonlinear saturation of the instability (⁠ $t≈6τ$⁠), this current disappears, from which we can conclude that the cause of the instability is the relative rotation of the electron and ion components of the plasma inside the injection region. In contrast to the modified Simon–Hoh instability, the drift difference here is associated not with the different influence of the radial electric field on the magnetized and non-magnetized plasma components (⁠ $Er=0$ inside the injection region), but exclusively with the ion and electron pressure gradients. We also note that the growth of the $m=3$ mode begins at the moment when the magnetic field in the injection region becomes so low that the local ion-cyclotron frequency $Ωi$ there begins to coincide with the drift frequency $kϕVdi$ calculated from the directional azimuthal velocity of ions $Vdi=Jϕi/ni$⁠. It is also clear from Fig. 7(f) that the deepening of the magnetic well is accompanied (with some delay) by a decrease in the velocity of the ions, so that the resonance condition $kϕVdi∼Ωi$⁠, taking into account the non-stationary character of equilibrium, is satisfied for the most unstable mode for a long time. The velocity of plasma rotation established inside the injection region at $t>7τ$ [see Fig. 7(e)] is much smaller than the thermal ion velocity, which allows us to consider $Πrr$ and $Πϕϕ$ as plasma pressure.

PIC modeling of the continuous injection of particles into a cylindrical plasma column confined by a uniform vacuum magnetic field showed that, due to the formation of an electric potential jump at the plasma boundary, the equilibrium current is mainly created by electrons, and the characteristic width of its profile $λFWHM∼32 ρe$ should be regarded as an electron scale, which is reproduced even in plasma with $ρi=ρe$⁠. Drift instabilities at harmonics of the ion-cyclotron frequency (⁠ $ω∼nΩvi$⁠) developing in the case of $mi≠me$ lead to the broadening of the region of current flow only at the rarefied plasma periphery having almost no effect on the width of this current profile at half maximum. The formation of such a thin, sub-ion transition layer leads to a distortion of the circular rotation of particles and produces a non-gyrotropic addition to pressure, which makes it possible to exceed the MHD limit $β=1$ by 15%. The growth of $β$ above the MHD limit is saturated due to the development of a low-frequency MHD flute-like instability with an azimuthal number $m=3$⁠. This unstable mode runs in azimuth at the speed of ions inside the injection region and resonates with the local ion-cyclotron frequency at a sufficiently large depth of the magnetic well. As a development of this work, it is planned to carry out three-dimensional simulations of mirror and cusp configurations of the vacuum magnetic field. These studies will allow us to answer the questions of whether the constraint on further growth of $β$ will be relaxed in fields with a favorable curvature of magnetic field lines due to the stabilization of the observed MHD instability and whether it is possible to create the potential jumps necessary to obtain a thin transition layer in a mirror trap by applying the potential difference to the isolated sectioned plates of the plasma receiver. The future 3D simulations accounting for realistic boundary conditions will make it possible to compare simulation results with data from forthcoming experiments at the CAT facility.

This work was supported by the Russian Science Foundation (Grant No. 21-72-10071).

The authors have no conflicts to disclose.

Igor Timofeev: Conceptualization (lead); Formal analysis (lead); Investigation (equal); Supervision (equal); Writing – original draft (equal). Vladislav Kurshakov: Investigation (equal); Visualization (lead); Writing – original draft (equal). Evgeny Berendeev: Investigation (equal); Software (lead); Validation (lead).

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