Plasma equilibrium in diamagnetic trap with neutral beam injection

Abstract

This paper presents a theoretical model of plasma equilibrium in the diamagnetic confinement mode in an axisymmetric mirror device with neutral beam injection. The hot ionic component is described within the framework of the kinetic theory, since the Larmor radius of the injected ions appears to be comparable to or even larger than the characteristic scale of the magnetic field inhomogeneity. The electron drag of the hot ions is taken into account, while the angular scattering of the hot ions due to Coulomb collisions is neglected. The background warm plasma, on the contrary, is considered to be in local thermal equilibrium, i.e. has a Maxwellian distribution function and is described in terms of magnetohydrodynamics. The density of the hot ions is assumed to be negligible compared with that of the warm plasma. Both the conventional gas-dynamic loss and the non-adiabatic loss specific to the diamagnetic confinement mode are taken into account. In this work, we do not consider the effects of the warm plasma rotation as well as the inhomogeneity of the electrostatic potential. A self-consistent theoretical model of the plasma equilibrium is constructed. In the case of the cylindrical bubble, this model is reduced to a simpler one. The numerical solutions in the limit of a thin transition layer of the diamagnetic bubble are found. Examples of the equilibria corresponding to the gas-dynamic multiple-mirror trap device are considered.

1 Introduction

Diamagnetic confinement, or the diamagnetic bubble (Beklemishev 2016; Khristo & Beklemishev 2019, 2022; Chernoshtanov 2020, 2022; Kotelnikov 2020; Soldatkina et al. 2023), is a new operating mode designed to significantly enhance confinement in linear systems (Dimov 2005; Steinhauer 2011b; Ivanov & Prikhodko 2017). The main idea of this regime is to form a high-pressure plasma bubble in the central part of a mirror device. Inside such a bubble, the plasma pressure $p$ reaches the equilibrium limit corresponding to $\beta = 8\pi p/B^2_{\mathrm{v}}\rightarrow 1$, where $B_{\mathrm{v}}$ is the vacuum magnetic field without plasma, and the magnetic field tends to zero due to the plasma diamagnetism. The lifetime of the particles in the diamagnetic trap, i.e. the mirror device operating in the diamagnetic confinement mode, is estimated by Beklemishev (2016) in the framework of magnetohydrodynamics (MHD) as follows:

(1.1)\begin{equation} \tau_{\mathrm{DC}}\sim\sqrt{\tau_{\mathrm{GD}}\tau_{{\perp}}}, \end{equation}

where $\tau _{\mathrm {GD}}$ is the lifetime in the gas-dynamic trap (GDT) (Ivanov & Prikhodko 2017), and $\tau _{\perp }$ is the transverse transport time. Since plasma transport across the magnetic field in axisymmetric traps is usually far below axial losses: $\tau _{\perp }\gg \tau _{\mathrm {GD}}$, a significantly enhanced overall confinement is expected in the diamagnetic trap: $\tau _{\mathrm {DC}}\sim \tau _{\mathrm {GD}}\sqrt {\tau _{\perp }/\tau _{\mathrm {GD}}}\gg \tau _{\mathrm {GD}}$. In this regard, there is a practical interest in the experimental and theoretical study of this mode. In particular, one of the goals of the new generation linear machine gas-dynamic multiple-mirror trap (GDMT) (Beklemishev et al. 2013; Bagryansky, Beklemishev & Postupaev 2019; Skovorodin et al. 2023) is to experimentally verify the concept of diamagnetic confinement. In addition, the study of the diamagnetic regime is planned at the compact axisymmetric toroid (CAT) (Bagryansky et al. 2016) currently operating at the Budker Institute of Nuclear Physics (Budker INP).

One can find a number of earlier theoretical works related to $\beta \sim 1$ plasma confinement in open traps (see Grad 1967; Newcomb 1981; Lansky 1993; Lotov 1996; Kotelnikov, Bagryansky & Prikhodko 2010; Kotelnikov 2011). Effective confinement of plasma with $\beta \sim 1$ in a gas-dynamic system was demonstrated experimentally in the 2MK-200 test facility by Zhitlukhin et al. (1984). High $\beta$ plasma confinement was also studied in the so-called magnetoelectrostatic traps (see Pastukhov 1978, 1980; Ioffe et al. 1981; Pastukhov 2021). Structures similar to the diamagnetic bubble, called magnetic holes, are often observed in space plasmas (see Kaufmann, Horng & Wolfe 1970; Turner et al. 1977; Tsurutani et al. 2011; Kuznetsov et al. 2015).

The idea of diamagnetic confinement was originally proposed by Beklemishev (2016). The possibility of a discharge transition into the diamagnetic confinement mode is briefly discussed, and a stationary MHD model of a diamagnetic bubble equilibrium is constructed in the cylindrical approximation. Further, this hydrodynamic equilibrium model was extended by Khristo & Beklemishev (2019, 2022) to the case of a non-paraxial axisymmetric trap. In particular, diamagnetic bubble equilibria in the GDMT are computed, and the effect of magnetic field corrugation on equilibrium in a diamagnetic trap is studied. Beklemishev (2016) and Khristo & Beklemishev (2022) also briefly discuss the possibility of MHD stabilization by a combination of the vortex confinement (see Soldatkina, Bagryansky & Solomakhin 2008; Beklemishev et al. 2010; Bagryansky et al. 2011) and the conducting wall (see Kaiser & Pearlstein 1985; Berk, Wong & Tsang 1987; Kotelnikov et al. 2022).

As already noted, the magnetic field inside the diamagnetic bubble is close to zero. In this case, the Larmor radius and the mean free path of high-energy particles can be comparable to or even exceed the characteristic scale of magnetic field inhomogeneity, which is beyond the scope of MHD theory. Therefore, there is a need for a detailed kinetic model of fast particles in the diamagnetic trap. Research on the kinetic theory of the diamagnetic regime is already underway at the Budker INP. Kotelnikov (2020) constructed a fully kinetic equilibrium model in a cylindrical geometry with a distribution function isotropic in the transverse plane inside the bubble. The collisionless dynamics of individual particles in a diamagnetic trap was studied by Chernoshtanov (2022). In addition, a particle-in-cell model for simulating high-pressure plasma in an open trap is currently being developed at the Budker INP by Kurshakov & Timofeev (2023). Work is also underway to create a code for numerical simulation of the diamagnetic regime (see Boronina et al. 2020; Efimova, Dudnikova & Vshivkov 2020; Chernoshtanov et al. 2023, 2024).

Modern experiments on linear devices, such as the already mentioned GDMT and CAT, commonly involve the injection of high-energy (of the order of several tens of electronvolts) neutral beams to heat the plasma (Ivanov & Prikhodko 2017; Belchenko et al. 2018). With this in mind, in the present paper, we aim to construct a theoretical model of diamagnetic bubble equilibrium in a GDT with neutral beam injection. The neutral beams are absorbed by the background plasma at a temperature much lower than the energy of the injected atoms. For this reason, to describe the equilibrium in such a system, we assume the plasma to consist of two fractions: the background warm plasma and the hot ions resulting from the neutral beam injection. The former is considered to be in local thermodynamic equilibrium and is described in terms of MHD; we take as a basis the hydrodynamic model constructed in the earlier works (Beklemishev 2016; Khristo & Beklemishev 2019, 2022). The hot ions, on the contrary, are assumed to have a non-Maxwellian distribution function and are described within the framework of the kinetic theory. Similar hybrid equilibria were considered earlier in application to field reversed configurations (FRCs) (Rostoker & Qerushi 2002; Qerushi & Rostoker 2002a, b, 2003; Steinhauer 2011a).

