Effects observed in numerical simulation of high-beta plasma with hot ions in an axisymmetric mirror machine

Abstract

We present the results of numerical simulation by two-dimensional hybrid particle-in-cell code of high-beta plasma with hot ions in an axisymmetric mirror machine. Two particular effects are discussed: the self-rotating of plasma with Maxwellian ions in regime of diamagnetic confinement and the excitation of axisymmetric magnetosonic waves in a high-beta plasma with sloshing ions.

1. Introduction

One of the advantages of axisymmetric mirror machines is the possibility to confine the plasma with high ratio of plasma pressure to magnetic field pressure. This allows thermonuclear systems with high plasma pressure and high density of thermonuclear power to be constructed. In the limiting case, when the plasma pressure is equal to the pressure of the magnetic field, a qualitative change of regime of the plasma confinement occurs, in particular, the transition to the regime of diamagnetic confinement takes place (Beklemishev 2016). In this regime the magnetic field is displaced from the region occupied by the plasma and the so-called diamagnetic bubble is formed. Analytical estimates based on magnetohydrodynamic (MHD) (Beklemishev 2016; Beklemishev & Khristo 2019) and kinetic approaches (Chernoshtanov 2022; Khristo & Beklemishev 2022) have shown that the transition to the regime of diamagnetic confinement allows us to increase the lifetime of particles in the trap by at least one order of magnitude.

The potential benefits of this regime have led to efforts to validate it experimentally. In particular, similar regimes are being investigated on a C2-W device of the Tri Alpha Energy company (Gota et al. 2019) and a high-beta plasma experiment is planned on a CAT device of the Budker Institute of Nuclear Physics (BINP) (Bagryansky et al. 2016). Moreover, the experiment with high-beta plasma on a 2XIIB device (Turner et al. 1989) and experiments in cusps (Haines 1977) should also be mentioned.

The non-trivial structure of the magnetic field and of the distribution function of plasma particles results in numerous difficulties in analytical consideration of the diamagnetic regime. The numerical simulation in turn can be useful to study in detail the different effects arising in this regime. In this article we discuss the plasma rotation and the excitation of axisymmetric oscillations observed in simulations with hybrid particle-in-cell (PIC) numerical code MTRAP (Boronina et al. 2019).

The article is organized as follows. In § 2 we briefly discuss the equations on which the numerical model is based. Special attention is paid on the treatment of ion–ion Coulomb collisions. In § 3 the simulation results, in particular, the plasma rotation and the associated effects are presented. In § 4 the axisymmetric oscillations observed in the simulation of a diamagnetic bubble in a mirror machine with skew neutral beam injection are discussed. A brief conclusion is given in § 5.

2. Description of the model

In this section, we describe the basic equations of our model, while the details of numerical realization can be found in Boronina et al. (2019). We consider a plasma which consists of electrons and one type of singly charged ions (in this paper the ions are protons). We treat all plasma ions kinetically and use the PIC method (see e.g. Hockney & Eastwood 2021) to solve the equation for the ion distribution function using the macroions of the same charge-to-mass ratio as the real ions in the plasma. Plasma electrons are instead treated as a charge-neutralizing massless fluid.¹

2.1. Basic equations

The equations of motion for the macroions read

(2.1a,b)\begin{equation} \frac{{\rm d}{\boldsymbol{r}}_i}{{\rm d}t}={\boldsymbol{v}}_i,\quad m_i\frac{{\rm d}{\boldsymbol{v}}_i}{{\rm d}t} =e{\boldsymbol{E}}(\boldsymbol{r},t)+\frac{e}{c}{\boldsymbol{v}}_i \times{\boldsymbol{B}(\boldsymbol{r},t)}+m_e\frac{{\boldsymbol{V}}_e(\boldsymbol{r},t) -{\boldsymbol{V}}_i(\boldsymbol{r},t)}{\tau_{ie}}, \end{equation}

