Wall stabilization of high-beta anisotropic plasmas in an axisymmetric mirror trap

The stabilization of the ‘rigid’ flute and ballooning modes both with and without the effect of additional MHD anchors in an axisymmetric mirror trap, with the help of a perfectly conducting lateral wall, is studied using a model pressure distribution of a plasma with a steep-angle neutral beam injection at the magnetic field minimum. The calculations were performed for an anisotropic plasma in a model that simulates the pressure distribution during the injection of beams of fast neutral atoms into the magnetic field minimum. It is assumed that the lateral wall is an axisymmetric shell surrounding the plasma that follows the shape of the magnetic field line and is placed at a certain distance to the plasma border. It has been found that for the effective stabilization of the modes by the lateral wall only the parameter beta (β, ratio of plasma pressure to the magnetic field pressure) must exceed some critical value β crit. When combined with the conducting end plates imitating the MHD end stabilizers, there are two critical beta values and two stability zones 0<β<βcrit1 and βcrit2<β<1 that can merge, making the entire range of allowable beta values 0<β<1 stable. The dependence of the critical betas on the degree of plasma anisotropy, the mirror ratio, and the width of the vacuum gap between the plasma and the lateral wall is studied. In contrast to the previous studies focusing on a plasma model with a sharp boundary, we calculated the stability zones for a number of diffuse radial pressure profiles and several axial magnetic field profiles.


Introduction
Continuing the study of ballooning instabilities in isotropic plasma started in our article [1], in this paper we present * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. the results of the study of the so-called wall stabilization of rigid flute and ballooning perturbations with azimuthal mode number m = 1 in a mirror trap (also called an open trap or a linear trap) with anisotropic plasma. In its essence, the wall stabilization of plasma with sufficiently high pressure is achieved due to excitation of image currents in the conducting lateral walls that can either be a part of the vacuum vessel or be specially designed conducting components surrounding the plasma. These currents are in the opposite direction to the diamagnetic current in the plasma, and the opposite currents, as is known, are pushed apart. Such pushing returns the pop-up 'tongues' of plasma back towards the trap axis.
In order to emphasize the significance of stabilization of a single m = 1, we should recall what constitutes a ballooning instability in an open axially symmetric trap. It is generally accepted to think of ballooning perturbations as small-scale deformations of plasma equilibrium with a large azimuthal mode number m ≫ 1. It has long been known that the finite Larmor radius (FLR) effects can stabilize small-scale flutetype perturbations if m > 1 [2]. Ballooning perturbations are, in fact, related to flute perturbations; therefore, it is commonly believed that they are also stabilized by the FLR effects.
In their pure form, flute perturbations can be detected in a low-pressure plasma (that is, in the limit of β → 0, where β is the ratio of plasma pressure to the magnetic field pressure). They represent a magnetic flux tube, which, together with the plasma captured in it, floats to the periphery of the plasma column without changing its shape and without distorting the magnetic field. The classical method for stabilizing flute disturbances is to attach the end MHD anchors (for example, of the cusp type) to the central section of on open trap, which 'clamp' the ends of the flux tube. The same 'clamping' effect is achieved by placing conductive plates at the ends of the plasma column. If the ends of the flux tube are 'frozen' into the end plates, the tube cannot float without bending and deforming the magnetic field. Thus, the flute perturbation is transformed into a balloon-type perturbation. Deformation of the magnetic field requires energy, which must be withdrawn from the thermal energy of the plasma. Therefore, ballooning instability is possible if plasma beta exceeds a certain threshold value; it is important to emphasize that this threshold value in axially symmetric open traps is very high, on the order of 60%-70% [3][4][5].
The FLR effects are not capable of stabilizing flute and ballooning perturbations with m = 1. However, they impose rigidity on such perturbations, so that the plasma density distribution does not change in each cross section of the plasma column. The plasma column bends, shifting from the trap axis to different distances in different sections. Stabilization of such a 'rigid' mode m = 1 would mean stabilization of all flute and ballooning perturbations, provided that the modes m > 1 are stabilized by the FLR effects.
Unlike many other studies on the stability of the rigid ballooning mode m = 1 [6][7][8][9][10][11][12][13][14][15][16], which were focused exclusively on a plasma model with a stepwise pressure profile along the radius, in our article [1] we studied the stability of the m = 1 mode in a plasma with a diffuse pressure profile, but at that time limited ourselves to the case of an isotropic plasma. Specifically, we calculated the so-called critical beta β crit , such that the rigid flute and ballooning modes m = 1 can be stabilized by the lateral conducting wall without using any other MHD stabilization methods if β > β crit . A second case was also analyzed when the lateral conductive wall was supplemented with conducting end plates installed in the magnetic mirror throats, which simulated the end MHD stabilizers, such as cusp. In this case, two stability zones were found: the first one at small beta, β < β crit1 , the second one at large beta, β > β crit2 , and these two zones can merge.
In this article, we study wall stabilization on one particular model of anisotropic plasma, which allows us to perform a significant part of the calculations in an analytical form. The key equation for solving the problem was derived by LoDestro [17] in 1986, but, it seems that neither she nor anyone else has ever used this equation. We found a few papers [18][19][20][21] that refer to LoDestro's paper but do not use her equation. LoDestro's work was probably in oblivion for many years due to an early termination of the TMX (Tandem Mirror Experiment) and MFTF-B (Mirror Fusion Test Facility B) projects in the USA in the same year [22]. However, the achievement of high electron temperature and high beta in the gas dynamic trap (GDT) at Budker Institute of Nuclear Physics (INP) in Novosibirsk [23][24][25][26][27][28][29][30][31][32], emergence of new ideas [33], and new projects [34][35][36][37][38] makes us revisit old results.
Another reason for the loss of interest in the LoDestro equation could be the insufficient power of computers of that era. It is possible to calculate the critical beta for a plasma with a sharp boundary without this equation, which was done by LoDestro's predecessors. As the plasma model becomes more complex, the demands on computer performance and resources increase rapidly. It takes less than an hour on a modern desktop PC with eight cores to recalculate all the graphs presented in our article [1] for the isotropic plasma model. But calculations for some anisotropic plasma models, not yet published, take weeks and even months of continuous work.
To avoid unnecessary repetitions, we will not analyze the content of the articles cited above, since their detailed review was included in our previous article [1], and we will immediately proceed to the description of the new calculations. In section 2, LoDestro's equation will be written and the necessary notation will be introduced. In section 3, the anisotropic pressure model used below will be formulated. The results of calculating some coefficients in the LoDestro equation for this model are presented in appendices A and B. Section 4 presents the results of calculating the critical beta in the limit when the conducting wall surrounding the plasma column almost conforms to the lateral boundary of the plasma column, but does not touch it. In this limit, LoDestro's ordinary differential equation of second order reduces to an integral along the z coordinate on the trap axis; the integral turns to zero at the critical beta value. Section 5 describes the solution of the LoDestro equation by the shooting method and presents the results of calculations for several model profiles of pressure and magnetic field. In section 6, the shooting method is again used to solve the LoDestro equation with other boundary conditions that imitate the effect of conducting end plates or MHD anchors installed in magnetic plugs. The final section 7 summarizes our results.