For simplicity, we consider the approximation of a long, axisymmetric diamagnetic trap. In this case, in the absence of dissipation and scattering, the azimuthal angular momentum and the total energy of a single charged particle are conserved. At the same time, since the magnetic field in the diamagnetic confinement regime has considerable gradients, the magnetic moment is not globally conserved. However, it was shown independently by Chernoshtanov (2020, 2022) and Kotelnikov (2020) that, in a diamagnetic trap, the adiabatic invariant $I_{r}=(2{\rm \pi} )^{-1}\oint p_{r}\,{\rm d}r$ can be conserved under specific conditions, where $p_{r}$ is the radial component of the particle momentum. In addition, it is known that in such axisymmetric systems there is a region in phase space where particles with a sufficiently large canonical angular momentum are absolutely confined (see Morozov & Solov'ev 1966; Lovelace, Larrabee & Fleischmann 1978; Larrabee, Lovelace & Fleischmann 1979; Hsiao & Miley 1985).

In the present paper, we assume the hot ions to be injected into the region of the phase space, where, first, the absolute confinement criterion is met and, second, the condition of the adiabaticity is violated. The former saves us from solving the complex problem of taking into account the hot ion loss. The latter leads to the distribution function of the hot ions being homogeneous on the hypersurface of the constant azimuthal angular momentum and total energy, and hence not depending on the adiabatic invariant $I_{r}$. Violation of the adiabaticity occurs typically due to the non-paraxial ends of the bubble configuration (Beklemishev 2016). We also assume the energy of the hot ions to be much higher than the temperature of the warm plasma. In this approximation, the hot ions are mainly slowing down on the electrons of the warm plasma and hardly collide with the ions. On this basis, when calculating the hot ion distribution function, we completely neglect the angular scattering due to Coulomb collisions and take into account only the weak drag force from the warm electrons.

The article is structured as follows: § 2 focuses on the detailed problem statement; § 3 defines the equations describing the equilibrium of the magnetic field; § 4 derives the equilibrium equations for the warm plasma. Section 5 finds the equilibrium distribution function of the hot ions; § 6 reduces the theoretical model to the case of the cylindrical diamagnetic bubble; § 7 is devoted to the solution of the equilibrium equations in the approximation of a thin transition layer at the bubble boundary; § 8 considers examples of equilibria corresponding to the diamagnetic confinement regime in the GDMT device; § 9 summarizes the main results of this paper and discusses the issues to be addressed in future work.

2 Basic assumptions of the theoretical model

Consider the stationary equilibrium of a diamagnetic bubble in an axisymmetric GDT. At the periphery of the bubble, the magnetic field is close to the vacuum magnetic field $B_{{\rm v}}$, i.e. the magnetic field without plasma. In the interior of the bubble, the magnetic field is vanishingly small $B\ll B_{{\rm v}}$, being almost completely expelled by diamagnetic plasma; this region we further refer to as the core of the diamagnetic bubble. We denote the radius of the bubble core as $r_{0}$, which, generally speaking, can vary along the trap: $r_{0}=r_{0}(z)$.¹ The core radius in the central section of the trap, $z=0$, we define as $a=r_{0}(0)$. The region at the boundary of the bubble, inside which the magnetic field changes from $B\simeq 0$ in the core to $B=B_{{\rm v}}$ at the periphery, we further refer to as the transition layer of the diamagnetic bubble.

Let the plasma consist of hot ions, resulting from the neutral beam injection, and background warm plasma. The warm plasma is assumed to be in local thermal equilibrium with the temperature $T=T(r,z)$ and described in terms of MHD. The hot ions, on the contrary, are expected to have a non-Maxwellian distribution function and are described within the framework of kinetic theory.

As was previously found by Khristo & Beklemishev (2022), the specific choice of the warm plasma transport model in the bubble core seems to have little effect on the equilibrium. For this reason, we further assume that, inside the bubble core, the magnetic field is identically zero $B\equiv 0$, and the warm plasma electrical conductivity $\sigma _{{\rm w}}$ and transverse diffusion coefficient $\mathcal {D}$ are extremely high: $\sigma _{{\rm w}}\rightarrow \infty$ and $\mathcal {D}\rightarrow \infty$. In what follows, we also consider the approximation of a long paraxial bubble with short non-paraxial ends. In addition, we do not take into account the rotation of the warm plasma and the effect of the electrostatic potential inhomogeneity. The warm plasma radial electric currents are also neglected. We understand that these issues are definitely important and thus should be addressed in future work.

In an axisymmetric system, the azimuthal canonical angular momentum $\mathcal {P}$ and the total energy $\mathcal {E}$ of a charged particle are conserved

(2.1a,b)\begin{equation} \mathcal{P}=m_{{\rm s}}rv_{\theta}+\frac{e_{{\rm s}}}{c}\frac{\psi}{2{\rm \pi}}=\mathrm{const}.,\quad \mathcal{E}=\frac{m_{{\rm s}}v^{2}}{2}=\mathrm{const}.,\end{equation}

where $m_{{\rm s}}$, $e_{{\rm s}}$ are the mass and the electric charge of a particle of species ${{\rm s}}$, respectively, $v_{\theta }$ is the azimuthal component of the particle velocity and $\psi =\int _{0}^{r}B_{z}(r',z)2{\rm \pi} r'\,{\rm d}r'$ is the flux of the longitudinal (‘poloidal’) component of the magnetic field.² Assuming the magnetic flux in the mirrors to be approximately equal to $\psi _{{\rm m}}\simeq B_{{\rm m}}{\rm \pi} r^{2}$, we arrive at the absolute confinement criterion in the form (see Morozov & Solov'ev 1966; Lovelace et al. 1978; Larrabee et al. 1979; Hsiao & Miley 1985):

(2.2)\begin{equation} {-}\mathcal{R}\varOmega_{{\rm s}}\mathcal{P}>\mathcal{E},\end{equation}

where $\varOmega _{{\rm s}}=e_{{\rm s}}B_{{\rm v}}/m_{{\rm s}}c$ is the cyclotron frequency in the vacuum magnetic field of the central section of the trap $B_{{\rm v}}$, $\mathcal {R}=B_{{\rm m}}/B_{{\rm v}}$ is the vacuum mirror ratio, $B_{{\rm m}}$ is the mirror magnetic field, and $c$ is the speed of light. Worth noting is that the particles are absolutely confined only if the direction of their rotation and the direction of the Larmor rotation are the same. In other words, positively charged particles with $\varOmega _{{\rm s}}>0$ are confined if $\mathcal {P}<0$ and negatively charged particles with $\varOmega _{{\rm s}}<0$ are confined if $\mathcal {P}>0$. In order to avoid solving the complex problem of taking into account the hot ion loss, in this paper, we consider the neutral beam being injected into the absolute confinement region (2.2). It is clear that, in practice, the injection should be carried out in this way to avoid undesirable loss of the hot ions.

It was discovered by Kotelnikov (2020) and Chernoshtanov (2020, 2022) that there exists the adiabatic invariant $I_{r}=(2{\rm \pi} )^{-1}\oint p_{r}\,{\rm d}r$ for particles in a diamagnetic trap. As additionally shown by Chernoshtanov (2020, 2022), $I_{r}$ is conserved for particles with not too high velocity $v_{\parallel }$ along the magnetic field, namely, the adiabaticity criterion can be approximately written as a limitation on the pitch angle $\xi$

(2.3)\begin{equation} \tan\xi=\frac{v_{{\perp}}}{\left|v_{{\parallel}}\right|}\gtrsim\left|\frac{{\rm d}r_{0}}{{\rm d}z}\right|_{\max},\end{equation}

where $v_{\perp }$ is the velocity transverse to the magnetic field, $|{\rm d}r_{0}/{\rm d}z|_{\max }$ is the maximum inclination of the field line corresponding to the bubble core boundary $r_{0}=r_{0}(z)$. The criterion (2.3) imposes considerable restrictions on the uniformity of the magnetic field. In addition to the requirement of sufficient global smoothness of the longitudinal profile of the field lines, the small-scale ripples of the magnetic field, resulting from the discreteness of the magnetic system and instabilities, should also be insignificant. Due to this, the region in the phase space where the criterion (2.3) is met usually proves to be relatively narrow. To simplify the theoretical model, we further assume that the condition (2.3) is violated and the adiabatic invariant $I_{r}$ is not conserved. As a result, the hot ion dynamics becomes chaotic, and the trajectories fill the hypersurfaces with constant angular momentum and energy $\mathcal {P}=\mathrm {const}.$, $\mathcal {E}=\mathrm {const}$. Based on this, we approximately consider the distribution function of the hot ions depending only on $\mathcal {P}$ and $\mathcal {E}$. The collisionless dynamics of the hot ions in a diamagnetic trap may resemble that in FRCs. In particular, Landsman, Cohen & Glasser (2004) provide a classification of particle trajectories in FRCs, where the regions of phase space corresponding to chaotic motion are observed as well.

As shown by Chernoshtanov (2020, 2022), in a diamagnetic trap there should also arise an additional axial loss, the so-called non-adiabatic loss. It results from the particles escaping the trap through a ‘leak’ in the phase space outside the region of the absolute confinement (2.2). The characteristic non-adiabatic loss time for plasma with a Maxwellian distribution function at the centre of the trap can be estimated as

(2.4)\begin{equation} \tau_{{\rm s}\parallel0}\sim\tau_{\mathrm{GD}\,{\rm s}}\frac{a}{\rho_{{\rm s}}}, \end{equation}

where $\tau _{\mathrm {GD}\,{\rm s}}=\mathcal {R}L/v_{T{{\rm s}}}$ is the gas-dynamic confinement time, $a$ and $L$ are the characteristic radial size and length of the diamagnetic bubble, $\rho _{{\rm s}}=v_{T{{\rm s}}}/|\varOmega _{{\rm s}}|$ is the characteristic Larmor radius and $v_{T{{\rm s}}}=\sqrt {T/m_{{\rm s}}}$ is the thermal velocity. Note that the time $\tau _{{\rm s}\parallel 0}$ does not depend on the particle mass since it is cancelled in $|\varOmega _{{\rm s}}|/v_{T{{\rm s}}}^{2}$. As can be seen, despite the additional non-adiabatic loss, confinement of the warm plasma in the diamagnetic trap can still be significantly enhanced compared with the ‘classical’ GDT, provided the ratio $\rho _{{\rm s}}/a$ is sufficiently small.

In this work, we assume the injection energy $\mathcal {E}_{\mathrm {NB}}$ to be large compared with the temperature of the warm plasma $T$. However, due to the significant difference in masses, the characteristic velocity of the injected ions $\sqrt {\mathcal {E}_{\mathrm {NB}}/m_{{\rm i}}}$ is still much less than the thermal velocity of the warm electrons $\sqrt {T/m_{{\rm e}}}$. Summarizing the above, we have

(2.5)\begin{equation} 1\ll\frac{\mathcal{E}_{\mathrm{NB}}}{T}\ll\frac{m_{{\rm i}}}{m_{{\rm e}}}, \end{equation}

where $m_{{\rm e}}$ and $m_{{\rm i}}$ are the electron and ion masses, respectively. We also consider the density³ of the hot ions to be negligible compared with that of the warm plasma, so the quasi-neutrality condition is satisfied for the latter: $n_{{\rm i}}\simeq n_{{\rm e}}/Z$, where $n_{{\rm e}}$ and $n_{{\rm i}}$ are the densities of warm electrons and ions, respectively, and $Z$ is the ion atomic number. In addition, we neglect collisions of the hot ions with ions and take into account only the drag force from the warm electrons. In other words, we use a simple linear kinetic equation for the hot ions, in which we keep only the terms related to the slowing down of the hot ions on the warm electrons, neglecting the angular scattering due to Coulomb collisions. This approximation apparently corresponds to the high-energy limit $\mathcal {E}\gg T$ of the injected ions.

Let us briefly summarize the formulation of the problem. We consider the stationary axisymmetric equilibrium of a diamagnetic trap in the paraxial approximation. The plasma consists of the warm plasma in thermal equilibrium and non-Maxwellian hot ions. The warm plasma is described in terms of MHD. Both gas-dynamic and non-adiabatic losses of the warm plasma are taken into account. The electrical conductivity and transverse diffusion coefficient of the warm plasma inside the bubble core are assumed to be extremely high. We neglect the rotation of the warm plasma and the inhomogeneity of the electrostatic potential. Hot ions are considered within the framework of kinetic theory. A linear kinetic equation is used, which takes into account only the electron drag of the hot ions. The distribution function of the hot ions is considered to depend only on azimuthal angular momentum and energy.

3 Magnetic field equilibrium distribution

The equilibrium stationary distribution of the magnetic field can be found from Maxwell's equation with a source

(3.1)\begin{equation} \boldsymbol{\nabla}\times\boldsymbol{B}=4{\rm \pi} c^{-1}\boldsymbol{J},\end{equation}

where $\boldsymbol {J}$ is the total electric current density. Since the radial and longitudinal currents, as well as the azimuthal magnetic field, are assumed to be zero, we have $\boldsymbol {J}=J_{\theta }\hat {\boldsymbol {\theta }}$ and $\boldsymbol {A}=A_{\theta }\hat {\boldsymbol {\theta }}$. In the axisymmetric case before us, we can also fix the vector potential gauge as follows:

(3.2)\begin{equation} \partial_{\theta}A_{\theta}=0\quad\Leftrightarrow\quad A_{\theta}=A_{\theta}\left(r,z\right). \end{equation}

Then the vector potential $A_{\theta }$ proves to be related to the axial magnetic flux $\psi$

(3.3a,b)\begin{equation} A_{\theta}=\frac{\psi}{2{\rm \pi} r},\quad\psi=\int _{0}^{r}B_{z}\left(r',z\right)2{\rm \pi} r'\,{\rm d}r'. \end{equation}

Therefore, (3.1) reduces to the following two-dimensional second-order differential equation for $\psi$:

(3.4)\begin{equation} r\partial_{r}\left(r^{-1}\partial_{r}\psi\right)+\partial_{z}^{2}\psi=-8{\rm \pi}^{2}c^{-1}rJ_{\theta},\end{equation}

which is the Grad–Shafranov equilibrium equation (Grad & Rubin 1958; Shafranov 1958) in the case of an axisymmetric mirror device. The electric current density $J_{\theta }$ on the right-hand side of (3.4) contains both the plasma diamagnetic current $J_{{\rm p}\theta }$ and the current in the external coils $J_{{\rm v}\theta }$, i.e. $J_{\theta }=J_{{\rm p}\theta }+J_{{\rm v}\theta }$. Having solved (3.4), one can evaluate the magnetic field $\boldsymbol {B}$ via the found magnetic flux $\psi$ as follows:

(3.5)\begin{equation} \boldsymbol{B}=\frac{1}{2{\rm \pi} r}\boldsymbol{\nabla}\psi\times\hat{\boldsymbol{\theta}}, \end{equation}

and then

(3.6)\begin{equation} B=\left|\boldsymbol{B}\right|=\frac{1}{2{\rm \pi} r}\sqrt{\left(\partial_{r}\psi\right)^{2}+\left(\partial_{z}\psi\right)^{2}}. \end{equation}

We consider the external boundary condition of (3.4) to be the regularity of the magnetic field at infinity

(3.7)\begin{equation} \left.B\right|_{r^{2}+z^{2}\rightarrow+\infty}\rightarrow B_{\infty},\quad0\leq B_{\infty}<+\infty.\end{equation}

This corresponds to a plasma with a free boundary, i.e. confined only due to the magnetic field of the external coils. At the same time, by definition, there is no magnetic field in the interior of the bubble core, i.e. $B\equiv 0$, and hence $\psi \equiv 0$, for $r\leq r_{0}(z)$, where $r_{0}=r_{0}(z)$ is the bubble core radius. Then one can set the following internal boundary conditions:

(3.8)\begin{gather} \left.B\right|_{r=r_{0}}=0, \end{gather}

(3.9)\begin{gather}\left.\psi\right|_{r=r_{0}}=0. \end{gather}

Relations (3.7) and (3.8) together define the boundary condition for (3.4). In turn, relation (3.9) determines the bubble core boundary $r_{0}=r_{0}(z)$.

The total plasma electric current density $J_{{\rm p}\theta }$ is the sum of the warm plasma current density $J_{{\rm w}\theta }$ and the current density of the hot ions $J_{{\rm h}\theta }$. The former is to be found by solving hydrodynamic equilibrium equations for the warm plasma; as a basis, we take the MHD models constructed in previous works (see Beklemishev 2016; Khristo & Beklemishev 2019, 2022). The latter is determined by the distribution function of the hot ions, which should be found from the solution of the corresponding kinetic equation. These two issues are dealt with in the following §§ 4 and 5, respectively.

4 Warm plasma equilibrium

As mentioned above, we consider the warm plasma to be in local thermal equilibrium. We also assume the hot ion density to be negligible compared with the density of the warm plasma. This allows us to consider the warm plasma as quasi-neutral: $n_{{\rm i}}\simeq n_{{\rm e}}/Z$, where $n_{{\rm e}}$ and $n_{{\rm i}}$ are the densities of the warm electrons and ions, respectively, and $Z$ is the ion atomic number. In addition, we completely neglect the warm plasma rotation and the electrostatic potential inhomogeneity. The radial and longitudinal electric currents, as well as the azimuthal magnetic field, are set equal to zero.

4.1 Force equilibrium

The electric current density of the warm plasma $\boldsymbol {J}_{{\rm w}}$, which appears in (3.4), can be found from the force balance equation for the warm plasma⁴

(4.1)\begin{equation} \boldsymbol{\nabla} p=c^{-1}\boldsymbol{J}_{{\rm w}}\times\boldsymbol{B},\end{equation}

where $p$ is the warm plasma pressure.⁵ It can be seen that, due to the stationary axisymmetric equilibrium being considered, the sum of the azimuthal friction forces acting on the entire warm plasma is equal to zero.

It immediately follows from (4.1) that the warm plasma pressure $p$ is the constant on the flux surfaces $\psi ={\rm const}.$ and, therefore, is a function of the magnetic flux $\psi$ only: $p=p(\psi )$. In addition, due to the high longitudinal electron thermal conductivity, the temperature of the warm plasma $T$ appears to be constant along the field lines. Together with axial symmetry, this results in the temperature also being a flux function: $T=T(\psi )$. Eventually, assuming the pressure, temperature and density⁶ of the warm plasma to be related by the equation of state

(4.2)\begin{equation} p=\left(1+Z\right)n_{{\rm i}}T,\end{equation}

we arrive at the same for the warm plasma density: $n_{{\rm i}}=n_{{\rm i}}(\psi )$.

As a result, (4.1) allows the warm plasma current density to be expressed in terms of the magnetic field and pressure gradient

(4.3)\begin{equation} J_{{\rm w}\theta}=2{\rm \pi} rc\frac{{\rm d}p}{{\rm d}\psi}.\end{equation}

The magnetic flux distribution $\psi =\psi (r,z)$ is determined by Maxwell's equations (3.4), and the pressure profile $p=p(\psi )$ should be found from the solution of the mass and energy conservation equations for the warm plasma.

4.2 Mass conservation

To obtain the warm plasma equilibrium equation, we consider the material balance in a flux tube $\psi =\mathrm {const}$.

(4.4)\begin{equation} 2\varPhi_{{\rm i}\parallel{\rm m}}+\varPhi_{{\rm i}\perp}-W_{{\rm i}}=0,\end{equation}

where $\varPhi _{{\rm i}\parallel {\rm m}}$ is the axial loss of the warm ions from the flux tube, determined by the longitudinal flow through a mirror throat, $\varPhi _{{\rm i}\perp }$ is the warm ion flow transverse to the flux surface $\psi =\mathrm {const}.$ and $W_{{\rm i}}=W_{{\rm i}}(\psi )$ is the total warm ion source inside the flux tube. In other words, $W_{{\rm i}}$ represents the number of warm ion–electron pairs supplied by an external source to the interior of the surface $\psi =\mathrm {const}.$ per unit time. Longitudinal and transverse ion fluxes are determined as follows:

(4.5a,b)\begin{equation} \varPhi_{{\rm i}\parallel{\rm m}}=\int_{S_{{\rm m}}\left(\psi\right)}n_{{\rm i}}\boldsymbol{u}_{{\rm i}}\boldsymbol{\cdot} {\rm d}\boldsymbol{S},\quad\varPhi_{{\rm i}\perp}=\int_{S_{{\perp}}\left(\psi\right)}n_{{\rm i}}\boldsymbol{u}_{{\rm i}}\boldsymbol{\cdot} {\rm d}\boldsymbol{S}, \end{equation}

where $\boldsymbol {u}_{{\rm i}}$ is the warm ion flow velocity. Integration is carried out over the cross-section of the magnetic surface in the mirror $S_{{\rm m}}(\psi )$ and the flux surface $S_{\perp }(\psi )$ between the mirrors, respectively.

Gas-dynamic confinement implies the axial mirror loss being described by an outflow through the nozzle

(4.6)\begin{equation} \varPhi_{{\rm i}\parallel{\rm m}}=\int_{S_{{\rm m}}\left(\psi\right)}n_{{\rm i}}u_{{\rm m}}\,{\rm d}S,\end{equation}

where $u_{{\rm m}}\sim \sqrt {T/m_{{\rm i}}}$ is the warm plasma flow velocity in the mirror throats. The ion flow velocity transverse to the flux surface $\psi =\mathrm {const}.$ can be found from the generalized Ohm's law, which can be written in the form

(4.7)\begin{equation} \frac{1}{c}\left[\boldsymbol{u}_{{\rm i}}\times\boldsymbol{B}\right]_{\theta}=\frac{J_{{\rm w}\theta}}{\sigma _{{\rm w}}},\end{equation}

where $\sigma _{{\rm w}}$ is the electrical conductivity of the warm plasma.⁷ The force balance (4.1) together with the azimuthal projection of Ohm's law (4.7) yields the warm plasma transverse flow velocity

(4.8)\begin{equation} \boldsymbol{u}_{{\rm i}\perp}=-\frac{c^{2}}{\sigma _{{\rm w}}B^{2}}\boldsymbol{\nabla} p.\end{equation}

Then the transverse ion flux is

(4.9a,b)\begin{equation} \varPhi_{{\rm i}\perp}=-\varLambda_{{\rm i}\perp}\frac{{\rm d}p}{{\rm d}\psi},\quad \varLambda_{{\rm i}\perp}=4{\rm \pi}^{2}c^{2}n_{{\rm i}}\int_{\gamma_{{\perp}}\left(\psi\right)}\frac{r^{2}\,{\rm d}l}{ \sigma _{{\rm w}}B},\end{equation}

where integration is carried out in $(r,z)$ space along the curve $\gamma _{\perp }(\psi )$ corresponding to the flux surface $\psi =\mathrm {const}$.

Finally, substituting (4.6) and (4.9a,b) into the material balance (4.4) and differentiating it with respect to $\psi$, we obtain the mass conservation equation for the warm plasma in the following form:

(4.10)\begin{equation} {-}\frac{{\rm d}}{{\rm d}\psi}\left(\varLambda_{{\rm i}\perp}\frac{{\rm d}p}{{\rm d}\psi}\right)+\frac{2u_{{\rm m}}}{B_{{\rm m}}}n_{{\rm i}}= \frac{{\rm d}W_{{\rm i}}}{{\rm d}\psi},\end{equation}

where $B_{{\rm m}}$ is the magnetic field in the mirror throats; here, we also take into account that ${\rm d}\psi =B_{{\rm m}}\,{\rm d}S_{{\rm m}}$.

4.3 Energy conservation

For a given distribution of the magnetic flux, the equilibrium state of the warm plasma is fully described by three parameters: density $n_{{\rm i}}$, temperature $T$ and pressure $p$. A closed system of equations can be obtained by supplementing (4.2) and (4.10) with another equation relating these three parameters. In the previous MHD models (see Beklemishev 2016; Khristo & Beklemishev 2019, 2022), the plasma temperature $T$ is simply assumed to be constant. However, introducing an equation describing the energy balance of the warm plasma would be more proper.

The law of thermodynamics for a moving warm plasma element reads

(4.11)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{3}{2}p\boldsymbol{u}_{{\rm i}}\right)= \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varkappa_{{\rm w}}\boldsymbol{\nabla} T\right)+\frac{J_{{\rm w}}^{2}}{\sigma _{{\rm w}}}+q_{{\rm h}}-p\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}_{{\rm i}},\end{equation}

where $-\varkappa _{{\rm w}}\boldsymbol {\nabla } T$ is the heat flux density due to thermal conductivity with a coefficient $\varkappa _{{\rm w}}$, $q_{{\rm h}}$ is the density of the heating power from the hot ions (see expression (5.18) in § 5). The left-hand side of the equation is the change in the internal energy of the moving element with time. The first three terms on the right-hand side describe the total heat release, and the last term is related to the mechanical work. Here, we consider the heating to be due to the neutral beam injection. If necessary, heating from other sources can also be taken into account by simply adding the corresponding terms to the right-hand side of the equation.

Equations (4.1) and (4.8) together yield the resistive heating power density

(4.12)\begin{equation} \frac{J_{{\rm w}}^{2}}{\sigma _{{\rm w}}}=-\boldsymbol{u}_{{\rm i}}\boldsymbol{\cdot}\boldsymbol{\nabla} p. \end{equation}