here $m_i, m_e$ are the ion and electron mass, $\boldsymbol {E}$ and $\boldsymbol {B}$ are the electric and magnetic fields, $\tau _{ie}$ is the electron-ion collision time (time of the electron drag), $\boldsymbol {r}_i, \boldsymbol {v}_i$ are position and velocity of each individual macroion and $\boldsymbol {V}_i(\boldsymbol {r},t)$ is the mean velocity of macroions, obtained (along with the ion charge density $n_i$) from their spatial distribution. We assume the plasma to be quasineutral, $n_e=n_i,$ and write the equation of motion for the massless electron fluid as

(2.2)\begin{equation} e\boldsymbol{E}(\boldsymbol{r},t)+\frac{e}{c}\boldsymbol{V}_e(\boldsymbol{r},t)\times\boldsymbol{B}(\boldsymbol{r},t)+\frac{\boldsymbol{\nabla} p_e}{n_e}+m_e\frac{\boldsymbol{V}_e(\boldsymbol{r},t)-\boldsymbol{V}_i(\boldsymbol{r},t)}{\tau_{ie}}=0, \end{equation}

here $p_e=n_eT_e$ is the scalar electron pressure. In the hybrid algorithm this equation is used to obtain the electric field $\boldsymbol {E}(\boldsymbol {r},t).$ The heat equation for the electron fluid reads

(2.3)\begin{equation} n_e\left(\frac{\partial T_e}{\partial t}+(\boldsymbol{V}_e(\boldsymbol{r},t)\boldsymbol{\cdot}\boldsymbol{\nabla})T_e\right)+(\gamma-1)\left(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{V}_e(\boldsymbol{r},t)\right)p_e=(\gamma-1)\left(Q_e-\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q}_e\right). \end{equation}

Here $\gamma =\tfrac {5}{3}$ is the adiabatic index of an atomic gas, $Q_e={J^2}/{\sigma }$ is the heat generated in electrons, with the electric conductivity $\sigma = {e^2 n_e\tau _{ie}}/{m_e}$ and the total current density $\boldsymbol {J}=en_e(\boldsymbol {V}_i-\boldsymbol {V}_e)$. The electronic heat flux reads $\boldsymbol {q}_e=\kappa \boldsymbol {\nabla } T_e$ with $\kappa$ being the thermal conductivity. In our hybrid model we consider low-frequency processes and do not take into account the displacement current. Therefore, the total current density $\boldsymbol {J}$ can be obtained from Ampere's law leading to the following equation for the electron fluid velocity:

(2.4)\begin{equation} \boldsymbol{V}_e(\boldsymbol{r},t)=\boldsymbol{V}_i(\boldsymbol{r},t)-\frac{c}{4{\rm \pi}}\frac{1}{e n_e(\boldsymbol{r},t)}\boldsymbol{\nabla}\times\boldsymbol{B}(\boldsymbol{r},t). \end{equation}

Finally, the magnetic field evolves according to Faraday's law

(2.5)\begin{equation} \frac{\partial\boldsymbol{B}(\boldsymbol{r},t)}{\partial t} ={-}c \boldsymbol{\nabla}\times\boldsymbol{E}(\boldsymbol{r},t). \end{equation}

2.2. Initial set-up

At $t = 0$ the cold uniform background plasma with density $n_0$ is located inside a cylindrical chamber of radius $R_0$ and length $L.$ A beam of neutral hydrogen plasma is injected into the chamber. The particular initial distributions of the beam ions are specified in §§ 3.1 and 4.1. As for the initial configuration of the magnetic field, for the simulations presented in § 4 the vacuum field is homogeneous, while in § 3 it is generated by a pair of identical current coils, located at both ends of the chamber. We assume that the current flows through the coils strictly in the azimuthal direction $\varphi,$ and that its amplitude does not depend on the azimuthal angle. In this case at $t = 0$ the whole system has an axial symmetry. The amplitude of the magnetic field in the centre of the trap at the point $r = 0, z = 0$ is equal to $B_0.$ The amplitude of the field at the points $r = 0, z = \pm L/2$ is defined by the given mirror ratio $R_v.$ In this paper we consider the reduced two-dimensional cylindrical problem, with $z$-axis directed along the symmetry axis of the trap, i.e. we assume that the field configuration and the evolution of the plasma do not depend on the azimuthal angle through the whole process. In cylindrical geometry the simulation box is a rectangle $[-L/2, L/2] \times [0,R_0].$ In order to find the initial distribution of the magnetic field inside the trap for the simulations, discussed in § 3, we introduce the vector potential $\boldsymbol {A}=(0,A_\varphi, 0)$ such that $\boldsymbol {\nabla }\times \boldsymbol {A}=\boldsymbol {B}=(B_r,0,B_z).$ The components of the magnetic field then read

