Hostname: page-component-669899f699-g7b4s Total loading time: 0 Render date: 2025-04-25T19:17:29.212Z Has data issue: false hasContentIssue false

On the stability of the m = 1 rigid ballooning mode in a mirror trap with high-beta sloshing ions

Published online by Cambridge University Press:  14 April 2025

Igor A. Kotelnikov*
Affiliation:
Novosibirsk State University, Novosibirsk, Russia Budker Institute of Nuclear Physics, Novosibirsk, Russia
*
Email address for correspondence: [email protected]

Abstract

Stability of the ‘rigid’ ($m = 1$) ballooning mode in a mirror axisymmetric trap is studied for the case of oblique neutral beam injection (NBI), which creates an anisotropic population of fast sloshing ions. Since small-scale modes with azimuthal numbers $m>1$ in long thin (paraxial) mirror traps are easily stabilized by finite-Larmor-radius (FLR) effects, suppression of the rigid ballooning and flute modes would mean stabilization of all magnetohydrodynamic (MHD) modes, with the exception of the mirror and firehose disturbances, which are intensively studied in geophysics, but have not yet been identified in mirror traps. Large-scale ballooning mode can, in principle, be suppressed either by the lateral perfectly conducting wall, or by the end MHD anchors such as the cusp, by biased limiters or by a combination of these two methods. The effects of the wall shape, vacuum gap width between the plasma column and the lateral wall, angle of oblique NBI, radial profile of the plasma pressure and axial profile of the vacuum magnetic field are studied. It is confirmed that the lateral conducting wall still creates the upper stability zone, where the ratio $\beta$ of the plasma pressure to the pressure of vacuum magnetic field exceeds the second critical value $\beta_{\text{cr2}}$, $\beta >\beta _{\text {cr2}}$. However, in many cases the upper zone is clamped from above by mirror instability. When the lateral wall is combined with end MHD anchors, a lower stability zone $\beta <\beta _{\text {cr}1}$ appears, where $\beta$ is below the first critical value $\beta_{\text{cr1}}$. These two zones can overlap in the case of a sufficiently smooth radial pressure profile, and/or a sufficiently low mirror ratio and/or a sufficiently narrow vacuum gap between the plasma column and the lateral wall. However, even in this case, the range of permissible values of beta is limited from above by the threshold of mirror instability $\beta_{\text{mm}}$, so that $\beta <\beta _{\text {mm}}<1$, in contrast to the case of transversal NBI, when neutral beams are injected perpendicularly to the magnetic field.

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

1. Introduction

In a series of articles, the stability of a so-called rigid ballooning mode with azimuthal number $m=1$ was studied for several model configurations of an axisymmetric mirror trap (also called open or linear traps) using a numerical solution of the LoDestro equation (LoDestro Reference LoDestro1986). The conditions for stabilization of these modes were found assuming both a lateral perfectly conducting cylindrical wall surrounding the plasma column and traditional magnetohydrodynamic (MHD) stabilizers attached to the central cell of such a trap. A few model radial and axial plasma pressure profiles and model axial vacuum magnetic field profiles we examined.

In quite realistic experimental conditions, oscillations of a ballooning perturbation with an azimuthal number $m\geqslant 2$ are easily stabilized by finite-Larmor-radius (FLR) effects, which, however, cannot stabilize the $m=1$ mode, although they can make it ‘rigid’ in a certain sense. Thus, stabilization of the rigid ballooning mode $m=1$ is, in practice, equivalent to the stabilization of all types of ballooning and flute disturbances.

In the first article (Kotelnikov et al. Reference Kotelnikov, Zeng, Prikhodko, Yakovlev, Zhang, Chen and Yu2022b), which was followed by the two more papers (Kotelnikov, Prikhodko & Yakovlev Reference Kotelnikov, Prikhodko and Yakovlev2023; Zeng & Kotelnikov Reference Zeng and Kotelnikov2024), a plasma model with isotropic pressure was implemented. It was assumed that the plasma column is surrounded by a cylindrical perfectly conducting chamber with a variable radius, which is proportional to the plasma radius with a constant coefficient. This shape was later called ‘proportional’ because it reproduced the shape of the plasma column on a larger scale.

The numerical code was then upgraded to implement a transverse neutral beam injection (transverse NBI) model, which simulates the anisotropic fast ion pressure that occurs when a beam of fast neutral atoms is injected into a relatively cold target plasma at right angles to the axis of trap at a magnetic field minimum (Kotelnikov et al. Reference Kotelnikov, Prikhodko and Yakovlev2023).

The next update to the numeric code added the ability to select a custom conductive chamber shape. This possibility was demonstrated assuming an example of a conducting chamber in the shape of a straight cylinder. Calculations published in Zeng & Kotelnikov (Reference Zeng and Kotelnikov2024) made it possible to compare the stabilizing effect of a proportional conducting chamber with the effect of such a ‘straightened’ chamber for a model of the transverse NBI.

The present paper reports the results of study for the case on oblique NBI. Oblique NBI yields a population of fast ions, which are called sloshing ions. They are characterized by a pressure distribution with a peak located somewhere between the median plane of the mirror trap and magnetic mirror plugs. It is significant that, in addition to the ballooning instability, mirror and firehose-type instabilities can also be raised under oblique injection. It will be shown below that, in the high mirror ratio trap, it is these instabilities that determine the ultimate pressure limit. All such limits are expressed through the dimensionless parameter $\beta$ (beta), which is defined by (3.8) as the ratio of the transverse plasma pressure at the point of intersection of the trap axis with the midplane to the vacuum field pressure in the same midplane.

The current article completes the series of three above cited publications that used model functions to approximate radial and axial profiles of the plasma pressure and vacuum magnetic fields. They were convenient for identifying general patterns in how the ballooning instability threshold and margin depend on the mirror ratio, width of the magnetic mirror plugs, shape and steepness of the radial and axial plasma pressure profile, size and shape of the vacuum gap between the lateral conducting wall and the lateral surface of the plasma column, as well as upon the stability margin created due to the end MHD anchors such as magnetic cusps. In addition, the step-by-step complication of model functions made it possible to break down the creation of a numerical code into stages, each of which gave a new meaningful result.

Each stage began with a code modernization, and each modernization was accompanied by a complication of the numerical code and an increase in the duration of calculations. If for an isotropic plasma the recalculation of the data necessary to reproduce all the graphs published so far took (after all subsequent improvements to implemented algorithms) no more than an hour on a desktop computer with a typical number of cores, then in the case of a transverse NBI and a straightened camera, this would take few weeks. The oblique NBI calculations presented in this paper initially took few months of intensive work and only the transfer of calculations to a multi-core cluster made it possible to return the duration of calculations to reasonable limits.

The reason for the slowdown in the numerical code is trivial: with each update, the number of integrations that could be performed analytically inevitably decreased.

In the current version, the developed numerical code consists of a set of executable modules written in the Wolfram Language© and combined into the PEK package, so named in memory of my grandfather, Pavel E. Kotelnikov. As a result of recent modernization, PEK has acquired the ability to use interpolated experimental profiles of plasma pressure, magnetic field and lateral wall shape in addition to the model profiles that have been used up to the present time.

To reduce inevitable repetition, the traditional review of the contents of many articles on the topic of interest is omitted and the reader is referred to the publications cited above. In addition to these, it is worth mentioning the work of Kesner et al. (Kesner Reference Kesner1985; Li, Kesner & Lane Reference Li, Kesner and Lane1985, Reference Li, Kesner and Lane1987a; Li, Kesner & LoDestro Reference Li, Kesner and LoDestro1987b) with which the results of this work will be compared. In particular, it is worth noting the statement of Li, Kesner and LoDestro in Li et al. (Reference Li, Kesner and LoDestro1987b). They wrote that, in a proportional chamber, the walls of which are as close as possible to the surface of the plasma column, anisotropic plasma is always stable near the threshold of mirror instability. Our calculations only partially confirm this statement.

Also, in order to shorten the introductory part of the article, the formulation of the LoDestro equation and boundary conditions has been moved to Appendix A. It has been supplemented by recent findings, which have made it possible to calculate some more integrals analytically. To understand the main part of the article, it is enough to explain that it deals with the Sturm–Liouville problem for the LoDestro equation in the formulation described in detail in §§ 2 and 5 of Kotelnikov et al. (Reference Kotelnikov, Prikhodko and Yakovlev2023).

Section 3 describes the anisotropic pressure distribution model in a mirror trap under both transverse and oblique NBI. The beta limits imposed by the thresholds of mirror and firehose instabilities are also calculated. Sections 4, 5 and 6 step-by-step present and discuss the results of calculations of the critical beta values $\beta _{\text {cr1}}$ and $\beta _{\text {cr2}}$, first for the case when there are no other means of MHD stabilization, except for the lateral conducting wall, next for the case when only end MHD anchors are mounted and finally for the case when infinitely strong MHD anchor stabilizers are used together in combination with the lateral wall. In addition, § 6 explains how the PEK package simulates MHD anchors with different stability margins.

2. Magnetic field model

In the paraxial approximation, the magnetic flux $\psi$ is related to the distance $r$ from the axis $z$ of an axially symmetric mirror trap by the equation

(2.1)\begin{equation} \frac{r^{2}}{2} = \int_{0}^{\psi} \frac{\mathrm{d}{\psi}}{B} . \end{equation}

The paraxial approximation means that the radius of curvature of the magnetic field lines significantly exceeds both the plasma radius and the distance between the magnetic mirrors. The equivalent condition for the applicability of such an approximation is formulated as the smallness of the plasma column radius compared with the distance between the magnetic mirrors.

The true magnetic field $B$ weakened by the diamagnetic effect is related to the vacuum magnetic field $B_{v}$ and transverse pressure $p_{\bot}$ by the transverse equilibrium equation (where the rationalized electromagnetic units are used, they are also known as Heaviside–Lorentz units)

(2.2)\begin{equation} B^{2} + 2p_{\mathrel\perp} = B_{v}^{2} . \end{equation}

In this form, it is approximately true in the paraxial approximation, when the radius of the plasma column $a=a(z)$ is everywhere small compared with the distance between the magnetic mirrors.

In what follows, dimensionless notation is used, such that $b=B/B_{R}$, $b_{v}=B_{v}/B_{R}$, with $B_{R}$ being the magnetic field at the turning point of fast ions, where the plasma transverse pressure drops to zero, $p_{\perp }=p_{\|}=0$, as explained in the next section. The radial coordinate $r$ is normalized in such a way that $\psi =1$ at the plasma lateral boundary $r=a(z)$. The three-parameter function

(2.3)\begin{equation} b_{v}(z)= \left[ 1 + (M-1)\sin^{q}({\rm \pi} z/2) \right]/R , \end{equation}

earlier used by Kotelnikov et al. (Reference Kotelnikov, Zeng, Prikhodko, Yakovlev, Zhang, Chen and Yu2022b, Reference Kotelnikov, Prikhodko and Yakovlev2023) and Zeng & Kotelnikov (Reference Zeng and Kotelnikov2024), approximates the axial profile of the vacuum magnetic field and which we used earlier. It assumes normalization of both the true magnetic field $B$ and the vacuum magnetic field $B_{v}$ by their common value at the tuning point $B_{R}$ of fast ions. In addition, we assume that the magnetic plugs with mirror ratio M are located in the $z=\pm 1$ planes, which means that longitudinal coordinate $z$ is normalized by the distance between the midplane and the mirror plug throat. The repeated choice of the function $b_{v}(z)=B_{v}(z)/B_{R}$ in the form (2.3) is due to the desire to simplify the comparison of new results with previous publications.

Parameter $q$ in (2.3) determines the steepness and width of the magnetic mirrors: larger $q$ corresponds to a smaller fraction of the mirror plug in the total length of the trap, as shown in figure 1.

Figure 1. Axial profile of the vacuum magnetic field (2.3) for mirror ratio $M=32$ and three values of the index $q$ indicated on the graphs as compared with the magnetic field in the gas-dynamic trap (GDT, Bagryanskij et al. Reference Bagryanskij, Ivanov, Klesov, Kotel'nikov, Krasnikov, Roslyakov, Tsidulko, Breun, Molvik and Casper1990) and the Wisconsin HTS axisymmetric mirror (WHAM, Endrizzi et al. Reference Endrizzi, Anderson, Brown, Egedal, Geiger, Harvey, Ialovega, Kirch, Peterson and Petrov2023) devices.

Parameter $R=B_{R}/\min (B_{v})$ has the meaning of the mirror ratio between the turning point $B=B_{v}=B_{R}$, where the fast ion pressure drops to zero, and the minimum value $\min (B_{v})$ of the magnetic vacuum field in the middle plane of the trap.

Parameter $M=\max (B_{v})/\min (B_{v})=B_{v}(\pm 1)/B_{v}(0)$ is equal to the mirror ratio in the traditional sense.

3. Anisotropic pressure models

In publications on the stability of the rigid ballooning mode, two analytic models of anisotropic pressure have previously been proposed. Kesner in his paper (Kesner Reference Kesner1985) indicates that in the first model the transverse pressure $p_{\mathrel \perp }$ in a non-uniform magnetic field $b\leq 1$ varies according to the law

(3.1)\begin{equation} p_{\mathrel\perp} = p(\psi) \left( 1-b^{2} \right), \end{equation}

while in the second model

(3.2)\begin{equation} p_{\mathrel\perp} = p(\psi)\, b^{2}\left( 1-b \right)^{n-1} , \end{equation}

where $p(\psi )$ is a function of magnetic flux $\psi$, $b=1$ is the normalized magnetic field at a ‘turning point’ where the fast ion pressure drops to zero and $n\geq 2$ is an integer index. Both models assume that, in the region $b>1$, the plasma pressure is negligible, $p_{\perp }=0$, but there is a relatively cold plasma with low pressure that extends up to the magnetic field peak $b_{\max }=\max (B_{v})/\min (B_{v})=M/R$ in a magnetic plug. Cold plasma provides electrical contact with perfectly conducting end plates which imitate the end MHD anchors.

The first model (3.1) was previously applied to the case of transverse injection of neutral atom beams (transverse NBI) by Kotelnikov et al. (Reference Kotelnikov, Prikhodko and Yakovlev2023). It describes the pressure distribution in an anisotropic plasma with a pressure peak near the magnetic field minimum at the midplane of the central cell of a mirror trap. The second model (3.2) is suitable for the case of an oblique NBI, when the pressure peak is located between the middle plane of the trap and the magnetic mirror. It is this model with indices $n=2$ and $n=3$ that is used in this article.

Kinetic theory proves that if the transverse pressure $p_{\mathrel \perp }$ is given as a function of $b$, then the longitudinal pressure $p_{\|}$ is uniquely determined using the parallel equilibrium equation with a properly posed boundary condition (Newcomb Reference Newcomb1981). The latter can be rewritten in terms of the partial derivative with respect to $b$ for a constant magnetic flux $\psi$ as

(3.3)\begin{equation} p_{\mathrel\perp} ={-} b^{2}\frac{\partial }{\partial b}\frac{p_{\|}}{b} . \end{equation}

It would be wrong therefor to choose at random a couple of functions for the transverse and longitudinal pressures, since they must correspond to a distribution function, which is a solution to the kinetic equation. In § 3.4 the distribution functions that lead to (3.1) and (3.2) are restored under some dedicated assumptions.

Another key result of the kinetic theory is the assertion that the function $p_{\mathrel \perp }/b^{2}$ always decreases as $B$ increases, i.e.

(3.4)\begin{equation} \frac{\partial }{\partial b}\frac{p_{\mathrel\perp}}{b^{2}}\leqslant0. \end{equation}

Obviously, both models of anisotropic pressure satisfy the condition (3.4). Note, however, that the inequality (3.4) was proved under the assumption that the class of distribution of fast ions $F(\varepsilon,\mu )$ as a function of energy $\varepsilon$ and magnetic moment $\mu$, is narrowed by the condition

(3.5)\begin{equation} \frac{\partial F}{\partial \varepsilon }\leq 0, \end{equation}

for all $\varepsilon$ and $\mu$, so that $F$ is everywhere monotone decreasing in $\varepsilon$. It can be shown, in fact, that monotonicity in $\varepsilon$ is required as a condition for local variational stability in the sense of the Grad energy criterion (Grad Reference Grad1967). It also underlies the derivation of the Kruskal–Oberman energy principle (Kruskal & Oberman Reference Kruskal and Oberman1958) and is considered by some authors as a necessary requirement for the correct formulation of the energy principle (Stupakov Reference Stupakov1987). Refusal of the condition (3.5) leads, in particular, to the appearance of new instabilities associated with resonances between longitudinal vibrations of ions and the wave (Skovorodin, Zaytsev & Beklemishev Reference Skovorodin, Zaytsev and Beklemishev2013).

It is known that anisotropic plasma can be subject to the mirror (sometimes called the diamagnetic) and firehose instabilities (sometimes called the garden hose instability). According to the fluid Chew–Goldberger–Low theory, the stability of firehose-type disturbances requires the inequality

(3.6)\begin{equation} p_{\|} \leqslant p_{\mathrel\perp}+b^{2}, \end{equation}

to hold for any $b$ (Kalsrud Reference Kalsrud1983; Ilgisonis Reference Ilgisonis1993a). Stability against the mirror mode implies that (Rudakov & Sagdeev Reference Rudakov and Sagdeev1961; Thompson Reference Thompson1964; Newcomb Reference Newcomb1981; Ilgisonis Reference Ilgisonis1993b; Southwood & Kivelson Reference Southwood and Kivelson1993)

(3.7)\begin{equation} \frac{\partial }{\partial b}\left( p_{\mathrel\perp}+\frac{b^{2}}{2}\right)>0. \end{equation}

Simple, intuitively watertight, ways to derive criteria (3.6) and (3.7) are collected in the student problems after lecture 25 in the author's 2nd volume of ‘Lectures on plasma physics’ (Kotelnikov Reference Kotelnikov2021). All specified criteria of stability are satisfied for the pressure model (3.1), which, we recall, corresponds to the case of transverse NBI. For the oblique NBI model (3.2), these conditions impose restrictions on the value of beta, which are specified below.

The firehose instability is essentially a modification of the Alfvén wave in an anisotropic plasma. In a homogeneous plasma, its increment is proportional to the longitudinal wavenumber $k_{\|}$. In the low-frequency limit, firehose perturbations are strongly elongated along magnetic field lines and represent bending or torsional vibrations of magnetic flux tubes. Inhomogeneity of the magnetic field changes the threshold for the firehose instability (Mirnov Reference Mirnov1986).

In a homogeneous magnetic field, the increment of mirror instability is also proportional to the longitudinal wavenumber $k_{\|}$ and, oddly enough, inversely proportional to the number of resonant ions with zero longitudinal velocity (Southwood & Kivelson Reference Southwood and Kivelson1993; Pokhotelov et al. Reference Pokhotelov, Treumann, Sagdeev, Balikhin, Onishchenko, Pavlenko and Sandberg2002). Near the threshold, mirror disturbances, like the flute ones, are strongly elongated along magnetic field lines.

The criteria (3.6) and (3.7) should be interpreted as formal and reference. For the pressure models described below, they first break down locally near two spatially separated points on the axis of the plasma column. It has not yet been studied whether the lateral conductive wall can influence the behaviour of the localized mirror and firehose disturbances.

In PEK's internal classification, the isotropic plasma variant is designated ‘A0’ and the anisotropic pressure (3.1), which is formed under transverse NBI, is designated ‘A1’. The oblique NBI modelled by (3.2) with an arbitrary index $n\geq 2$ is denoted ‘An.’ The two cases $n=2$ and $n=3$ discussed in detail in this article are abbreviated ‘A2’ and ‘A3.’

Method of the MHD stabilization is specified by an additional abbreviation. Variants with conducting lateral wall stabilization are marked with the shorthand ‘Lw’ (for lateral wall), and variants with joint lateral wall and conducting end plates installed in the throat of the magnetic plug are marked with the shorthand ‘Cw’ (for combined wall). The design of the lateral conducting wall in the form of a proportional chamber is labelled ‘Pr’ (proportional). Thus, the label ‘A3-LwPr’ in a figure caption should be deciphered as the variant of plasma with anisotropic pressure of the ‘A3’ type in a proportional conducting chamber without the end MHD anchor. Continuing the line of Zeng & Kotelnikov (Reference Zeng and Kotelnikov2024), this paper compares, when it is suitable, the stabilizing effect of a proportional chamber with a straightened chamber, which is assigned the abbreviation ‘St’ (straightened).

In what follows, the parameter beta is defined as the maximum of the ratio $2p_{\mathrel \perp }/b_{v}^{2}$

(3.8)\begin{equation} \beta=\max(2p_{\mathrel\perp}/b_{v}^{2}) . \end{equation}

In case the radial pressure profile is peaked at the trap axis, and the inequality (3.4) holds, the maximum is reached on the trap axis (where $\psi =0$) at the vacuum field minimum (where $\min (B_{v})=1/R$). Combining definition (3.8) and the dimensionless version of (2.2) readily yields true magnetic field at the minimum

(3.9)\begin{equation} b_{\min} = \sqrt{1-\beta }/R . \end{equation}

As shown in Kotelnikov et al. (Reference Kotelnikov, Zeng, Prikhodko, Yakovlev, Zhang, Chen and Yu2022b, Reference Kotelnikov, Prikhodko and Yakovlev2023), there can exist from none to two stable zones on the stability maps depending on the availability of lateral wall and MHD anchors. Using notations introduced in Zeng & Kotelnikov (Reference Zeng and Kotelnikov2024), one can specify than the Lw configuration can provide an upper stability zone at a sufficiently large beta that exceeds the second critical value, $\beta > \beta _{\text {cr}2}$. In the Cw, Bw (blind wall) and Rw (ring wall) configurations, a lower zone $\beta < \beta _{\text {cr}1}$ can appear in addition to the upper one. These two zones can merge when the vacuum gap between the plasma column and the inner surface of the conducting shell is sufficiently narrow.

3.1. Pressure model A1

The A1 pressure model is described in detail in Zeng & Kotelnikov (Reference Zeng and Kotelnikov2024). It is intended for modelling the transverse NBI. To facilitate comparison with other pressure models, the main properties of the A1 model are briefly summarized below, starting with the definition

(3.10)\begin{equation} \left. \begin{aligned} p_{\mathrel\perp} & = p_{0} f_k(\psi )\, \left(1-b^{2}\right),\\ p_{\|} & = p_{0} f_k(\psi )\, \left(1-b\right)^2. \end{aligned} \right\} \end{equation}