Substituting it into (4.11) we arrive at

(4.13)\begin{equation} \tfrac{5}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\,p\boldsymbol{u}_{{\rm i}}\right)- \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\varkappa_{{\rm w}}\boldsymbol{\nabla} T\right)=q_{{\rm h}}. \end{equation}

Integrating this equation over the volume of the flux tube $\psi =\mathrm {const}.$ between the mirrors results in

(4.14)\begin{equation} 2\varPhi_{{\rm E}\parallel}+\varPhi_{{\rm E}\perp}=Q_{{\rm h}},\end{equation}

where $\varPhi _{{\rm E}\parallel }$ is the axial energy loss from the flux tube

(4.15)\begin{equation} \varPhi_{{\rm E}\perp}=\int_{S_{{\perp}}\left(\psi\right)}\frac{5}{2}p\boldsymbol{u}_{{\rm i}}\boldsymbol{\cdot} {\rm d}\boldsymbol{S}-\int_{S_{{\perp}}\left(\psi\right)}\varkappa_{{\rm w}}\boldsymbol{\nabla} T\boldsymbol{\cdot} {\rm d}\boldsymbol{S} \end{equation}

is the energy flow transverse to the flux surface $\psi =\mathrm {const}.$, and $Q_{{\rm h}}$ is the total heating power from the hot ions inside the flux tube.

The axial energy loss in a GDT can be described as follows:

(4.16)\begin{equation} \varPhi_{{\rm E}\parallel}=\alpha_{{\rm E}}\int_{S_{{\rm m}}\left(\psi\right)}Tn_{{\rm i}}u_{{\rm m}}\,{\rm d}S. \end{equation}

In other words, an ion–electron pair escaping the trap carries away energy $\alpha _{{\rm E}}T$ through the mirror (see Ryutov 2005; Skovorodin 2019). The factor $\alpha _{{\rm E}}$ depends on the nature of the plasma outflow through the mirror; for the GDT facility in Budker INP, this parameter lies in the range of $6\div 8$ (see Soldatkina et al. 2020). Further, taking into account (4.8), the transverse energy flow is written in the form

(4.17)\begin{equation} \varPhi_{{\rm E}\perp}=\frac{5}{2}\frac{p}{n_{{\rm i}}}\varPhi_{{\rm i}\perp}-\varLambda_{{\rm E}\perp}\frac{{\rm d}T}{{\rm d}\psi}, \end{equation}

where

(4.18)\begin{equation} \varLambda_{{\rm E}\perp}=4{\rm \pi}^{2}\int_{\gamma_{{\perp}}\left(\psi\right)}\varkappa_{{\rm w}}Br^{2}\,{\rm d}l. \end{equation}

As a result, after taking the derivative of the energy balance (4.14) with respect to $\psi$, we arrive at the energy conservation equation for the warm plasma in the following form:

(4.19)\begin{equation} {-}\frac{{\rm d}}{{\rm d}\psi}\left(\frac{5}{2}\frac{p}{n_{{\rm i}}}\varLambda_{{\rm i}\perp} \frac{{\rm d}p}{{\rm d}\psi}\right)-\frac{{\rm d}}{{\rm d}\psi} \left(\varLambda_{{\rm E}\perp}\frac{{\rm d}T}{{\rm d}\psi}\right)+ \frac{2\alpha_{{\rm E}}u_{{\rm m}}}{B_{{\rm m}}}n_{{\rm i}}T=\frac{{\rm d}Q_{{\rm h}}}{{\rm d}\psi}.\end{equation}

4.4 Bubble core equilibrium

Strictly speaking, (4.3), (4.10) and (4.19) cannot be used to describe the equilibrium of the warm plasma in the bubble core, where $\psi \rightarrow 0$. In addition, the magnetized plasma approximation obviously ceases to be applicable there. For this reason, the equilibrium of the warm plasma inside the bubble core should be considered separately.

Since the magnetic field inside the bubble core is close to zero, the warm plasma transport in the core should increase significantly compared with the external transverse diffusion. Therefore, we assume the pressure, the temperature and the density of the warm plasma to be uniform inside the entire bubble core. In addition, we further suppose the warm plasma in the bubble core to be perfectly conducting, which results in the magnetic field there being constant and identically equal to zero

(4.20)\begin{equation} B\equiv0,\quad r< r_{0}\left(z\right). \end{equation}

This is possible only if there is no plasma current inside the core, i.e. the electric current of the hot ions is completely compensated by the inductive current of the warm plasma

(4.21)\begin{equation} J_{{\rm w}\theta}=-J_{{\rm h}\theta},\quad r< r_{0}\left(z\right). \end{equation}

It is shown by Chernoshtanov (2020, 2022) that the so-called non-adiabatic loss should arise in the diamagnetic trap. In the case of a sufficiently collisional warm plasma with an isotropic distribution function, there are particles outside the absolute confinement region (2.2). These particles may reach the mirrors and escape the trap directly from the bubble core. The thermal equilibrium distribution of the warm ions in a diamagnetic trap can be approximately represented as follows (see Chernoshtanov 2020, 2022):

(4.22)\begin{equation} f_{{\rm i}}\left(\mathcal{E},\mathcal{P}\right)=n_{{\rm i}0}\left(\frac{m_{{\rm i}}}{2{\rm \pi} T_{0}}\right)^{{3}/{2}}{\rm e}^{-{\mathcal{E}}/{T_{0}}}\varTheta\left(a-\frac{\left|\mathcal{P}\right|}{\sqrt{2m_{{\rm i}} \mathcal{E}}}\right),\end{equation}

where $n_{{\rm i}0}$ and $T_{0}$ are the warm plasma density and temperature in the bubble core, $\varTheta (x)$ is the Heaviside step function, $a=\max [r_{0}(z)]=r_{0}(0)$ is the maximum radius of the core and $\mathcal {P}$ and $\mathcal {E}$ are the canonical angular momentum and the total energy, respectively, defined as in (2.1a,b). Such a distribution function corresponds to the collisional plasma when the characteristic Maxwellization time is much less than the confinement time. It can be obtained by direct averaging of the ‘ordinary’ Maxwellian distribution over the hypersurface of the constant $\mathcal {P}$ and $\mathcal {E}$ (see the definition of averaging (5.13) in § 5). The total axial loss from the core of the bubble is determined by the flow through the mirrors

(4.23)\begin{equation} \varPhi_{{\rm i}\parallel0}\simeq\frac{3}{4}n_{{\rm i}0}v_{T{\rm i}0}\frac{2{\rm \pi} a\rho_{{\rm i}0}}{\mathcal{R}},\end{equation}

where $v_{T{\rm i}0}=\sqrt {T_{0}/m_{{\rm i}}}$ is the warm ion thermal velocity, $\rho _{{\rm i}0}=v_{T{\rm i}0}m_{{\rm i}}c/eZB_{{\rm v}}$ is the characteristic warm ion Larmor radius in the vacuum field $B_{{\rm v}}$, $\mathcal {R}=B_{{\rm m}}/B_{{\rm v}}$ is the vacuum mirror ratio and $B_{{\rm m}}\simeq \mathrm {const}.$ is the mirror magnetic field. The expression (4.23) is obtained in the limit $a\gg \rho _{{\rm i}0}/\mathcal {R}$; for details, refer to Appendix A.

The non-adiabatic loss (4.23) is to be taken into account when setting the internal boundary condition for the mass conservation equation (4.10). Namely, the total transverse ion flow through the boundary of the bubble core is the total warm ion source inside the core $W_{{\rm i}0}$ minus the non-adiabatic loss through the mirrors $2\varPhi _{{\rm i}\parallel 0}$

(4.24)\begin{equation} \left.\left(-\varLambda_{{\rm i}\perp}\frac{{\rm d}p}{{\rm d}\psi}\right)\right|_{\psi\rightarrow0}=W_{{\rm i}0}- 2\varPhi_{{\rm i}\parallel0}.\end{equation}

The corresponding axial energy loss is defined as follows:

(4.25)\begin{equation} \varPhi_{{\rm E}\parallel0}=\alpha_{{\rm E}0}T_{0}\varPhi_{{\rm i}\parallel0}, \end{equation}

where $\alpha _{{\rm E}0}$ is a coefficient similar to the previously introduced $\alpha _{{\rm E}}$ for gas-dynamic energy loss. Then, the internal boundary condition for the energy conservation equation (4.19) is

(4.26)\begin{equation} \left.\left(-\frac{5}{2}\frac{p}{n_{{\rm i}}}\varLambda_{{\rm i}\perp}\frac{{\rm d}p}{{\rm d}\psi}- \varLambda_{{\rm E}\perp}\frac{{\rm d}T}{{\rm d}\psi}\right)\right|_{\psi\rightarrow0}=Q_{{\rm h}0}-2\alpha_{{\rm E}0}T_{0} \varPhi_{{\rm i}\parallel0},\end{equation}

where $Q_{{\rm h}0}$ is the total heat release from the hot ions inside the bubble core.

To complete the formulation of the boundary value problem, the external boundary conditions are to be set. As such, for example, the conditions

(4.27a,b)\begin{equation} \left.n_{{\rm i}}\right|_{r=a_{\mathrm{lim}}}=0,\quad\left.T\right|_{r=a_{\mathrm{lim}}}=0\end{equation}

can be chosen, which correspond to an external limiter located at the radius $a_{\mathrm {lim}}$. If the particle source is localized inside the bubble core, as we assume in the present paper, due to the large axial loss, the actual boundary of the warm plasma typically does not reach the limiter. This means that, in this case, the equilibrium appears to be weakly dependent on the external boundary conditions (4.27a,b).

5 Hot ion equilibrium

Consider the dynamics of a single ion with charge $Ze$ and mass $m_{{\rm i}}$ in a given axisymmetric magnetic field, which is determined by the magnetic flux distribution $\psi =\psi (r,z)$ (see § 3). Then the equations of motion have the form

(5.1)\begin{equation} \left.\begin{gathered} \dot{r} =\frac{P_{r}}{m_{{\rm i}}}, \quad \dot{P}_{r} =\frac{P_{\theta}-\varPsi}{m_{{\rm i}}r^{2}} \left(\frac{P_{\theta}-\varPsi}{r}+\partial_{r}\varPsi\right),\\ \dot{\theta} =\frac{P_{\theta}-\varPsi}{m_{{\rm i}}r^{2}}, \quad \dot{P}_{\theta} =0,\\ \dot{z} =\frac{P_{z}}{m_{{\rm i}}}, \quad \dot{P}_{z} =\frac{P_{\theta}-\varPsi}{m_{{\rm i}}r^{2}}\partial_{z}\varPsi, \end{gathered}\right\} \end{equation}

where $(r,\theta,z,P_{r},P_{\theta },P_{z})$ is a set of canonically conjugate Hamiltonian variables, overdot indicates a time derivative and, for brevity, the normalized magnetic flux $\varPsi =Ze\psi /2{\rm \pi} c$ is introduced.

There are two global regular integrals of motion: the total energy

(5.2)\begin{equation} \mathcal{E}=\frac{P_{r}^{2}}{2m_{{\rm i}}}+\frac{\left(P_{\theta}-\varPsi\right)^{2}}{2m_{{\rm i}}r^{2}}+\frac{P_{z}^{2}}{2m_{{\rm i}}}, \end{equation}

and the angular momentum

(5.3)\begin{equation} \mathcal{P}=P_{\theta}=m_{{\rm i}}r^{2}\dot{\theta}+\varPsi. \end{equation}

Excluding the cyclic variable $\theta$, for a given $\mathcal {P}$, the number of degrees of freedom is reduced to two, and the ion dynamics is equivalent to the motion of a particle with energy $\mathcal {E}$ in two-dimensional $(r,z)$ space in the effective potential

(5.4)\begin{equation} \varphi_{\mathrm{eff}}=\frac{\left(\mathcal{P}-\varPsi\right)^{2}}{2m_{{\rm i}}r^{2}}.\end{equation}

Therefore, the conservation of $\mathcal {P}$ and $\mathcal {E}$ leads to the ion moving in a bounded area: $\mathcal {E}-\varphi _{\mathrm {eff}}>0.$

For given $\mathcal {E}$ and $\mathcal {P}$ the equality $\mathcal {E}=\varphi _{\mathrm {eff}}$ determines a curve in $(r,z)$ space corresponding to the boundary of the region of ion motion. Then, taking into account that the magnetic flux in the mirrors is approximately $\psi _{{\rm m}}=B_{{\rm m}}{\rm \pi} r^{2}$, we can find the radial boundaries in the mirror throats from the equation

(5.5)\begin{equation} \left(\frac{\mathcal{R}\varOmega_{{\rm i}}}{2}\right)^{2}\left(m_{{\rm i}}r^{2}\right)^{2}-2 \left(\frac{\mathcal{R}\varOmega_{{\rm i}}}{2}\mathcal{P}+\mathcal{E}\right)\left(m_{{\rm i}}r^{2}\right)+\mathcal{P}^{2}=0, \end{equation}

where $\varOmega _{{\rm i}}=ZeB_{{\rm v}}/m_{{\rm i}}c$ is the ion cyclotron frequency in the vacuum magnetic field $B_{{\rm v}}$. If this equation has no real roots, then the region of ion motion does not reach the mirrors, and the ion is confined absolutely. This brings us to the absolute confinement criterion (2.2) for the ions

(5.6)\begin{equation} \left(\frac{\mathcal{R}\varOmega_{{\rm i}}}{2}\mathcal{P}+\mathcal{E}\right)^{2}-\left(\frac{\mathcal{R} \varOmega_{{\rm i}}}{2}\right)^{2}\mathcal{P}^{2}<0,\quad\Rightarrow\quad-\mathcal{R} \varOmega_{{\rm i}}\mathcal{P}>\mathcal{E}. \end{equation}

Otherwise, an ion outside the absolute confinement region may escape the trap in a finite time.

If there is a third conserved quantity in addition to $\mathcal {E}$ and $\mathcal {P}$, the system is integrable, and the dynamics of the ion is regular. In particular, it was found by Chernoshtanov (2020, 2022) that, for an ion with not too high longitudinal velocity, satisfying the condition (2.3), there is an adiabatic invariant $I_{r}=\oint p_{r}\,\textrm {d}r$. Therefore, the trajectory of such an ion is fully described by the given $\mathcal {E}$, $\mathcal {P}$ and $I_{r}$. Otherwise, if the criterion (2.3) is violated, there is no adiabatic invariant $I_{r}$ and the dynamics of the ion becomes chaotic. This means that the ion trajectory ergodically fills a finite volume of phase space on the invariant hypersurface of constant $\mathcal {E}$ and $\mathcal {P}$ (Sagdeev, Usikov & Zaslavsky 1988; Lichtenberg & Lieberman 1992).

The chaotic behaviour of ions can be qualitatively explained by collisionless scattering on the longitudinal inhomogeneities of the magnetic field on the bubble boundary. ‘Reflecting’ from the magnetic field at the boundary of the bubble leads to a change in the pitch angle of an ion by approximately $\Delta \xi \sim \arctan |\textrm {d}r_{0}/\textrm {d}z|$. Therefore, the adiabaticity criterion (2.3) is met when the maximum scattering angle is not too large: $(\Delta \xi )_{\max }/\xi \lesssim 1$, while large-angle scattering: $(\Delta \xi )_{\max }/\xi \gtrsim 1$, on the considerable field inhomogeneities: $|\textrm {d}r_{0}/\textrm {d}z|\sim 1$, may result in violation of adiabaticity and hence to dynamic chaos. With this in mind, we can also suppose that the characteristic time of ergodization appears to be of the order of several free pass times between the large-scale inhomogeneities: $|\textrm {d}r_{0}/\textrm {d}z|\sim 1$.

It is useful to point out the following fundamental feature of the dynamics of individual particles in a diamagnetic trap. Particles with a sufficiently large canonical angular momentum $\mathcal {P}$ satisfying the criterion (2.2) remain absolutely confined even in the absence of adiabaticity, when condition (2.3) is violated. This is different from ‘classical’ mirror machines, where non-adiabatic particles are lost over a time scale of the order of several bounce periods.