(2.6a,b)\begin{equation} B_r={-}\frac{\partial A_\varphi}{\partial z},\quad B_z=\frac{1}{r}\frac{\partial (rA_\varphi)}{\partial r}. \end{equation}

The initial distribution of the magnetic field is found from the numerical solution of the elliptic equation for $A_\varphi,$ which reads

(2.7)\begin{equation} \frac{\partial}{\partial r}\left[\frac{1}{r}\frac{\partial (rA_\varphi)}{\partial r}\right]+\frac{\partial^2 A_\varphi}{\partial z^2}={-}\frac{4{\rm \pi}}{c}j_\varphi, \end{equation}

with $j_\varphi$ being the current density in the coils. The details of the numerical algorithm can be found in Liseykina, Vshivkov & Kholiyarov (2024).

2.3. Simulation procedure and boundary conditions

The simulation procedure is as follows. (i) Setting the initial homogeneous distribution of the cold background plasma and the trap magnetic field (homogeneous in § 4 or produced by a pair of current coils in § 3): $n_i=n_0; \boldsymbol {V}_i=0; \boldsymbol {V}_e=0; T_e=0; p_e=0, \boldsymbol {B}=(B_r,0,B_z).$ The electric field at $t=0$ is obviously set to zero. (ii) Injection of beam ions. (iii) Update of the ion positions and velocities, (2.1a,b) and calculation of $n_i$ and $\boldsymbol {V}_i$ from the spatial distribution of ions. (iv) Calculation of the electron fluid velocity $\boldsymbol {V}_e$ using (2.4). (v) Update of the electric field $\boldsymbol {E},$ (2.2). (vi) Update of the magnetic field $\boldsymbol {B}$, (2.5). (vii) Update of the electron temperature $T_e,$ (2.3). (viii) Repeat steps (ii)–(vii) until necessary. Note that although the density is calculated with second-order accuracy in space and the velocity with second-order accuracy in time, the overall numerical scheme is first order in time and space. The boundary conditions on the axis $r=0$ for vector functions $\boldsymbol {F}=\boldsymbol {E},\boldsymbol {B},\boldsymbol {V}$ are $F_r=F_\varphi =0, \partial F_z/\partial r=0,$ while for scalar functions they read $\partial n/\partial r= \partial T_e/\partial r=0.$ Particles can leave the computational region through lateral boundaries at $z=\pm L/2.$ Radial component of the electric field $E_r$ and the derivative $\partial E_\varphi /\partial z$ are zero at $z=\pm L/2$.

2.4. Normalization and typical simulation parameters

In the simulations all quantities are normalized to the ion cyclotron frequency $\varOmega _0=eB_0/(m_ic)$, the Alfvén velocity $v_A=B_0/(4{\rm \pi} m_in_0)^{1/2}$ and the ion inertial length $c/\omega _{pi}=(m_ic^2/(4{\rm \pi} n_0e^2))^{1/2}.$ We choose the following normalization parameters: $B_0=2\,{\rm kG}$ and $n_0=10^{13}\,{\rm cm}^{-3}$, so that $\varOmega _0=1.9\times 10^7\,{\rm s}^{-1}$, $v_A=1.38\times 10^8\,{\rm cm}\,{\rm s}^{-1}$ and $c/\omega _{pi}=7.19\,{\rm cm}$.

