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 (Reference Liseykina, Vshivkov and Kholiyarov2024).

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 Reference Takizuka and Abe1977) 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. Reference Vshivkov, Soloviev, Chernoshtanov and Efimova2021). 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. (Reference Vshivkov, Soloviev, Chernoshtanov and Efimova2021).

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,Footnote ² $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 Reference Beklemishev2016), where ion losses are also caused by the electron drag in the transition layer.Footnote ³

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.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$.