In the presence of non-paraxial regions, where $|\textrm {d}r_{0}/\textrm {d}z|\sim 1$, which seems to be typical for the diamagnetic bubble equilibrium (see Kotelnikov et al. 2010; Kotelnikov 2011; Beklemishev 2016; Khristo & Beklemishev 2019, 2022), the adiabaticity criterion (2.3) seems to be violated for the majority of the injected ions. In particular, figure 1 illustrates an example of the ion Poincaré map in the central plane, $z=0$, for various initial conditions and fixed energy $\mathcal {E}$ and angular momentum $\mathcal {P}$. The magnetic field distribution is taken from the MHD simulations of the diamagnetic bubble equilibrium in GDMT (see Khristo & Beklemishev 2022). Vacuum magnetic field in the central plane is $B_{{\rm v}}\simeq 1.5\ \mathrm {T}$, the ion Larmor radius in the vacuum magnetic field is $\rho _{\mathcal {E}}=\sqrt {2m_{{\rm i}}\mathcal {E}}c/ZeB_{{\rm v}}\simeq 1.6\ \mathrm {cm}$ and the minimum distance from the trap axis that the ion can approach is $r_{\min }=|\mathcal {P}|/\sqrt {2m_{{\rm i}}\mathcal {E}}\simeq 3.5\ \mathrm {cm}$. It can be seen that a considerable region, corresponding to transverse momentum in the range $P_{\perp }=\sqrt {2m_{{\rm i}}\mathcal {E}-P_{z}^{2}}\lesssim 0.6\div 0.7$, is filled with chaotic trajectories.

Figure 1. An example of the ion Poincaré map in the central plane, $z=0$, for various initial conditions and fixed $\mathcal {E}$ and $\mathcal {P}$. The magnetic field distribution is taken from MHD simulations for the GDMT configuration (see Khristo & Beklemishev 2022). Vacuum magnetic field: $B_{{\rm v}}\simeq 1.5\ \mathrm {T}$; ion Larmor radius in the vacuum field: $\rho _{\mathcal {E}}=\sqrt {2m_{{\rm i}}\mathcal {E}}c/ZeB_{{\rm v}}\simeq 1.6\ \mathrm {cm}$; minimum distance from the trap axis that the ion can approach: $r_{\min }=|\mathcal {P}|/\sqrt {2m_{{\rm i}}\mathcal {E}}\simeq 3.5\ \mathrm {cm}$.

Strictly speaking, the ion dynamics should be considered separately for the chaotic and regular trajectories. In the chaotic region, the distribution function can be considered depending only on the integrals of motion: energy $\mathcal {E}$ and angular momentum $\mathcal {P}$. At the same time, in the region of regular motion, these are also supplemented by the dependence on a third conserved quantity, for instance, the adiabatic invariant $I_{r}$. On top of that, required is also a correct description of the transition region, separating chaotic and regular trajectories. This issue appears to be quite complex and needs to be addressed in an individual paper. For this reason, in the present article, we limit ourselves to considering the simplest case. Namely, we further assume that the region of the hot ion regular motion has a negligible measure, while the dynamics of the hot ions is globally chaotic and their trajectories ergodically fill the hypersurfaces of constants $\mathcal {E}$ and $\mathcal {P}$. Since the chaotic behaviour appears to result from the collisionless scattering on the non-paraxial end regions, $|\textrm {d}r_{0}/\textrm {d}z|\sim 1$, we also assume that the characteristic ergodization time is of the order of several periods of longitudinal oscillations $\tau _{{\rm b}}\sim L/\sqrt {\mathcal {E}/m_{{\rm i}}}$.

5.1 Hot ion distribution function

As noted above, we consider the hot ion density to be small enough to neglect their collisions with each other. In addition, we neglect the angular scattering of the hot ions due to Coulomb collisions⁸ and take into account only the warm electron drag force

(5.7)\begin{equation} \boldsymbol{F}_{{\rm s}}\simeq-\nu_{\mathrm{ie}}m_{{\rm i}}\boldsymbol{v},\end{equation}

where $\boldsymbol {v}$ is the hot ion velocity,

(5.8)\begin{equation} \nu_{\mathrm{ie}}=1.6\times10^{-9}\mu_{{\rm i}}^{-1}Z^{3}\varLambda_{\mathrm{ie}}\frac{n_{{\rm i}}}{{\rm cm}^{-3}}\left(\frac{T}{{\rm eV}}\right)^{-{3}/{2}}{\rm s}^{-1},\end{equation}

is the inverse ion slowing down time (Trubnikov 1965), $\mu _{{\rm i}}=m_{{\rm i}}/m_{{\rm p}}$ is the ion mass $m_{{\rm i}}$ normalized to proton mass $m_{{\rm p}}$, $\varLambda _{\mathrm {ie}}$ is the ion–electron Coulomb logarithm. In what follows, we assume that the neutral beam is injected into the absolute confinement region (2.2); in this case, the loss of the hot ions can be ignored. Therefore, in the stationary case under consideration, the kinetic equation for the distribution function of the hot ions $f_{{\rm h}}$ has the following invariant form:

(5.9)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{X}}\boldsymbol{\cdot}\left(\dot{\boldsymbol{X}}f_{{\rm h}}\right)=g_{{\rm h}},\end{equation}

where $\boldsymbol {X}$ is the six-dimensional vector of generalized phase variables, $\dot {\boldsymbol {X}}$ is the corresponding phase velocity vector, $\boldsymbol {\nabla }_{\boldsymbol {X}}$ is the nabla differential operator in the $\boldsymbol {X}$-space and $g_{{\rm h}}=g_{{\rm h}}(\boldsymbol {X})$ is the phase density of the hot ion source intensity, the explicit form of which is determined by the injection. Specifying the variables as $\boldsymbol {X}=(r,\theta,z,P_{r},P_{\theta },P_{z})$, we obtain the equations of motion in the form

(5.10)\begin{equation} \left.\begin{gathered} \dot{r} =\frac{P_{r}}{m_{{\rm i}}}, \quad \dot{P}_{r} =\frac{P_{\theta}-\varPsi}{m_{{\rm i}}r^{2}} \left(\frac{P_{\theta}-\varPsi}{r}+\partial_{r}\varPsi\right)-\nu_{\mathrm{ie}}P_{r},\\ \dot{\theta} =\frac{P_{\theta}-\varPsi}{m_{{\rm i}}r^{2}}, \quad \dot{P}_{\theta} =- \nu_{\mathrm{ie}}\left(P_{\theta}-\varPsi\right),\\ \dot{z} =\frac{P_{z}}{m_{{\rm i}}}, \quad \dot{P}_{z} =\frac{P_{\theta}-\varPsi}{m_{{\rm i}}r^{2}}\partial_{z}\varPsi-\nu_{\mathrm{ie}}P_{z}. \end{gathered}\right\} \end{equation}

Since we assume axial symmetry, the distribution function $f_{{\rm h}}$ and the magnetic flux $\varPsi$ do not depend on the cyclic variable $\theta$.

As mentioned above, we consider that the dynamics of ions is globally chaotic, and the ergodization time is of the order of several periods of longitudinal oscillations of the ion $\tau _{{\rm b}}\sim L/\sqrt {\mathcal {E}/m_{{\rm i}}}$. Since the characteristic slowing down time $\nu _{\mathrm {ie}}^{-1}$ significantly exceeds the characteristic bounce time $\tau _{{\rm b}}$ (typically $\nu _{\mathrm {ie}}\tau _{{\rm b}}\sim 10^{-3}\div 10^{-4}$), the drag force (5.7) can be considered weak on time scales of the period of longitudinal oscillations $\tau _{{\rm b}}$. By means of the conventionally used Krylov–Bogoliubov–Mitropolsky averaging approach (Bogoliubov & Mitropolskii 1961), the slow phase space diffusion associated with the drag force (5.7) can be separated from the fast longitudinal oscillations. Based on this, we replace the exact kinetic equation (5.9) by an averaged one, assuming that the averaging should be performed over the $\mathcal {E},\mathcal {P}=\textrm {const}.$ surface in the phase space.

Instead of $(P_{r},P_{\theta },P_{z})$, it is convenient to use the new variables $(\mathcal {E},\mathcal {P},\alpha )$, where $\alpha \in [0,2{\rm \pi} )$ is defined by

(5.11a,b)\begin{equation} P_{r}=\sqrt{2m_{{\rm i}}\left(\mathcal{E}-\varphi_{\mathrm{eff}}\right)}\sin\alpha,\quad P_{z}=\sqrt{2m_{{\rm i}}\left(\mathcal{E}-\varphi_{\mathrm{eff}}\right)}\cos\alpha. \end{equation}

Hence, we have

(5.12)\begin{equation} \frac{\partial\left(P_{r},P_{\theta},P_{z}\right)}{\partial\left(\mathcal{E},\mathcal{P},\alpha\right)}=m_{{\rm i}}, \end{equation}

and we should also keep in mind the additional condition $\mathcal {E}>\varphi _{\mathrm {eff}}$, which defines the region of permissible values of the new variables. Then averaging of some quantity $\mathcal {Q}$ over the surface $\mathcal {E},\mathcal {P}=\textrm {const}.$ is performed as follows:

(5.13)\begin{equation} \left\langle \mathcal{Q}\right\rangle _{\varGamma}\overset{\mathrm{def}}=\varGamma^{-1}\idotsint_{\mathcal{E}>\varphi_{\mathrm{eff}}} \mathcal{Q}m_{\mathrm{i}}\,{\rm d}r\,{\rm d}\theta \,{\rm d}z\,{\rm d}\alpha,\end{equation}

where

(5.14)\begin{equation} \varGamma=\idotsint _{\mathcal{E}>\varphi_{\mathrm{eff}}}m_{{\rm i}}\,{\rm d}r\,{\rm d}\theta \,{\rm d}z\,{\rm d}\alpha=4{\rm \pi}^{2}m_{{\rm i}}\iint _{\mathcal{E}>\varphi_{\mathrm{eff}}}\,{\rm d}r\,{\rm d}z, \end{equation}

is the phase volume of the hypersurface $\mathcal {E},\mathcal {P}=\mathrm {const}$. Finally, at the leading order with respect to $\nu _{\mathrm {ie}}\tau _{{\rm b}}\ll 1$, when the distribution function $f_{{\rm h}}$ approximately depends only on $\mathcal {E}$ and $\mathcal {P}$, averaging the perturbed system yields

(5.15)\begin{equation} \left.\begin{gathered} \varGamma^{-1}\partial_{\mathcal{E}}\left(\varGamma\left\langle \dot{\mathcal{E}}\right\rangle _{\varGamma}f_{{\rm h}} \right)+\varGamma^{-1}\partial_{\mathcal{P}}\left(\varGamma\left\langle \dot{\mathcal{P}}\right\rangle _{\varGamma}f_{{\rm h}} \right)=\left\langle g_{{\rm h}}\right\rangle _{\varGamma},\\ \left\langle \dot{\mathcal{E}}\right\rangle _{\varGamma}=-2\left\langle \nu_{\mathrm{ie}}\right\rangle _{\varGamma} \mathcal{E},\quad\left\langle \dot{\mathcal{P}}\right\rangle _{\varGamma}=-\left\langle \nu_{\mathrm{ie}}\right\rangle _{\varGamma}\mathcal{P}+\left\langle \nu_{\mathrm{ie}}\varPsi\right\rangle _{\varGamma}.\\ \left\langle \dot{r}\right\rangle _{\varGamma}=0,\quad\left\langle \dot{z}\right\rangle _{\varGamma}=0,\quad\left\langle \dot{\alpha}\right\rangle _{\varGamma}=0. \end{gathered}\right\}\end{equation}

Solving the kinetic equation (5.15) one can find the averaged distribution function of the hot ions $f_{{\rm h}}\simeq f_{{\rm h}}(\mathcal {E},\mathcal {P})$. This in turn enables the density of some quantity $\mathcal {Q}$ to be obtained

(5.16)\begin{equation} \left\langle \mathcal{Q}\right\rangle _{{\rm h}}\overset{\mathrm{def}}=\frac{1}{r}\iiint\mathcal{Q}f_{{\rm h}}\,{\rm d}P_{r}\,{\rm d}P_{\theta}\, {\rm d}P_{z}=\frac{m_{{\rm i}}}{r}\int_{0}^{2{\rm \pi}}{\rm d}\alpha\int _{-\infty}^{+\infty}{\rm d}\mathcal{P}\int_{\varphi_{\mathrm{eff}}}^{+\infty}{\rm d}\mathcal{E} \mathcal{Q}f_{{\rm h}}. \end{equation}

In particular, the azimuthal electric current density of the hot ions is given by

(5.17)\begin{equation} J_{{\rm h}\theta}=\left\langle Zev_{\theta}\right\rangle _{{\rm h}}=\frac{2{\rm \pi} Ze}{r^{2}}\int _{-\infty}^{+\infty}{\rm d}\mathcal{P}\int_{\varphi_{\mathrm{eff}}}^{+\infty}{\rm d}\mathcal{E}\left(\mathcal{P}-\varPsi\right)f_{{\rm h}}, \end{equation}

the power density of heating the warm plasma by the hot ions is defined as

(5.18)\begin{equation} q_{{\rm h}}=\left\langle 2\left\langle \nu_{\mathrm{ie}}\right\rangle _{\varGamma}\mathcal{E}\right\rangle _{{\rm h}}=\frac{4{\rm \pi} m_{{\rm i}}}{r}\int _{-\infty}^{+\infty}{\rm d}\mathcal{P}\int _{\varphi_{\mathrm{eff}}}^{+\infty}{\rm d}\mathcal{E}\left\langle \nu_{\mathrm{ie}}\right\rangle _{\varGamma}\mathcal{E}f_{{\rm h}}.\end{equation}

6 Cylindrical bubble model

The complete system of equations describing the equilibrium of the diamagnetic bubble with neutral beam injection consists of the equation for the magnetic field (3.4) and the equations of mass (4.10) and energy (4.19) conservation for the warm plasma. This system should be supplemented with the equation of state for the warm plasma (4.2) and the expression for the warm plasma current density (4.3), as well as the expressions for the hot ion current density (5.17) and heating power density (5.18). Finally, the formulation of the problem is completed by setting the boundary conditions for the magnetic field: (3.7), (3.8) and (3.9), the latter of which determines the boundary of the bubble core, and also the internal (4.24), (4.26) and external (4.27a,b) conditions for the warm plasma equilibrium equations.

Just as it was done by Beklemishev (2016) for the case of MHD equilibrium, we further consider a simplified model of a cylindrical diamagnetic bubble. In the central part of such a bubble, there is a cylindrical core of length $L$ and radius $a$ with zero magnetic field: $\psi \equiv 0$; outside the core, for $r>a$, the magnetic field lines are straight: $\psi =\psi (r)$. At the ends of the central cylindrical part, there are non-paraxial regions gradually turning into the mirror throats with the magnetic field $B_{{\rm m}}\simeq \mathrm {const}$. In addition, the source of the warm plasma is further considered to be entirely contained inside the bubble core. In this case, the density of the warm plasma inside the core is expected to be much higher than outside. Since the drag force (5.7) is proportional to the density of the warm plasma, the hot ions seem to slow down and release the energy mainly inside the bubble core as well.

6.1 Magnetic field distribution

The equilibrium equation for the magnetic field (3.4) can only be considered in the central cylindrical region, and the magnetic flux in the mirrors is considered to be given: $\psi _{{\rm m}}=B_{{\rm m}}{\rm \pi} r^{2}$. Then, after substituting the current density of the warm plasma (4.3), (3.4) predictably reduces to the paraxial equilibrium with the Lorentz force from the hot ions

(6.1)\begin{gather} \frac{{\rm d}}{{\rm d}\psi}\left(\frac{B^{2}}{8{\rm \pi}}+p\right)=-\frac{J_{{\rm h}\theta}}{2{\rm \pi} rc}, \end{gather}

(6.2)\begin{gather}\frac{{\rm d}r}{{\rm d}\psi}=\frac{1}{2{\rm \pi} r}\frac{1}{B}. \end{gather}

In this equations we use the flux coordinate $\psi$, so $r=r(\psi )$ is the inverse function of $\psi =\psi (r)$ and is essentially the radius of the magnetic surface that corresponds to the magnetic flux $\psi$. In addition, the magnetic field $B$ and the current density of the hot ions $J_{{\rm h}\theta }$ should also be considered here as functions of $\psi$. The corresponding boundary conditions (3.7), (3.8) and (3.9) can be reduced to

(6.3a,b)\begin{gather} \left.B\right|_{\psi=0}=0,\quad\left.r\right|_{\psi=0}=a, \end{gather}

(6.4)\begin{gather}\left.B\right|_{\psi\rightarrow+\infty}=B_{{\rm v}}, \end{gather}