The numerical parameters are the following: time step $10^{-4}\varOmega _0^{-1}$; grid size $[102\times 42]$. In the quasistationary regime the number of macroparticles in calculation area is $2\times 10^4$ for hot ions and $1.5\times 10^4$ for background ions. Typical simulation until $t\simeq 200 \varOmega _0^{-1}$ on one processor Intel Xeon Platinum 8268 (2.9 GHz) takes four hours.

2.5. Simulation of ion–ion collisions

The accurate treatment of the ion–ion collisions with the conservation of the full momentum and energy of ions is important, for example, for the simulation of the outflow of Maxwellian plasma. When applying the frequently used Takizuka–Abe method (Takizuka & Abe 1977) to model pair collisions, the amount of computation grows quadratically with increasing number of macroparticles. Therefore, we applied the null-collision Monte Carlo algorithm (Vshivkov et al. 2021). The idea of this method is following. We calculate the density of ions $n_i$, their mean velocity $\boldsymbol {V}_i$ and the root mean square velocity spread $\Delta v_i=(\sum (\boldsymbol {v}_i-\boldsymbol { V}_i)^2)^{1/2}$ in each numerical cell. Then each ion is scattered in the same way as in case of scattering on plasma with Maxwellian ions with density $n_i$, mean velocity $\boldsymbol {V}_i$ and temperature $T_i=m_i({\Delta v_i^2}/{2})$. To ensure the conservation of the full energy and momentum of ions the procedure of renormalization is applied. Details of the particular realization can be found in Vshivkov et al. (2021).

3. Simulation of the outflow of Maxwellian plasma

There exist semianalytical MHD (Beklemishev 2016) and kinetic (Chernoshtanov 2022) models describing gas-dynamic confinement of Maxwellian plasma in a long quasicylindrical diamagnetic bubble. The MHD models are applied when the mean free path of the ions is smaller than the width of the transition layer, while the kinetic model is applicable when the ion free path is shorter than the distance between magnetic mirrors $L.$ Both models predict that the confinement time $\tau$ depends linearly on the width of the transition layer $\lambda.$ In particular,

(3.1a,b)\begin{equation} \tau=\tau_{gdt}\frac{a}{\lambda},\quad\tau_{gdt}=\frac{R_vL}{v_t}, \end{equation}

where $\tau _{gdt}$ is the time of gas-dynamic outflow from the vacuum magnetic field with the mirror ratio $R_v$, $v_t$ is the thermal velocity of ions and $a$ is the bubble radius. The width of the transition layer of the diamagnetic bubble $\lambda$ depends either on the plasma resistivity in the MHD model or on Larmor radius of ions $v_t/\varOmega _0$ and the electrostatic potential in the transition layer in the kinetic model.

It should be noted that the plasma is assumed to be non-rotating in both models. On the other hand, there are different factors resulting in plasma rotation, such as the neutral beam injection, the radial electric field induced by the sectioned plasma plates (Ivanov & Prikhodko 2017) and limiters, the ambipolar plasma potential, as well as different kinetic effects. Moreover, in the diamagnetic trap and in the field reversal configurations (FRCs) the important effect is plasma spin-up (Steinhauer 2011). This effect arises because of the differences in the confinement of particles with different azimuthal momentum $p_\varphi.$ Namely, particles in axisymmetric systems with high-beta plasma can move chaotically due to the non-conservation of magnetic momentum $\mu$. Chaotically moving particles with $\varOmega _0p_\varphi >0$ predominantly leave the trap, while the particles with $\varOmega _0p_\varphi <0$ can be confined in the trap due to the so-called absolute confinement (Steinhauer 2011). Therefore, in the diamagnetic confinement regime plasma can spin-up in a direction which coincides with the direction of the cyclotron rotation of ions. This rotation can sufficiently influence the plasma confinement.

3.1. Formulation of the problem