The dimensionless function $f_{k}(\psi )$ defines the radial pressure profile and is normalized so that $f_{k}(0)=1$. It was chosen in Kotelnikov et al. (Reference Kotelnikov, Zeng, Prikhodko, Yakovlev, Zhang, Chen and Yu2022b) so that the integrals over $\psi$ in the coefficients of the LoDestro equation could be calculated analytically, at least for the case of an isotropic plasma. For integer values of index $k$

(3.11)\begin{equation} f_{k}(\psi) = \begin{cases} 1 - \psi^{k}, & \text{if } 0\leq \psi \leq 1, \\ 0, & \text{otherwise}, \end{cases} \end{equation}

and for $k=\infty$

(3.12)\begin{equation} f_{\infty}(\psi) = H(1-\psi) , \end{equation}

where $H(x)=0$ for $x<0$ and $H(x)=1$ for $x>0$. In the case of oblique NBI, the integrals over $\psi$ in the coefficients of the LoDestro equations can be calculated only numerically, and therefore there is no need to restrict the class of functions $f_{k}(\psi )$, except for the comparison of the new and earlier results. Function $f_{1}$ describes the most smooth radial pressure profile. Maximal index $k=\infty$ corresponds to a stepwise profile with a sharp boundary.

The equation of transverse equilibrium (2.2) is readily resolved regarding the true magnetic field as

(3.13)\begin{equation} b(z,\psi) = \sqrt{ \frac{b_{v}^{2}(z)-2p_{0}f_{k}(\psi)}{1-2p_{0}f_{k}(\psi)}}. \end{equation}

The plasma's beta defined by (3.8) is related to parameter $p_{0}$ in (3.10) by

(3.14)\begin{equation} \beta = \frac{2 \left(R^2-1\right)p_{0}}{1 - 2p_{0}} . \end{equation}

Parameter $p_{0}$ can vary within the range

(3.15)\begin{equation} 0< p_{0}<1/2R^{2}, \end{equation}

and $\beta \to 1$ for $p_{0}\to 1/2R^{2}$. For $\beta >1$, transverse equilibrium near the median plane of the trap is impossible, at least in the paraxial approximation and within the MHD approach.Footnote 1

Figure 2(a) shows the axial profile of the transverse pressure for the model A1 in magnetic field (2.3). The pressure peak in this model is located at the midplane of the trap and narrows as the local mirror ratio $R$ at the turning point of fast ions decreases.

Figure 2. Axial profiles of the transverse pressure in models A1 (a), A2 (b) and A3 (c) described respectively by (3.10), (3.18) and (3.28) in magnetic field (2.3) at the axis of the trap for mirror ratio $M=8$, index $q=4$ and various mirror ratios $R$ at the turning points where the pressure of hot plasma component drops to zero. The value $\beta =0.3$ for this figure is chosen so as not to exceed the threshold for excitation of the mirror and firehose instabilities for all values of the parameter $R$ indicated in the figure legend.

It seems intuitively obvious that the degree of anisotropy also increases as $R$ decreases. However, it is not so easy to give a completely satisfactory definition of the degree of anisotropy suitable for all pressure models to be considered in this paper. If the degree of anisotropy would be related to the pressure in the minimum magnetic field as

(3.16)\begin{equation} A_{1} = \frac{p_{\mathrel\perp}-p_{\|}}{p_{\mathrel\perp}+p_{\|}} , \end{equation}

then

(3.17)\begin{equation} A_{1} = {\sqrt{1-\beta}}/{R}, \end{equation}

which it is indeed inversely proportional to the parameter $R$.

It can be proved that the A1 pressure model is stable against the mirror and firehose modes. For example, putting the first of (3.10) into (3.7) yields the inequality $b\,(1-2p_{0})>0$, which is incompatible with condition $0< p_{0}<1/2R^{2}$ within the range $0< b<1$.

3.2. Pressure model A2

The model (3.2) with index $n=2$ reduces to

(3.18)\begin{equation} \left. \begin{aligned} p_{\mathrel\perp} & = p_{0} f_k(\psi )\, b^2\left(1-b\right),\\ p_{\|} & = \frac{1}{2} p_{0} f_k(\psi )\, b\left(1-b\right)^2, \end{aligned} \right\} \end{equation}

in the region $b<1$.

In case of pressure model A1 (3.1), (2.2) reduces to a second-order polynomial equation with respect to $b$. For the model A2 the same equation is of the third order, and the model A3 at $n=3$ yields a fourth-order algebraic equation. As is well known, the solution to such equations can always be expressed in terms of radicals, which saves us from the need to solve these expressions numerically. However, the resulting formulas for a polynomial of the third (and even more so, fourth) order are extremely cumbersome. Therefore, the result of solving (2.2) is written below in a more compact form, which is introduced by Wolfram Mathematica ©