where $B_{{\rm v}}=\mathrm {const}.$ is the vacuum magnetic field at the central section of the trap.

6.2 Warm plasma equilibrium

Warm plasma equilibrium equations (4.10) and (4.19) mainly remain the same except for the transport coefficients $\varLambda _{{\rm i}\perp }$ and $\varLambda _{{\rm E}\perp }$, which are reduced to a simpler form in the cylindrical bubble approximation

(6.5)\begin{gather} -\frac{{\rm d}}{{\rm d}\psi}\left(\varLambda_{{\rm i}\perp}\frac{{\rm d}p}{{\rm d}\psi}\right)+ \frac{2u_{{\rm m}}}{B_{{\rm m}}}n_{{\rm i}}=0, \end{gather}

(6.6)\begin{gather}-\frac{{\rm d}}{{\rm d}\psi}\left(\frac{5}{2}\frac{p}{n_{{\rm i}}}\varLambda_{{\rm i}\perp}\frac{{\rm d}p}{{\rm d}\psi}\right)- \frac{{\rm d}}{{\rm d}\psi}\left(\varLambda_{{\rm E}\perp}\frac{{\rm d}T}{{\rm d}\psi}\right)+ \frac{2\alpha_{{\rm E}}u_{{\rm m}}}{B_{{\rm m}}}n_{{\rm i}}T=0, \end{gather}

(6.7a,b)\begin{gather}\varLambda_{{\rm i}\perp}\simeq4{\rm \pi}^{2}c^{2}n_{{\rm i}}\frac{r^{2}L}{\sigma _{{\rm w}}B}, \quad\varLambda_{{\rm E}\perp}\simeq4{\rm \pi}^{2}\varkappa_{{\rm w}}Br^{2}L. \end{gather}

Here, we also take into account that, due to the hot ions mainly slowing down inside the bubble core, the external hot ion heating power release is negligible, i.e. $Q_{{\rm h}}\equiv 0$ and $W_{{\rm i}}\equiv 0$ for $\psi >0$. The boundary conditions (4.24), (4.26) and (4.27a,b), in turn, remain exactly the same, except that $Q_{{\rm h}0}$ and $W_{{\rm i}0}$ now have the meaning of the total absorbed injection power and the total warm plasma source, respectively.

6.3 Hot ion equilibrium

Consider the injection of a monoenergetic neutral beam with the injection energy $\mathcal {E}_{\mathrm {NB}}$ and the total absorbed injection power $Q_{{\rm h}0}$. The injection is carried out at the angle $\xi _{\mathrm {NB}}\in (0,{\rm \pi} /2)$ to the axis of the trap, and the distance from the beam to the trap axis is finite and equal to $r_{\mathrm {NB}}< a$. In this case, the source in the kinetic equation (5.15) has the form⁹

(6.8)\begin{equation} \left\langle g_{{\rm h}}\right\rangle _{\varGamma}=\frac{1}{\varGamma}\frac{Q_{{\rm h}0}}{\mathcal{E}_{\mathrm{NB}}} \delta\left(\mathcal{E}-\mathcal{E}_{\mathrm{NB}}\right)\delta\left(\mathcal{P}-\mathcal{P}_{\mathrm{NB}}\right),\end{equation}

where $\delta (x)$ is the Dirac delta function

(6.9)\begin{equation} \mathcal{P}_{\mathrm{NB}}=-\sqrt{2m_{{\rm i}}\mathcal{E}_{\mathrm{NB}}}r_{\mathrm{NB}}\sin\xi_{\mathrm{NB}}<0,\end{equation}

is the angular momentum of the injected ions. Due to the beam slowing down mainly in the core of the bubble, averaging in (5.15) yields

(6.10a,b)\begin{equation} \left\langle \nu_{\mathrm{ie}}\right\rangle _{\varGamma}\simeq\nu_{\mathrm{ie}0},\quad\left\langle \nu_{\mathrm{ie}}\varPsi\right\rangle_{\varGamma}\simeq0, \end{equation}

where $\nu _{\mathrm {ie}0}=\nu _{\mathrm {ie}}|_{T=T_{0},n_{{\rm i}}=n_{{\rm i}0}}$. The resulting kinetic equation

(6.11)\begin{equation} {-}2\nu_{\mathrm{ie}0}\partial_{\mathcal{E}}\left(\mathcal{E}\varGamma f_{{\rm h}}\right)-\nu_{\mathrm{ie}0}\partial_{\mathcal{P}}\left(\mathcal{P}\varGamma f_{{\rm h}}\right)\simeq\frac{Q_{{\rm h}0}}{\mathcal{E}_{\mathrm{NB}}}\delta \left(\mathcal{E}-\mathcal{E}_{\mathrm{NB}}\right)\delta\left(\mathcal{P}-\mathcal{P}_{\mathrm{NB}}\right), \end{equation}

is satisfied by

(6.12a,b)\begin{equation} f_{{\rm h}}\simeq\frac{1}{2\nu_{\mathrm{ie}0}}\frac{1}{\varGamma}\frac{Q_{{\rm h}0}}{\mathcal{E}_{\mathrm{NB}}} \frac{1}{\mathcal{E}^{{3}/{2}}}\delta\left(\frac{\mathcal{P}}{\sqrt{\mathcal{E}}}- \frac{\mathcal{P}_{\mathrm{NB}}}{\sqrt{\mathcal{E}_{\mathrm{NB}}}}\right),\quad \mathcal{E}<\mathcal{E}_{\mathrm{NB}},\end{equation}

where the phase volume $\varGamma$ is approximately equal to

(6.13)\begin{equation} \varGamma=4{\rm \pi}^{2}m_{{\rm i}}L\int _{\mathcal{E}>\varphi_{\mathrm{eff}}}\,{\rm d}r, \end{equation}

and $\varphi _{\mathrm {eff}}$ is the effective potential (5.4) in the central plane.

The remaining integral in $\varGamma$ can be represented in a simpler form. To do this, we further assume that, first, there is no reversed field and, second, the magnetic field either increases with the radius or decreases not faster than $r^{-1}$. In other words, the following conditions are met:

(6.14a,b)\begin{equation} B\geq0,\quad\frac{r}{B}\frac{{\rm d}B}{{\rm d}r}\geq -1. \end{equation}

It can be shown that, in this case, the effective potential $\varphi _{\mathrm {eff}}$ has only one minimum. Consequently, the equation $\mathcal {E}=\varphi _{\mathrm {eff}}$ has no more than two roots: $r_{\min }=r_{\min }(\mathcal {E},\mathcal {P})$ and $r_{\max }=r_{\max }(\mathcal {E},\mathcal {P})$, which are the radial boundaries of the integration domain $\mathcal {E}>\varphi _{\mathrm {eff}}$. These roots are essentially the minimum and maximum distance from the axis available to a hot ion with energy $\mathcal {E}$ and angular momentum $\mathcal {P}$. In what follows, $r_{\min }$ and $r_{\max }$ are shortly referred to as the minimum radius and the maximum radius, respectively, and we also define the corresponding minimum and maximum radii for the injected ions: $\bar {r}_{\min }=r_{\min }(\mathcal {E}_{\mathrm {NB}},\mathcal {P}_{\mathrm {NB}})$ and $\bar {r}_{\max }=r_{\max }(\mathcal {E}_{\mathrm {NB}},\mathcal {P}_{\mathrm {NB}})$. Eventually, the phase volume can be represented in the following form:

(6.15)\begin{equation} \varGamma=4{\rm \pi}^{2}m_{{\rm i}}L\left(r_{\max}-r_{\min}\right). \end{equation}

It is worth noting that the minimum radius for the ions passing through the interior of the bubble core can be found explicitly. Indeed, since $\varPsi \equiv 0$ in the bubble core, for $r_{\min }< a$ we have

(6.16)\begin{equation} \mathcal{E}=\varphi_{\mathrm{eff}}\left(r_{\min}\right)=\frac{\left|\mathcal{P}\right|^{2}}{2m_{{\rm i}} r_{\min}^{2}},\quad\Rightarrow\quad r_{\min}=\frac{\left|\mathcal{P}\right|}{\sqrt{2m_{{\rm i}}\mathcal{E}}}. \end{equation}

Note that the solution (6.12a,b) also implies that, in the approximation considered, the hot ion energy and angular momentum are related as follows: $\mathcal {P}\sqrt {\mathcal {E}_{\mathrm {NB}}}=\mathcal {P}_{\mathrm {NB}}\sqrt {\mathcal {E}}$. This means that the minimum radius for the ions injected in the bubble core ($r_{\mathrm {NB}}< a$) remains constant

(6.17)\begin{equation} r_{\min}=\bar{r}_{\min}=\frac{\left|\mathcal{P}_{\mathrm{NB}}\right|}{\sqrt{2m_{{\rm i}}\mathcal{E}_{\mathrm{NB}}}}. \end{equation}

Given (6.9), the minimum radius can also be expressed in terms of the injection radius $r_{\mathrm {NB}}$ and injection angle $\xi _{\mathrm {NB}}$: $\bar {r}_{\min }=r_{\mathrm {NB}}\sin \xi _{\mathrm {NB}}$. In addition, it follows that the sign of the angular momentum does not change during the ion slowing down, and the ions initially injected into the absolute confinement region (2.2) always stay in it

(6.18)\begin{equation} {-}\mathcal{R}\varOmega_{{\rm i}}\mathcal{P}_{\mathrm{NB}}>\mathcal{E}_{\mathrm{NB}},\quad\Rightarrow \quad-\mathcal{R}\varOmega_{{\rm i}}\mathcal{P}>\sqrt{\mathcal{E}\mathcal{E}_{\mathrm{NB}}}\geq \mathcal{E}. \end{equation}

Therefore, the axial loss of the hot ions can indeed be neglected.¹⁰

Eventually, using the resulting distribution function (6.12a,b) and the definition (5.17), we obtain the azimuthal diamagnetic current density of the hot ions

(6.19)\begin{align} \frac{J_{{\rm h}\theta}}{2{\rm \pi} rc}& =-\frac{\bar{r}_{\min}}{a}\frac{\varPi_{{\rm h}}}{\psi_{{\rm h}}} \frac{a^{3}}{r^{3}}\varTheta\left(\frac{r-\bar{r}_{\min}}{a-\bar{r}_{\min}}- \frac{\psi}{\psi_{{\rm h}}}\right)\nonumber\\ & \quad \times\int_{(({a-\bar{r}_{\min}})/({r-\bar{r}_{\min}}))({\psi}/{\psi_{{\rm h}}})}^{1}\frac{a-\bar{r}_{\min}}{ r_{\max}^{*}\left(\eta\right)-\bar{r}_{\min}}\left(1+\frac{a-\bar{r}_{\min}}{\bar{r}_{\min}} \frac{\psi}{\psi_{{\rm h}}}\frac{1}{\eta}\right){\rm d}\eta, \end{align}

where $\psi _{{\rm h}}=2{\rm \pi} B_{{\rm v}}(a-\bar {r}_{\min })\rho _{\mathrm {NB}}$ is the characteristic magnetic flux induced by the hot ion current, $\rho _{\mathrm {NB}}=\sqrt {2m_{{\rm i}}\mathcal {E}_{\mathrm {NB}}}c/ZeB_{{\rm v}}$ is the characteristic Larmor radius of injected ions, the coefficient

(6.20)\begin{equation} \varPi_{{\rm h}}=\frac{Q_{{\rm h}0}}{\nu_{\mathrm{ie}0}{\rm \pi} a^{2}L}, \end{equation}

is the characteristic energy density of the hot ions and

(6.21)\begin{equation} r_{\max}^{*}\left(\eta\right)\overset{\mathrm{def}}=r_{\max}\left(\mathcal{E}_{\mathrm{NB}}\eta^{2},\mathcal{P}_{\mathrm{NB}}\eta\right), \end{equation}

is the maximum radius of the hot ions with the distribution (6.12a,b), which is to be found from the equation $\mathcal {E}=\varphi _{\mathrm {eff}}$ reduced to explicit form

(6.22)\begin{equation} r_{\max}^{*}\left(\eta\right)=\left.r\left(\psi\right)\right|_{\psi=\psi_{{\rm h}}\eta (({r_{\max}^{*}\left(\eta\right)-\bar{r}_{\min}})/({a-\bar{r}_{\min}}))}.\end{equation}

In particular, the outer boundary is $\bar {r}_{\max }=r_{\max }^{*}(1)$.

7 Thin transition layer limit

The resulting complete system of equilibrium equations for a cylindrical bubble (6.1), (6.2), (6.5) and (6.6) accompanied by the expression for the hot ion current (6.19) proves to be essentially nonlinear. An exact equilibrium solution can only be obtained numerically. Detailed numerical simulations and analysis of the corresponding numerical equilibria are planned to be provided in future work. In the present paper, we focus on considering a limiting case that allows a significant simplification of the equations.

As mentioned at the beginning of § 2, at the boundary of the diamagnetic bubble, right beyond the core, there is a transition layer inside which the magnetic field changes from $B=0$ to $B=B_{{\rm v}}$. The total thickness of this layer – the radial size of the region such that $0< B< B_{{\rm v}}$ – is further denoted by $\lambda$. The inner boundary of the layer $r=a$ corresponds to the radial boundary of the bubble core, where the magnetic field is zero: $B|_{r\leq a}=0$. The outer boundary $r=a+\lambda$ is the radius at which the diamagnetic current density of the plasma vanishes and the magnetic field reaches the vacuum value: $B|_{r\geq a+\lambda }=B_{{\rm v}}$. In other words, the outer boundary of the transition layer represents the boundary of the plasma. In the same way, we define separately the boundaries for the warm plasma $r=a+\lambda _{{\rm w}}$ and hot ion ions $r=a+\lambda _{{\rm h}}$, beyond which the corresponding diamagnetic currents vanish; $\lambda _{{\rm w}}$ and $\lambda _{{\rm h}}$ are further naturally called the thicknesses of the transition layer for the warm plasma and the hot ions, respectively. It is clear that the total thickness of the transition layer is determined by the largest one: $\lambda =\max \{ \lambda _{{\rm w}},\lambda _{{\rm h}}\}$.

The transition layer thickness for the warm plasma $\lambda _{{\rm w}}$ is determined by the characteristic scale of the warm plasma resistive transverse diffusion across the magnetic field, which is normally quite small. In particular, the MHD equilibrium model (Beklemishev 2016) shows that the thickness of the transition layer proves to be approximately equal to $\lambda _{{\rm w}\,\mathrm {MHD}}=7\lambda _{\mathrm {GD}}$, where

(7.1)\begin{equation} \lambda_{\mathrm{GD}}=\sqrt{\frac{c^{2}}{4{\rm \pi}\sigma _{{\rm w}}}\tau_{\mathrm{GD}}},\end{equation}

is the characteristic thickness of magnetic field diffusion into the plasma and $\tau _{\mathrm {GD}}=\mathcal {R}L/2u_{{\rm m}}$ is the gas-dynamic lifetime. In the case of ‘classical’ Spitzer conductivity $\sigma _{{\rm w}}=Z^{2}e^{2}n_{{\rm i}}/m_{{\rm i}}\nu _{\mathrm {ie}}$, for the typical parameters: $T\sim 1\ \mathrm {keV}$, $\mathcal {R}\sim 15$, $L\sim 500\ \mathrm {cm}$, the quantity $\lambda _{\mathrm {GD}}$ is of the order of tenths of a centimetre. However, in the presence of the hot ion component, MHD is not applicable, and this estimate ceases to be valid.

For the distribution (6.12a,b), the outer boundary, at which the current of the hot ions (6.19) vanishes, corresponds to the maximum radius of the injected ions $\bar {r}_{\max }$. Then the thickness of the hot ion transition layer in this case is equal to $\lambda _{{\rm h}}=\bar {r}_{\max }-a$. On the other hand, the hot ion transition layer is essentially the boundary region inside which the hot ions are ‘reflected’ from the external vacuum magnetic field outside the bubble core, and its thickness seems to be proportional to the Larmor radius of the injected ions $\lambda _{{\rm h}}\sim \rho _{\mathrm {NB}}$. Under typical conditions, the hot ion Larmor radius $\rho _{\mathrm {NB}}$ is of the order of centimetres and normally proves to be much greater than the warm plasma transverse diffusion scale $\lambda _{\mathrm {GD}}$. Thus, the total transition layer thickness $\lambda =\max \{ \lambda _{{\rm w}},\lambda _{{\rm h}}\}$ could be expected mainly determined by the hot ion transition layer thickness: $\lambda _{{\rm h}}\gg \lambda _{{\rm w}}$.