We consider an axisymmetric mirror machine, the magnitude of the magnetic field in the minimum is $B_0$, the mirror ratio is $R_v$ and distance between the mirrors is $L$. Ions with a Maxwell velocity distribution $f_i=\exp ({-v^2/v_t^2})/({\rm \pi} ^{3/2}v_t^3)$ arise in the region bounded by the fixed field line, namely in the region with $\varPsi < B_0\psi _0,$ where $\varPsi =\int _0^r r'B_z(r')\,{\rm d}r'$ is the magnetic flux. For a short time (20 cyclotron periods), Maxwell ions are injected into the trap, resulting in the formation of a diamagnetic bubble. Subsequently the injection rate is decreased by several times and the transition to a quasistationary state, when all variables are almost time-independent, is investigated. We choose the following basic simulation parameters: the distance between mirrors $L=10c/\omega _{pi}$; the thermal velocity of ions is $v_t=0.2 v_A$; the mirror ratio of the vacuum field is $R_v=2,$ $1/\tau _{ie}=0.03\varOmega _0$. The injection rate corresponds to $q=2.6\times 10^{23}$ real protons injected per second for $t<20/\varOmega _0$ and it is reduced down to $q=3.6\times 10^{22}$ for $t>20/\varOmega _0.$

To simplify the interpretation of the simulation results, we exclude most of sources of plasma rotation, such as off-axis neutral beam injection and the radial electric field generated by external plasma plates. The main source of plasma rotation in our case is the predominant loss of ions with $\varOmega _0 p_\varphi >0.$

The parameters were chosen so that the characteristic values of plasma density and thermal ion velocity were of the order of (see below) $10^{13}\,{\rm cm}$ and $10^7\,{\rm cm}\,{\rm s}^{-1}$ (ion temperature of the order of 100 eV). These values correspond to the mean free path of ions with respect to the ion-to-ion Coulomb collisions of the order of 10 cm. It is smaller than the distance between the mirrors (71.9 cm), so the plasma is gas-dynamically confined. Increasing the frequency of Coulomb collisions slightly broadens the plasma radial distribution due to the growth of radial diffusion.

3.2. Simulation results

Let us discuss the time-dependence of the plasma parameters. The dependence of the full number $N$ of the injected ions in the trap, the transverse $W_\perp ={\sum }_im_i(v_r^2+v_\varphi ^2)/2$ and the longitudinal $W_{\|}={\sum }_im_iv_z^2/2$ parts of the kinetic energy of ions over time are shown in figure 1. After decreasing the injection rate the kinetic energy of ions transits to a quasistationary state within the time period which is of the order of time of flight from one mirror to another $L/v_t=50\varOmega _0^{-1}$. Due to Coulomb collisions the transverse and the longitudinal pressures are approximately equal,² $W_\perp /2\approx W_\|$.

Figure 1. Dependence of full number of ions $N$ (a) and full transverse $W_\perp$ (solid line) and longitudinal $W_\|$ (dashed line) kinetic energy of ions (b) in time. Parameters: $R_v=2$; $\psi _0=0.08(c/\omega _{pi})^2$; $1/\tau _{ie}=0.03\varOmega _0$; injection rate is $q=2.6\times 10^{23}\,\mathrm {particles}\,\mathrm {s}^{-1}$ at $\varOmega _0 t<20$ and $q=3.6\times 10^{22}\,\mathrm {particles}\,\mathrm {s}^{-1}$ at $\varOmega _0 t>20$. The energy is measured in units of $m_iv_A^2,$ the injection rate and the number of ions both refer to the number and injection rate of real protons.

The decrease of the magnetic field inside the diamagnetic bubble is due to the diamagnetic current flowing in the transition layer of the bubble. This means that the average azimuthal velocities of electrons and ions are different and that electron drag slows down the ions in the transition layer. This friction force causes ions to drift in the radial direction outwards of the bubble. This drift, combined with Coulomb scattering, leads to broadening of the transition layer and diffusion of ions outwards of the bubble. A part of these ions is cooled by the electron drag, so they leave the trap slowly. The existence of this population of cold ions at the periphery of the plasma explains the slow growth of the total number of ions in the trap $N,$ see figure 1.