LoDestro equation
The LoDestro equation is a second-order ordinary differential equation for the function: which depends on one coordinate z along the trap axis and is expressed in terms of the variable radius of the plasma column boundary a = a(z), vacuum magnetic field B v = B v (z), and virtual small displacement ξ n = ξ n (z) of the plasma column from the axis. It is obtained on the assumption that • FLR effects are strong enough, m = 1; • the paraxial (long-thin) approximation applies; • the pressure tensor depends only on two functions: The LoDestro equation does not take into account the resistance of plasma and conductive walls. It is also possible that the assumption about the dominance of the FLR effects is violated near the boundary of the plasma column, which may result in an incorrect account of plasma drifts and rotation. Plasma rotation is taken into account by Beklemishev's theory of vortex confinement [39], but it is limited to the β = 0 approximation. Plasma resistance was taken into account by Kang et al [40] in an isotropic plasma model with a sharp boundary.
In its final form, the LoDestro equation reads: where the derivative d/dz in the first two lines acts on all factors to the right of it, the prime ( ′ ) is a shortcut for d/dz, and ω is the oscillation frequency. Other notations are defined as follows: Parameter ψ a , used here has the meaning of the reduced (i.e. divided by 2π) magnetic flux through the plasma cross section π a 2 . It is related to the plasma radius a = a(z) by equation (3). Equation (4) relates the radial coordinate r and the magnetic flux ψ through a ring of radius r in the z plane. The magnetic field B = B(ψ, z), weakened by the plasma diamagnetism, in the paraxial (long-thin) approximation (i.e. with a small curvature of the field lines) is related to the vacuum magnetic field B v = B v (z) by the transverse equilibrium equation (5).
Function Λ = Λ(z) is expressed in terms of the actual radius of the plasma/vacuum boundary a = a(z) and the radius of the conducting shell r w = r w (z), which surrounds the plasma column. In the remainder of the paper we assume that Λ is a constant. The case Λ(z) = Λ 0 = const will be referred to below as the proportional chamber. The shape of the conducting walls of such a chamber on larger scale repeats the shape of the plasma column, that is, r w /a = const ⩾ 1. Such a chamber is difficult to manufacture, but the assumption Λ(z) = const greatly simplifies the solution of the LoDestro equation. The larger the value of Λ, the closer the conducting lateral wall is to the plasma boundary. The limit Λ → ∞ corresponds to the case when the conducting lateral wall is as close as possible to the plasma boundary, repeating its shape, but not touching the plasma. The limit Λ → 1 means that the lateral conducting wall is removed to infinity.
Kinetic theory predicts (see, for example, [41]) that the transverse and longitudinal plasma pressures can be considered as functions of B and ψ, i.e. p ⊥ = p ⊥ (B, ψ), p ∥ = p ∥ (B, ψ). In equation (2), one must assume that the magnetic field B is already expressed in terms of ψ and z, and therefore p ⊥ = p ⊥ (ψ, z), p ∥ = p ∥ (ψ, z). The angle brackets in equation (2) denote the mean value of an arbitrary function of ψ and z over the plasma cross section. In particular, the average value ⟨ρ⟩ of the density ρ = ρ(ψ, z) is calculated using a formula similar to (8).
The boundary conditions for equation (2) and similar equations in the study of ballooning instability are traditionally set at the ends of the plasma column at the magnetic field maxima B v = B max , where B ′ v = 0 and p ⊥ = p ∥ = 0. In accordance with the geometry of actually existing open traps, it is usually assumed that the magnetic field is symmetrical with respect to the median plane z = 0, and the magnetic mirrors (i.e. field maxima) are located at z = ±L.
Traditionally, two types of boundary conditions are considered. In the presence of conducting end plates directly in magnetic mirrors, it is required that the boundary condition: be satisfied at z = ±L. A similar boundary condition is usually used in studying the stability of small-scale ballooning disturbances, thereby modeling the presence of a stabilizing cell behind a magnetic mirror (see, for example, [5]). If the plasma ends are electrically insulating, the boundary condition: is applied. Basically, it implies that other methods of MHD stabilization in addition to stabilization by a conducting lateral wall are not used. It is this boundary condition (11) that was used earlier in the works on the stability of the m = 1 ballooning mode. The LoDestro equation (2) with boundary conditions (10) or (11) represents a standard Sturm-Liouville problem. At the first glance, it may seem that the solution of such a problem is not difficult to find. However, the LoDestro equation has the peculiarity that its coefficients can be singular. In the anisotropic pressure model, which will be formulated in the next section, the singularity appears near the minimum of the magnetic field in the limit β → 1. There are indications that our predecessors were aware of the singularity problem. For example, in Kesner's article [13], the graphs in figure 2 break off at β ≈ 0.9. Unfortunately, he did not leave us recipes for dealing with this singularity.