Further, in this work, we consider the limit of the thin transition layer $\lambda \ll a$, which allows the model equations to be greatly simplified. Since $\lambda$ normally decreases with increasing magnetic field, this approximation seems to correspond to the case of a sufficiently strong vacuum field $B_{{\rm v}}$.

7.1 Equilibrium magnetic field distribution

The magnetic field equilibrium is determined by (6.1) and (6.2) with the boundary conditions (6.3a,b) and (6.4). When solving these equations, we consider the pressure profile of the warm plasma $p=p(\psi )$ given, implying that it can be found from the corresponding equilibrium equations (6.5) and (6.6). At the same time, the diamagnetic current of the hot ions is determined by the expression (6.19).

It is further convenient to use the following dimensionless quantities:

(7.2a-d)\begin{equation} \phi=\frac{\psi}{\psi_{{\rm h}}},\quad x=\frac{r}{a},\quad R=\frac{B}{B_{{\rm v}}},\quad\beta_{{\rm w}}=\frac{8{\rm \pi} p}{B_{{\rm v}}^{2}}, \end{equation}

where $\beta _{{\rm w}}$ is the warm plasma pressure normalized to the energy density of the vacuum magnetic field $B_{{\rm v}}^{2}/8{\rm \pi}$, also referred to simply as the relative pressure of the warm plasma. According to this, we also define the dimensionless minimum and maximum radii as follows:

(7.3a,b)\begin{equation} \bar{x}_{\min}=\frac{\bar{r}_{\min}}{a},\quad x_{\max}^{*}\left(\eta\right)=\frac{r_{\max}^{*}\left(\eta\right)}{a}, \end{equation}

where the latter is determined by the equation

(7.4)\begin{equation} x_{\max}^{*}\left(\eta\right)=\left.x\left(\phi\right)\right|_{\phi=\eta (({x_{\max}^{*}\left(\eta\right)-\bar{x}_{\min}})/{(1-\bar{x}_{\min}}))},\end{equation}

which results from the corresponding normalization of (6.22). Then, (6.1) and (6.2) with substituted hot ion current (6.19) take the following dimensionless form:

(7.5)\begin{gather} \frac{{\rm d}}{{\rm d}\phi}\left(R^{2}+\beta_{{\rm w}}\right)\nonumber\\ =\frac{\varLambda_{{\rm h}}}{x^{3}}\varTheta\left(\frac{x-\bar{x}_{\min}}{1-\bar{x}_{\min}}-\phi\right) \int_{(({1-\bar{x}_{\min}})/({x-\bar{x}_{\min}}))\phi}^{1} \frac{1-\bar{x}_{\min}}{x_{\max}^{*}\left(\eta\right)-\bar{x}_{\min}}\left(1+ \frac{1-\bar{x}_{\min}}{\bar{x}_{\min}}\frac{\phi}{\eta}\right){\rm d}\eta, \end{gather}

(7.6)\begin{gather} \frac{{\rm d}x}{{\rm d}\phi}=\frac{\epsilon}{x}\frac{1}{R},\end{gather}

where

(7.7a,b)\begin{equation} \varLambda_{{\rm h}}=\bar{x}_{\min}\frac{8{\rm \pi}\varPi_{{\rm h}}}{B_{{\rm v}}^{2}},\quad \epsilon=\left(1-\bar{x}_{\min}\right)\frac{\rho_{\mathrm{NB}}}{a}. \end{equation}

Normalizing the boundary conditions (6.3a,b) and (6.4) yields

(7.8a,b)\begin{gather} \left.R\right|_{\phi=0}=0,\quad\left.x\right|_{\phi=0}=1, \end{gather}

(7.9)\begin{gather}\left.R\right|_{\phi\rightarrow+\infty}=1. \end{gather}

As mentioned above, the hot ion transition layer thickness $\lambda _{{\rm h}}$ appears to be of the order of the Larmor radius of the injected ions $\rho _{\mathrm {NB}}$, which means that, in the thin transition layer limit before us $\lambda _{{\rm h}}\ll a$, we should also expect the ratio $\rho _{\mathrm {NB}}/a$ to be small. Given this, it seems appropriate to apply the method of successive approximations by considering the coefficient $\epsilon \propto \rho _{\mathrm {NB}}/a$ on the right-hand side of (7.6) as a small expansion parameter. Namely, we further assume that the solution of the system (7.5) and (7.6) can be represented in the form of an asymptotic expansion.

For a given magnetic field profile $R=R(\phi )$, by integrating (7.6) and taking into account the boundary condition (7.8a,b), the radius as a function of the magnetic flux can be explicitly expressed

(7.10)\begin{equation} x\left(\phi\right)=\sqrt{1+2\epsilon\int_{0}^{\phi}\frac{{\rm d}\phi^{\prime}}{R\left(\phi^{\prime}\right)}}.\end{equation}

At the leading order of the approximation $\epsilon \ll 1$, when small corrections are completely neglected, the expression (7.10) reduces to

(7.11)\begin{equation} x\left(\phi\right)\simeq1+o\left(\epsilon^{0}\right).\end{equation}

Substituting (7.11) into (7.4) yields the corresponding approximation for the maximum radius

(7.12)\begin{equation} x_{\max}^{*}\left(\eta\right)\simeq1+o\left(\epsilon^{0}\right).\end{equation}

The obtained formulas (7.11) and (7.4) correspond to $r\simeq a$. In other words, on the scale of the transition layer, the radius $r$ can be considered almost constant and equal to the bubble radius $a$, while the magnetic flux $\psi$ varies greatly. At a qualitative level, this is associated with the high density of the diamagnetic current inside the layer.

Taking into account (7.11) and (7.12), the integral in the right-hand side of (7.5) is evaluated explicitly

(7.13)\begin{align} \frac{{\rm d}}{{\rm d}\phi}\left(R^{2}+\beta_{{\rm w}}\right)& \simeq \varLambda_{{\rm h}}\varTheta\left(1-\phi\right)\int _{\phi}^{1}\left[1+\left(\bar{x}_{\min}^{-1}-1\right) \frac{\phi}{\eta}\right]{\rm d}\eta+o\left(\epsilon^{0}\right)\nonumber\\ & =\varLambda_{{\rm h}}\left[1-\phi-\left(\bar{x}_{\min}^{-1}-1\right)\phi\ln\phi\right] \varTheta\left(1-\phi\right)+o\left(\epsilon^{0}\right). \end{align}

As can be seen, the outer boundary of the hot ions at this order of approximation corresponds to $\phi =1$. Given the boundary conditions (7.8a,b) and (7.9), the resulting equation yields

(7.14)\begin{equation} R\left(\phi\right)\simeq R^{\left(0\right)}\left(\phi\right)+o\left(\epsilon^{0}\right),\end{equation}

where

(7.15)\begin{align} R^{\left(0\right)}\left(\phi\right)=\begin{cases} \sqrt{\beta_{{\rm w}0}-\beta_{{\rm w}}\left(\phi\right)+\varLambda_{{\rm h}}\phi\left[1- \dfrac{\phi}{2}-\left(\bar{x}_{\min}^{-1}-1\right)\dfrac{\phi}{2}\left(\ln\phi-\dfrac{1}{2}\right)\right]}, & \phi<1,\\[10pt] {}\sqrt{1-\beta_{{\rm w}}\left(\phi\right)}, & \phi\geq1, \end{cases} \end{align}

and $\beta _{{\rm w}0}=\beta _{{\rm w}}(0)=8{\rm \pi} p_{0}/B_{{\rm v}}^{2}$ is the relative pressure of the warm plasma inside the bubble core. In the absence of surface currents, the solution (7.14) should be continuous at the boundary $\phi =1$, from which we also obtain the following condition:

(7.16)\begin{equation} 1\simeq\beta_{{\rm w}0}+\beta_{{\rm h}0}+o\left(\epsilon^{0}\right).\end{equation}

Here, the quantity

(7.17)\begin{equation} \beta_{{\rm h}0}=\frac{1}{4}\left(\bar{x}_{\min}^{-1}+1\right)\varLambda_{{\rm h}}= \left(1+\bar{x}_{\min}\right)\frac{2{\rm \pi}\varPi_{{\rm h}}}{B_{{\rm v}}^{2}}, \end{equation}

can be interpreted as the characteristic relative energy density of the hot ions in the bubble. The formula (7.16) essentially expresses the balance between the plasma pressure inside the bubble and the pressure of the external magnetic field.

By applying the successive approximations to find higher orders of the expansion, the solution can be further refined. In particular, taking into account first-order corrections in the expression (7.10) yields

(7.18)\begin{equation} x\left(\phi\right)\simeq1+\epsilon\int_{0}^{\phi}\frac{{\rm d}\phi^{\prime}}{R^{\left(0\right)}\left(\phi^{\prime}\right)}+o\left(\epsilon^{1}\right). \end{equation}

Using the obtained dependence, the corresponding correction to the maximum radius is also found from (7.4)

(7.19)\begin{equation} x_{\max}^{*}\left(\eta\right)\simeq1+\epsilon\int_{0}^{\eta}\frac{{\rm d}\phi}{R^{\left(0\right)} \left(\phi\right)}+o\left(\epsilon^{1}\right). \end{equation}

For $\eta =1$ this formula determines the thickness of the hot ion transition layer

(7.20)\begin{equation} \frac{\lambda_{{\rm h}}}{a}=x_{\max}^{*}\left(1\right)-1\simeq\epsilon\int_{0}^{1}\frac{{\rm d}\phi}{R^{ (0)}\left(\phi\right)}+o\left(\epsilon^{1}\right).\end{equation}

Further, in this paper, we consider the equilibrium configuration of the magnetic field to be approximately defined by (7.11) and (7.14) at $\epsilon \rightarrow 0$. In other words, we take into account only the leading order of the asymptotic expansion. The applicability criterion for this approximation is determined by the limit

(7.21)\begin{equation} \frac{1-x_{\max}^{*}\left(1\right)}{1-\bar{x}_{\min}}=\frac{\lambda_{{\rm h}}}{a-\bar{r}_{\min}}\ll1, \end{equation}

which is actually assumed in the derivation of (7.14). After substituting (7.20), we arrive at the following condition:

(7.22)\begin{equation} \int _{0}^{1}\frac{{\rm d}\phi}{R^{\left(0\right)}\left(\phi\right)}\ll\frac{a}{\rho_{\mathrm{NB}}}. \end{equation}

An upper bound for the integral involved can be obtained as follows:

(7.23)\begin{equation} \int_{0}^{1}\frac{{\rm d}\phi}{R^{\left(0\right)}\left(\phi\right)}\leq \frac{1}{\sqrt{1-\beta_{{\rm w}0}}} \mathcal{W}\left(\frac{a}{\bar{r}_{\min}}\right), \end{equation}

where

(7.24)\begin{equation} \mathcal{W}\left(z\right)\overset{\mathrm{def}}=\frac{\sqrt{z+1}}{2}\int_{0}^{1} \dfrac{{\rm d}\phi}{\sqrt{\phi\left[1-\dfrac{\phi}{2}-\left(z-1\right) \dfrac{\phi}{2}\left(\ln\phi-\dfrac{1}{2}\right)\right]}}.\end{equation}

As can be seen from figure 2, for all reasonable $a/\bar {r}_{\min }$ the function $\mathcal {W}$ is of the order of $2\div 3$.

Figure 2. Function $\mathcal {W}=\mathcal {W}(z)$ defined by the expression (7.24).

Before moving on to considering the equilibrium of the warm plasma, it is also useful to examine the expression (7.16) in greater detail. For given parameters of the warm plasma inside the bubble core: density $n_{{\rm i}0}$ and temperature $T_{0}$, which are to be found from the solution of the equilibrium equations (6.5) and (6.6), the expression (7.16) defines the relation between the radius of the core $a$ and the fixed external parameters: absorbed heating power $Q_{{\rm h}0}$, vacuum magnetic field $B_{{\rm v}}$ and geometric parameters $L$ and $\bar {r}_{\min }$. After expressing the warm plasma density $n_{{\rm i}0}$ from the equation of state (4.2), the formula (7.16) is reduced to a cubic equation for the normalized core radius $k=a/\bar {r}_{\min }=\bar {x}_{\min }^{-1}$

(7.25)\begin{equation} k^{3}-\left(4\tau_{{\rm s}}\nu_{\mathrm{ie}0}\right)^{-1}\left(k+1\right)\simeq0,\end{equation}

where

(7.26)\begin{equation} \tau_{{\rm s}}=\beta_{{\rm h}0}\frac{B_{{\rm v}}^{2}}{8{\rm \pi}}\frac{{\rm \pi}\bar{r}_{\min}^{2}L}{Q_{{\rm h}0}}. \end{equation}

Provided that the core radius $a$ exceeds the minimum radius $\bar {r}_{\min }$, i.e. $k\geq 1$, this equation has exactly one real root for a fixed $\beta _{{\rm w}0}\in [0,1]$

(7.27)\begin{equation} k = \begin{cases} \dfrac{3}{\mathfrak{D}}\cos\left(\dfrac{1}{3}\arccos\mathfrak{D}\right), & 0<\mathfrak{D}\leq1,\\[10pt] \dfrac{3}{\mathfrak{D}}\cosh\left(\dfrac{1}{3}\mathrm{arccosh}\,\mathfrak{D}\right), & 1<\mathfrak{D}, \end{cases} \end{equation}

where $\mathfrak {D}=\sqrt {27\tau _{{\rm s}}\nu _{\mathrm {ie}0}}$ and $\mathrm {arccosh}\,x=\ln (x+\sqrt {x^{2}-1})$. This means that, for a given temperature in the core $T_{0}$, this root specifies the functional dependency of the core radius on the pressure in the core: $k=k(\beta _{{\rm w}0})$.

7.2 Warm plasma equilibrium

The equilibrium of the warm plasma is described by (6.5) and (6.6) with the boundary conditions (4.24), (4.26) and (4.27a,b). In what follows, the electrical and thermal conductivity of the warm plasma are considered classical (Braginskii 1965)

(7.28)\begin{equation} \left.\begin{gathered} \sigma _{{\rm w}}=\frac{Z^{2}e^{2}n_{{\rm i}}}{m_{{\rm i}}\nu_{\mathrm{ie}}}=8.7\times10^{13}Z^{-1}\varLambda_{\mathrm{ie}}^{-1} \left(\frac{T}{\mathrm{eV}}\right)^{{3}/{2}}\ \mathrm{s}^{-1},\\ \varkappa_{{\rm w}}=10^{-3}\left(\mu_{{\rm i}}^{{1}/{2}}+0.077Z\right)Z^{2}\varLambda_{\mathrm{ie}} \left(\frac{B}{\mathrm{G}}\right)^{-2}\left(\frac{n_{{\rm i}}}{{\rm cm}^{-3}}\right)^{2}\left(\frac{T}{{\rm eV}}\right)^{-{1}/{2}}{\rm cm}^{-1}\ {\rm s}^{-1}. \end{gathered}\right\} \end{equation}

In addition, the warm plasma flow velocity in the mirror is assumed to be

(7.29)\begin{equation} u_{{\rm m}}=\sqrt{\frac{2}{\rm \pi}}v_{T{\rm i}}=\sqrt{\frac{2T}{{\rm \pi} m_{{\rm i}}}}, \end{equation}

which corresponds to the gas-dynamic loss in the case of short mirrors and a filled loss cone (Ivanov & Prikhodko 2017). In what follows, we also consider $r\simeq a$, assuming the leading order of the thin transition layer approximation $\lambda _{{\rm w}}\ll a$.

Let us show that in this case there exists an equilibrium of the warm plasma such that the temperature can be considered slowly varying on the scale of the warm plasma inhomogeneity

(7.30)\begin{equation} \frac{n_{{\rm i}}}{T}\frac{{\rm d}T}{{\rm d}n_{{\rm i}}}\simeq\frac{p}{T}\frac{{\rm d}T}{{\rm d}p}=\gamma_{T}\ll1.\end{equation}

Subtracting (6.5) multiplied by $\alpha _{{\rm E}}T$ from (6.6) yields

(7.31)\begin{equation} \alpha_{{\rm E}}T{\rm d}\left(\varLambda_{{\rm i}\perp}\frac{{\rm d}p}{{\rm d}\psi}\right)- {\rm d}\left(\frac{5\left(1+Z\right)}{2}T\varLambda_{{\rm i}\perp}\frac{{\rm d}p}{{\rm d}\psi}\right)- {\rm d}\left(\varLambda_{{\rm E}\perp}\frac{{\rm d}T}{{\rm d}\psi}\right)=0. \end{equation}