(3.19)\begin{equation} b = \operatorname{Root}\left[ 2 {\#}^3 p+{\#}^2 ({-}2 p-1) +b_v^2 \ \& , \ 2 \right] . \end{equation}

Literally, it means the second root of the cubic equation

(3.20)\begin{equation} 2 {b}^3 p+{b}^2 ({-}2 p-1) +b_v^2 =0 \end{equation}

represented as a ‘pure function’ by the ampersand $\& $ in the first argument of the $\operatorname {Root}$ built-in utility. Parameter $p$ in this formula stands for $p_{0}f_{k}(\psi )$. The $\operatorname {Root}$ utility arranges the roots of the polynomials in an order known only to it, but in such a way that the real roots come first.

Relative plasma pressure, i.e. the parameter beta, is defined by (3.8) as the maximum of the ratio $2p_{\mathrel \perp }/b_{v}^{2}$. According to (3.4), this ratio is a decreasing function of $b$ so that the maximum is reached on the axis of the trap (where $\psi =0$) at the minimum of the vacuum field (where $b_{v}=1/R$). Hence

(3.21)\begin{equation} p_{0} = \frac{\beta}{2 (1-\beta) \left(1-\sqrt{1-\beta }/R\right)} . \end{equation}

Inversion of (3.21) yields

(3.22)\begin{align} \beta & = \operatorname{Root} \left[ 4 {\#}^3 p_{0}^2 + {\#}^2 \left( 4 p_{0}^2R^2-12 p_{0}^2+4 p_{0} R^2 + R^2 \right) \right.\nonumber\\ & \quad \left. -\,{\#} \left(8 p_{0}^2 R^2-12 p_{0}^2+4 p_{0} R^2\right) - 4 p_{0}^2 R^2+4 p_{0}^2 \& ,\ 2 \right] . \end{align}

As follows from (3.21), parameter $p_{0}$ tends to infinity as $\beta \to 1$. In fact, the mirror mode stability condition (3.7) imposes a more stringent condition

(3.23)\begin{equation} 0< p_{0}<1. \end{equation}

If $p_{0}>1$, there is no continuous solution to (2.2). This can be verified by plotting the left-hand side of the equation as a function of $b$. For $p_{0}>1$, the maximum of this function on the interval $0< b<1$ exceeds $1$, while the left-hand side of the equation is less than $1$. This means that paraxial equilibrium is impossible above the mirror instability threshold. Near this threshold,  unexpected phenomena such as equilibrium hysteresis appear. Plasma can reside in two different states, between which a transition is possible by relieving excess pressure (Kotelnikov Reference Kotelnikov2011). Signs of such a transition were found in experiments at the GDT facility (Kotelnikov, Bagryansky & Prikhodko Reference Kotelnikov, Bagryansky and Prikhodko2010). Similar structures are observed in the solar wind, see e.g. Winterhalter et al. (Reference Winterhalter, Neugebauer, Goldstein, Smith, Bame and Balogh1994), Pantellini (Reference Pantellini1998) and Jiang et al. (Reference Jiang, Verscharen, Li, Wang and Klein2022).

In practice, it turned out that the upper limit of the interval (3.23) is sometimes inaccessible to the shooting method when solving the LoDestro equation (A2). The solution of this equation with initial conditions (A13) and (A14) on the left boundary $z=0$ of the domain of definition of the Sturm–Liouville problem (see Appendix A) sometimes went to infinity, before reaching the right boundary $z=z_{E}$. Most often this happened if the parameter $p_{0}$ was close enough to $p_{0}=1$.

The threshold value of parameter $p_{0}$, above which the mirror instability is excited, will be further written as a function

(3.24)\begin{equation} p_{\text{mm}}(R) = 1, \end{equation}

where the subscript ‘mm’ stands for ‘mirror mode’. The beta limit is obtained from (3.22) by substituting $p_{0}=1$

(3.25)\begin{equation} \beta_{\text{mm}}(R) = \operatorname{Root} \left[ 4 {\#}^3 + {\#}^2 \left(9 R^2-12\right) + {\#} \left(12-12 R^2\right)+4 R^2-4 \ \& ,\ 2 \right] . \end{equation}

The firehose instability threshold is described by more cumbersome formulas

(3.26)\begin{gather} p_{\text{fh}}(R) = \operatorname{Root} \left[ {\#}^3 \left( 4 R^4 + 4 R^2\right) + {\#}^2 \left(R^4-30 R^2-27\right) - 24 {\#} R^2 -4 R^2 \ \& , \ j \right], \end{gather}
(3.27)\begin{gather}\beta_{\text{fh}}(R) = \operatorname{Root} \left[ {\#}^3 + {\#}^2\left(R^2-9\right) + 24 {\#} - 16 \ \& , \ j \right], \end{gather}

where $j=1$ if $1< R<\sqrt {9+6 \sqrt {3}}$ and $j=3$ if $R>\sqrt {9+6 \sqrt {3}}$; the subscript ‘fh’ is shorthand for ‘firehose’ mode.

The stability zone bounded by the conditions $\beta <\beta _{\textrm {mm}}(R)$ and $\beta <\beta _{\textrm {fh}}(R)$ is shown in figure 3(a). Its width shrinks to zero both in the $R\to 1$ limit and in the $R\to \infty$ limit. The maximum width corresponds to the point of intersection of the curves $\beta =\beta _{\textrm {mm}}(R)$ and $\beta =\beta _{\textrm {fh}}(R)$, where $\beta = 0.6202$, $R=3.358$. At the same time, it is worth noting that, at the maximum, these functions are equal to $\beta _{\textrm {mm}}(\infty )=2/3$ and $\beta _{\textrm {fh}}(1)=0.9126$.

Figure 3. Threshold of mirror and firehose instabilities in an anisotropic plasma with oblique NBI within the framework of models A2 (a) and A3 (b). Unstable areas are shaded. The same shading scheme is used below in the stability maps without further reminder.

Limitation of the stability zone in an anisotropic plasma by the threshold of mirror instability was taken into account earlier by Kesner et al. (Kesner Reference Kesner1985; Li et al. Reference Li, Kesner and LoDestro1987b), but the firehose instability did not attract the attention of these authors.

Exceeding the threshold of mirror instability within the framework of the applicability of the LoDestro equation is impossible, since this is equivalent to a violation of the equilibrium configuration of the plasma in the paraxial approximation, which is indirectly confirmed by experiments in a gas-dynamic trap (Kotelnikov et al. Reference Kotelnikov, Bagryansky and Prikhodko2010). The PEK code refused to continue computing if $p_{0}>p_{\textrm {mm}}$. However, with a few tweaks, it works above the firehose instability threshold. In some cases, it finds a solution in the interval $p_{\textrm {fh}}(R)< p_{0}< p_{\textrm {mm}}(R)$.

There are a number of theoretical papers (Kotelnikov et al. Reference Kotelnikov, Bagryansky and Prikhodko2010; Kotelnikov Reference Kotelnikov2011; Beklemishev Reference Beklemishev2016) that predict the stabilization of mirror instability in a non-uniform magnetic field. According to these works, when the mirror instability threshold is exceeded, a region appears in the plasma near the trap axis, where something like a magnetic field jump is formed, and non-paraxiality effects smooth out this jump. It is not yet clear at the moment whether the LoDestro equation can be modified to extend its range of applicability to equilibria with narrow non-paraxial jumps inside the plasma column.

The condition (3.6) guarantees the stability of small-scale perturbations of the firehose type, which, if the parameter $p_{0}$ slightly exceeds $p_{\text{fh}}(R)$, can become unstable near the axis of the plasma column. It is expected that such perturbations are not critical for the rigid ballooning mode.

Figure 4 confirms the above statement that the relation $\beta (z) = 2p_{\mathrel \perp }/B_{v}^{2}$, which has the meaning of a local beta, has a global maximum in the median plane of the trap. Moreover, $\beta (z)$ decreases monotonically in a monotonically increasing magnetic field. However, kinetic calculations of the distribution function of sloshing ions using various numerical codes, in particular the DOL code (Yurov, Prikhodko & Tsidulko Reference Yurov, Prikhodko and Tsidulko2016), show that, if the condition (3.5) is violated, the $\beta (z)$ dependence may turn out to be non-monotonic and a local peak appears on the $\beta (z)$ graph near fast ion turning point.

Figure 4. Axial profile of the local beta $\beta (z) = 2p_{\mathrel \perp }/B_{v}^{2}$ in the magnetic field (2.3) for mirror ratio $M=8$, index $q=4$ and various mirror ratios $R$ at the turning points where the pressure of the hot plasma component drops to zero: (a) transversal NBI, pressure model A1; (b) oblique NBI, model A2; (c) oblique NBI, model A3. The values of the $\beta$ parameter for each $R$ value indicated on the graphs are chosen to be equal to the smallest of the two stability thresholds for the mirror and firehose instabilities.

3.3. Pressure model A3

In the model A3, the pressure functions have the form

(3.28)\begin{equation} \left. \begin{aligned} p_{\mathrel\perp} & = p_{0} f_k(\psi )\, b^2\left(1-b\right)^{2},\\ p_{\|} & = \tfrac{1}{3}\, p_{0} f_k(\psi )\, b\left(1-b\right)^3 . \end{aligned} \right\} \end{equation}

In at least one respect, the A3 model is more realistic than the other anisotropic pressure models discussed so far. The point is that in models A1 and A2 the derivative ${\partial p_{\mathrel \perp }}/{\partial b}$ suffers a discontinuity at $b=1$. It should be expected that the actual pressure distribution in a mirror trap should not have such a discontinuity, since it will inevitably be smoothed out by Coulomb collisions of plasma particles. At the discontinuity point, the derivative ${\partial b}/{\partial z}$ tends to infinity, as on the threshold of mirror instability. In the A3 model, the derivative ${\partial b}/{\partial z}$ also tends to infinity as it approaches the mirror instability threshold, but this happens at some distance from the turning point $b=1$, namely, at $b=3/4$.

Comparing figures 2(a), 2(b) and 2(c), which show the axial transverse pressure profiles for models A1, A2 and A3, respectively, one can see that, in the latter panel, the pressure profile is smoother near the turning point. Another difference is that, for the same values of the parameter $R$, the maximum pressure in the right panel is shifted closer to the centre of the trap compared with the central panel.

It was noted in Kotelnikov et al. (Reference Kotelnikov, Prikhodko and Yakovlev2023) and Zeng & Kotelnikov (Reference Zeng and Kotelnikov2024) that the parameter $R$ characterizes the degree of plasma anisotropy in the A1 model. According to (3.17), the degree of anisotropy is higher for smaller $R$. In the case of oblique NBI (3.2) this connection does not seem obvious. Kesner (Reference Kesner1985) defines the degree of anisotropy as the ratio of the transverse pressure to the longitudinal pressure at the magnetic field minimum

(3.29)\begin{equation} A_{{K}} = {p_{\mathrel\perp}}/{p_{\|}}, \quad \text{at}\ b=b_{\min}=\sqrt{1-\beta }/R. \end{equation}

Examining the quantity $p_{\mathrel \perp }/p_{\|}=nb/(1-b)$ at an arbitrary $b$ reveals it monotonically increases up to infinity as it approaches the turning point $b=1$. The point $b=b_{\min }$ is distinguished only by the fact that the local beta maximum is located there. However, in the model (3.2) the transverse pressure peak is in the field $b=2/(n+1)$, where the ratio ${p_{\mathrel \perp } }/{p_{\|}}=2n/(n-1)$ does not depend on $R$. This and other similar facts indicate that model (3.2) with any $n$ actually describes the pressure distribution with an approximately fixed degree of anisotropy.

In § 3.4 it is shown that the parameter $R$ is related to the angle of injection of neutral atoms. Note, however, that it is possible to invent such a definition of the degree of anisotropy that yields the same result as (3.16) for the A1 model, namely

(3.30)\begin{equation} A_{n} = \frac{p_{\mathrel\perp}}{p_{\mathrel\perp}+n p_{\|}} = \frac{\sqrt{1-\beta }}{R} . \end{equation}

In the A3 model, the true magnetic field in the plasma is expressed in terms of the real root of a fourth-degree polynomial equation

(3.31)\begin{equation} b = \operatorname{Root}\left[ 2 {\#}^4 p - 4 {\#}^3 p + {\#}^2 (2 p+1) -b_{v}^2 \ \& , \ 2 \right] . \end{equation}

Parameter $\beta$ is still defined by (3.8), but (3.21) should be replaced by

(3.32)\begin{equation} p_{0} = \frac{\beta}{2 (1-\beta) \left(1-\sqrt{1-\beta }/R\right)^2} . \end{equation}

The inverse expression is noticeably more complex than (3.22)

(3.33)\begin{align} \beta & = \operatorname{Root}\left[ 4 \#^4 p_{0}^2 + \#^3 \left(8 p_{0}^2 R^2-16 p_{0}^2-4 p_{0} R^2\right) \right.\nonumber\\ & \quad \left. +\, \#^2\left(4 p_{0}^2 R^4-24 p_{0}^2 R^2+24 p_{0}^2+4 p_{0} R^4+8p_{0} R^2+R^4\right) \right.\nonumber\\ & \quad \left. +\, \# \left({-}8 p_{0}^2 R^4+24 p_{0}^2 R^2-16 p_{0}^2-4 p_{0} R^4-4 p_{0} R^2\right) \right.\nonumber\\ & \quad \left. +\, 4 p_{0}^2 R^4-8 p_{0}^2 R^2+4 p_{0}^2\ \& ,\ 1\right]. \end{align}

The threshold value of $p_{0}$, above which the mirror instability is excited, is now greater than prescribed by (3.24)

(3.34)\begin{equation} p_{\text{mm}}(R) = 4 . \end{equation}

The beta threshold is computed as a root of (3.33) at $p=4$

(3.35)\begin{align} \beta_{\text{mm}}(R) & = \operatorname{Root}\left[ 64 {\#}^4 + {\#}^3 \left(112 R^2-256\right) + {\#}^2 \left(81 R^4-352 R^2+384\right) \right.\nonumber\\ & \quad \left. +\, {\#} \left({-}144 R^4+368 R^2-256\right) + 64 R^4-128 R^2+64 \ \& , \ 1 \right] . \end{align}

The firehose instability threshold is described by equally cumbersome formulas

(3.36)\begin{align} p_{\text{fh}}(R) & = \operatorname{Root}\left[ {\#}^3 \left(9R^6+18R^4+9R^2\right) \right.\nonumber\\ & \quad \left.+\, {\#}^2 \left(2R^6-150R^4-264R^2-128\right) + {\#} \left(96R^2-117R^4\right) - 18R^4 \ \& , \ j \right] , \end{align}

where $j=1$ if $R < \sqrt {\tfrac {1}{2} (59+11\sqrt {33})}$, and $j=3$ if $R> \sqrt {\tfrac {1}{2} (59+11\sqrt {33})}$. But the firehose threshold over beta is written almost as simply as (3.27)

(3.37)\begin{equation} \beta_{\text{fh}}(R) = \operatorname{Root}\left[ 4 {\#}^3+{\#}^2 \left(R^2-28\right) +60 {\#}-36 \ \& , \ j \right] , \end{equation}

where $j=1$ if $1 R < R_{\ast }=\tfrac {1}{2}\operatorname {Root}[{\#}^4-236 {\#}^2-2048\ \& , \ 2]$, and $j=3$ if $R> R_{\ast }$.

The stability zone bounded by the conditions $\beta <\beta _{\textrm {mm}}(R)$ and $\beta <\beta_{\text{fh}}(R)$ is shown in figure 3(b). Its width shrinks to zero both in the $R\to 1$ limit and in the $R\to \infty$ limit. The maximum width corresponds to the point of intersection of the curves $\beta =\beta _{\textrm {mm}}(R)$ and $\beta =\beta_{\text{fh}}(R)$, where $\beta =0.8386$, $R=2.071$. At the same time, it is worth noting that the maxima of these functions are respectively $\beta _{\textrm {mm}}(\infty )=8/9$ and $\beta_{\text{fh}}(1)=0.9468$. Comparing (a) and (b) in figure 3 one can see that, in model A3, this zone is noticeably larger than in model A2. As can be assumed by looking at the graphs in figure 5 in the next section, this fact is associated with the observation that the distribution function in model A3 is slightly closer to isotropic than in model A2.

Figure 5. Angular distribution of the fast ions in the A1, A2, A3 and A8 models for a set of parameters $R$ indicated in the graphs.

3.4. Angle distribution of fast ions

For an arbitrary index $n$, the oblique injection model reads

(3.38)\begin{gather} p_{\mathrel\perp} = p_{0} f_k(\psi )\, b^2\left(1-b\right)^{n-1}, \end{gather}
(3.39)\begin{gather}p_{\|} = p_{0} f_k(\psi )\, \frac{b}{n}\left(1-b\right)^{n}, \end{gather}
(3.40)\begin{gather}p_{0} = \frac{\beta }{ 2\left(1-\beta\right) \left(1-\sqrt{1-\beta}/R \right)^{n-1}} . \end{gather}

The threshold value of parameter $p_{0}$ with respect to the excitation of the mirror instability is given by

(3.41)\begin{equation} p_{\text{mm}}(R) = \left( \frac{n-2}{n+1} \right)^{2-n} . \end{equation}

It is not possible to write formulas for $\beta _{\textrm {mm}}(R)$, $p_{\text{fh}}(R)$, $\beta_{\text{fh}}(R)$ in compact form.

Let us try to find an answer to the question: How realistic are the analytical models of anisotropic pressure? Is it possible to specify the distribution function of fast ions, which forms the pressure profiles that are given by (3.1) and (3.2)?

It is known that if the distribution of fast ions is given by a function $F(\varepsilon,\mu )$ of two variables, the energy $\varepsilon$ and the magnetic moment $\mu$, then the dependence of the density $N$, transverse $p_{\mathrel \perp }$ and longitudinal pressure $p_{\|}$ on the magnetic field can be calculated through double integrals of the distribution function over the variables $\varepsilon$ and $\mu$ with different weight factors

(3.42)\begin{gather} N(B) = \frac{2 \sqrt{2} {\rm \pi}B}{m^{3/2}} \int_{0}^{\infty} \left(\int_{\epsilon/B_{R}}^{\epsilon/B} \frac{F}{\sqrt{\epsilon -\mu B}} \, \mathrm{d}\mu \right) \, \mathrm{d}\epsilon, \end{gather}
(3.43)\begin{gather}p_{\mathrel\perp}(B) = \frac{2 \sqrt{2} {\rm \pi}B^{2}}{m^{3/2}} \int_{0}^{\infty} \left(\int_{\epsilon/B_{R}}^{\epsilon/B} \frac{\mu F}{\sqrt{\epsilon -\mu B}} \, \mathrm{d}\mu \right) \, \mathrm{d}\epsilon, \end{gather}
(3.44)\begin{gather}p_{\|}(B) = \frac{4 \sqrt{2} {\rm \pi}B}{m^{3/2}} \int_{0}^{\infty} \left(\int_{\epsilon/B_{R}}^{\epsilon/B} \sqrt{\epsilon -\mu B} \, F\, \mathrm{d}\mu \right) \, \mathrm{d}\epsilon . \end{gather}

The inverse problem of recovering the distribution function $F$ depending on two variables from a given dependence of $p_{\mathrel \perp }(B)$ on one variable obviously does not have a unique solution. The situation will not change if $p_{\|}(B)$ is added to $p_{\mathrel \perp }(B)$, since these two functions are related by the longitudinal equilibrium equation. The function $N(B)$ is not yet known to us.

The situation changes radically if the distribution function can be rewritten in the form of $F=F(\varepsilon,\mu /\varepsilon )$, for example, $F=g(\varepsilon )\,F_{n}(\mu B_{R}/\varepsilon )$. Then the integrals are separated. By replacing $\mu =\varepsilon \sin ^{2}\theta /B_{\min }$, $B=B_{R}b$, $B_{R}=B_{\min }R$, $\sin ^{2}\theta =X/R$ in (3.43), one obtains

(3.45)\begin{equation} p_{\mathrel\perp}(b) = \frac{2\sqrt{2}{\rm \pi} b^{2}}{m^{3/2}} \left( \int_{0}^{\infty} \varepsilon^{3/2}g(\epsilon ) \, {\rm d}\epsilon \right) \left( \int_{1}^{1/b} \frac{X F_{n}(X)}{\sqrt{1 -b X}} \, {\rm d}X \right) . \end{equation}

The integral in the first pair of brackets gives a constant independent of $b$. Omitting this and other constants leads to the Volterra integral equation of the first kind (see, e.g. Polyanin & Manzhirov Reference Polyanin and Manzhirov2008, chap. 10) for determining the function $F_{n}(X)$

(3.46)\begin{equation} (1-b)^{n-1} = \int_{1}^{1/b} \frac{X F_{n}(X)}{\sqrt{1-b X}} \, \mathrm{d}{X} .\end{equation}

It is solved using the Laplace transform and the convolution theorem. To bring the last equation to the canonical formulation of the convolution theorem, one needs to make a few more changes of variables: $b=1/t$; $t=1+\tau$, $X=1+x$. This yields the equation

(3.47)\begin{equation} h_{n}(\tau) = \frac{\tau^{n-1}}{\left(1+\tau\right)^{n-1/2}} = \int_{0}^{\tau} K(\tau-x)\,Q_{n}(x) \mathrm{d}{x} , \end{equation}

where $K(\tau )=1/\sqrt {\tau }$ is the kernel of the integral equation, and $Q_{n}(x)=(x+1)F_{n}(x+1)$. By the convolution theorem, the Laplace image of the desired function $Q_{n}(x)$ is equal to the Laplace image of the function $h_{n}(\tau )$ on the left-hand side of (3.47) divided by the Laplace image of $K (\tau )$. After performing the inverse Laplace transform, it is possible to restore the function $F_{n}(X)$. The result takes an unexpectedly simple form

(3.48)\begin{equation} F_{n}(X) = \frac{\varGamma(n)}{\sqrt{\rm \pi}\varGamma(n-1/2)} \frac{(X-1)^{n-3/2}}{X^{n+1}} . \end{equation}

Recall that $X=R\sin ^{2}\theta$, where $\theta$ is the pitch angle in the median plane. Function (3.48) passes through a maximum at $X=2(n+1)/5$. In terms of the pitch angle, the last relation means that

(3.49)\begin{equation} \sin^{2}\theta = 2(n+1)/5R . \end{equation}

Thus, parameter $R$ can be interpreted as a measure of the slope of the NBI.

Examples of the angular distribution of charged particles are shown in figure 5 for the A1, A2, A3 and A8 pressure models. Comparing these figures, it can be concluded that the width of the angular distribution increases, and the degree of anisotropy in its intuitive interpretation decreases as the index $n$ gets larger.

Knowing the function $F_{n}(X)$, one can find the plasma density distribution

(3.50)\begin{gather} N = n_0 f_k(\psi) \frac{b}{2n}\left(1-b\right)^{n-1} \left(1+(2n-1)\,b\right) , \end{gather}
(3.51)\begin{gather}p_{0}/n_{0} = \left.\int_{0}^{\infty }\varepsilon^{3/2}g(\varepsilon )\,\mathrm{d}\varepsilon \right/ \int_{0}^{\infty }\varepsilon^{1/2}g(\varepsilon )\,\mathrm{d}\varepsilon . \end{gather}

In a similar way, it is possible to restore the angular part of the distribution function in model A1

(3.52)\begin{equation} F_{1}(X) = \frac{ X^2+3 \sqrt{X-1}\, X^2 \tan ^{{-}1} \left(\sqrt{X-1}\right) -X-2 }{2{\rm \pi} \sqrt{X-1}\, X^2} , \end{equation}

and then the dependence of density on the magnetic field

(3.53)\begin{equation} N = n_0 f_k(\psi)\, \frac{1}{2}\, (1-b) (b+3) . \end{equation}

4. MHD stabilization by lateral wall

It is advisable to divide the study of MHD stabilization of plasma in an open trap using a perfectly conducting lateral wall and end MHD stabilizers into three parts. To begin with, this section presents calculations for the case when only the lateral wall is used. The next § 5 details stabilization by the MHD end stabilizers. Finally, in § 6 maps of stability zones under the combined action of the lateral wall and end stabilizers are presented.

If the lateral conducting wall is the only means of suppressing MHD instabilities, any solution of the Sturm–Liouville problem for the LoDestro equation gives the second critical beta, $\beta _{\textrm {cr}2}$, which determines the lower margin of the upper stability zone $\beta >\beta _{\textrm {cr}2}$. Results of such solutions are reported below in this section for the magnetic field (2.3) with mirror ratios selected from the set $M \in \{16,8,4,2\}$ for some most informative combinations of parameters $k\in \{1,2,4,\infty \}$, $q \in \{2,4,8\}$ and $R$. As can be seen in figure 1, the real magnetic field in the GDT and WHAM devices is better approximated by large values of $M$ and $q$. However, these values are not very convenient for identifying trends caused by changes in the proportion of magnetic mirrors in the total length of the magnetic trap. In addition, these values are reserved for comparison with upcoming calculations with the real magnetic field in the GDT.

Parameter $R$ varies from $R=1.1$ to $R=M$, taking discrete values from some predefined set, usually $R \in \{1.1,1.2,1.5,2,4,8,16\}$, but the maps in the plane $\{R,\beta \}$ are based on a more dense set.

The width of the vacuum gap between the lateral conducting wall and the lateral surface of the plasma column is determined by the ratio of the radius of lateral conducting wall $r_{w}(z)$ to the radius of the plasma column $a(z)$. The ratio $r_{w}(z)/a(z)$ enters the LoDestro equation (A2) through the function

(4.1)\begin{equation} \varLambda(z) = \frac { r_{w}^{2}(z) + a^{2}(z) }{ r_{w}^{2}(z) - a^{2}(z) } . \end{equation}

The PEK code accepts as an input parameter the ratio of the wall radius $r_{w0}=r_{w}(0)$ to the radius of the plasma column $a_{0}=a(0)$ in the median plane $z=0$ of an imaginary trap, where the magnetic field is minimal. Calculation was made for a discrete set of $r_{w0}/a_{0}$ values from $\sqrt {501/499}$ to $\infty$, which corresponds to $\lambda _{0}=\varLambda (0)$ values from $500$ to $1$.

The results of calculations are presented below for both the proportional lateral wall shape (labelled by Pr) and the straightened one (St). In the first case, the ratio $r_{w}(z)/a(z)=r_{w0}/a_{0}=\operatorname {const.}$ is the same in all cross-sections $z=\operatorname {const.}$ of an imaginary trap, so the function (4.1) is also constant, $\varLambda (z)\equiv \varLambda _{0}$. The radius of the plasma column in all cases decreases monotonically from the maximum value $a_{0}$ in the median plane $z=0$ through $a_{R}=\sqrt {2}$ at the turning point to the minimum value in the magnetic mirrors $a_{\min }=a(\pm 1)=\sqrt {2R/M}$. Note that $a_{0}$ in the theoretical treatment depends on parameter $R$, radial index $k$ and the value of $p_{0}$ to be given as an input parameter at the code start. On the contrary, in an experimental environment both $r_{w0}$ and $a_{0}$ can be considered as given quantities.

The width of the vacuum gap in the proportional conducting shell decreases monotonically. On the contrary, in a straightened chamber, the radius of the conducting wall is fixed, $r_{z}=r_{w}(0)$, so the vacuum gap increases monotonically as an observation point moves from the centre of the trap to the mirror plugs. Thus, for the same ratio $r_{w0}/a_{0}$ the average width of the vacuum gap is larger in case of straightened lateral wall. It therefore turns out that the stable zones in the stability maps are always smaller in the latter case (Zeng & Kotelnikov Reference Zeng and Kotelnikov2024).

4.1. Minimal vacuum gap

There is hardly any need to prove that the stabilizing effect of the lateral conducting wall is maximum when it is located closest to the surface of the plasma column. Therefore, as a first study of a new configuration, it is reasonable to carry out calculations of the critical beta with the conducting wall located as close as possible to the surface of the plasma column. For the pressure model A1, corresponding to the transverse NBI, such calculations are given in Kotelnikov et al. (Reference Kotelnikov, Prikhodko and Yakovlev2023) and Zeng & Kotelnikov (Reference Zeng and Kotelnikov2024). Figures 6 and 7 are compiled to compare the stabilizing effect of the lateral wall for the three pressure models: A1, A2, A3. They show a series of graphs for the case $\varLambda _{0}=500$ and illustrate the dependence of $\beta _{\textrm {cr}2}$ on $R$ and the lateral wall shape. Graphs in odd and even rows are drawn for proportional and straightened conducting walls, respectively.

Figure 6. Second critical beta vs $R$ in the limit $\varLambda \to \infty$ for three pressure models: A1 (a,d,g,j), A2 (b,e,h,k) and A3 (c,f,i,l). Stability zone of the ballooning rigid mode for a radial pressure profile with a given index $k$ is located above the margin curve, coloured as indicated in the legend under (jl).

Figure 7. Second critical beta vs $R$ in the limit $\varLambda \to \infty$ for three pressure models: A1 (a,d,g,j), A2 (b,e,h,k) and A3 (c,f,i,l). Stability zone of the ballooning rigid mode for a radial pressure profile with a given index $k$ is located above the margin curve, coloured as indicated in the legend under (jl).

Both figures confirm the previously discovered strong dependence of the critical beta $\beta _{\textrm {cr}2}$ on the shape of the radial pressure profile, represented by the index $k$. The region of stability with respect to the flute and ballooning disturbances in these figures is located above the marginal beta curve $\beta _{\textrm {cr}2}(R)$, the colour of which is associated with the index $k$ according to the legend under the bottom row of graphs in each figure. The areas of mirror and firehose instabilities are hatched by the scheme presented in figure 3. Minimal area of stability in respect to ballooning modes is shaded with the colour of the most upper curve. Most often, this is the blue colour corresponding to the most smooth radial profile $k=1$, but in figure 7(b,h) this blue zone is invisible since the corresponding blue curve is too short because it crosses the lower margin of the mirror mode. The PEK code quits above the mirror threshold where paraxial equilibrium is not possible.

Previously, it was found for model A1 (Kotelnikov et al. Reference Kotelnikov, Prikhodko and Yakovlev2023; Zeng & Kotelnikov Reference Zeng and Kotelnikov2024) that, when the end MHD stabilizers are switched off, the critical value $\beta _{\textrm {cr}2}$ depends only very weakly on both the mirror ratio $M$ and index $q$, which controls the axial profile of the magnetic field. For this reason, in figure 6, illustrating the dependence on the parameter $M$, only graphs with index $q=4$ are kept, and the set of values of $M$ is reduced to two, $M=4$ and $M=16$. Similarly, in figure 7, illustrating the dependence on $q$, only graphs with the mirror ratio $M=8$ and a pair of index values $q$ are kept, namely: $q=2$ and $q=8$. In general, all these graphs only confirm the conclusion about the weak influence of the magnetic field on the stabilizing properties of the lateral conducting wall within the same model of anisotropic pressure.

However, the differences between the three studied anisotropic pressure models A1, A2 and A3 are visible to the naked eye. These conclusions partially coincide with those of Kesner (Kesner Reference Kesner1985). He stated that the critical beta was insensitive to the axial profile of the magnetic field, but did not specify the value of the index $n$ in (3.2) used in his calculations.

In model A2, the stabilization boundary of ballooning instability is closely adjacent to the threshold of mirror instability, so that the region of joint stability $\beta _{\textrm {cr}2} < \beta < \beta _{\textrm {mm}}$ turns out to be very narrow. In model A3 this region is noticeably wider, especially in the region $R\lesssim 2$. In both models, the zone of stabilization of ballooning modes partially overlaps the region of firehose instability at $R\gtrsim 4$.

Thus, it can be concluded that the stability properties of an anisotropic plasma are very sensitive to details of the axial pressure profile.

In the most general terms, the last statement can hardly be disputed. At the same time, it should be noted that model A2 apparently does not adequately describe a real experiment. As follows from the discussion of (3.28), the point of origin of mirror instability in this model coincides with the turning point $b=1$ of fast ions, since the derivative ${\partial p_{\mathrel \perp }}/{\partial b}$ at this point experiences a jump as a result of mathematical idealization. A quick look at figure 5 shows that both the increase in anisotropy and the decrease in NBI tilt angle make it more difficult to stabilize an anisotropic plasma.

Comparison of graphs in odd and even rows of figures 6 and 7 reveals that the shape of the conducting chamber significantly deforms the shape of the $\beta _{\textrm {cr}2}(R)$ curves in the case of oblique NBI simulated by the models A2 and A3. For the model A1, Zeng & Kotelnikov (Reference Zeng and Kotelnikov2024) came to the opposite conclusion that the effect of the shape of the conducting chamber is insignificant. The apparent contradiction can be explained by the fact that the width of the gap between the conducting wall and the plasma at the location of the pressure peak plays a decisive role. In model A1 (transverse NBI), the pressure peak is located in the midplane of the trap. For the same value of the parameter $\varLambda _{0}$, the vacuum gap in the midplane is the same for the proportional and straightened chambers, so the difference in the gap width in other sections of the trap is not very significant. In the case of the A2 and A3 models (oblique NBI), the peak is located outside the midplane. The size of the gap at the location of the peak in the strengthened chamber will be larger than in the proportional chamber; accordingly, the stabilizing effect of the conducting chamber will be less.

Comparison of graphs in successive columns of figures 6 and 7 shows that the critical beta $\beta _{\textrm {cr}2}$ is definitely smaller in the case of oblique injection (which is good), but the stability zone disappears at large values of parameter $R$, since the plot of $\beta _{\textrm {cr}2}(R)$ crosses the mirror instability threshold $\beta _{\textrm {mm}}(R)$ at approximately $R\approx 3.5\unicode{x2013}4$. Considering that (3.49) relates the parameter $R$ to the injection angle $\theta _{\textrm {inj}}$, one obtains $\theta _{\textrm {inj}}\approx 40^{\circ }$ in model A2 and $\theta _{\textrm {inj}}\approx 47^{\circ }$ in model A3. In a straightened chamber, the $\beta _{\textrm {cr}2}$ curve intersects the $\beta _{\textrm {mm}}$ curve at a lower value of the parameter $R\approx 2$, which corresponds to $\theta _{\textrm {inj }}\approx 56^{\circ }$ in model A2 and $\theta _{\textrm {inj}}\approx 70^{\circ }$ in model A3.

The second observation is that the straightened chamber stabilizes the smoothest radial pressure profile with index $k=1$ noticeably worse in the sense that the range of values of the $R$ parameter at which this profile can be stabilized is noticeably narrower compared with the proportional chamber. In § 6, it is shown that, in the presence of the end MHD stabilizers, on the contrary, the smoothest radial profile is the easiest to stabilize because the unstable zones on the stability maps in figures 11–18 are the smallest for the $k=1$ radial profile.

It is useful to compare the above depicted calculations with the results of predecessors, in particular, with the works of Kesner and his co-authors (Kesner Reference Kesner1985; Li et al. Reference Li, Kesner and Lane1985, Reference Li, Kesner and Lane1987a, Reference Li, Kesner and LoDestrob). In all these papers, the authors analyse the plasma model with a sharp boundary, which corresponds to the infinite index $k=\infty$. Direct comparison is possible only with the first of the listed publications, and then only with certain reservations. In other works, the differences in the formulation of the problem are too great.

In his first work (Kesner Reference Kesner1985), Kesner approximated the vacuum magnetic field with a parabola. In the model of magnetic field (2.3), a parabola can approximate the field near the midplane in the case of $q=2$. Therefore, results can be compared in cases where the pressure peak is concentrated near the median plane, that is, at small values of the parameter $R\in \{1.1,1.2,1.3\}$. The necessary data can be extracted fromFootnote 2 figure 2 of Kesner (Reference Kesner1985), which shows the calculation results (in our terms) for the model A1-LwPr and zero vacuum gap. Comparing these data with figure 2(a,d,g,j) of Kotelnikov et al. (Reference Kotelnikov, Zeng, Prikhodko, Yakovlev, Zhang, Chen and Yu2022b), one can see quite satisfactory agreement.

It is more difficult to compare the results of calculations in the case of oblique injection. Figure 3 of Kesner (Reference Kesner1985) shows the graph of $\beta _{\textrm {cr}2}$ depending on Kesner's anisotropy degree (3.29) in the range from $A_{{K}}=1$ to $A_{{K}} =6$. The explanation of that figure indicates that the plasma occupies the entire length of the trap up to the magnetic mirrors, i.e. $R=M$, but the magnetic field was still approximated by a parabola. The values of the parameters $M$, $R$, $n$ are not specified in the article, however, taking into account the definition (3.29) of the anisotropy parameter in that article, one can assume that there should be $n\gg M$ in order to have such values of $A_{{K}}$. Formally, this means that comparison with my calculations for $n=2$ (i.e. model A2) and $n=3$ (model A3) is impossible.

An alternative possibility to ‘organize’ a large parameter $A_{{K}}$ involves the simultaneous limits $R\to 1$ and $\beta \to 0$. Parameter $A_{{K}}$ may well fall into the range $1\unicode{x2013} 6$ for $n=2\unicode{x2013} 3$ if $\beta \ll 1$ and $R=1.1\unicode{x2013} 2$. In order to draw figure 8 in $\{A_{{K}}, \beta \}$ coordinates, suitable for comparison with figure 3 in Kesner (Reference Kesner1985), $\beta _{\textrm {cr}2}$ was calculated for a set of discrete values of $R$ with four values of the mirror ratio $M\in \{2, 4, 8, 16\}$. Comparing figure 8 and figure 3 of Kesner (Reference Kesner1985), one can see some qualitative discrepancies. The most evident of them is that, in figure 8(b,c), the width of the stable zone between the $\beta _{\textrm {cr}2}$ and $\beta _{\textrm {mm}}$ curves decreases as $A_{{K}}$ becomes smaller whereas figure 8 demonstrates an opposite tendency.

Figure 8. Second critical beta $\beta _{\textrm {cr}2}$ vs Kesner's degree of anisotropy $A_{{K}}$ defined by (3.29) as in Kesner (Reference Kesner1985) for three models of anisotropic pressure and four mirror ratios, $M\in \{4,8,16\}$. Dashed curve shows the mirror mode thresholds. Compare with figure 3 in Kesner (Reference Kesner1985).

It does not follow from this fact that Kesner's results are wrong. In fact, we do not quite understand his prerequisites because some important details are missed. In addition, the choice of parameters $A_{{K}}$ and $\beta$ for figure 3 of Kesner (Reference Kesner1985) axes is unsuccessful, since parameter $A_{{K}}$ itself depends on $\beta$.

4.2. Effect of the vacuum gap

The width of the vacuum gap between the lateral surface of the plasma column and the inner surface of the conducting chamber of a selected type is determined by the parameter

(4.2)\begin{equation} r_{w0}/a_{0}=\sqrt{(\varLambda_{0}+1)/(\varLambda_{0}-1)} . \end{equation}

In § 4.1 it was shown that, in the limit $r_{w0}/a_{0}\to 1$, the critical beta remains essentially unchanged when both the mirror ratio $M$ and the index $q$ change. This conclusion holds for any value of $r_{w0}/a_{0}$, so it is sufficient to show only one (say, the average) value of the index and one value of the mirror ratio. In addition, the range of changes in the parameter $R$ is limited below by the interval $R\in [1.1\ldots 4]$, since at $R\gtrsim 4$ the margin of stability with respect to ballooning oscillations lies inside the firehose instability zone.

Figures 9 and 10 contain series of graphs of $\beta _{\textrm {cr}2}$ vs ratio $r_{w0}/a_{0}$ at fixed mirror ratio $M = 8$ and fixed axial profile index $q=4$. Figure 9 is compiled for a proportional chamber LwPr, and figure 10 repeats the same graphs for the LwSt configuration.

Figure 9. Stability maps for the LwPr configuration and three pressure models: A1 (a,d,g,j,m), A2 (b,e,h,k,n) and A3 (c,f,i,l,o). The second critical beta is drawn as a function of $r_{w}/a$ for the model magnetic field (2.3) with index $q=4$, set of mirror ratios $R\in \{1.1,1.2,1.5,2,4,8\}$ at the turning point and fixed mirror ratio $M=8$. The stable zone for a radial pressure profile with an index $k$ is located above the curve $\beta _{\textrm {cr}2}$ coloured according to the legend under (mo). Compare with figure 10.

Figure 10. Stability maps for the LwSt configuration and three pressure models: A1 (a,d,g,j,m), A2 (b,e,h,k,n) and A3 (c,f,i,l,o). The second critical beta is drawn as a function of $r_{w}/a$ for the model magnetic field (2.3) with index $q=4$, set of mirror ratios $R\in \{1.1,1.2,1.5,2,4,8\}$ at the turning point and fixed mirror ratio $M=8$. The stable zone for a radial pressure profile with an index $k$ is located above the curve $\beta _{\textrm {cr}2}$ coloured according to the legend under (mo). Compare with figure 9.

Comparison of graphs within any row demonstrates a strong dependence of $\beta _{\textrm {cr}2}$ on the pressure model. It can be seen that, for $r_{w}/a\approx 2$, the lowest critical beta is much smaller in the case of the oblique NBI imitated by model A2 (middle column) compared with the transverse NBI described by model A1 (left column). In other words, with an oblique NBI it is easier to achieve a stable plasma confinement mode. The effect is even more significant in model A3 (right column). However, it is very difficult to stabilize a plasma with a maximally smooth radial pressure profile $k=1$ during oblique injection. In model A2, the blue curve corresponding to the index $k=1$ is completely absent from most graphs. In addition, in model A2, the boundary of the stability zone $\beta =\beta _{\textrm {cr}2}$ rests on the threshold of mirror instability $\beta =\beta _{\textrm {mm}}$ even for a not very large ratio $r_{w}/a$. As is noted above, model A2 is not realistic enough. Model A3 predicts a wider range of $r_{w}/a$ values at which stabilization by a properly designed conducting chamber is possible.

5. MHD stabilization by end anchors

A proven method for suppressing MHD instabilities in an axially symmetric mirror trap is to attach end MHD anchors to the central mirror cell on the side of each of the two magnetic mirrors (Ryutov et al. Reference Ryutov, Berk, Cohen, Molvik and Simonen2011). The PEK package simulates MHD end stabilizers using boundary conditions of the form (A11) in the $z=z_{E}$ plane, where an imaginary conducting end plate is located. The further such a plate is placed from the centre of the trap, the less its stabilizing effect (Zeng & Kotelnikov Reference Zeng and Kotelnikov2024).

If $z_{E}>1$, the end plate is located behind the magnetic mirror and may not be virtual but real, as in some experiments at the GDT installation (Soldatkina et al. Reference Soldatkina, Anikeev, Bagryansky, Korzhavina, Maximov, Savkin, Yakovlev, Yushmanov and Dunaevsky2017, Reference Soldatkina, Maximov, Prikhodko, Savkin, Skovorodin, Yakovlev and Bagryansky2020). It is interesting and to some extent unexpected that such a plate, contrary to well-founded fears, did not lead to degradation of plasma parameters, even if it was placed relatively close to the neck of the magnetic plug. Variant $z_{E}>1$ simulates a terminal MHD stabilizer of the cusp type (Taylor Reference Taylor1963; Logan Reference Logan1980, Reference Logan1981; Baker et al. Reference Baker, Garner, Parks, Sleeper, Okamura, Adati, Aoki, Fujita, Hidekuma and Hattori1984; Anikeev et al. Reference Anikeev, Bagryansky, Deichuli, Ivanov, Karpushov, Maximov, Pod'minogin, Stupishin and Tsidulko1997; Li et al. Reference Li, Zhu, Ren, Ying, Yang and Sun2023).

If $z_{E}<1$, the imaginary end plate is located in front of the magnetic plug. In the problem of stabilizing the rigid ballooning mode, we can assume that its role is played by a limiter in the form of a ring, which is used in the vortex confinement method (Bagryansky et al. Reference Bagryansky, Lizunov, Zuev, Kolesnikov and Solomachin2003; Bagryansky, Beklemishev & Soldatkina Reference Bagryansky, Beklemishev and Soldatkina2007; Beklemishev Reference Beklemishev2008; Bagryansky et al. Reference Bagryansky, Anikeev, Beklemishev, Donin, Ivanov, Korzhavina, Kovalenko, Kruglyakov, Lizunov and Maximov2011).

Calculation of the stability margin that is created by one or another design of the end MHD anchor is a separate task. It is beyond the scope of this article (see, for example, Nagornyi & Stupakov Reference Nagornyi and Stupakov1984; Kuz'min & Lysyanskij Reference Kuz'min and Lysyanskij1990). The end plate, located behind the magnetic plug, imitates an MHD anchor, which has a smaller stability margin. On the contrary, the end plate installed in front of the magnetic plug imitates an MHD anchor, which has a larger stability margin.

According to the historically established classification inside the PEK package, the variants $z_{E}=1$, $z_{E}>1$ and $z_{E}<1$ are designated as Cw (combined wall), Bw (blind wall) and Rw (ring wall), respectively, although, in the current version of the PEK package, the calculation of all three options is carried out by a common module. To limit the number of figures, only the case $z_{E}=1$ (Cw) is presented below, when an imaginary conducting plate is located in the neck of the magnetic plug.

The absence of a lateral conducting shell corresponds to the limit $r_{w}/a_{0}\to \infty$ (that is, $\varLambda \to 1$). In this limit, the shape of the conducting shell obviously does not matter and, in this sense, the options CwPr, CwSt (as well as BwPr, BwSt, RwPr, RwSt) are equivalent. The PEK package in this limit finds at most one root, but now, in contrast to § 4, this is the upper margin $\beta _{\textrm {cr}1}$ of the lower stability zone.

The series of figures 11–13 illustrates the dependence of $\beta _ {\textrm {cr1}}$ on the parameter $R$ for three models of anisotropic pressure A1, A2, A3, three values of the mirror ratio $M\in \{24,16,8\}$ and three magnetic field profiles, corresponding to index values $q\in \{2,4,8\}$. Analysis of the graphs in these figures shows that, when using MHD stabilizers, all previously made conclusions for the case of lateral wall stabilization have to be reversed.

Figure 11. First critical beta $\beta _{\textrm {cr}1}$ for the CwPr configuration and three pressure models: A1 (ac), A2 (df) and A3 (gi) at mirror ratio $M=15$ vs parameter $R$ in the limit $\varLambda \to 1$. The stability zone of rigid ballooning perturbation for a radial pressure profile with a given index $k$ is located below the curve, coloured as indicated in the legend under (gi).

Figure 12. Same maps as in figure 11 but for mirror ratio $M=8$.

Figure 13. Same maps as in figure 11 but for mirror ratio $M=4$. The rows A2-CwPr and A3-CwPr are not shown as they do not have an unstable zone for rigid ballooning modes.

The stability zone for the radial profile with index $k$ is located below the curve $\beta =\beta _{\textrm {cr1}}(R)$, coloured as indicated in the legend under the bottom row in each figure. Absence of a curve for a certain range of $R$ means that the corresponding radial profile $k$ is stable over the entire interval of beta values allowed by the inequalities $0<\beta < 1$ (in the case of transverse NBI as described by the A1 model) or $0<\beta <\beta _{\textrm {mm}}$ below the threshold of mirror instability (in the case of oblique NBI simulated by the models A2, A3, …). The areas where the mirror and firehose modes are unstable are shaded in blue and brown colours, respectively. Computation was performed for a wider set of mirror ratios $M$ then mentioned above but no $\beta _{\textrm {cr}1}$ was found for $M\lesssim 4$ .

Comparison of figures 11–13 reveals a strong dependence of the ballooning instability threshold on the mirror ratio $M$. This drastically differentiates the lower zone of stability from the upper one which exists exclusively due to the lateral conducting wall, as discussed in § 4. The lower stability zone tends to expand as the mirror ratio becomes smaller. This fact is quite consistent with the intuitive expectations that a decrease in $M$ improves the contact of the end MHD stabilizer with the central section of the open trap.

Comparison of the first, second and third rows in any of the three figures (more precisely: comparison of graphs in the same column, but in different rows) demonstrates that there is a qualitative and quantitative difference between the anisotropic pressure models. It is especially noticeable when passing from transverse injection (model A1) to oblique injection (models A2 and A3). The difference between models A2 and A3 is limited in quantitative indicators, although it is quite noticeable. It can also be noted that the influence of the radial profile on the value of $\beta _{\textrm {cr1}}$ in model A1 is more noticeable than in models A2 and A3. This statement follows from the fact that the distance between curves of different colours in the graphs in the first row of each of figures 11–13 is significantly larger than in the second and third rows. The absence of a blue, yellow or green curves means that the corresponding fairly smooth radial profiles are stable with respect to ballooning perturbations of plasma equilibrium. Again, this effect is more pronounced in the case of transverse NBI, whereas, in the case of oblique injection, the firehose instability threshold $\beta_{\text{fh}}$ in terms of beta value is often lower than the ballooning instability threshold $\beta _{\textrm {cr1} }$.

Comparison of graphs within every row of any of the three figures shows that the shape of the axial profile of the magnetic field has a significant effect on the value of $\beta _{\textrm {cr}1}$. This fact is most noticeable in the case of transverse injection (model A1, first row of each figure) and also significantly distinguishes stabilization using end MHD stabilizers from stabilization using a lateral wall. The general tendency is that the narrowing and steepening of magnetic mirrors expands the lower stability zone.

With transverse injection, the entire range $0<\beta <1$ can be made stable for any radial pressure profile if $M$, $R$ are sufficiently low. With oblique injection, the $\beta =1$ limit on transverse equilibrium is unattainable due to the development of mirror instability. However, the last statement was proven only within the paraxial approximation. A number of publications give reason to assume that the threshold of shear instability can be exceeded in non-paraxial open traps (Lansky Reference Lansky1993; Lotov Reference Lotov1996; Kotelnikov et al. Reference Kotelnikov, Bagryansky and Prikhodko2010; Kotelnikov Reference Kotelnikov2011; Beklemishev Reference Beklemishev2016).

Comparing figures 6 and 7 from § 4.1 with figures 11–13 it is easy to notice that, as $R$ increases, the boundaries of the lower and upper stability zones shift in opposite directions, and in such a way that both zones contract. Note also that in all the figures the blue curves are located above the yellow ones, the yellow ones are above the green ones and the green ones are above the red ones. In other words, the boundary of both the upper stability zone and the boundary of the upper stability zone shift downward as the radial pressure profile steepens. But if for the upper zone this order means that it becomes wider as the radial profile steepens, then the expansion of the lower zone, on the contrary, occurs as the radial profile is smoothed.

6. Combined MHD stabilization

Taking into account the existence of upper and lower stability zones when the two stabilization methods are applied separately, which were successively described in §§ 4 and 5, it is easy to believe that, with the simultaneous application of these two stabilization methods, both stability zones will be preserved. The first indication of the presence of two zones of stability can be found in the works of D'Ippolito and Hafizi (D'Ippolito & Hafizi Reference D'Ippolito and Hafizi1981) and D'Ippolito and Myra (D'Ippolito & Myra Reference D'Ippolito and Myra1984). After the author's work, it was shown earlier (Kotelnikov, Lizunov & Zeng Reference Kotelnikov, Lizunov and Zeng2022a; Kotelnikov et al. Reference Kotelnikov, Zeng, Prikhodko, Yakovlev, Zhang, Chen and Yu2022b, Reference Kotelnikov, Prikhodko and Yakovlev2023; Zeng & Kotelnikov Reference Zeng and Kotelnikov2024) that the simultaneous existence of two stability zones has become a proven fact using the examples of isotropic plasma and anisotropic plasma in the special case of transverse NBI. The lower zone exists at low plasma pressure, at $0<\beta <\beta _{\textrm {cr}1}$, and the second, at high pressure, at $\beta _{\textrm {cr}2}<\beta <$1 . These two zones merge for larger $\varLambda$, providing overall stability for any beta in the range $0<\beta <1$. It will be shown below that the last statement should be corrected in the case of inclined NBI.

When the lateral wall and end stabilizers act together, PEK often does not find solutions for the set of parameters discussed in the previous sections. This could mean either that the ballooning perturbations are stable over the entire range of available beta values below the mirror instability threshold, or that there is an error in the program. Therefore, calculations in the Cw, Bw and Rw modes usually began with the case $R=M$, when the instability zone has maximum dimensions for a fixed set of parameter values $k$, $q$, $M$, $z_{E}$. The remaining part of the article also begins with the case $R=M$. In order not to overload the article with an excessive number of figures, the results of calculations in the Bw and Rw configurations are omitted below. The main conclusion from the analysis of these configurations is quite banal. The Rw configuration is more stable than the Cw configuration, and the Bw configuration, on the contrary, is less stable in the sense that the instability zone is larger than in the Cw configuration.

6.1. Wide pressure peak, $M=R$

The series figures 14–16 visualizes the calculation results for the special case $R=M$ when the hot plasma component occupies the entire mirror trap from plug to plug. Initially, this case was highlighted among many others in Kotelnikov et al. (Reference Kotelnikov, Chen, Bagryansky, Sudnikov, Zeng, Yakovlev, Wang, Ivanov and Wu2020) in order to compare the stability of anisotropic plasma, which is formed under transverse NBI, with the stability of isotropic plasma. In the A1 model, the $R=M$ variant corresponds to the minimum anisotropy. However, in models A2 and A3, parameter $R$ cannot serve as a measure of anisotropy, as shown in § 3.4, but it is related to the injection angle $\theta _{\textrm {inj}}$ by (3.49). Too large values of $R$, $R\gg 2(n+1)/5$ correspond to too small angles of injection $\theta _{\textrm {inj}}$, which can hardly be of interest from a practical point of view. Taking into account this explanation, the range of the parameter $R$ is further reduced compared with the range that was adopted in previously published analogues for the transverse NBI model, in particular in figure 9 of Kotelnikov et al. (Reference Kotelnikov, Prikhodko and Yakovlev2023).

Figure 14. Stability maps for model magnetic field (2.3) and anisotropic plasma pressure models (3.10) (af), (3.18) (gl) and (3.28) (mr) at combined MHD stabilization by lateral wall and end MHD anchors, $q\in \{2, 4, 8\}$, $M=R=16$. The unstable zone is located between the lower $\beta _{\textrm {cr}1}(r_{w}/a_{0})$ and upper branches $\beta _{\textrm {cr}2}(r_{w}/a_{0})$ of every curve. Correspondence of the index $k$ to the colour of the curves is shown at the bottom of the figure. Shaded common zone of stability lies to the left of the curve for the most steep radial pressure profile ($k = \infty$). (a,d,g,j,m,p) Show the maps for a ‘parabolic’ magnetic field with the index $q=2$, (c,f,i,l,o,r) show the maps for the ‘quasi-flat hole’ magnetic field with $q=8$. Panels (ac,gi,mo) and (df,jl,pr) show maps for the CwPr and CwSt configurations, respectively.

Figure 15. Same as in figure 14 but for $M=R=8$.

Figure 16. Same as in figure 14 but for $M=R=4$. Rows for the A2 and A3 models are dropped since the rigid ballooning mode is stable in the entire region below the mirror instability threshold.

In the first two rows, figures 14–16 present the stability maps for the A1 model, corresponding to the transverse NBI. In figures 14 and 15, composed respectively for $M=16$ and $M=8$, they are followed by two rows of graphs for the oblique NBI model A2 and the last two rows contain maps for oblique NBI model A3. Figure 16 has only two tows since the plasma produced by oblique NBI is stable against ballooning perturbation in the entire range of $\beta$ below the mirror instability threshold.

The difference between the first two rows and the four following ones is so obvious that it almost could not be commented on. With oblique injection, there is practically no upper stability zone $\beta _{\textrm {cr}2}$, which was present both during transverse NBI and in isotropic plasma, as shown in figure 10 in Kotelnikov et al. (Reference Kotelnikov, Prikhodko and Yakovlev2023). It can be said that, with oblique NBI, the upper stability zone is absorbed into the region of mirror mode instability.

One can also notice a significant decrease in the value of $\beta _{\textrm {cr}1}$ for all radial pressure profiles. This effect is accompanied by the appearance of an instability zone for smooth radial profiles (with indices $k=1$, $k=2$) in the range of small values of parameter $r_{w}/a$ at which such zones were not present in the A1 model. The brown hatching from the third to sixth rows means that, in the A2 and A3 models, the margin of the stability zone of the rigid ballooning mode passes entirely inside the zone of firehose MHD instability. It is important to emphasize that the threshold of firehose instability is determined by the criterion (3.6), which is obtained for a homogeneous magnetic field. In an inhomogeneous magnetic field, the threshold should be slightly lower, as shown in V. Mirnov's dissertation (Mirnov Reference Mirnov1986). We are not aware of experiments with plasma in open traps in which manifestations of hose instability would be reliably recorded. Therefore, it would be premature to assert that a transition beyond the formal boundary of firehose instability is not possible.

As for the threshold of mirror instability, the PEK code is currently unable to perform calculations in the area above its threshold.

As can be seen from the analysis of figures 14–16, with joint stabilization by the lateral wall and end MHD stabilizers, the instability zone $\beta _{\textrm {cr}1} <\beta < \beta _{\textrm {cr}2}$ decreases as $q$ increases. In other words, open traps with short magnetic plugs are more stable than traps in which the magnetic field increases smoothly as you move from the centre of the trap to the magnetic plugs. This observation is consistent with the long-standing calculations of V. Mirnov and O. Bushkova. They calculated the threshold of ballooning instability with respect to small-scale disturbances (Bushkova & Mirnov Reference Bushkova and Mirnov1986), ignoring the effects of the finite Larmor radius. Later, the results of those calculations were confirmed in the work (Kotelnikov et al. Reference Kotelnikov, Lizunov and Zeng2022a).

With fixed parameters $q$, $M$, $R$, the instability zone is maximal for the steepest radial pressure profile ($k=\infty$) and may be absent altogether for smooth profiles ($k=1$, $k=2$). Without stabilization by the end MHD anchors, the opposite situation occurs: the stability zone is smaller and may be absent altogether for smooth profiles.

6.2. Effect of pressure peak width, $1< R< M$

The case $R=M$, considered in § 6.1, corresponds to the widest pressure profile along the trap axis. Recall that the parameter $R$ in the magnetic field model (2.3) makes sense of a plug ratio in the section of an open trap where the pressure of hot ions decreases to almost zero. Thus, parameter $R$ is a characteristic of the extent of the plasma pressure peak localization region, which is universally suitable for all A1, A2, A3, $\ldots$ anisotropic pressure models studied so far. In addition, in the A1 model, which simulates transverse NBI, parameter $R$ characterizes the degree of plasma anisotropy, whereas in the case of oblique NBI, described by the A2 and A3 models, it is associated with the angle of inclination of the NBI to the direction of the magnetic field.

To study the effect of the extent of the region occupied by fast ions, which are formed during the injection of beams of neutral atoms, in real experiments at the GDT facility at the Budker Institute of Nuclear Physics, the magnetic field is reprofiled in the axial direction, forming a short well (Shmigelsky et al. Reference Shmigelsky, Lizunov, Meyster, Pinzhenin, Solomakhin and Yakovlev2024). A similar operation in the numerical study with the model field (2.3) in use involves increasing parameters $q$ and $M$. As can be seen from figure 1, short magnetic plugs as in the GDT correspond to higher values of $q$ and $M$ than those adopted in the current and earlier works. Corresponding calculations for model magnetic field (2.3) are planned for publication together with calculations for the real magnetic field in GDT and ALIANCE facilities.

Maps of stability for the case $q=4$, $M=16$ and five values of the parameter $R\in {1.2, 1.5, 2, 4, 8}$ are presented in figure 17 for the proportional lateral wall chamber and figure 18 for the straightened chamber. The first thing that catches eye when analysing these figures is the absence of an unstable zone at small values of the parameter $R$ for all three models of anisotropic pressure A1, A2, A3, and the marginal value of $R$ for transverse NBI is noticeably smaller than for oblique NBI. From this point of view, oblique NBI appears to be preferable to transverse NBI.

Figure 17. Stability maps vs ratio $r_{w}/a$ for the A1, A2 and A3 pressure models in the CwPr configuration simulating combined stabilization by the proportional lateral conducting chamber and end MHD anchors; $q=4$, $M=16$, $R\in \{1.2, 1.5, 2, 4, 8\}$. The instability zone is located between $\beta _{\textrm {cr}1}(r_{w}/a)$ (lower branch of the marginal curve) and $\beta _{\textrm {cr}2}(r_{w}/a)$ (upper branch of the same colour); in the case of inclined injection, which corresponds to the A2 and A3 models, the upper branch is completely or partially absorbed by the region of mirror instability; the stability zone is shaded for a plasma with a sharp boundary ($k=\infty$), for which it has the minimum dimensions.

Figure 18. Same as in figure 17 but for CwSt configuration simulating combined stabilization by the straightened lateral conducting wall and the end MHD anchors.

The second conclusion that arises when comparing graphs in the same row  in these figures is that transverse NBI, with other things being equal, allows us to achieve higher values of relative pressure $\beta$ than oblique NBI. This circumstance is due to the fact that the limiting beta in the case of oblique injection is limited by the mirror instability, not the ballooning one. From this point of view, transverse NBI appears to be preferable to oblique NBI.

7. Conclusions

This article completes a series of three publications (Kotelnikov et al. Reference Kotelnikov, Zeng, Prikhodko, Yakovlev, Zhang, Chen and Yu2022b, Reference Kotelnikov, Prikhodko and Yakovlev2023; Zeng & Kotelnikov Reference Zeng and Kotelnikov2024) in which the stability of a rigid ballooning mode with an azimuthal number $m=1$ was studied using a model magnetic field profile (2.3) along the axis of a linear open trap and model pressure profiles of the hot plasma component along the radius (3.11). The use of model profiles at the stage of development, debugging and testing of the PEK software package was a natural and reasonable solution from many points of view. In particular, as the models become more complex, the calculation time increases by tens, hundreds and thousands of times. In addition, the use of model functions, depending on a small set of parameters, to approximate the vacuum magnetic field and the radial distribution of plasma pressure, made it possible to obtain relatively simple recommendations for achieving maximum beta in open traps. The main results can be formulated in several statements:

  1. (i) When only a lateral conducting wall is used for MHD stabilization, stable plasma confinement occurs in the upper zone $\beta > \beta _{\textrm {cr}2}$, which expands with a decrease in the mirror ratio $R$ at the turning point of fast ions produced by NBI. A steep radial pressure profile is the most stable in this case. The dependence on the magnetic field profile is minimal for transverse NBI, but significant for oblique NBI. The second marginal beta $\beta _{\textrm {cr}2}$ is only slightly affected by the mirror ratio $M$ and axial profile of the vacuum magnetic field.

  2. (ii) When only end MHD anchors are used, stable confinement occurs in the lower zone $\beta < \beta _{\textrm {cr}1}$, which expands with a decrease in $M$ and/or $R$. A smooth radial pressure profile and steep axial magnetic field profile are the most stable. The first marginal beta $\beta _{\textrm {cr}1}$ is strongly affected by the mirror ratio $M$ and the axial profile of the vacuum magnetic field.

  3. (iii) In the case of oblique NBI at a small angle of injection ($R\gtrsim 4$), beta is limited by the threshold of mirror instability rather than by the ballooning instability.

  4. (iv) Wall stabilization of ballooning perturbations is especially effective in combination with end MHD anchors, which make it possible to suppress all MHD oscillations in the range $0<\beta <1$ for transverse NBI and $0<\beta <\beta _{\textrm {mm}}(R)$ for oblique NBI.

  5. (v) The size of the stable zones in all cases essentially depends on the radial and axial profiles of the plasma pressure, width and shape of the vacuum gap, plasma anisotropy, angle of injection and stability margin of MHD anchors. The axial profile of the magnetic field is significant in all cases when the end MHD stabilizers are switched on.

Further development of the PEK software package should include a series of minor and major upgrades. The first step in this queue will be to connect the modules for calculating the actual magnetic field in the GDT, CAT and ALIANCE facilities. At the program level, this step has already been completed. The other most anticipated major modernization should be the inclusion of dissipative processes, both inside the plasma and in the conductors surrounding it.

List of notations

$M$

Mirror ratio at the neck of the magnetic plug

$R$

Mirror ratio at the turning point of the sloshing ions

$L$

Mirror ratio at the location of the limiting ring or conducting end wall simulating the MHD anchor

$\psi$

Normalized magnetic flux with $\psi =1$ at the lateral boundary of the plasma column

$z$

Normalized coordinate along the axes of the plasma column with $z=0$ in the midplane of the trap and $z=\pm 1$ on the magnetic plugs

$r$

Normalized radial coordinate $r=0$ on the axes of the trap and $r=a$ on the lateral surface of the plasma column

$a$

Normalized radius of the plasma column as a function $z$

$a_{0}$

Value of $a(z)$ in the trap midplane

$r_{w}$

Radius of the conducting cylinder surrounding the plasma column

$r_{w0}$

Value of $r_{w}$ in the trap midplane

$\varLambda$

Dimensionless function of the ratio $r_{w}(z)/a(z)$, included in the LoDestro equation

$\varLambda _{0}$

Value of $\varLambda (z)$ in the trap midplane

$b_{v}$

Normalized vacuum magnetic field as a function of $z$ with $b_{v}=1$ at turning points of sloshing ions and $b_{v}=1/R$ in the mirror trap midplane

$b$

Normalized true magnetic field with $b=1$ at the turning points of sloshing ions

$p_{\mathrel \perp }$

Normalized transverse plasma pressure as a function of $\psi$ and $z$

$p_{\|}$

Normalized longitudinal plasma pressure as a function of $\psi$ and $z$

$p_{0}$

Value of $p_{\mathrel \perp }$ on the axis of the plasma column in the midplane of the trap

$F$

Distribution function of sloshing ions as a function of energy $\varepsilon$ and magnetic moment $\mu$

$f_{k}$

Dimensionless function $\psi$ associated with the radial pressure profile with $f_{k}=1$ on the axes of the plasma column and $f_{k}=1$ on the lateral surface

$k$

Index of function $f_{k}$, where $k=\infty$ corresponds to the step profile pressure with a sharp boundary, and $k=1$ corresponds to a smooth quasi-parabolic profile

$q$

Index of the axial profile of the vacuum magnetic field, where $q=2$ corresponds to a quasi-parabolic profile, and larger $q$ values simulate a profile with shorter magnetic plugs

Lw

Lateral wall configuration without MHD end anchors

Cw

Combined wall configuration with lateral wall and MHD end anchors, simulated by conducting end plates located at the magnetic plugs

Bw

Blind wall configuration similar Cw to with weaker end MHD stabilizers imitated by the conducting end plates located behind magnetic plugs

Rw

Ring wall configuration with annular limiters located beyond the turning points and before magnetic plugs imitated by a conducting end plate located in front of the magnetic plug

Pr

Proportional shape of the wall of the conducting cylinder surrounding the plasma column

St

Straightened shape of the wall of the conducting cylinder surrounding the plasma column

$\beta$

The beta parameter is defined as the maximum of the ratio $2p_{\mathrel \perp }/b_{v}^{2}$

$\beta _{\textrm {cr}1}$

The first critical beta value defining the upper boundary of the lower stability zone against rigid ballooning modes due to MHD end anchors

$\beta _{\textrm {cr}2}$

Second critical beta value defining the lower boundary of the upper stability zone against rigid ballooning modes due to the lateral wall

$\beta _{\textrm {mm}}$

Beta value at the upper boundary of the stability zone against mirror modes

$\beta_{\text{fh}}$

Beta value at the upper boundary of the stability zone against hose modes

$\theta$

Pitch angle

$\theta _{\textrm {inj}}$

Inclination angle of NBI

Acknowledgements

The author thanks P. Bagryansky, C. Forest, A. Ivanov, V. Prikhodko, S. Putvinsky, E. Shmigelsky, D. Yakovlev, P. Yushmanov and Q. Zeng for numerous discussions and valuable comments at different stages of development of the PEK package. Special thanks to A. Beklemishev for discussing the conditions for the applicability of the LoDestro equation, V. Ilgisonis for discussing the firehose instability threshold, A. Milstein for discussing the properties of second-order ordinary differential equations with singular coefficients and V. Mirnov for answering many of my questions.

Editor C. Forest thanks the referees for their advice in evaluating this article.

Funding

This work was supported by the Russian Science Foundation under the Grant 24-12-00309 awarded to Novosibirsk State University. This work is also a part of the state assignment of the Russian Federation for the Budker Institute of Nuclear Physics.

Declaration of interests

The author reports no conflict of interest.

Data availability statement

The data that support the findings of this study are not openly available.

Appendix A. LoDestro equation

The LoDestro equation is a second-order ordinary differential equation for the function

(A1)\begin{equation} \phi(z) = a(z) B_{v}(z) \xi_{n}(z) , \end{equation}

which depends on single coordinate $z$ along the trap axis and is expressed in terms of the variable radius of the plasma column boundary $a=a(z)$, the vacuum magnetic field $B_{v}=B_{v}(z)$ and the virtual small displacement $\xi _{n}=\xi _{n}(z)$ of the plasma column from the axis. In its final form, the LoDestro equation reads

(A2)\begin{align} 0 & = \frac{\mathrm{d}{}}{\mathrm{d}{z}} \left[ \varLambda + 1 - \frac{2\left\langle \bar{p}\right\rangle}{B_{v}^{2}} \right] \frac{\mathrm{d}{\phi}}{\mathrm{d}{z}}\nonumber\\ & \quad +\phi \left[ - \frac{\mathrm{d}{}}{\mathrm{d}{z}}\left( \frac{B_{v}'}{B_{v}} + \frac{2a'}{a} \right) \left( 1 - \frac{\left\langle \bar{p}\right\rangle}{B_{v}^{2}} \right) + \frac{\omega^{2}\left\langle \rho\right\rangle}{B_{v}^{2}} \right.\nonumber\\ & \quad \left. -\, \frac{2\left\langle \bar{p}\right\rangle}{B_{v}^{2}}\frac{a_{v}''}{a_{v}} - \frac{1}{2}\left( \frac{B_{v}'}{B_{v}} + \frac{2a'}{a} \right)^{2} \left( 1 - \frac{\left\langle \bar{p}\right\rangle}{B_{v}^{2}} \right) \right] , \end{align}

where the derivative ${\mathrm {d}{}}/{\mathrm {d}{z}}$ in the first two lines acts on all factors to the right of it, and the prime ($'$) is a shorthand for ${\mathrm {d}{}}/{\mathrm {d}{z}}$. Other notations are defined as follows:

(A3)\begin{gather} \frac{a^{2}}{2} = \int_{0}^{1} \frac{\mathrm{d}{\psi}}{B}, \end{gather}
(A4)\begin{gather}\frac{r^{2}}{2} = \int_{0}^{\psi} \frac{\mathrm{d}{\psi}}{B}, \end{gather}
(A5)\begin{gather}B^{2} = B_{v}^{2} -2p_{\mathrel\perp}, \end{gather}
(A6)\begin{gather}a_{v}(z) = \sqrt{\frac{2}{B_{v}(z)}}, \end{gather}
(A7)\begin{gather}\bar{p} = \frac{p_{\mathrel\perp} + p_{\|}}{2}, \end{gather}
(A8)\begin{gather}\left\langle \bar{p}\right\rangle=\frac{2}{a^{2}}\int_{0}^{1} \frac{\mathrm{d}{\psi}}{B}\,\bar{p}, \end{gather}
(A9)\begin{gather}\varLambda = \frac { r_{w}^{2} + a^{2} }{ r_{w}^{2} - a^{2} } . \end{gather}

Equation (A4) relates the radial coordinate $r$ and the magnetic flux $\psi$ through a ring of radius $r$ in the $z$ plane. The magnetic field $B=B(\psi,z)$, weakened by the plasma diamagnetism, in the paraxial (long–thin) approximation (which assumes small curvature of field lines) is related to the vacuum magnetic field $B_{v}=B_ {v }(z)$ by the equation of transverse equilibrium equation (A5). Kinetic theory predicts (see, for example, Newcomb Reference Newcomb1981) that the transverse and longitudinal plasma pressures can be considered as functions of the magnetic field $B$ and magnetic flux $\psi$, i.e. $p_{\mathrel \perp }=p_{\mathrel \perp }(B ,\psi )$, $p_{\|}=p_{\|}(B,\psi )$. In (A2), one must assume that the magnetic field $B$ is already expressed in terms of $\psi$ and $z$, and therefore $p_{\mathrel \perp }=p_{\mathrel \perp }(\psi,z)$, $p_{\|}=p_{\|}(\psi,z)$. The angled brackets in (A2) denote the mean value of an arbitrary function of $\psi$ and $z$ over the plasma cross-section. In particular, the average value $\left \langle \rho \right \rangle$ of the density $\rho =\rho (\psi,z)$ is calculated using a formula similar to (A8), and $\omega$ is the angular frequency of oscillations.

Parameter $\varLambda$ is, generally speaking, a function of the $z$ coordinate. It implicitly depends on the plasma parameters and magnetic field through the dependence of the plasma column radius $a=a(z)$ on them. In the special case of a proportional chamber, when $r_{w}(z)/a(z)=\operatorname {const.}$, the function $\varLambda (z)$ becomes constant, which simplifies the equation somewhat. Namely, for the sake of such simplification, in most papers by other authors they assume that $\varLambda =\operatorname {const.}$.

Traditionally, two types of boundary conditions are considered. In the presence of conducting end plates located in the magnetic mirrors at $z=\pm 1$, the boundary condition

(A10)\begin{equation} \phi(\pm1) = 0 \end{equation}

should be chosen. In a more general case, when the conducting end plate is installed somewhere in the behind-the-mirror region, namely, in the plane with coordinates $z=\pm z_{E}$, the zero boundary condition must obviously be assigned to this plane

(A11)\begin{equation} \phi({\pm} z_{E}) = 0 . \end{equation}

By solving the LoDestro equation with the boundary condition (A11), it is possible to simulate the effect of the end MHD anchors with different stability margins.

If the plasma ends are electrically isolated, the boundary condition

(A12)\begin{equation} \phi'({\pm} z_{E}) = 0 \end{equation}

is applied at $z=\pm z_{E}$. As a rule, it implies that other methods of MHD stabilization in addition to stabilization by a conducting lateral wall are not used.

An obvious fit to the LoDestro equation with boundary conditions (A11) or (A12) is the trivial solution $\phi \equiv 0$. To eliminate the trivial solution, we impose a normalization in the form of one more condition

(A13)\begin{equation} \phi(0)=1. \end{equation}

Taking into account the symmetry of the magnetic field in actually existing open traps with respect to the median plane $z=0$, it suffices to find a solution to the LoDestro equation at half the distance between the magnetic mirrors, for example, in the interval $0< z<1$. Due to the same symmetry, the desired function $\phi (z)$ must be even, therefore

(A14)\begin{equation} \phi'(0)=0. \end{equation}

It is convenient to search for a solution to the LoDestro equation by choosing the boundary conditions (A13) and (A14). In theory, a second-order linear ordinary differential condition with two boundary conditions must always have a solution. However, the third boundary condition (A11) or (A12) can only be satisfied for a certain combination of parameters. If the parameters of the plasma, magnetic field and geometry of the lateral conducting wall are given, the third boundary condition should be considered as a nonlinear equation for the square $\omega ^{2}$ of frequency $\omega$. If the root of such an equation is positive, the MHD oscillations with azimuthal number $m=1$ are stable; if $\omega ^{2}<0$, then instability takes place. On the margin of the stability zone $\omega ^{2}=0$. In this case, the solution of the boundary-value problem (A2), (A13), (A14) with the additional boundary condition (A11) or (A12) gives the critical value of beta, respectively both $\beta _{\textrm {cr}1}$ and $\beta _{\textrm {cr}2}$ or only $\beta _{\textrm {cr}2}$.

As was mentioned by Kotelnikov et al. (Reference Kotelnikov, Prikhodko and Yakovlev2023), the LoDestro equation (A2) with boundary conditions (A11) or (A12) constitutes the standard Sturm–Liouville problem. At first glance, it may seem that the solution of such a problem is rather standard. However, the equation has the peculiarity that its coefficients could be singular. In the anisotropic pressure model A1 relevant to the transverse NBI, the singularity appears near the minimum of the magnetic field in the limit $\beta \to 1$. In the case of oblique NBI models A2, A3, …, some coefficients of the LoDestro equations can have a singularity near point $z$ where the condition (3.7) of the mirror stability breaks at the plasma column axis.

For example, in the A3 model, the singularity occurs at $p_{0} \to p_{\textrm {mm}}=4$ in a magnetic field $b=3/4$, which corresponds to $b_{v}=\sqrt {27/32}$. Specifically, the coefficient at the function $\phi (z)$ in (A2) has a second-order pole, $-{217 b_{v}'(z)^2}/{36 (p-p_{\textrm {mm}})^{2}}$, and the coefficient at the first derivative $\phi '(z)$ has a first-order pole, $11b_{v}'(z)/9\sqrt {6}(p-p_{\textrm {mm}})$. For the magnetic field model (2.3), the coordinate of the singular point can be found explicitly

(A15)\begin{equation} z_{\text{mm}} = \frac{2}{{\rm \pi} } \sin^{{-}1}\left( \left( \frac{3 \sqrt{6} R-8}{8M-8} \right)^{{1}/{q}} \right) . \end{equation}

Reproducing the method of § 5 of Kotelnikov et al. (Reference Kotelnikov, Prikhodko and Yakovlev2023), we denote by $z_{R}$ the coordinate of the turning point on the $z$ axis, where $b=b_{v}=1$ and $p_{\mathrel \perp } =p_{\|}=0$. For the magnetic field model (2.3) we have

(A16)\begin{equation} z_{R} = \frac{2}{{\rm \pi} }\, \arcsin\left( \left( \frac{R-1}{M-1} \right)^{{1}/{q}} \right) . \end{equation}

In the adopted models of an anisotropic plasma, its pressure is zero in the region $z_{R}< z<1$ between the turning point $b=b_{v}=1$ and the magnetic mirror throat $b=b_{v}=M/R$. At zero pressure, the LoDestro equation (A2) takes an extremely simple form

(A17)\begin{equation} 0 = \frac{\mathrm{d}{}}{\mathrm{d}{z}} \left[ \left( \varLambda + 1 \right) \frac{\mathrm{d}{\phi}}{\mathrm{d}{z}} \right] . \end{equation}

Its solution is the equality

(A18)\begin{equation} \left[ \varLambda(z) + 1 \right] \phi'(z) =\operatorname{const.}, \end{equation}

where the constant on its right-hand side can be found from the boundary condition at $z=z_{E}$.

In the case of a plasma with electrically isolated ends, the boundary condition (A12) at $z=z_{E}$ should be applied, so that $\phi '(z_{E})=0$. Hence, the constant on the right-hand side of (A18) is also zero. Since the factor $\varLambda (z) + 1$ is greater than zero everywhere, we conclude that $\phi '(z)=0$ in the entire region $z_{R} < z < z_{E}$. Thus, it is sufficient to find a numerical solution of the original equation (A2) in the region $0 < z < z_{R}$, but one must be careful when setting the boundary conditions for $z=z_{R}$.

It should be taken into account that the derivative $\phi '(z)$ undergoes a jump at $z=z_{R}$. Indeed, integrating (A2) over an infinitesimal neighbourhood of the point $z_{R}$ from $z_{R}^{-}$ to $z_{R}^{+}$ yields the equation

(A19)\begin{equation} \left[ \varLambda(z_{R}) + 1 \right] \left[ \phi'({z_{R}^{+}}) - \phi'({z_{R}^{-}}) \right] = \left[ Q({z_{R}^{+}}) - Q({z_{R}^{-}}) \right] \phi(z_{R}) , \end{equation}

in which we took into account that $\varLambda$, $\phi$ and $\left \langle \bar {p}\right \rangle$ are continuous at the point $z=z_{R}$, in contrast to the derivative $\phi '(z)$ and the coefficient

(A20)\begin{equation} Q(z) = \frac{B_{v}'}{B_{v}} + \frac{2a'}{a} = \frac{2a'}{a} - \frac{2a_{v}'}{a_{v}} . \end{equation}

The jump in the coefficient is due to the fact that, for $b=b_{v}=1$, the derivative of functions such as (3.18) jumps (only in the case of pressure models A1 and A2). Since $\phi '({z_{R}^{+}})=0$ and $Q({z_{R}^{+}})=0$, (A19) yields the value $\phi '(z_{R}^{-})$ that the derivative of $\phi '(z)$ must have at the point $z_{R}^{-}$ on the right boundary of the interval $0< z< z_{R}$ from its inner side

(A21)\begin{equation} \phi'({z_{R}^{-}}) = \frac{ Q({z_{R}^{-}}) }{ \varLambda(z_{R}) + 1 }\, \phi(z_{R}) . \end{equation}

When solving (A2) on the interval $0< z< z_{R}^{-}$, the boundary condition (A21) should be used instead of (A12). Note that the $z_{E}$ coordinate is not included in (A21).

For all pressure models parameter $Q({z_{R}^{-}})$ can be found in analytic form. Namely, for model A1

(A22)\begin{gather} Q_{\infty}^{-} ={-} \frac{2 p}{(1-2 p)} b_{v}'(1) , \end{gather}
(A23)\begin{gather} Q_{1}^{-} = \left( 1 + \frac{\log (1-2 p)}{2p} \right) b_{v}'(1), \end{gather}
(A24)\begin{gather} Q_{2}^{-} = \frac{1}{2} \left(\frac{\sqrt{2} \sqrt{\dfrac{1}{p}-2} \sin^{{-}1}\left(\sqrt{2} \sqrt{p}\right)}{2 p-1}+2\right) b_{v}'(1) , \end{gather}
(A25)\begin{gather}\begin{aligned} Q_{4}^{-} & = \frac{1}{2} \left( \left(\frac{2p}{2 p-1}\right)^{3/4} \left( \tan^{{-}1} \left( \sqrt[4]{\frac{2p}{2 p-1}} \right) + \tanh^{{-}1} \left( \sqrt[4]{\frac{2p}{2p-1}} \right) \right) \right.\nonumber\\ & \quad \left. +\, \left\{\vphantom{\left(1+\frac{1+i}{\sqrt[4]{1-\dfrac{1}{2 p}}} \right)} 8\sqrt[4]{p}+\sqrt[4]{2-4 p} \log\left( \sqrt{1-2 p} + \sqrt{2p} - 2^{3/4}\sqrt[4]{(1-2 p) p } \right) \right. \right.\nonumber\\ & \quad \left. \left. -\, \sqrt[4]{2-4 p} \log\left( \sqrt{1-2 p} + \sqrt{2p} + 2^{3/4}\sqrt[4]{(1-2p)p} \right) \right. \right.\nonumber\\ & \quad \left. \left. -\, 2 ({-}1)^{3/4} \sqrt[4]{4p-2} \tan^{{-}1}\left( 1-\frac{1+i}{\sqrt[4]{1-\dfrac{1}{2 p}}} \right) \right. \right.\nonumber\\ & \quad \left. \left.\left. +\, 2({-}1)^{3/4} \sqrt[4]{4 p-2} \tan^{{-}1}\left( 1+\frac{1+i}{\sqrt[4]{1-\dfrac{1}{2 p}}} \right) \right\} \right/ \left\{ 4\sqrt[4]{p} \right\} \right) b_{v}'(1);\end{aligned} \end{gather}

for model A2

(A26)\begin{gather} Q_{\infty}^{-} = \frac{p}{p-1} b_{v}'(1), \end{gather}
(A27)\begin{gather}Q_{1}^{-} = \frac{p+\log (1-p)}{p} b_{v}'(1), \end{gather}
(A28)\begin{gather}Q_{2}^{-} = \frac{1}{2} \left( 2-\frac{2 \sin^{{-}1}\left(\sqrt{p}\right)}{\sqrt{(1-p) p}} \right) b_{v}'(1) , \end{gather}
(A29)\begin{gather}\begin{aligned} Q_{4}^{-} & = \frac{1}{2} \left\{\vphantom{\left(\sqrt[4]{\frac{p}{p-1}}\right)} \left[\vphantom{\left(\sqrt[4]{\frac{p}{p-1}}\right)} 8 \sqrt[4]{p} + \sqrt{2} \sqrt[4]{1-p} \log \left( \sqrt{1-p}+\sqrt{p}-\sqrt{2} \sqrt[4]{(1-p) p} \right) \right. \right.\nonumber\\ & \quad \left. \left. -\, \sqrt{2} \sqrt[4]{1-p} \log\left( \sqrt{1-p}+\sqrt{p}+\sqrt{2} \sqrt[4]{(1-p) p} \right) \right. \right.\nonumber\\ & \quad \left. \left. +\, 2 \sqrt{2} \sqrt[4]{1-p} \tan^{{-}1}\left( 1-\sqrt{2}\sqrt[4]{\frac{p}{1-p}} \right) \right. \right.\nonumber\\ & \quad \left.\left. \left. -\, 2 \sqrt{2} \sqrt[4]{1-p} \tan^{{-}1}\left( \sqrt{2} \sqrt[4]{-\frac{p}{p-1}}+1 \right) \right] \right/ { \left[ 4\sqrt[4]{p} \right] } \right.\nonumber\\ & \quad \left. +\, \left( \frac{p}{p-1}\right)^{3/4} \left( \tan^{{-}1}\left(\sqrt[4]{\frac{p}{p-1}}\right) + \tanh^{{-}1}\left(\sqrt[4]{\frac{p}{p-1}}\right) \right) \right\} b_{v}'(1),\end{aligned} \end{gather}

where

(A30)\begin{equation} b_{v}'(1) = \frac{{\rm \pi} q (R-1) }{2R}\, \sqrt{\left(\frac{M-1}{R-1}\right)^{2/q}-1} . \end{equation}

For model A3 and above

(A31)\begin{equation} Q_{\infty}^{-} = Q_{1}^{-} = Q_{2}^{-} = Q_{4}^{-} = 0 . \end{equation}

Let us move on to solving (A2) with boundary conditions (A13) and (A14) for $z=0$ and (A11) for $z=z_{E}$. As shown above, in the region $z_{R}< z< z_{E}$ the desired solution satisfies (A18) but the constant in this equation is now equal to $[ \varLambda (z_{R}) + 1 ] \phi '(z_{R}^{+})$, rather than zero. Hence, the derivative $\phi '(z)$ at $z=z_{R}^{+}$ is equal to

(A32)\begin{equation} \phi'(z_{R}^{+}) = \frac{- \phi(z_{R})} { \left[ \varLambda(z_{R}) + 1 \right] \displaystyle\int_{z_{R}}^{z_{E}}\dfrac{\mathrm{d}{z}}{\varLambda(z) + 1} } . \end{equation}

Substituting $\phi (z_{R}^{+})$ into (A19) and taking into account that $Q({z_{R}^{+}})=0$ yields the derivative $\phi '({z_{R}^{-}})$ on the right boundary of the interval $0< z< z_{R}$ from its inner side

(A33)\begin{equation} \phi'({z_{R}^{-}}) = \left[ \frac{Q({z_{R}^{-}})}{ \varLambda(z_{R}) + 1 \vphantom{\int_{z_{R}}^{z_{E}}} } - \frac{1} { \left[ \varLambda(z_{R}) + 1 \right] \displaystyle\int_{z_{R}}^{z_{E}}\dfrac{\mathrm{d}{z}}{\varLambda(z) + 1} } \right] \phi(z_{R}) . \end{equation}

This is the boundary condition that should be used instead of (A11) in the problem of ballooning instability with a combination of lateral wall stabilization and stabilization by conducting end plates. In the case of $\varLambda (z)=\operatorname {const.}$ under consideration in this paper it reduces to

(A34)\begin{equation} \phi'({z_{R}^{-}}) = \left[ \frac{ Q({z_{R}^{-}}) }{ \varLambda + 1 } - \frac{ 1 }{ z_{E} - z_{R} } \right] \phi(z_{R}) . \end{equation}

If a conducting limiter is installed before the magnetic plug at mirror ratio $L$, $R< L< M$, as in the RwPr or RwSt configurations, then the limiter coordinate $z_{E}$ is calculated by a formula similar to (A16), in which $R\to L$ should be replaced. If the conducting end plate is placed behind the magnetic plug as in BwSt or BwPr configuration, then

(A35)\begin{equation} z_{E} = 2 - \frac{2}{{\rm \pi} }\, \arcsin\left( \left( \frac{L-1}{M-1} \right)^{{1}/{q}} \right) . \end{equation}

The integrals in (A32) and (A33) were calculated analytically. For example, for the magnetic field model (3.28) with the index $q=2$, a relatively short expression is found

(A36)\begin{gather} \int_{z_{R}}^{1} \frac{\varLambda(z_{R}) + 1}{\varLambda(z) + 1} \mathrm{d}{z} = \frac{ 1 }{ {\rm \pi}\,\sqrt{M} \left(r_w^2-2\right) } \times \left[ -{\rm \pi}\left(\sqrt{M}\, r_w^2 \left(z_s-1\right)+M+1\right) + {\rm \pi}\,(M-1) \cos\left({\rm \pi} z_s\right) \right.\nonumber\\ \left. +\, 2 \left( M+1 - (M-1) \cos\left({\rm \pi} z_s\right) \right) \tan ^{{-}1}\left( \sqrt{M} \tan \left(\frac{{\rm \pi} z_s}{2}\right) \right) \right] , \end{gather}

where

(A37)\begin{equation} r_{w}^{2} = a_{0}^{2} \frac{\varLambda_{0}+1}{\varLambda_{0}-1}>2 \end{equation}

denotes the dimensionless square of the radius of the conducting wall in the model of a straightened cylinder. For magnetic field profiles with indices $q=4$ and $q=8$, much more cumbersome formulas are obtained, but they still speed up the calculation of integrals by several times as compared with numerical integration.

Footnotes

1 Equilibrium with $\beta >1$ is still possible in systems where the ionic Larmor radius is large in a certain sense, see Kotelnikov (Reference Kotelnikov2020), Kurshakov & Timofeev (Reference Kurshakov and Timofeev2023), Timofeev, Kurshakov & Berendeev (Reference Timofeev, Kurshakov and Berendeev2024) and Khristo & Beklemishev (Reference Khristo and Beklemishev2025).

2 A reference to a figure in another article is hereinafter supplemented with a superscript indicating the reference number in the bibliography.

References

Anikeev, A.V., Bagryansky, P.A., Deichuli, P.P., Ivanov, A.A., Karpushov, A.N., Maximov, V.V., Pod'minogin, A.A., Stupishin, N.V. & Tsidulko, Y.A. 1997 Observation of magnetohydrodynamic stability limit in a cusp-anchored gas-dynamic trap. Phys. Plasmas 4 (2), 347354.CrossRefGoogle Scholar
Bagryanskij, P.A., Ivanov, A.A., Klesov, V.V., Kotel'nikov, I.A., Krasnikov, Y.I., Roslyakov, G.V., Tsidulko, Y.A., Breun, R.A., Molvik, A.W. & Casper, T.A. 1990 Experimental MHD stability limit in the gas dynamic trap. In Plasma Physics and Controlled Nuclear Fusion Research, vol. 2, pp. 655–662. IAEA.Google Scholar
Bagryansky, P.A., Anikeev, A.V., Beklemishev, A.D., Donin, A.S., Ivanov, A.A., Korzhavina, M.S., Kovalenko, Y.V., Kruglyakov, E.P., Lizunov, A.A., Maximov, V.V., et al. 2011 Confinement of hot ion plasma with $\beta =0.6$ in the gas dynamic trap. Fusion Sci. Technol. 59 (1 T), 3135.CrossRefGoogle Scholar
Bagryansky, P.A., Beklemishev, A.D. & Soldatkina, E.I. 2007 Influence of radial electric field on high-beta plasma confinement in the gas dynamic trap. Fusion Sci. Technol. 51 (2 T), 340342.CrossRefGoogle Scholar
Bagryansky, P.A., Lizunov, A.A., Zuev, A.A., Kolesnikov, E.Y. & Solomachin, A.L. 2003 Experiments with controllable application of radial electric fields in GDT central cell. Fusion Sci. Technol. 43 (1 T), 152156.CrossRefGoogle Scholar
Baker, D.R., Garner, H.R., Parks, P.B., Sleeper, A.M., Okamura, S., Adati, K., Aoki, T., Fujita, H., Hidekuma, S., Hattori, K., et al. 1984 Stability studies of a hollow plasma in the double cusp experiment. Phys. Fluids 27 (11), 27232729.CrossRefGoogle Scholar
Beklemishev, A.D. 2008 Shear-flow effects in open traps. AIP Conf. Proc. 1069 (1), 314.CrossRefGoogle Scholar
Beklemishev, A.D. 2016 Diamagnetic “bubble” equilibria in linear traps. Phys. Plasmas 23 (8), 082506.CrossRefGoogle Scholar
Bushkova, O.A. & Mirnov, V.V. 1986 Influence of the configuration of the magnetic field on the MHD stability of the gas-dynamic trap. VANT [Questions Atomic Sci. Tech., Ser. Nucl. Fus.] (2), 1924, (in Russian).Google Scholar
D'Ippolito, D.A. & Hafizi, B. 1981 Low-$m$ ballooning stability of an axisymmetric sharp-boundary tandem mirror. Phys. Fluids 24 (12), 22742279.CrossRefGoogle Scholar
D'Ippolito, D.A. & Myra, J.R. 1984 Stability of mirrors with inverted pressure profiles. Phys. Fluids 27 (9), 22562263.CrossRefGoogle Scholar
Endrizzi, D., Anderson, J., Brown, M., Egedal, J., Geiger, B., Harvey, R., Ialovega, M., Kirch, J., Peterson, E., Petrov, Y., et al. 2023 Physics basis for the wisconsin hts axisymmetric mirror (WHAM). J. Plasma Phys. 89 (5), 975890501.CrossRefGoogle Scholar
Grad, H. 1967 Toroidal containment of a plasma. Phys. Fluids 10 (1), 137154.CrossRefGoogle Scholar
Ilgisonis, V.I. 1993 a Guiding-center theory for three-dimensional collisionless finite Larmor radius plasmas. Phys. Fluids B: Plasma Phys. 5 (7), 23872397.CrossRefGoogle Scholar
Ilgisonis, V.I. 1993 b Guiding-center theory for three-dimensional collisionless finite larmor radius plasmas. Phys. Fluids B: Plasma Phys. 5 (7), 23872397.CrossRefGoogle Scholar
Jiang, W., Verscharen, D., Li, H., Wang, C. & Klein, K.G. 2022 Whistler waves as a signature of converging magnetic holes in space plasmas. Astrophys. J. 935 (2), 169.CrossRefGoogle Scholar
Kalsrud, R. 1983 Magnetohydrodynamic Description of Plasma, vol. 1, p. 115. North-Holland.Google Scholar
Kesner, J. 1985 Axisymmetric, wall-stabilized tandem mirrors. Nucl. Fusion 25 (3), 275282.CrossRefGoogle Scholar
Khristo, M.S. & Beklemishev, A.D. 2025 Plasma equilibrium in diamagnetic trap with neutral beam injection. J Plasma Phys 91 (1), E3.CrossRefGoogle Scholar
Kotelnikov, I.A. 2011 Equilibrium of a high-$\beta$ plasma with sloshing ions above the mirror instability threshold. Fusion Sci. Technol. 59 (1 T), 4750.CrossRefGoogle Scholar
Kotelnikov, I. 2020 On the structure of the boundary layer in a beklemishev diamagnetic bubble. Plasma Phys. Control. Fusion 62 (7), 075002.CrossRefGoogle Scholar
Kotelnikov, I.A. 2021 Magnetic Hydrodynamics, 4th edn., Lectures on Plasma Physics, vol. 2. Lan, (in Russian).Google Scholar
Kotelnikov, I.A., Bagryansky, P.A. & Prikhodko, V.V. 2010 Formation of a magnetic hole above the mirror-instability threshold in a plasma with sloshing ions. Phys. Rev. E 81, 067402.CrossRefGoogle Scholar
Kotelnikov, I., Chen, Z., Bagryansky, P., Sudnikov, A., Zeng, Q., Yakovlev, D., Wang, F., Ivanov, A. & Wu, Y. 2020 Summary of the 2nd international workshop on gas-dynamic trap based fusion neutron source (GDT-FNS). Nucl. Fusion 60 (6), 067001.CrossRefGoogle Scholar
Kotelnikov, I., Lizunov, A. & Zeng, Q. 2022 a On the stability of small-scale ballooning modes in axisymmetric mirror traps. Plasma Sci. Technol. 24 (1), 015102.CrossRefGoogle Scholar
Kotelnikov, I., Prikhodko, V. & Yakovlev, D. 2023 Wall stabilization of high-beta anisotropic plasmas in an axisymmetric mirror trap. Nucl. Fusion 63 (6), 066027.CrossRefGoogle Scholar
Kotelnikov, I., Zeng, Q., Prikhodko, V., Yakovlev, D., Zhang, K., Chen, Z. & Yu, J. 2022 b Wall stabilization of the rigid ballooning $m = 1$ mode in a long-thin mirror trap. Nucl. Fusion 62 (9), 096025.CrossRefGoogle Scholar
Kruskal, M.D. & Oberman, C.R. 1958 On the stability of plasma in static equilibrium. Phys. Fluids 1 (4), 275280.CrossRefGoogle Scholar
Kurshakov, V.A. & Timofeev, I.V. 2023 The role of electron current in high-$\beta$ plasma equilibria. Phys. Plasmas 30 (9), 092513.CrossRefGoogle Scholar
Kuz'min, S.V. & Lysyanskij, P.B. 1990 MHD stability margin of a cusp and nonparaxial magnetic mirror against rigid mode. Fiz. Plazmy 16 (8), 10011010, (in Russian).Google Scholar
Lansky, I.M. 1993 On the paraxial equilibrium of the finite $\beta$ plasma in open magnetic configuration. Tech. Rep. BudkerINP 93-96. Budker Institute of Nuclear Physics.Google Scholar
Li, X.-Z., Kesner, J. & Lane, B. 1985 Conducting-wall and pressure profile effect on MHD stabilization of axisymmetric mirror. Nucl. Fusion 25 (8), 907917.CrossRefGoogle Scholar
Li, X.Z., Kesner, J. & Lane, B. 1987 a MHD stabilization of a high beta mirror plasma partially enclosed by a conducting wall. Nucl. Fusion 27 (1), 101107.CrossRefGoogle Scholar
Li, Z.-Z., Kesner, J. & LoDestro, L. 1987 b Wall stabilized high beta mirror plasma in a rippled magnetic field. Nucl. Fusion 27 (8), 12591266.CrossRefGoogle Scholar
Li, Q., Zhu, G., Ren, B., Ying, J., Yang, Z. & Sun, X. 2023 Experimental studies of cusp stabilization in Keda Mirror with AXisymmetricity (KMAX). Plasma Sci. Technol. 25 (2), 025102.CrossRefGoogle Scholar
LoDestro, L.L. 1986 The rigid ballooning mode in finite-beta axisymmetric plasmas with diffuse pressure profiles. Phys. Fluids 29 (7), 23292332.CrossRefGoogle Scholar
Logan, B.G. 1980 An axisymmetric high-beta tandem-mirror reactor. Comments Plasma Phys. Control. Fusion 5 (6), 271274.Google Scholar
Logan, B.G. 1981 Improved axisymmetric-cusp plugs for tandem mirror reactors. Comments Plasma Phys. Control. Fusion (United Kingdom) 6 (6), 199207.Google Scholar
Lotov, K.V. 1996 Spontaneous formation of zero magnetic field region near the axis of a high-$\beta$ mirror device. Phys. Plasmas 3 (4), 14721473.CrossRefGoogle Scholar
Mirnov, V.V. 1986 Theory of open traps with a short free path of ions. PhD thesis, Budker INP, Novosibirsk.Google Scholar
Nagornyi, V. & Stupakov, G. 1984 Flute instability of plasma in a cusp. Sov. J. Plasma Phys. 10, 275.Google Scholar
Newcomb, W.A. 1981 Equilibrium and stability of collisionless systems in the paraxial limit. J. Plasma Phys. 26 (3), 529584.CrossRefGoogle Scholar
Pantellini, F.G.E. 1998 A model of the formation of stable nonpropagating magnetic structures in the solar wind based on the nonlinear mirror instability. J. Geophys. Res.: Space Phys. 103 (A3), 47894798.CrossRefGoogle Scholar
Pokhotelov, O.A., Treumann, R.A., Sagdeev, R.Z., Balikhin, M.A., Onishchenko, O.G., Pavlenko, V.P. & Sandberg, I. 2002 Linear theory of the mirror instability in non-Maxwellian space plasmas. J. Geophys. Res.: Space Phys. 107 (A10), SMP 18-1SMP 18–11.CrossRefGoogle Scholar
Polyanin, A.D. & Manzhirov, A.V. 2008 Handbook of Integral Equations. CRC Press.CrossRefGoogle Scholar
Rudakov, L.I. & Sagdeev, R.Z. 1961 On the instability of a nonuniform rarefied plasma in a strong magnetic field. Dokl. Akad. Nauk 138 (3), 581583.Google Scholar
Ryutov, D.D., Berk, H.L., Cohen, B.I., Molvik, A.W. & Simonen, T.C. 2011 Magneto-hydrodynamically stable axisymmetric mirrors. Phys. Plasmas 18 (9), 092301.CrossRefGoogle Scholar
Shmigelsky, E.A., Lizunov, A.A., Meyster, A.K., Pinzhenin, E.I., Solomakhin, A.L. & Yakovlev, D.V. 2024 Recent results on high-$\beta$ plasma confinement studies in the gas dynamic trap. J. Plasma Phys. 90 (2), 975900206.CrossRefGoogle Scholar
Skovorodin, D.I., Zaytsev, K.V. & Beklemishev, A.D. 2013 Global sound modes in mirror traps with anisotropic pressure. Phys. Plasmas 20 (10), 102123.CrossRefGoogle Scholar
Soldatkina, E., Anikeev, M., Bagryansky, P., Korzhavina, M., Maximov, V., Savkin, V., Yakovlev, D., Yushmanov, P. & Dunaevsky, A. 2017 Influence of the magnetic field expansion on the core plasma in an axisymmetric mirror trap. Phys. Plasmas 24 (2), 022505.CrossRefGoogle Scholar
Soldatkina, E., Maximov, V., Prikhodko, V., Savkin, V., Skovorodin, D., Yakovlev, D. & Bagryansky, P. 2020 Measurements of axial energy loss from magnetic mirror trap. Nucl. Fusion 60 (8), 086009.CrossRefGoogle Scholar
Southwood, D.J. & Kivelson, M.G. 1993 Mirror instability. 1. Physical mechanism of linear instability. J. Geophys. Res.: Space Phys. 98 (A6), 91819187.CrossRefGoogle Scholar
Stupakov, G.V. 1987 Possibility of MHD-stable plasma confinement in an axisymmetric mirror system. Preprint 1987-014. Budker INP, Novosibirsk, in Russian.Google Scholar
Taylor, J.B. 1963 Some stable plasma equilibria in combined mirror-cusp fields. Phys. Fluids 6 (11), 15291536.CrossRefGoogle Scholar
Thompson, W.B. 1964 An Introduction to Plasma Physics, 2nd edn. Pergamon Press.Google Scholar
Timofeev, I.V., Kurshakov, V.A. & Berendeev, E.A. 2024 Formation of cylindrical plasma equilibria with $\beta > 1$. Phys. Plasmas 31 (8), 082512.CrossRefGoogle Scholar
Winterhalter, D., Neugebauer, M., Goldstein, B.E., Smith, E.J., Bame, S.J. & Balogh, A. 1994 Ulysses field and plasma observations of magnetic holes in the solar wind and their relation to mirror-mode structures. J. Geophys. Res.: Space Phys. 99 (A12), 2337123381.CrossRefGoogle Scholar
Yurov, D.V., Prikhodko, V.V. & Tsidulko, Y.A. 2016 Nonstationary model of an axisymmetric mirror trap with nonequilibrium plasma. Plasma Phys. Rep. 42, 210225.CrossRefGoogle Scholar
Zeng, Q. & Kotelnikov, I. 2024 Influence of the shape of a conducting chamber on the stability of rigid ballooning modes in a mirror trap. Plasma Phys. Control. Fusion 66 (7), 075020.CrossRefGoogle Scholar
Figure 0

Figure 1. Axial profile of the vacuum magnetic field (2.3) for mirror ratio $M=32$ and three values of the index $q$ indicated on the graphs as compared with the magnetic field in the gas-dynamic trap (GDT, Bagryanskij et al.1990) and the Wisconsin HTS axisymmetric mirror (WHAM, Endrizzi et al.2023) devices.

Figure 1

Figure 2. Axial profiles of the transverse pressure in models A1 (a), A2 (b) and A3 (c) described respectively by (3.10), (3.18) and (3.28) in magnetic field (2.3) at the axis of the trap for mirror ratio $M=8$, index $q=4$ and various mirror ratios $R$ at the turning points where the pressure of hot plasma component drops to zero. The value $\beta =0.3$ for this figure is chosen so as not to exceed the threshold for excitation of the mirror and firehose instabilities for all values of the parameter $R$ indicated in the figure legend.

Figure 2

Figure 3. Threshold of mirror and firehose instabilities in an anisotropic plasma with oblique NBI within the framework of models A2 (a) and A3 (b). Unstable areas are shaded. The same shading scheme is used below in the stability maps without further reminder.

Figure 3

Figure 4. Axial profile of the local beta $\beta (z) = 2p_{\mathrel \perp }/B_{v}^{2}$ in the magnetic field (2.3) for mirror ratio $M=8$, index $q=4$ and various mirror ratios $R$ at the turning points where the pressure of the hot plasma component drops to zero: (a) transversal NBI, pressure model A1; (b) oblique NBI, model A2; (c) oblique NBI, model A3. The values of the $\beta$ parameter for each $R$ value indicated on the graphs are chosen to be equal to the smallest of the two stability thresholds for the mirror and firehose instabilities.

Figure 4

Figure 5. Angular distribution of the fast ions in the A1, A2, A3 and A8 models for a set of parameters $R$ indicated in the graphs.

Figure 5

Figure 6. Second critical beta vs $R$ in the limit $\varLambda \to \infty$ for three pressure models: A1 (a,d,g,j), A2 (b,e,h,k) and A3 (c,f,i,l). Stability zone of the ballooning rigid mode for a radial pressure profile with a given index $k$ is located above the margin curve, coloured as indicated in the legend under (jl).

Figure 6

Figure 7. Second critical beta vs $R$ in the limit $\varLambda \to \infty$ for three pressure models: A1 (a,d,g,j), A2 (b,e,h,k) and A3 (c,f,i,l). Stability zone of the ballooning rigid mode for a radial pressure profile with a given index $k$ is located above the margin curve, coloured as indicated in the legend under (jl).

Figure 7

Figure 8. Second critical beta $\beta _{\textrm {cr}2}$ vs Kesner's degree of anisotropy $A_{{K}}$ defined by (3.29) as in Kesner (1985) for three models of anisotropic pressure and four mirror ratios, $M\in \{4,8,16\}$. Dashed curve shows the mirror mode thresholds. Compare with figure 3 in Kesner (1985).

Figure 8

Figure 9. Stability maps for the LwPr configuration and three pressure models: A1 (a,d,g,j,m), A2 (b,e,h,k,n) and A3 (c,f,i,l,o). The second critical beta is drawn as a function of $r_{w}/a$ for the model magnetic field (2.3) with index $q=4$, set of mirror ratios $R\in \{1.1,1.2,1.5,2,4,8\}$ at the turning point and fixed mirror ratio $M=8$. The stable zone for a radial pressure profile with an index $k$ is located above the curve $\beta _{\textrm {cr}2}$ coloured according to the legend under (mo). Compare with figure 10.

Figure 9

Figure 10. Stability maps for the LwSt configuration and three pressure models: A1 (a,d,g,j,m), A2 (b,e,h,k,n) and A3 (c,f,i,l,o). The second critical beta is drawn as a function of $r_{w}/a$ for the model magnetic field (2.3) with index $q=4$, set of mirror ratios $R\in \{1.1,1.2,1.5,2,4,8\}$ at the turning point and fixed mirror ratio $M=8$. The stable zone for a radial pressure profile with an index $k$ is located above the curve $\beta _{\textrm {cr}2}$ coloured according to the legend under (mo). Compare with figure 9.

Figure 10

Figure 11. First critical beta $\beta _{\textrm {cr}1}$ for the CwPr configuration and three pressure models: A1 (ac), A2 (df) and A3 (gi) at mirror ratio $M=15$ vs parameter $R$ in the limit $\varLambda \to 1$. The stability zone of rigid ballooning perturbation for a radial pressure profile with a given index $k$ is located below the curve, coloured as indicated in the legend under (gi).

Figure 11

Figure 12. Same maps as in figure 11 but for mirror ratio $M=8$.

Figure 12

Figure 13. Same maps as in figure 11 but for mirror ratio $M=4$. The rows A2-CwPr and A3-CwPr are not shown as they do not have an unstable zone for rigid ballooning modes.

Figure 13

Figure 14. Stability maps for model magnetic field (2.3) and anisotropic plasma pressure models (3.10) (af), (3.18) (gl) and (3.28) (mr) at combined MHD stabilization by lateral wall and end MHD anchors, $q\in \{2, 4, 8\}$, $M=R=16$. The unstable zone is located between the lower $\beta _{\textrm {cr}1}(r_{w}/a_{0})$ and upper branches $\beta _{\textrm {cr}2}(r_{w}/a_{0})$ of every curve. Correspondence of the index $k$ to the colour of the curves is shown at the bottom of the figure. Shaded common zone of stability lies to the left of the curve for the most steep radial pressure profile ($k = \infty$). (a,d,g,j,m,p) Show the maps for a ‘parabolic’ magnetic field with the index $q=2$, (c,f,i,l,o,r) show the maps for the ‘quasi-flat hole’ magnetic field with $q=8$. Panels (ac,gi,mo) and (df,jl,pr) show maps for the CwPr and CwSt configurations, respectively.

Figure 14

Figure 15. Same as in figure 14 but for $M=R=8$.

Figure 15

Figure 16. Same as in figure 14 but for $M=R=4$. Rows for the A2 and A3 models are dropped since the rigid ballooning mode is stable in the entire region below the mirror instability threshold.

Figure 16

Figure 17. Stability maps vs ratio $r_{w}/a$ for the A1, A2 and A3 pressure models in the CwPr configuration simulating combined stabilization by the proportional lateral conducting chamber and end MHD anchors; $q=4$, $M=16$, $R\in \{1.2, 1.5, 2, 4, 8\}$. The instability zone is located between $\beta _{\textrm {cr}1}(r_{w}/a)$ (lower branch of the marginal curve) and $\beta _{\textrm {cr}2}(r_{w}/a)$ (upper branch of the same colour); in the case of inclined injection, which corresponds to the A2 and A3 models, the upper branch is completely or partially absorbed by the region of mirror instability; the stability zone is shaded for a plasma with a sharp boundary ($k=\infty$), for which it has the minimum dimensions.

Figure 17

Figure 18. Same as in figure 17 but for CwSt configuration simulating combined stabilization by the straightened lateral conducting wall and the end MHD anchors.