Anisotropic pressure
In publications on the stability of the rigid ballooning mode, two models of anisotropic pressure have previously been used. Kesner in his paper [13] indicates that in the first model the transverse pressure in a nonuniform magnetic field B ⩽ B max varies according to the law: while in the second model: The first model approximately describes the pressure distribution in an open trap, which arises when beams of neutral atoms (NB) are injected into plasmas at a right angle to the trap axis at the magnetic field minimum. This is the so-called normal NB injection (NBI). The second model corresponds to the oblique NBI, which forms a population of so-called sloshing ions.
In the present paper, we will restrict ourselves to the first model for three reasons. First, it is relevant to the Compact Axisymmetric Toroid (CAT) at Budker INP [28,29,42] as well as the new Wisconsin HTS Axisymmetric Mirror (WHAM++) under development at the University of Wisconsin [37]. Second, in this model, the coefficients in the LoDestro equation can be calculated without using numerical integration. Finally, in the second model, the ultimate beta can actually be limited by either mirror or the firehose instabilities.
In our internal classification, the isotropic plasma variant is designated 'A0', the anisotropic pressure (12), which is formed during normal NBI, is designated 'A1'. The oblique injection simulated by the formula (13) with an arbitrary index n is denoted by 'An'. The two cases n = 2 and n = 3 will be dealt with in our next article and are abbreviated 'A2' and 'A3'. Variations with conducting lateral wall stabilization will be marked with the letters 'LW' (shortcut for lateral wall), and variations with combined lateral wall stabilization and conductive end plates installed in the throat of the magnetic plug will be marked with the letters 'CW' (shortcut for combined wall). The design of the lateral conductive wall in the form of a proportional chamber will be labeled 'Pr' (proportional).
Thus, the label 'A1-LWPr' in a figure corresponds to the variant of plasma stabilization with anisotropic pressure of the 'A1' type in a proportional conducting chamber without the effect of the end MHD stabilizers. The influence of the shape of the lateral wall on the stability of the rigid ballooning mode will be studied in a separate article. In particular, there will be a comparison of the stabilizing effect of a proportional chamber with a straight chamber, which is assigned the abbreviation 'St' (straight).
In kinetic theory [41], it is proved that if either of the two pressures p ⊥ and p ∥ is given as a function of B, then the other is uniquely determined using the parallel equilibrium equation. The latter can be rewritten in terms of the partial derivative with respect to B for constant magnetic flux ψ as: Another key result of the kinetic theory is the assertion that the function p ⊥ /B 2 always decreases as B increases, i.e.