Taking into account the condition (7.30), the left-hand side of the resulting expression is approximately reduced to an exact differential

(7.32)\begin{equation} {\rm d}\left(\left(\alpha_{{\rm E}}-\frac{5\left(1+Z\right)}{2}\right)T\varLambda_{{\rm i}\perp} \frac{{\rm d}p}{{\rm d}\psi}\right)-{\rm d}\left(\varLambda_{{\rm E}\perp}\frac{{\rm d}T}{{\rm d}\psi}\right)\simeq0, \end{equation}

which further allows the equation to be explicitly integrated. The constant of integration should be set equal to zero, since at the boundary of the warm plasma, where $n_{{\rm i}}\rightarrow 0$, the transverse fluxes of mass and energy should vanish. Finally, we find the slope factor

(7.33)\begin{align} \gamma_{T}=\frac{p}{T}\frac{{\rm d}T}{{\rm d}p}\simeq\left(\alpha_{{\rm E}}- \frac{5\left(1+Z\right)}{2}\right)\frac{\varLambda_{{\rm i}\perp}}{\varLambda_{{\rm E}\perp}} p\simeq\frac{1}{60}\left(\alpha_{{\rm E}}- \frac{5\left(1+Z\right)}{2}\right)\frac{1+Z}{\left(\mu_{{\rm i}}^{{1}/{2}}+0.077Z\right)Z}.\end{align}

For hydrogen plasma $Z=1$, $\mu _{{\rm i}}=1$ and $\alpha _{{\rm E}}=8$, which is typical for GDT (Ivanov & Prikhodko 2017; Skovorodin 2019; Soldatkina et al. 2020), the slope factor proves to be quite small: $\gamma _{T}\sim 10^{-1}\ll 1$.

When the condition (7.30) is met, the equilibrium equations (6.5) and (6.6) are equivalent and reduced to

(7.34a,b)\begin{equation} \frac{{\rm d}\mathcal{F}}{{\rm d}\chi}+2\beta_{{\rm w}}\simeq0,\quad \mathcal{F}=-\frac{\beta_{{\rm w}}}{R}\frac{{\rm d}\beta_{{\rm w}}}{{\rm d}\chi},\end{equation}

where we use the dimensionless magnetic flux $\chi =\psi /\psi _{\mathrm {GD}}$ normalized to $\psi _{\mathrm {GD}}=2{\rm \pi} a\lambda _{\mathrm {GD}}B_{{\rm v}}$, and $\lambda _{\mathrm {GD}}$ is defined by (7.1). We also introduce the quantity $\mathcal {F}$, which has the meaning of the warm plasma transverse flow. The magnetic field distribution is given by (7.14) taking into account the corresponding renormalization of the magnetic flux: $\phi =\chi /\chi _{{\rm h}}$, where $\chi _{{\rm h}}=\psi _{{\rm h}}/\psi _{\mathrm {GD}}=(1-k^{-1})\rho _{\mathrm {NB}}/\lambda _{\mathrm {GD}}$ corresponds to the boundary of the hot ions $\phi =1$. When solving (7.34a,b), according to the approximation (7.30), the temperature should be considered constant and equal to $T_{0}$. Further, a weak temperature dependence can be found from the equality (7.33)

(7.35)\begin{equation} \frac{T}{T_{0}}=\left(\frac{p}{p_{0}}\right)^{\gamma_{T}}.\end{equation}

Internal boundary conditions (4.24) and (4.26) also merge into one

(7.36)\begin{gather} \left.\mathcal{F}\right|_{\chi\rightarrow0}\simeq\mathcal{S}, \end{gather}

(7.37)\begin{gather}\mathcal{S}=\beta_{{\rm w}0}\frac{3}{2}\sqrt{\frac{\rm \pi}{2}}\frac{\rho_{{\rm i}0}}{ \lambda_{\mathrm{GD}}}\left(\frac{W_{{\rm i}0}}{2\varPhi_{{\rm i}\parallel0}}-1\right), \end{gather}

with an additional condition

(7.38)\begin{equation} \alpha_{{\rm E}}T_{0}\left(W_{{\rm i}0}-2\varPhi_{{\rm i}\parallel0}\right)=Q_{{\rm h}0}- 2\alpha_{{\rm E}0}T_{0}\varPhi_{{\rm i}\parallel0},\end{equation}

which defines the relationship between the temperature and the pressure of the warm plasma in the bubble core. For simplicity, we further assume $\alpha _{{\rm E}0}=\alpha _{{\rm E}}$, then the temperature of the warm plasma $T_{0}$ can be explicitly found from the equality (7.38). It proves to be independent of the pressure and is determined by the ratio of the sources

(7.39)\begin{equation} T_{0}\simeq\frac{Q_{{\rm h}0}}{\alpha_{{\rm E}}W_{{\rm i}0}}.\end{equation}

External boundary conditions (4.27a,b) are reduced to the following:

(7.40)\begin{equation} \left.\beta_{{\rm w}}\right|_{\chi=\chi_{{\rm w}}}=0,\quad\Leftrightarrow\quad\left. \mathcal{F}\right|_{\chi=\chi_{{\rm w}}}=0,\end{equation}

where $\chi =\chi _{{\rm w}}$ corresponds to the outer boundary of the warm plasma $r=a+\lambda _{{\rm w}}$. The second condition is essentially the vanishing of the transverse flow at the boundary, which is consistent with the pressure gradient and the current density being finite. It is also worth noting that the position of the boundary $\chi _{{\rm w}}$ is not specified and should be found self-consistently from the solution of the equilibrium equation (7.34a,b). This means that the two conditions (7.40) are related by $\chi _{{\rm w}}$ and are formally reducible to one. In other words, one of the expressions (7.40) can be considered as a boundary condition, and the other as a definition of the boundary position $\chi =\chi _{{\rm w}}$.

7.3 Numerical solution

The boundary value problem (7.34a,b), (7.36) and (7.38) described above turns out to be rather complex, and to find its exact solution we use numerical methods.

As mentioned above, provided that the temperature is given by the expression (7.39), (7.25) relates the radius of the bubble core and the pressure of warm plasma in the core: $k=k(\beta _{{\rm w}0})$. This means that $\beta _{{\rm w}0}$ remains the only free (i.e. unknown yet) parameter that determines the equilibrium. In that regard, the boundary value problem before us is convenient to reformulate as the following Cauchy problem. The equilibrium equations (7.34a,b) can be considered as the system of ordinary differential equations for the functions $\beta _{{\rm w}}=\beta _{{\rm w}}(\chi )$ and $\mathcal {F}=\mathcal {F}(\chi )$ with the initial conditions (7.36) and $\beta _{{\rm w}}|_{\chi =0}=\beta _{{\rm w}0}$. At the same time, the quantity $\beta _{{\rm w}0}$ should be treated as a variable parameter, which is found from the external boundary condition (7.40).

To solve the formulated problem, we apply an approach similar to the shooting method. For a given parameter $\beta _{{\rm w}0}$, the equations (7.34a,b) are integrated by means of the Runge–Kutta methods, starting from the corresponding initial conditions at $\chi =0$. Computing continues until either $\beta _{{\rm w}}=\beta _{{\rm w}}(\chi ;\beta _{{\rm w}0})$ or $\mathcal {F}=\mathcal {F}(\chi ;\beta _{{\rm w}0})$ hits zero at some point $\chi =\chi _{{\rm w}}(\beta _{{\rm w}0})$.¹¹ Thus, by varying the parameter $\beta _{{\rm w}0}$, we can formally obtain the functional dependence $\chi _{{\rm w}}=\chi _{{\rm w}}(\beta _{{\rm w}0})$. On the other hand, for a given relation $\chi _{{\rm w}}=\chi _{{\rm w}}(\beta _{{\rm w}0})$, the conditions (7.40) are actually reduced to an equation for $\beta _{{\rm w}0}$ on the finite interval $0\leq \beta _{{\rm w}0}\leq 1$. The root of this equation, which can be found using standard numerical methods, is a true value of the parameter $\beta _{{\rm w}0}$ corresponding to the real equilibrium profile $\beta _{{\rm w}}=\beta _{{\rm w}}(\chi )$.

Figure 3 shows the numerical solution of the warm plasma equilibrium equation, found using the procedure described above. There are also plotted the corresponding analytical asymptotics in the vicinity of the bubble core and the outer boundary of the warm plasma, respectively,

(7.41)\begin{align} \beta_{{\rm w}}& \sim\beta_{{\rm w}0}-\frac{4\mathcal{S}\mathfrak{C}_{{\rm h}}}{3\beta_{{\rm w}0}}\chi^{{3}/{2}}-\frac{\mathcal{S}^{2}}{6\beta_{{\rm w}0}^{2}}\chi^{2}\nonumber\\ & \quad +\frac{\mathfrak{C}_{{\rm h}}}{5}\left\{ 8-\frac{\mathcal{S}^{3}}{36\mathfrak{C}_{{\rm h}}^{2}\beta_{{\rm w}0}^{3}}+\frac{\mathcal{S}}{\beta_{{\rm w}0}\chi_{{\rm h}}} \left[1+\left(k-1\right)\left(\ln\left(\frac{\chi}{\chi_{{\rm h}}}\right)-\frac{9}{10}\right)\right]\right\} \chi^{{5}/{2}}\nonumber\\ & \quad +o\left(\chi^{{5}/{2}}\right),\quad\chi\rightarrow0, \end{align}

(7.42a,b)\begin{align} \beta_{{\rm w}}& \sim\tfrac{1}{3}R\left(\chi_{{\rm w}}\right)\left(\chi_{{\rm w}}-\chi\right)^{2}+o \left[\left(\chi_{{\rm w}}-\chi\right)^{2}\right],\quad\chi\rightarrow\chi_{{\rm w}},\end{align}

where

(7.43)\begin{equation} \mathfrak{C}_{{\rm h}}=\sqrt{\frac{1}{\chi_{{\rm h}}}\frac{\beta_{{\rm h}0}}{k+1}}. \end{equation}

The second term in the expansion (7.41) includes the contribution from the hot ions – it corresponds to the pressure profile of the warm plasma in the magnetic field induced by the hot ion current only. The third term, in turn, takes into account the diamagnetic current of the warm plasma. The non-polynomial logarithmic contribution of the hot ion current appears in the fourth term of the expansion.

Figure 3. Profile of the warm plasma pressure depending on the magnetic flux: the numerical solution (black solid curve), left (blue dashed curve) and right (red dash-dotted curve) asymptotics given by the expressions (7.41) and (7.42a,b), respectively. Simulation corresponds to the parameters of the GDMT10 regime (see table 1).

Table 1. Parameters of the numerical simulations. The thicknesses of transition layers $\lambda _{{\rm h}}$ and $\lambda _{{\rm w}}$ are the widths of the current profiles shown in figure 5. The warm plasma thermodynamic parameters: relative pressure $\beta _{{\rm w}0}$, density $n_{{\rm i}0}$ and temperature $T_{0}$ are the maxima of the corresponding quantities shown in figure 4. The parameter $\eta _{\parallel 0}$ is found from numerical simulations.

$^{a}$Thickness of the warm plasma transition layer according to MHD equilibrium (Beklemishev 2016): $\lambda _{{\rm w}\,\mathrm {MHD}}=7\lambda _{\mathrm {GD}}$. See the beginning of § 7.

$^{b}$Proportion of the non-adiabatic loss (4.23) in the total axial loss of the warm plasma: $\eta _{\parallel 0}=2\varPhi _{{\rm i}\parallel 0}/W_{{\rm i}0}$.

Figure 4. Radial profiles of the warm plasma (a) relative pressure, (b) density and (c) temperature outside the bubble core $r\ge a$ for the regimes GDMT05 (blue dashed curves), GDMT10 (black solid curves) and GDMT20 (red dash-dotted curves) specified in table 1.

8 Diamagnetic confinement in GDMT

As an application of the constructed theoretical model, we further investigate the equilibrium of the diamagnetic bubble corresponding to the design parameters of the GDMT device (Skovorodin et al. 2023). Hydrogen plasma is assumed: $Z=\mu _{{\rm i}}=1$; the axial loss constants correspond to GDT (Ivanov & Prikhodko 2017; Skovorodin 2019; Soldatkina et al. 2020): $\alpha _{{\rm E}0}=\alpha _{{\rm E}}=8$; the Coulomb logarithm is set equal to $\varLambda _{\mathrm {ie}}=15$. Energy, angle and impact parameter of the injection are $\mathcal {E}_{\mathrm {NB}}=30\ \mathrm {keV}$, $\xi _{\mathrm {NB}}={\rm \pi} /4$ and $r_{\mathrm {NB}}=5\ \mathrm {cm}$, respectively, which corresponds to the minimum radius of the hot ions equal to $\bar {r}_{\min }\simeq 3.54\ \mathrm {cm}$. The bubble length is estimated from the MHD equilibrium simulations of GDMT (Khristo & Beklemishev 2022) and is equal to $L=500\ \mathrm {cm}$. The total warm plasma source intensity is fixed at $W_{{\rm i}0}=5\times 10^{21}\ \mathrm {s}^{-1}$. The magnetic field in the mirrors is $B_{{\rm m}}=200\ \mathrm {kG}$.

We consider the conventional case of the vacuum magnetic field equal to $B_{{\rm v}}=10\ \mathrm {kG}$, as well as the regimes with halved and doubled fields: $B_{{\rm v}}=5\ \mathrm {kG}$ and $B_{{\rm v}}=20\ \mathrm {kG}$. In all the simulations, the radius of the bubble core is fixed at $a=20\ \mathrm {cm}$, which is the characteristic expected plasma radius in the GDMT. Then the required total absorbed injection powers prove to be $Q_{{\rm h}0}\simeq 6.11\ \mathrm {MW}$, $9.75\ \mathrm {MW}$ and $14.93\ \mathrm {MW}$ for $B_{{\rm v}}=5$, $10$ and $20\ \mathrm {kG}$, respectively. For convenience, the simulation parameters are listed in table 1, where the considered regimes are briefly called ‘GDMT05’, ‘GDMT10’ and ‘GDMT20’.

The results of the simulations are shown in figures 4 and 5. Figure 4 shows the radial profiles of the warm plasma relative pressure, density and temperature outside¹² the bubble core $r\geq a$. The radial profiles of the magnetic field, along with the current densities of the warm plasma and the hot ions, outside the core $r\geq a$, are presented in figure 5.

Figure 5. Equilibrium radial profiles of the magnetic field (black solid curves), the current densities of the warm plasma (blue dashed curves) and the hot ions (red dash-dotted curves) outside the bubble core $r\ge a$. Simulation parameters are presented in table 1: (a) GDMT05, (b) GDMT10 and (c) GDMT20.

8.1 Analysis of the solutions

1. (i) For all solutions the relative energy density of the hot ions is significantly greater than the warm plasma relative pressure (see table 1 and figure 4)

   (8.1)\begin{equation} \beta_{{\rm w}0}\ll\beta_{{\rm h}0}\sim1.\end{equation}
   Thus, it turns out that the energy content of the plasma as well as the diamagnetic current almost entirely correspond to the hot ions, while the contribution of the warm plasma is negligible. Since the warm plasma density in this case should be relatively small, the drag force (5.7) acting on the injected ions appears to be weak as well. Therefore, this results in the energy being accumulated by the hot component rather than being transferred to the warm plasma, which is actually consistent with (8.1).

2. (ii) In table 1, listed are the simulated values of the transition layer thickness for the warm plasma $\lambda _{{\rm w}}$ and the hot ions $\lambda _{{\rm h}}$, determined from the corresponding current profiles plotted in figure 5. As expected, the total transition layer thickness $\lambda =\max \{ \lambda _{{\rm w}},\lambda _{{\rm h}}\}$ is manly determined by the hot component, since the transition layer for the warm plasma $\lambda _{{\rm w}}$ proves to be much thinner than that for the hot ions $\lambda _{{\rm h}}$. The thickness of the transition layer for the hot ions naturally turns out to be of the order of the Larmor radius of the injected ions in the vacuum magnetic field: $\lambda _{{\rm h}}\sim \rho _{\mathrm {NB}}$. At the same time, the warm plasma transition layer thickness proves to be surprisingly close to the estimate made for MHD equilibrium (Beklemishev 2016): $\lambda _{{\rm w}\,\mathrm {MHD}}=7\lambda _{\mathrm {GD}}$. This effect is probably due to a combination of the following two factors. On the one hand, the warm plasma transition layer should apparently be wider in the presence of injected ions, since the characteristic scale of the warm plasma transverse diffusion proves to be greater in the magnetic field weakened by the diamagnetism of the hot component. On the other hand, for a fixed warm plasma source, the lower maximum relative pressure of the warm plasma (8.1) seems to correspond to a radial pressure drop in a more narrow region. A more clear explanation of this phenomenon and a proper quantitative estimate of the thickness $\lambda _{{\rm w}}$ should be made based on analysis the warm plasma equilibrium equations, which is planned to be addressed in future work.