Snapshots of the magnetic field amplitude and magnetic field lines at different time moments are shown in figure 2. These pictures demonstrate a typical structure of a diamagnetic bubble with a zero magnetic field in the central region and a vacuum magnetic field in the periphery.

Figure 2. The distribution of the amplitude of the magnetic field $|B(z,r)|$ inside the trap (colours) and the field lines of the magnetic field at $t=0$ (a) and $t=20\varOmega _0^{-1}$ (b). The field lines are labelled by the values of $rA_\varphi.$ The parameters of the simulations are the same as in figure 1.

When a particle moves in an axisymmetric electromagnetic field its azimuthal momentum $mrv_\varphi +e\varPsi (r)/c$ is conserved. Hence, to study plasma rotation it is convenient to plot the time-dependence of full azimuthal momentum ${\sum }_{j=e,i}(m_jrv_\varphi +e_j\varPsi (r)/c).$ Here summation is made for both species of plasma particles (ions and electrons). Note, that ${\sum }_j(e_j/c)\varPsi (r)\approx 0$ because of the plasma quasineutrality and the contribution of electrons to the sum ${\sum }_{j=e,i} m_jrv_\varphi$ is small compared with the contribution of ions due to $m_e/m_i \ll 1.$ In the following the reduced azimuthal momentum $P_\varphi ={\sum }_{j=i} m_jrv_\varphi$ is investigated.

The time-dependencies of the root mean square velocity spread $\bar {v}=\sqrt {(W_\perp +W_\|)/(3N)}$ and the mean azimuthal momentum of ions $\bar P_\varphi =N^{-1}{\sum }_im_irv_\varphi$ are shown in figure 3. Initially the mean velocity of ions decreases due to the electron drag and the longitudinal losses of most energetic ions. Besides, at $\varOmega _0 t<20,$ a part of the kinetic energy of ions is spent on the displacement of the magnetic field. The mean azimuthal momentum increases (in absolute value) initially because of the longitudinal losses of ions with $\varOmega _0 p_\varphi >0.$ Such spin-up is stopped later because of the electron drag in the transition layer. The mean azimuthal velocity of ions can be estimated as the ratio of $\bar P_\varphi$ to the bubble radius $a\approx 0.5 c/\omega _{pi}.$ The mean azimuthal velocity is approximately $0.03 v_A,$ which is two times smaller than the mean root square velocity spread of ions in a stationary state $\bar v\approx 0.07 v_A$.

Figure 3. The time dependence of the root mean square velocity of all ions $\bar v$ (a) and mean azimuthal momentum $\bar P_\varphi$ (b). The parameters are the same as in figure 1. The velocity is measured in units of $v_A$, the azimuthal momentum in units of $m_iv_Ac/\omega _{pi}$.

In the following we estimate the maximum of the total azimuthal momentum. The distribution function of the injected ions is approximately equal to

(3.2)\begin{equation} f\approx H(m^2v^2a^2-p_\varphi^2)\exp({-v^2/v_t^2})/({\rm \pi}^{3/2}v_t^3), \end{equation}

with $H(x)$ being the Heaviside function and $a$ being the radius of the diamagnetic bubble. If all ions with $\varOmega _0 p_\varphi >0$ leave the trap, the mean azimuthal momentum will be $\bar {P}_\varphi /m_i\sim -v_ta\approx 0.08 v_A c/\omega _{pi}.$ This value is several times bigger that the observed mean azimuthal momentum.

To verify the plausible assumption that plasma rotation is restricted by the electron drag in the transition layer, we performed a set of simulations with different electron–ion collision times, $\tau _{ie}\in [25,100]\varOmega _0^{-1}.$ The temporal evolution of $W_\perp$ plotted in figure 4 for three different values of $\tau _{ie}$ indicates that the full kinetic energy of ions in a quasistationary state depends on $\tau _{ie}$ as $W\propto \tau _{ie}^{1/2}.$ Such dependence agrees with the predictions of the MHD model (Beklemishev 2016), where ion losses are also caused by the electron drag in the transition layer.³