∂ ∂B
Using equation (14), one can rewrite this last result as: First, for stability against the so-called 'firehose mode', it is required that: Second, stability against the so-called 'mirror mode' implies that: All these equations are satisfied by the functions 1 They describe the plasma pressure profile with a peak near the magnetic field minimum in the median plane of an open trap. Both functions p ⊥ and p ∥ simultaneously vanish at B = B s . We assume that B s does not exceed the magnetic field B max in the magnetic mirrors, bearing in mind that there is some cold plasma in the region B s < B < B max , but its pressure is 1 The inequalities (15)- (17) are satisfied if p 0 > 0. The inequality in equation (18) formally leads to the condition 2p 0 < B 2 s , while the transverse equilibrium equation (5) has a solution in the entire region Bv < Bs under the obviously more stringent condition 2p 0 < min( negligible. Such a profile approximately describes the expected pressure distribution in the CAT device, where the field B s approximately corresponds to the stop point of sloshing ions. For the pressure profile (19), equation (5) can be solved and the true magnetic field B weakened by the diamagnetic effect can be explicitly expressed in terms of the vacuum field B v : In the second case (13), equation (5) can be solved analytically only if n = 2 or n = 3, but in any case, numerical integration is necessary when calculating functions such as a and ⟨p⟩, while for the functions in equations (17)- (19) these integrals can be calculated analytically.
The pressure functions expressed in terms of B v will be denoted by the capital letter P, Turning to dimensionless variables, we will take the magnetic field at the turning point B = B v = B s as a unit of measurement, and the vacuum mirror ratio B s / min(B v ) at the turning point will be further denoted R. Then, keeping the same notation for dimensionless quantities as for dimensional ones, we have: The unit of length along the radius r is further fixed by the fact that we take the magnetic flux ψ a through the plasma section as a unit of measurement, assuming further that ψ a = 1. As for the coordinate z along the trap axis, it will be further normalized to the distance L between the median plane z = 0 and the magnetic mirror. With this normalization, it turns out that the plugs are located in the planes z = ±1. Below we represent the dependence of pressure on the magnetic flux ψ as: where the dimensionless function f k (ψ) is defined as: for integer values of index k, and for k = ∞ is expressed in terms of a θ-function such that θ(x) = 0 for x < 0 and θ(x) = 1 for x > 0: Parameter beta, β, is defined as the maximum of the ratio The maximum is reached on the trap axis (where ψ = 0) at the minimum of the vacuum field (where B v = 1/R), so that in dimensionless notation: Parameter p 0 can vary within the range 0 < p 0 < 1/2R 2 , and p 0 → 1/2R 2 at β → 1. At β > 1, plasma equilibrium is impossible, since equation (5) does not have a continuous solution.
As in our first article [1], the subsequent calculations were performed in the Wolfram Mathematica © for four values of the index k ∈ {1, 2, 4, ∞}. The f 1 function describes the smoothest pressure profile. For β/B 2 v ≪ 1, it approximately gives the parabolic dependence of the pressure p on the coordinate r. The larger the index k, the more table-like distribution p near the axis of the plasma column is and the steeper it is near the boundary of the column. Index k = ∞ corresponds to a pressure profile in the form of a step with a sharp boundary.
For the above radial pressure profiles, Wolfram Mathematica © was able to calculate the integrals in equations (3) and (8) analytically (with the help of some custom subroutines as described in appendix D). The results of the calculations are summarized in appendices A and B.
An important characteristic of the plasma is the degree of its anisotropy. If the latter is formally defined as: relating it to the plasma pressure at the minimum of the magnetic field on the trap axis, we obtain: This shows that the anisotropy is greater for smaller R, reaching a maximum at R → 1. In connection with this parameter R will sometimes be called the anisotropy parameter. The limit R → ∞ corresponds to the case of an isotropic plasma. Calculating this limit in the formulas collected in appendices A and B, one can re-derive the formulas for isotropic plasma published earlier in [1]. However, in this article we restrict ourselves to the case of R ⩽ M, since the plasma pressure should tend to zero at the mirror throat, where dimensionless is the mirror ratio. We also note that the parameter R characterizes the spatial width of the plasma pressure distribution along the magnetic field lines. Larger values of R correspond to wider the axial pressure profiles.

Thin vacuum gap
As the lateral conducting wall approaches the plasma/vacuum boundary, where it produces its maximum stabilizing effect, parameter Λ tends to infinity. In the limit Λ → ∞, it is possible to make analytic progress in solving equation (2), and, hence, in assessing the effects of a diffuse profile and anisotropy on the rigid ballooning mode stability. Stabilization of the m = 1 mode by a conducting wall in the Λ → ∞ limit was previously studied by Kesner [13] and Li et al [14] in the case of a plasma with a sharp-boundary radial profile. In particular, apart from isotropic plasma Kesner analyzed the anisotropic pressure model, equation (13), and concluded that the critical beta becomes smaller as the degree of anisotropy increases.
As shown in [1], in the limit Λ → ∞, the Sturm-Liouville problem, equations (2) and (11), reduces to solving the integral equation: It allows one to calculate the squared oscillation frequency ω 2 if the radial profile of pressure p = (p ⊥ + p ∥ )/2, density ρ, and vacuum magnetic field B v are given. At the margins of the stable regime, the oscillation frequency is equal to zero, ω 2 = 0. In the stability region ω 2 > 0, and instability takes place if ω 2 < 0.
In a plasma with anisotropic pressure, the condition in equation (17) together with equation (5) guarantees that p < B 2 v /2. Hence, the multiplier (1 − p/B 2 v ) > 0 and the second term in the square brackets are always positive. Therefore, the rigid ballooning mode in a plasma with anisotropic pressure can also be stabilized by a conducting wall located close enough to the lateral surface of the plasma column.
Further in this section, we present the results of calculation of the critical value of beta in anisotropic plasma in the limit Λ = ∞. The calculation method is described in detail in our previous article using the example of isotropic plasma [1]. Here, we restrict ourselves to indicating that the critical value β crit of the parameter β, corresponding to the marginal stability ω 2 = 0, is to be determined as a root of the equation: where W(β) denotes the integral on the right-hand-side of equation (34). Some peculiarities of searching for roots of this strongly nonlinear equation in the Wolfram Mathematica © are described in appendix D.
As in [1], we have performed calculations for two vacuum magnetic field models. However, below we present the results only for the second model, since the first model revealed approximately the same trends as the second one.
In the second model, the vacuum magnetic field was given by a three-parameter family of functions: The calculations were carried out for discrete values of the parameter R in the range from R = 1.1 to R = M. As notedthe equation (33), the closer R is to unity, the greater is the degree of anisotropy.
Results of calculations are summarized in a series of graphs in figure 2. Within each graph, it is not difficult to detect a trend towards a decrease in the critical value of beta with an increase in the steepness of the radial pressure profile as the index k increases from k = 1 to k = ∞ for a fixed pair of parameters q and R.
The second trend is that critical beta β crit rapidly approaches zero at R → 1, but in our calculations we did not take R values less than 1.1, believing that they are unlikely to be achieved experimentally. At R = 1.1, critical beta ranged from 0.531 for the smoothest radial profile k = 1 to 0.178 for a plasma with a sharp boundary at k = ∞. For sufficiently large values of R ⩾ 3-5, critical beta approaches unity, and is closer to unity for smaller values of k. It is significant that with an increase in R, the stability zone disappears for smooth radial profiles (primarily for k = 1).
The dependence on the index q, which characterizes the width and steepness of the magnetic mirrors (the larger q, the narrower the mirror width), does not seem to be very significant, except that profiles with an increased q are more prone to the disappearance of the stability zone as the parameter R increases. In other words, stabilization of the rigid ballooning mode is more problematic in traps with short and steep magnetic mirrors.
We also noticed a moderate dependence of the critical beta on the mirror ratio M. Decreasing M from 16 to 4 reduces the value of β crit to the second or third decimal place.
Although the function in equation (36) is smooth everywhere, on the profile of the plasma boundary a(z) near the median plane z = 0, a 'swell' is formed in the form of a 'thorn' with a large curvature at a single point. An example of such a 'thorn' for q = 2 is shown in figure 3(a). At q = 8, the 'thorn' expands, forming the Beklemishev's diamagnetic 'bubble', as in figure 3(c). The graphs of the plasma boundary a(z) in figure 3 are plotted for different values of k and different values Critical beta in the model magnetic field, equation (36), versus the anisotropy parameter R for different mirror ratios M and different indices q in the limit Λ → ∞. The stability zone for the radial profile with index k is located above the corresponding curve. Correspondence of the colors and markers of the curves on the graphs to the index k is shown below the bottom row of graphs.  of β crit corresponding to them 2 , but with the same value of the magnetic flux ψ = ψ a = 1 taken up by the plasma. Interestingly, such plots of a(z) almost coincide, although the values of β crit for radial pressure profiles with different k differ quite significantly. A similar effect was reported in [1] for isotropic plasma. Note in passing that this effect simplifies the manufacturing of a proportional chamber.
The formation of a region of large curvature at β → 1 leads to the violation of the paraxial approximation, which was assumed in the derivation of the LoDestro equation. The formal condition for the applicability of paraxial approximation can be written as: where the functions a 2 k for the pressure profiles with different indices k are defined by equation (A1) in appendix A. According to their definition, equation (3), functions a 2 k are proportional to the magnetic flux ψ a trapped in the plasma, but when writing equation (A1) parameter ψ a was set to unity, so it was added explicitly to the condition in equation (37). Formulating the paraxiality condition, we assume that ψ a = a 2 v0 B v0 /2, where a v0 is + radius of the plasma column in the median plane at vacuum magnetic field.
An analysis of the graphs on the left-hand side of the condition in equation (37) for the magnetic field, equation (36), shows that it is most difficult to fulfill this condition in the case of q = 2, when the maximum curvature is reached at the median plane z = 0. In the other two variants q = 4 and q = 8, the first two derivatives of the vacuum field B v (z) are equal to zero at z = 0, so the curvature peak is formed at some distance from the median plane, where the vacuum field is larger and the local beta value is smaller. For the field of the type q = 2 near the median plane, we approximately have: where l determines the scale of the change in the vacuum field in this region. Substituting equation (38) into equation (A1) and calculating the left-hand side of equation (37) at z = 0, we obtain the desired constraint on the paraxiality parameter. For the smoothest pressure profile k = 1, this method yields: For the other three cases (k ∈ {2, 4, ∞}), respectively, we have: This shows that the most severe restriction on the value of beta takes place for the steepest pressure profile (k = ∞). From equation (39d) we obtain the condition: From the formal point of view of a refined mathematician, it can be true even for β → 1 if the parameter a v0 /l tends to zero even faster, but in real open traps the parameter a v0 /l, although small, has a finite value. Accordingly, our results will be unreliable if β is too close to unity. The displacement profiles of the plasma column ξ n (z) are shown in figure 4 for all radial pressure profiles k at the same value β = 0.9. We emphasize that the displacement is not constant, although, as shown in [1], ϕ(z) = const for Λ → ∞. At the critical values of beta indicated in figure 3, the displacement profiles would be practically the same for all k, since the profiles of the plasma boundary a(z) in figure 3 practically coincide.
Next, in section 5, we describe the results of solving the Sturm-Liouville problem for the LoDestro equation (2) with a finite value of Λ. We use the Λ → ∞ limit to test convergence of this solution.

Finite vacuum gap
Taking into account the symmetry of the magnetic field with respect to the median plane z = 0, it suffices to find a solution to equation (2) at a half-interval 0 < z < 1. Due to the symmetry, the desired function ϕ(z) must be even, hence: As for the point z = 1, according to equation (11), the function must satisfy the condition which corresponds to the case when the plasma ends are insulated from the conducting elements of the structure of the vacuum chamber. An obvious fit to the LoDestro equation with the boundary conditions in equations (41) and (42) is the trivial solution ϕ ≡ 0. To eliminate the trivial solution, we impose the normalization condition: It is possible to simultaneously satisfy the three boundary conditions, equations (41)-(43), when solving a second-order ordinary differential equation only for a certain combination of the parameters Λ 0 = Λ(0), β, and ω. If a pair of parameters Λ 0 (inverse width of the vacuum gap between the plasma column and the conducting wall) and β (plasma pressure) is given, for solution of the Sturm-Liouville problem in equation (2), equations (41)- (43) give the eigenfrequency ω of the rigid ballooning mode. We will be interested in solving this problem for the case of the proportional chamber when the ratio r w /a is constant (so that the function Λ(z) is also a constant independent of z), and the oscillation frequency is equal to zero, ω = 0. This solution gives the critical value of beta, β crit , as a function of Λ 0 . The application of the shooting method for solving such a problem in the Wolfram Mathematica © system was described earlier in [1] using the example of an isotropic plasma. Let us focus on its features in application to the anisotropic model.
In what follows, we denote by z s the coordinate of the turning point on the z axis, where B = B v = B s = 1, and p ⊥ = p ∥ = 0. For the second magnetic field model, equation (36), we have: In our model of anisotropic plasma, its pressure is zero in the region z s < z < 1 between the stop point B = B v = 1 and the magnetic mirror throat B = B v = M/R. At zero pressure, the LoDestro equation (2) takes an extremely simple form: Its solution is the equality: where we allow for a moment that Λ may be a function of z, leaving the possibility for a future analysis of a conducting wall with a shape different from a proportional one. The constant on the right-hand side can be found from the boundary condition at z = 1, where the derivative ϕ ′ (1) is equal to zero. From that it follows that the constant on the right-hand side of equation (46) is also zero. Since the factor Λ(z) + 1 is greater than zero everywhere, we conclude that ϕ ′ (z) = 0 in the entire region z s < z < 1. Thus, it suffices to find the numerical solution of the original equation (2) in the region 0 < z < z s . It should be taken into account that the derivative ϕ ′ (z) undergoes a jump at z = z s . Indeed, integrating equation (2) over an infinitesimal neighborhood of the point z s from z − s to z + s , we obtain the equation: in which we took into account that Λ, ϕ, and ⟨p⟩ are continuous at the point z = z s , in contrast to the derivative ϕ ′ (z) and the coefficient: The jump in the coefficient is due to the fact that for B = B v = 1 the derivative of the functions in equations (19) and (21) jumps. Since ϕ ′ (z + s ) = 0 and Q(z + s ) = 0, from equation (47) we find the value the derivative of ϕ ′ (z) must have at the point z − s on the right boundary of the interval 0 < z < z s from its inner side: When solving equation (2) in the interval 0 < z < z − s , the boundary condition from equation (49), should be used instead of equation (42).
The implementation of the shooting method in the Wolfram Mathematica © system for solving the Sturm-Liouville problem formulated above was described earlier in [1] using an isotropic plasma as an example. Let us briefly repeat this description, emphasizing the features characteristic of an anisotropic plasma.
First, the built-in ParametricNDSolveValue utility finds a solution to the ordinary differential equation (2) in the Wolfram Mathematica © system with boundary conditions from equations (41) and (43). This utility returns a reference pf to an interpolation function of the z coordinate, which also depends on the free parameters β and Λ 0 . Other parameters (M, R, q, ω) were given as numbers.
At the second step, we have to satisfy the third boundary condition, equation (49). In the Wolfram Mathematica © , function arguments are written in square brackets, so pf [β, Λ 0 ][z] denotes the solution at z for specified numerical values of the parameters β and Λ 0 . Accordingly, the analogue of the derivative ϕ ′ (z) is written as pf [β, Λ 0 ] ′ [z]. So to solve equation (49) and, thus, find the critical beta it is enough to pass the equation: to the built-in FindRoot utility or to the RootSearch custom utility developed by Ersek [43]. Because both of these utilities miss roots from time to time, we have developed our own XRS package (see appendix D). Our calculations were done for the magnetic field, equation (36), with mirror ratios M ∈ {4, 8, 16}, those combinations of parameters k ∈ {1, 2, 4, ∞}, q ∈ {2, 4, 8}, which are listed in section 4, and for discrete values of Λ(z) = Λ 0 ∈ {1, 1.001, 1.002, 1.003, . . . , 400, 450, 500}. For Λ 0 = 500, the critical beta value we calculated differed from the value found in previous section for Λ 0 = ∞ only in the third decimal place, while in the isotropic plasma model the difference was observed only in the fifth decimal place [1]. Parameter R again varied from R = 1.1 to R = M.
Critical beta values for the case Λ = 1, when the conducting lateral wall is removed to infinity, have not been found. However, such values were found for Λ → 1+. The discussion of the stability zone at values of Λ close to unity is relegated to appendix C, since it is of academic rather than practical interest because the stable zone β crit < β < 1 is extremely narrow is this case. Figures 5-7 show plots of β crit versus ratio r w /a = (Λ 0 + 1)/(Λ 0 − 1) respectively for mirror ratios 4, 8, and16 at various combinations of k, q, and R. Comparison of graphs in consecutive rows of each figure demonstrates strong dependence of critical betas on parameter R, which characterizes the degree of anisotropy (the smaller the R, the stronger the anisotropy, the lower the β crit ) and the spatial width of the pressure distribution (the larger the R, the wider the distribution). Comparison of graphs within each row confirms a weak tendency noted in section 4 to increase critical beta as the magnetic mirrors steepen with increasing index q. Comparison of figures 5-7 for different M demonstrates some decrease in the stability zone as the mirror ratio increases.
We also see that the most smooth radial profile (k = 1) is unstable for all ratios r w /a if R ⩾ 2 and q = 4 or q = 8. For the case of isotropic plasma [1], the stability zone disappeared also for the next steepest profile (k = 2) and sufficiently steep axial profiles of the magnetic field, but in the case of anisotropic plasma this profile can always be made stable by decreasing the vacuum gap (i.e. decreasing r w /a).
For the minimum anisotropy allowed in our calculations, that is, for R = M, the model of normal injection of neutral beams (12) and the model of isotropic plasma give close results. This fact is confirmed by figure 8. The small difference becomes more noticeable for smaller M and larger q.
The main trends identified by our calculations are listed below: • If the stability zone can in principle exist for a given set of parameters M, R, and q (i.e. if a numerical solution of equation (35) was found in section 4 for this set at Λ 0 = ∞), then it emerges if Λ 0 exceeds some minimum value Λ min > 1. • This minimum value is smaller for steeper radial pressure profile (larger k), smoother axial profile of the vacuum magnetic field (smaller q), smaller mirror ratio M, and greater anisotropy (smaller R). The stability zone narrows with increasing Λ min and may disappear altogether for smooth pressure profiles (k = 1). • Critical betas for the case Λ 0 = 1, when the conducting lateral wall is removed (r w /a = ∞), have not been found. • The narrower the vacuum gap between the plasma and the lateral conducting wall (the larger Λ 0 ), the wider the stability zone (the smaller β crit ). • The minimum β crit found for the studied set of radial and axial profiles and parameter of anisotropy R = 1.1 is about 18%, which is much lower than the value 70% reported for isotropic plasmas. We expect β crit → 0 at R → 1, but such a limit does not seem achievable in a real experiment.

Combined wall stabilization
Finally, we repeat calculations of the previous section with a minor adjustment. We take the boundary condition from equation (10), which means that the plasma is frozen into the conducting end plates, instead of the boundary condition, equation (11), describing insulating ends of a mirror trap. In our recent paper [1], we reported two critical beta values for the case of isotropic plasmas, β crit1 and β crit2 : one at low beta due to the balancing of the ponderomotive force with the curvature drive, and one at high β due to the proximity of the conducting lateral wall which enables magnetic line bending to balance the curvature drive. Corresponding to two critical values of beta, there are two zones of stability. The first zone exists at low plasma pressure, at 0 < β < β crit1 , and the second one exists at high pressure, at β crit2 < β <1. These two zones merge at larger Λ 0 providing overall stability at any beta. In this section we extend these results to the case of anisotropic plasma produced by normal NB injection. The main difference from the isotropic model is that the degree of anisotropy provides an additional 'knob' for controlling experimental parameters.
Let's move on to solving equation (2) with the boundary conditions from equation (10) for z = 1, and equations (41) and (43) for z = 0. As shown in section 5, in the region z s < z < 1 the desired solution satisfies equation (46). Noticing that the constant in this equation is now equal to [Λ(z s ) + 1] ϕ ′ (z + s ), rather than zero as in section 5, we find that the derivative ϕ ′ (z) at z = z + s is equal to: (51)    Substituting ϕ(z + s ) into equation (47) and taking into account that Q(z + s ) = 0, we find the derivative ϕ ′ (z − s ) on the right boundary of the interval 0 < z < z s from its inner side: This is the boundary condition that should be used instead of equation (49) in the problem of ballooning instability with a combination of wall stabilization and stabilization by conducting end plates. In the case of Λ(z) = const under consideration here it reduces to: Following the same scheme as in [1], we performed a series of calculations for the vacuum magnetic field, equation ( A series of graphs in figure 9 illustrates the results of calculations at the minimum degree of anisotropy, i.e. at R = M, when the instability zone has the maximum dimensions. The zone of instability in this and subsequent figures for each radial pressure profile with a given index k ∈ {1, 2, 4, ∞} lies between the lower and upper curves of the corresponding color. They represent β crit1 and β crit2 , respectively. If there is only one curve, as in figure 9(a) for the profile k = 1, we interpret it as β crit1 . If there are no curves of a given color on a graph, we suppose that the corresponding pressure profile is stable in the entire range of 0 < β < 1 and in the entire range of values of r w /a, which is shown in the graph. We have made the last two statements, analyzing the evolution of the instability zone as the parameters of the problem, such as r w /a, R, and M, change continuously. Strictly speaking, we should have verified these conclusions by calculating the sign of ω 2 . However, our numerical code does not currently provide this option.
In the range r w /a to the left of the left edge of the displayed interval (i.e. in the region where the ratio r w /a is relatively small), all radial profiles are stable over the entire interval 0 < β < 1. It is easy to see that smooth pressure profiles (k = 1, k = 2) are more stable than steep profiles (k = 4, k = ∞). The same trend also takes place for small-scale ballooning perturbations when using end MHD stabilizers [5].
In the opposite case r w /a → ∞, the upper stability zone becomes extremely thin (or even disappears for smooth radial profiles) but it is located in the immediate vicinity of the β = 1 boundary, where the paraxial approximation does not work. At the same time, the lower zone remains quite wide. This fact is in good agreement with the results of calculating the smallscale ballooning instability threshold [3][4][5], as we emphasized in the introductory section 1.
Surprise is caused by another fact, namely: the disappearance of the instability zone at r w /a → ∞ as the radial pressure profile is smoothed out. For the sharpest profile k = ∞, the instability zone is present on all graphs in figure 9; for the k = 4 profile it is visible in 7 out of 11 graphs, for the profile k = 2 only in 2, and for k = 1 only in 1. We interpret the absence of a gap between the lower and upper stability zones for a given ratio r w /a as the stability of rigid ballooning perturbations over the entire interval 0 < β < 1. But then the conclusion is inevitable that a sufficiently smooth radial pressure profile can be stable in the absence of a lateral wall for any value of β.
For comparison, figure 10 shows the results of calculations for isotropic plasma with the same combinations of indices q and mirror ratios M as in figure 9. As can be seen from these two figures 9 and 10, the dimensions of the instability zones for the same mirror ratio are relatively close, but in the case of an anisotropic plasma they are still smaller, which becomes more noticeable as the mirror ratio decreases. Both in the case of anisotropic and isotropic plasma models, the dimensions of the instability zone β crit1 < β < β crit2 (which is not shaded) are maximum for the smoothest magnetic field profile with index q = 2. They are minimal at q = 8.
In contrast, as can be seen from the analysis of figures 5-7 in section 5, with only wall stabilization implemented without end plates installed, the instability zone at β < β crit gets slightly larger as q increases. In the case of combined stabilization, for a fixed set of parameters q, M, and R the instability zone is the widest for the steepest radial pressure profile (k = ∞) and may be completely absent for smooth profiles (k = 1, k = 2). Without stabilization by the conducting end plates, the opposite situation takes place: the stability zone is smaller and may be completely absent for smooth profiles. It means that the effect of end plates dominates. Another important effect of end MHD stabilizers is a much stronger dependence on the mirror ratio M and the axial profile of the magnetic field characterized by the index q.
A series of graphs in figure 11 illustrates the dependence of the stability/instability zones on the anisotropy parameter R at a fixed index q = 2 and a fixed mirror ratio M = 16. The unshaded zone of instability noticeably decreases with increasing anisotropy (with decreasing parameter R) and may even disappear, as in figures 11(a) and (b) for all smooth radial profiles (k ∈ {1, 2, 4}).
A series of graphs in figure 12 shows the dependence of β crit1 on the anisotropy parameter R for different values of the mirror ratio M and index q characterizing the width of the magnetic mirrors in the case when the lateral conducting wall is removed (formally has an infinite radius r w ). All these graphs confirm the tendency noted above to expand the lower stability zone with smoothing of the radial pressure profile up to the occupation of the entire interval 0 < β < 1, which is admissible under the condition of transverse equilibrium of the plasma column in the mirror trap. Also noteworthy is the fact that the stability zone noticeably increases as the mirror ratio decreases. Note that the absence of graphs for the case M = 2, q ∈ {4, 8} in figure 12 means that there is no instability zone even for the sharpest radial profile k = ∞.
Summing up everything said in this section, we can state the following: • Two stability zones are found for moderate values of parameter Λ and a sufficiently large mirror ratio M. The lower zone β < β crit1 exists even for Λ = 1. • With other things being equal, the instability zone is the widest for the steepest pressure profile (k = ∞) and might not exist at all for smooth radial profiles (k = 1 or k = 2). Recall that in section 5 it is for the stepwise profile that the instability zone has the minimum dimensions. • For a fixed value of the parameter Λ, the stability zones expand and can merge with a decrease in the mirror ratio M and/or a smoothing of the radial pressure profile (with a decrease in k). • The zone of instability decreases and may even disappear as the plasma anisotropy increases (as the parameter R decreases). The dimensions of the unstable zone are maximum at the minimum studied degree of anisotropy R = M, considered for a given mirror ratio M. In this limit, the dimensions of the unstable zone are close to those found for an isotropic plasma. • If an instability zone exists between two stability zones for some combinations of the parameters k, q, M, and R, then it disappears if Λ 0 > Λ crit > 1. The value of Λ crit is smaller for lower values of k and M, and is also smaller for larger q.
The largest value Λ crit = 1.52 (r w /a = 2.201 40) in our calculations was found at q = 2, M = R =24, k = ∞. • With a not too large mirror ratio, and/or a sufficiently high plasma anisotropy, and/or a sufficiently steep magnetic field, and/or a sufficiently smooth pressure profile, and/or sufficiently close lateral wall, the rigid ballooning mode m = 1 can be stabilized at any feasible value of beta. Figure 12. Boundary of lower stability zone β crit1 for model magnetic field, equation (36), and anisotropic plasma model, equation (12), simulating normal injection of neutral beams, with stabilization only by end conducting plates without lateral conductive wall. The stability zone for the radial profile with index k is located below the corresponding curve. The absence of a curve means that the corresponding profile is stable over the entire interval 0 < β < 1.

Conclusions
In the present work, we show the feasibility of stabilization of the m = 1 rigid flute and ballooning modes in an axially symmetric mirror trap by a conducting axially symmetric shell that follows the shape of the plasma border. The stabilizing effect is calculated using an anisotropic plasma model, equation (12), which simulates normal (steep-angle) injection of fast neutral beams. Unlike most predecessors, who studied not quite realistic stepwise radial plasma pressure profiles, we considered a series of diffuse pressure profiles, equation (28), with different degrees of steepness, as well as several axial profiles of vacuum magnetic field given by the function in equation (36) with varying mirror ratio and axial gradient near the mirrors. Stabilization by a conducting wall without any additional means of MHD stabilization is achieved at a sufficiently high plasma pressure, when the parameter β (the dimensionless ratio of the plasma pressure to the magnetic field pressure) exceeds a certain critical value β crit . Therefore, our goal was to calculate this critical value and study its dependence on the degree of plasma anisotropy, the shape of the radial pressure profile, the axial profile of the magnetic field, the mirror ratio, and the size of the vacuum gap between the plasma and the conducting wall. For calculations, we developed a numerical code written in Wolfram language, which solved equation (2), previously derived by LoDestro.
Our calculations showed that the stability zone expands significantly due to a decrease in the critical beta β crit as the degree of plasma anisotropy increases. The mirror ratio and axial profile of the magnetic field have relatively smaller effect on the value of β crit . The dependence of β crit on the radial profile and the gap width between the plasma column and the lateral conducting wall are more significant.
From a practical point of view, a noticeable decrease in the value of β crit is achieved if the radius of the conducting wall r w exceeds the plasma radius a by no more than 2 times, r w /a ⩽ 2. The influence of the conducting wall practically disappears (in the sense that β crit → 1) if r w /a ⩾ 4. For effective wall stabilization, the ratio r w /a must be less than the indicated limits. On the other hand, too small a vacuum gap between the plasma and the wall can lead to direct losses of fast ions as well as gas desorption from the walls, resulting in degradation of the plasma parameters. For this reason, plasma stabilization using only the lateral wall might seem difficult to implement, but the parameter range near r w /a = 2 seems quite comfortable from the point of view of many experimenters.
Even more promising is the stabilization of the rigid ballooning mode with a combination of a conducting lateral wall and conducting end plates, which imitate the attachment of the end MHD anchors to the central cell of a mirror trap. In contrast to the pure wall stabilization, mirror ratio and magnetic field profile have strong effect on the combined stabilization.
Our calculations have shown the great efficiency of this method of stabilization. We found the existence of two stability zones; one at low beta due to the balancing of the lateral wall effect with the curvature drive, and one at high β due to the proximity of the conducting wall which enables magnetic line bending to balance the curvature drive. Corresponding to two critical values of beta, there are two zones of stability. The first zone exists at low plasma pressure, at 0 < β < β crit1 , and the second one exists at high pressure, at β crit2 < β <1. These two zones merge, making the entire range of allowable beta values 0 < β < 1 stable, as the mirror ratio decreases, as the vacuum gap between the plasma and the lateral wall decreases, or as the plasma pressure anisotropy increases.
Here F(ϕ |m) and E(ϕ |m) are the elliptic integral of the first and the second kinds, respectively., and F 1 (a; b 1 , b 2 ; c; x, y) is the Appell hypergeometric function of two variables.

Appendix C. Minimal Λ
The existence of a stability zone, even in the β → 1 limit, with a very large gap between the conducting wall and the plasma, was previously discovered for the case of an isotropic plasma [1], and then it was surprising. For Λ 0 = 1.01, the radius of the conducting wall r w is approximately 14 times larger than the plasma radius a. In the general case, the minimum Λ min for which the critical beta values were found was the smaller, the steeper the radial pressure profile (the larger the parameter k). It also decreased as the axial profile of the vacuum magnetic field was smoothed (with a decrease in the index q) and the mirror ratio M decreased. In the previous article [1], we did not set out to accurately calculate Λ min , but simply chose the minimum values of Λ 0 from the available list of discrete values for which we calculated β crit .
In the present work, we tried to calculate Λ min more accurately. The difficulty of such calculations is that the coefficients of the LoDestro equation are singular at the point z = 0 in the limit β → 1. On the other hand, we showed above that the LoDestro equation is inapplicable in this limit, since the paraxial approximation is violated at β → 1. Nevertheless, we searched for the roots of the dispersion equation (49) with respect to Λ 0 for four sufficiently large beta values, β ∈ {0.999, 0.9999, 0.999 99, 0.999 999}. The values of Λ min calculated in this way are given in table 1 for β = 0.999.
A cursory analysis of these tables shows that Λ min noticeably increases as the parameters q and R increase, that is, both as the magnetic mirrors steepen and as the anisotropy decreases. In contrast, Λ min decreases with an increase in the parameter k, that is, with a steepening of the radial pressure profile, and this decrease is especially noticeable at lower anisotropy. The same tendency towards a decrease in Λ min is observed as the mirror ratio K decreases.

Appendix D. Search for the roots of a nonlinear equation
The search for the roots of equations (35), (49) and (53) was initially carried out using the FindRoot utility built into the Wolfram Mathematica © library. FindRoot searches for a root near an initial guess β start passed to it. Success or failure in finding the root with this utility depends very much on luck in choosing β start . Therefore, as an additional means of searching for roots, the RootSearch package was included, which was developed by Ersek [43]. This package contains a utility of the same name that searches for all roots within a given interval.
In some intermediate version of our numerical code, a root found by the RootSearch utility was passed as β start to the FindRoot utility to recheck the result of calculation of β crit . In rare cases, when only one of the two utilities found a solution to equations (35), (49) or (53), the code was analyzed in order to improve it. In particular, the formulas that Wolfram Mathematica © obtained when calculating the integrals a 2 and ⟨p⟩ were improved. For example, in the formula (A1c) for a 2 4 , the number of Appell hypergeometric functions was reduced from four to one, and in the formula (B1c) for ⟨p⟩ 4 , from six to two. A side result of such code optimization was the 'expulsion' of complex numbers in the intermediate calculations of the integral in equation (35). Due to rounding errors for such numbers, it sometimes happened that the integral in equation (35) took on a complex value with a small but finite imaginary part. In such cases, FindRoot skipped the root of equation (35). In cases where both FindRoot and RootSearch did not find a solution to equation (35), it was considered that the solution did not exist.
Unfortunately, after all this effort, the FindRoot and RootSearch bundle of utilities kept missing roots of equations (35), (49), and (53). The problem of losing roots was solved by applying the ideas formulated in the post of the user with the nickname matheorem in the mathematica forum on the stackexchange.com portal. Using his code sample [44] and some of the code from Ersek's RootSearch [43] package, we wrote a XRS package that replaced the FindRoot and RootSearch utilities.