3. (iii) The equilibria presented in figures 4 and 5 are obtained within the thin transition layer limit $\lambda \ll a-\bar {r}_{\min }$ (see § 7). For all the solutions found, the expansion parameter $\lambda /(a-\bar {r}_{\min })$ proves to be less than unity. However, in the regime GDMT05 with relatively weak vacuum magnetic field (figure 5a), the expansion parameter is not too small: $\lambda /(a-\bar {r}_{\min })\simeq 7/20$. A proper description of such equilibria should be obtained by accurate numerical modelling that takes into account the finite thickness of the transition layer.

4. (iv) To simplify the system of the warm plasma equilibrium equations presented in § 7.2, we assume the temperature of the warm plasma to be approximately constant: $T\simeq \mathrm {const}$. As a result, the temperature profiles determined by the formula (7.35) actually prove to be relatively flat, but in addition, the corresponding pressure and density profiles turn out to be identically shaped (see figure 4). It is worth noting, however, that the approximation of a slowly varying temperature $T\simeq \mathrm {const}.$ is valid in the particular case of the hydrogen plasma with the energy loss factor $\alpha _{{\rm E}}$ corresponding to the GDT (Ivanov & Prikhodko 2017; Skovorodin 2019; Soldatkina et al. 2020), i.e. $Z=\mu _{{\rm i}}=1$ and $\alpha _{{\rm E}}=8$. At the same time, the nature of energy loss in the diamagnetic confinement mode may differ significantly from the gas-dynamic regime. In particular, the presence of the non-adiabatic loss (Chernoshtanov 2020, 2022) could greatly affect both longitudinal and transverse equilibrium.

5. (v) In § 6.3, we assume that the injected hot ions release the energy mainly inside the bubble core. As already mentioned, the warm plasma temperature radial profile turns out to be almost flat, while the density of the warm plasma $n_{{\rm i}}$ sharply drops just beyond the bubble core on the scale of $\lambda _{{\rm w}}$ (figure 4). In the thin transition layer approximation $\lambda _{{\rm w}}\leq \lambda \ll a-\bar {r}_{\min }$, this results in the ion–electron drag $\nu _{\mathrm {ie}}\propto n_{{\rm i}}T^{-{3}/{2}}$ indeed acting mainly inside the core. However, this appears to take place only when the source of the warm plasma is entirely contained inside the bubble core. Otherwise, the maximum density of the warm plasma could be located outside the core, which, apparently, could significantly increase the energy loss.

6. (vi) The simulations also yield the parameter $\eta _{\parallel 0}=2\varPhi _{{\rm i}\parallel 0}/W_{{\rm i}0}$ (see table 1), which represents the total proportion of the non-adiabatic loss (4.23) in the total warm plasma loss. It can be seen that the warm plasma loss turns out to be almost entirely non-adiabatic for all the considered regimes

   (8.2)\begin{equation} W_{{\rm i}0}\sim2\varPhi_{{\rm i}\parallel0},\end{equation}
   while the loss at the periphery of the bubble core $2\varPhi _{{\rm i}\parallel \mathrm {out}}=W_{{\rm i}0}-2\varPhi _{{\rm i}\parallel 0}$ proves to be negligible: $2\varPhi _{{\rm i}\parallel \mathrm {out}}\ll W_{{\rm i}0}$. One can also observe that the parameter $\eta _{\parallel 0}$ grows with increasing vacuum magnetic field. The reason for this is probably related to the high energy density of hot ions (8.1). The strong diamagnetism induced by the hot ions extending beyond the bubble core leads to a considerable weakening of the magnetic field $B$ in the vicinity of the bubble core $r\simeq a$. Therefore, the effective mirror ratio $\mathfrak {R}=B_{{\rm m}}/B$ inside the thin transition layer $a< r< a+\lambda _{{\rm w}}$ is greatly increased compared with the vacuum value $\mathcal {R}=B_{{\rm m}}/B_{{\rm v}}$. This should eventually result in the axial loss at the periphery $\varPhi _{{\rm i}\parallel \mathrm {out}}\sim n_{{\rm i}}u_{{\rm m}}S_{{\rm m}}$ being suppressed due to a decrease in the cross-section of the warm plasma flow in the mirrors $S_{{\rm m}}\sim a\lambda _{{\rm w}}/\mathfrak {R}$. At the same time, the non-adiabatic loss $\varPhi _{{\rm i}\parallel 0}\sim n_{{\rm i}0}u_{{\rm m}0}a\rho _{{\rm i}0}/\mathcal {R}$, being determined by the vacuum mirror ratio, $\mathcal {R}$, is not affected by the diamagnetism of the hot ions.

7. (vii) The warm plasma density $n_{{\rm i}0}$ turns out to be considerably lower in the regimes with a higher temperature of the warm plasma $T_{0}$ (see figure 4 and table 1). This appears to be due to the improved confinement resulting from the suppression of the axial loss in the regimes with a lower warm plasma temperature. When the non-adiabatic loss dominates, the warm plasma density can be approximately estimated from (8.2)

   (8.3)\begin{equation} n_{{\rm i}0}\sim\frac{W_{{\rm i}0}\tau_{{\rm i}\parallel0}}{{\rm \pi} a^{2}L}\propto\frac{W_{{\rm i}0}B_{{\rm m}}}{a}\frac{1}{T_{0}},\end{equation}
   where
   (8.4)\begin{equation} \tau_{{\rm i}\parallel0}\sim\frac{a}{\rho_{{\rm i}0}}\tau_{\mathrm{GD}}\propto aLB_{{\rm m}}\frac{1}{T_{0}}, \end{equation}
   is the characteristic time of the non-adiabatic loss. It turns out that the efficiency of the axial confinement does not depend explicitly on the mirror ratio $\mathcal {R}$, as for a ‘classical’ GDT $\tau _{\mathrm {GD}}\propto \mathcal {R}$, but it is enhanced with increasing the absolute value of the mirror magnetic field $B_{{\rm m}}$. In addition, the warm plasma density $n_{{\rm i}0}$ and the confinement time $\tau _{{\rm i}\parallel 0}$ indeed prove to be greater at a lower warm plasma temperature $T_{0}$. The inverse dependence of $n_{{\rm i}0}$ on $T_{0}$ is the result of the non-adiabatic loss (4.23) being proportional to $\rho _{{\rm i}0}v_{T{\rm i}0}\propto T_{0}$.

8. (viii) It can be observed from the simulations (see figure 4 and table 1) that the warm plasma temperature $T_{0}$ proves to be higher in the regimes with a stronger vacuum magnetic field $B_{{\rm v}}$. Indeed, sustaining the diamagnetic confinement equilibrium with stronger external field $B_{{\rm v}}$ (at fixed mirror field $B_{{\rm m}}$, length $L$ and radius $a$ of the bubble core) requires greater absorbed injection power $Q_{{\rm h}0}$, while the increased heating power $Q_{{\rm h}0}$, in turn, should result in a corresponding rise in the equilibrium temperature of the warm plasma $T_{0}$. This relation can be clarified by considering the force balance in the case of beam plasma (8.1)

   (8.5)\begin{equation} \varPi_{{\rm h}}\sim\varPi_{{\rm M}}, \end{equation}
   where
   (8.6a,b)\begin{equation} \varPi_{{\rm h}}\sim-\frac{Q_{{\rm h}0}}{\left\langle \dot{\mathcal{E}}\right\rangle _{\varGamma}V}\sim\frac{Q_{{\rm h}0}}{\nu_{\mathrm{ie}0}{\rm \pi} a^{2}L},\quad\varPi_{{\rm M}}\sim\frac{B_{{\rm v}}^{2}}{8{\rm \pi}}, \end{equation}
   are the characteristic energy densities of the hot ions and the vacuum magnetic field, respectively. This expression roughly corresponds to (7.25) in the limit: $\beta _{{\rm h}0}\rightarrow 1$, $\beta _{{\rm w}0}\rightarrow 0$ and $\bar {r}_{\min }/a\rightarrow 0$. Taking into account the estimate (8.3), we obtain the following relation:
   (8.7)\begin{equation} B_{{\rm v}}^{2}\sim\frac{8Q_{{\rm h}0}}{\nu_{\mathrm{ie}0}a^{2}L}\propto\frac{1}{aLB_{{\rm m}}}T_{0}^{{7}/{2}}.\end{equation}
   Therefore, the temperature of warm plasma indeed proves to increase with the vacuum magnetic field as $T_{0}\propto B_{{\rm v}}^{{4}/{7}}$.

9. (ix) Having fixed the mirror field $B_{{\rm m}}$, the warm plasma source $W_{{\rm i}0}$, the radius $a$ and the length $L$ of the bubble core, from the considered above qualitative estimates (8.3) and (8.7), we get

   (8.8a,b)\begin{equation} n_{{\rm i}0}\propto B_{{\rm v}}^{-{4}/{7}},\quad T_{0}\propto B_{{\rm v}}^{{4}/{7}}. \end{equation}
   Given the error due to the corresponding inaccuracy of the estimate (8.2), these relations agree well with the results of the simulations shown in figure 4 and table 1.

10. (x) In the present paper, we assume the classical model of transverse transport for the warm plasma, which corresponds to the Spitzer conductivity (7.28). At the same time, large gradients of equilibrium parameters, such as the magnetic field and the warm plasma density, inside the transition layer of the diamagnetic bubble could probably lead to non-classical transport, which corresponds to a much greater diffusion coefficient. In turn, a significant modification of the transverse transport could apparently have a qualitative impact on the equilibrium of the warm plasma. However, the issue concerning non-classical diffusion in the diamagnetic confinement mode remains poorly studied so far.

11. (xi) The model of the warm plasma equilibrium, constructed in § 4 within the framework of MHD, assumes isotropic pressure and gas-dynamic axial loss. This is known to correspond to the collisional regime with a filled loss cone. In the case of GDT (Ivanov & Prikhodko 2017), such a regime is realized when the gas-dynamic time $\tau _{\mathrm {GD}}\sim \mathcal {R}L/v_{T{\rm i}}$ significantly exceeds the kinetic time $\tau _{\mathrm {kin}}\sim \tau _{\mathrm {ii}}\ln \mathcal {R}$, where $\tau _{\mathrm {ii}}$ is the mean free time of ion–ion Coulomb collisions (Trubnikov 1965). However, when considering a diamagnetic trap, one should also take into account the effective increase in the mirror ratio $\mathfrak {R}\geq \mathcal {R}$, which results from the magnetic field weakening induced by the strong diamagnetic current of the high-pressure plasma. In addition, the exotic conditions of the diamagnetic confinement mode are likely to lead to anomalous collisionality with an effective mean free time $\tau _{\mathrm {eff}}$, which typically proves to be considerably shorter than the ‘classical’ time $\tau _{\mathrm {ii}}$. Thus, the warm plasma loss cone in the transition layer of a diamagnetic bubble can be considered filled when

    (8.9)\begin{equation} \frac{L\mathfrak{R}}{v_{T{\rm i}}}\gg\tau_{\mathrm{eff}}\ln\mathfrak{R}. \end{equation}
    This condition seems to be always satisfied at least in some neighbourhood of the bubble core $r\gtrsim a$, where the magnetic field approaches zero $B\rightarrow 0$ and, accordingly, the effective mirror ratio is greatly increased $\mathfrak {R}\rightarrow +\infty$. At the periphery of the warm plasma $r\lesssim a+\lambda _{{\rm w}}$, the loss cone may turn out to be only partially filled ($L\mathfrak {R}/v_{T{\rm i}}\sim \tau _{\mathrm {eff}}\ln \mathfrak {R}$) or even empty ($L\mathfrak {R}/v_{T{\rm i}}\ll \tau _{\mathrm {eff}}\ln \mathfrak{R}$), depending on the specific value of the anomalous mean free time $\tau _{\mathrm {eff}}$.

9 Summary

In the present work, we have constructed a theoretical model of plasma equilibrium in the diamagnetic confinement regime in an axisymmetric mirror device with neutral beam injection. To describe the equilibrium of the background warm plasma, we use the MHD equations of mass and energy conservation, as well as the force balance equation. Transverse transport is considered to be due to resistive diffusion, and the axial loss includes both the ‘classical’ gas-dynamic loss and the non-adiabatic loss (Chernoshtanov 2020, 2022) specific to the diamagnetic bubble. This model does not take into account the effects of the warm plasma rotation and inhomogeneity of electrostatic potential. The equilibrium of the hot ions resulting from the neutral beam injection is described by the distribution function, which is found from the corresponding kinetic equation. It takes into account the warm electron drag force acting on the hot ions, while the angular scattering due to Coulomb collisions is neglected. The chaotic nature of the dynamics of ions in the diamagnetic bubble is taken into consideration as well. Solving the equilibrium equations yields the equilibrium profiles of the warm plasma density, temperature and pressure; the hot ion current density is obtained from the corresponding distribution function; the equilibrium magnetic field is determined by Maxwell's equations.

Applying the constructed theoretical model, we have considered the case of the cylindrical core of the diamagnetic bubble. In this case, the equilibrium model is reduced to a simpler one-dimensional problem. To further simplify the equations, we have assumed the approximation of a thin transition layer. This allows the hot ion current density and the magnetic field to be explicitly expressed in terms of the magnetic flux. Finally, it has been found that the radial profile of the warm plasma typically turns out to be almost isothermal, which enables the system of two equilibrium equations for the warm plasma to be approximately reduced to a single one. For the resulting equation of the simplified equilibrium model, we have constructed a numerical solution algorithm that includes a variation of the shooting method in combination with the Runge–Kutta schemes.

The equilibria of the diamagnetic bubble have been found for the parameters of the GDMT device (Skovorodin et al. 2023). Three cases corresponding to different vacuum magnetic fields in the central plane, $B_{{\rm v}}=5,10,20\ \mathrm {kG}$, have been considered. All the other external parameters are fixed except the absorbed injection power, which is adjusted so that the radius of the bubble core remains the same $a=20\ \mathrm {cm}$. Equilibrium profiles of the warm plasma density, temperature and pressure have been found. In addition, the radial distributions of the magnetic field and the diamagnetic current densities of the warm plasma and the hot ions have been constructed.

It has been found that, for all the solutions obtained, the main contribution to the plasma energy content and the diamagnetism comes from the hot ions. This corresponds to a negligible relative pressure of the warm plasma $\beta _{{\rm w}0}\ll \beta _{{\rm h}0}\sim 1$. In addition, the non-adiabatic loss turns out to be dominant in all considered regimes, and its fraction in total loss is greater for stronger vacuum magnetic field. This seems to be due to the warm plasma in the transition layer being confined in the magnetic field weakened by the hot ion diamagnetism, which should lead to an increased mirror ratio and improved axial confinement. The transition layer of the diamagnetic bubble turns out to be quite thin in the regimes with the stronger vacuum fields, $B_{{\rm v}}=10,20\ \mathrm {kG}$. However, in the case of the weak external field, $B_{{\rm v}}=5\ \mathrm {kG}$, the thickness of the transition layer proves to be not too small, which may correspond to a reduced accuracy of the approximate solution. Qualitative estimates of the warm plasma density and the temperature of the warm plasma have also been obtained.

9.1 Future work

This work should be considered as a basis for further expansion of the theoretical model constructed here. In the near future, the equilibrium model presented in this paper will be used for a detailed investigation of the diamagnetic confinement regime in GDMT. In particular, it is planned to study the dependence of the bubble equilibrium on external conditions, such as heating power, vacuum magnetic configuration, warm plasma source, etc. In addition, the finite absorption efficiency of the injected neutral beam can be taken into account. The constructed model can also be refined by eliminating the thin transition layer approximation; this, however, would require more complex numerical simulations. In perspective, it is also planned to include in the model the effects of the warm plasma rotation and the inhomogeneity of the electrostatic potential. In the long term, collisional angular scattering of the hot ions can also be taken into account.