Figure 4. Dependence of full transverse kinetic energy $W_\perp$ of ions in time for $1/\tau _{ie}=0.02\varOmega _0$ (cyan line, divided by $\sqrt {3/2}$), $1/\tau _{ie}=0.03\varOmega _0$ (black line) and $1/\tau _{ie}=0.04\varOmega _0$ (blue line, multiplied to $\sqrt {4/3}$). The parameters are the same as in figure 1. The energy is measured in units of $m_iv_A^2$.

4. Axisymmetric waves in a plasma with the population of sloshing ions

Skew neutral beam injection is used in linear systems for the plasma confinement, for example in a Gas Dynamic Trap (GDT) device (Ivanov & Prikhodko 2017), a C2-W device (Gota et al. 2019) and a GOL-NB device (Postupaev et al. 2017), and was used in a TMX-U device (Simonen et al. 1983), as a method of heating and maintaining the plasma in the trap. Such injection results in the formation of a population of fast sloshing ions oscillating in the longitudinal direction between the magnetic mirrors. In the case of small angular spread of neutral beams, the density of the fast ions is peaked near the turning points. These effects can be used to construct the compact source of thermonuclear neutrons with peaked density of neutron flow (Ivanov & Prikhodko 2017).

Preliminary numerical simulations of a mirror machine with neutral beam injection show that axisymmetric oscillations with the frequency close to the bounce-frequency of sloshing ions between the magnetic mirrors are excited. These oscillations are similar to the global sound modes observed in a GDT device (Skovorodin, Zaytsev & Beklemishev 2013). To study the spatial structure of these oscillations, a set of simulations of plasma with sloshing ions was performed. To simplify the study of the wave structure and to accelerate the attainment of the quasistationary state, we model the plasma dynamics in a homogeneous magnetic field.

4.1. Formulation of the problem

A region with a uniform magnetic field with amplitude $B_0$ (mirror ratio $R_v=1$) and cold background plasma of density $n_0$ is filled by fast ions. The distribution function of fast ions of the source is

(4.1) \begin{align} & S(r,z,v_r,v_\varphi,v_z,t)=F(\boldsymbol{v})\cdot H(\Delta r^2-(r-a)^2) \nonumber\\ & \qquad {\cdot}\left.\left[H((z-z_\pm){\mbox{sign}}(v_z))+\delta(z-z_\pm)v_zt\right]\right/\left(L/|v_z|+t\right), \end{align}

where $L$ is length of the simulation region and $z_\pm =0$ if $v_z>0$ and $z_\pm =L$ if $v_z<0$. This source allows us to generate a population of fast ions with density which does not depend on the longitudinal coordinate. This results in a population of fast ions, localized near the surface of a cylinder with radius $a$. We choose $F(\boldsymbol {v})$ in the following form:

(4.2)\begin{equation} F(\boldsymbol{v})=\delta(v_r)\delta(v_\varphi-v_0)(\exp({-(v_z-v_{\|0})^2/w^2})+\exp({-(v_z+v_{\|0})^2/w^2}) \end{equation}

with $v_0\approx -\varOmega _0a$. Such velocity distribution models the distributions of fast ions in a mirror machine with skew off-axis neutral beam injection.

We choose the following parameters: injection rate, corresponding to the injection rate of real protons, $q=1.2\times 10^{21}\,{\rm particles}\,{\rm s}^{-1}$, Larmor radius of fast ions $a=1.0\times c/\omega _{pi}$ (so that $v_0=1.0\times v_A$), radial spread $\Delta r=0.1\times a$, the mean longitudinal velocity of fast ions $v_{\|0}=1.0\times v_A$ and the longitudinal velocity spread $w=0.5\times v_A$ or $w=1.0\times v_A.$

4.2. Numerical results

The typical time-dependence of the plasma parameters is shown in figure 5. Oscillations of the magnetic field amplitude are clearly visible. The period of oscillations is approximately $10/\varOmega _0$ which is close to the time of an overflow of fast ions $L/v_{\|0}$. The amplitude of oscillations decreases with increasing of the spread of the longitudinal velocity $w$.

Figure 5. The time-dependence of the full number of ions (a) and $z$-component of the magnetic field at $z=0, r=0.$ The parameters are: injection rate is $q=1.2\times 10^{21}\,\mathrm {particles}\,\mathrm {s}^{-1}$; $a=1.0\times c/\omega _{pi}$, $v_{\|0}=1.0\times v_A$; $w=0.5\times v_A$ (solid line) and $w=1.0\times v_A$ (dashed line). The magnetic field is measured in units of $B_0,$ the injection rate and the number of particles are both given (recalculated) for the number and injection rate of the real protons.

The typical time dependence of the azimuthal component of the electric field is shown in figure 6. The mean value of the field is not zero because of the monotonically decreasing magnetic field $B_z,$ which is displaced due to the growing diamagnetism of the plasma. The amplitude of the electric field oscillates with a frequency equal to the frequency of magnetic field oscillations.

Figure 6. The time-dependence of the azimuthal component of the electric field at $z=0$ and $r=1.0\times c/\omega _{pi}$. The parameters are the same as in figure 5.

The dependence of the magnetic field on the longitudinal coordinate (figure 7a) demonstrates the standing wave with wavelength equal to length of the simulation box. The standing wave is formed due to boundary conditions ($E_r=0$ and $\partial _zE_\varphi =0$) on the left and right boundaries of the computation region. It is convenient to show the radial structure of the wave through the azimuthal component of the electric field $E_\varphi,$ which is equal to zero if the wave is absent. The electric field is peaked at $r\approx a$, where the density of fast ions reaches a maximum. The phase velocity of the wave is close to the mean longitudinal velocity of fast ions $v_{\|0}$. The azimuthal electric field does not depend on time for a population of ions with the longitudinal velocity $v_\|=\omega /k_\|$, where $\omega$ is wave frequency and $k_\|$ is the longitudinal wavenumber. These ions can effectively transfer energy to the wave through the Landau mechanism. The increase of the longitudinal velocity spread $w$ reduces the fraction of ions which can interact resonantly with the wave, so that the amplitude of the wave decreases.

Figure 7. The dependence of $z$-component of the magnetic field on longitudinal coordinate $z$ (a) at $r=0$ and of the azimuthal component of electric field on radial coordinate $r$ (b) at $z=0$ at $\varOmega _0t=80$. The parameters are the same as in figure 5, $w=1.0\times v_A$.

5. Conclusion

In this paper we present the results of numerical simulation by a two-dimensional hybrid PIC MTRAP code of high-beta plasma with hot ions in an axisymmetric mirror machine. The injection of a plasma beam results in the displacement of the magnetic field from the region occupied by plasma and in the formation of the so-called diamagnetic bubble. The difference in confinement of ions with opposite sign of the azimuthal component of angular momentum in the diamagnetic bubble leads to plasma spin-up. Such a plasma spin-up was observed when modelling the outflow of a plasma with Maxwell ions from a diamagnetic bubble. The simulations show that the average azimuthal ion velocity is limited due to electron drag in the transition layer at the bubble boundary. A population of fast sloshing ions can efficiently transfer energy to axisymmetric magnetosonic waves via the Landau mechanism if the phase velocity of the oscillations is close to the mean longitudinal velocity of the ions. Such waves were observed in numerical simulations of the flow of plasma with sloshing ions in a homogeneous magnetic field. The most unstable are the oscillations with the maximum wavelength, which is limited by the length of the simulation box. The radial distribution of the azimuthal component of the oscillation electric field has a peak in the region where the density of fast ions reaches a maximum. Such oscillations are similar to the global sound modes which were observed earlier in the GDT device (Skovorodin et al. 2013).