MHD-PIC Supercomputer Simulation of Plasma Injection into Open Magnetic Trap

This paper presents a two-dimensional hybrid Magneto-Hydro-Dynamical-Particle-in-Cell numerical model to study the interaction of a beam of plasma, injected into an axisymmetric magnetic trap, with the background trap plasma. We apply a kinetic description for the positively charged ions, treat electrons as a massless charge neutralizing ﬂuid and assume that the direct coupling be-tween ions and electrons is due to the anomalous scattering on the ﬂuctuations of electromagnetic ﬁelds only. The model adequately describes nonlinear nonstationary evolution of the plasma and of the magnetic ﬁeld and allows to follow this evolution for large simulation times in a wide range of the initial magnetic ﬁeld and plasma parameters. We show in particular, that the continuous injection of plasma beam leads to the displacement of the magnetic ﬁeld and to the formation and growth of extended region with low amplitude of the ﬁeld. The numerical results demonstrate the accumulation and capture of the plasma in the magnetic cavity region.


Introduction
Open (open-ended) magnetic trap is one possible solution to the problem of plasma confinement and heating in laboratory experiments [11].The idea of such a facility was proposed in the 50 s by G.I. Budker [3] and R.F.Post [10] independently, and it was tested for the first time then.It uses the gradient of the magnetic field to trap particles, whereby the confinement of the particle is due to the adiabatic invariance of its magnetic moment, which takes place when the Larmor radius of the particle is small compared to the scale of change of the magnetic field.Although open traps can be operated under steady-state conditions and are attractive from the engineering point of view, they have an important drawback, namely the relatively high losses of plasma along the magnetic field lines.Several proposals of improved traps, that are largely free from this drawback, have been made in the last decades.The diamagnetic confinement is one of those.The main idea behind diamagnetic confinement is to suppress longitudinal losses from an axisymmetric open trap by creating a configuration with an extremely high plasma pressure equal to the magnetic field pressure.Diamagnetic confinement was proposed and theoretically justified in [1], where the magneto-hydro-dynamical (MHD) approximation was used to describe the equilibrium of a plasma with β = 1, with β being the ratio of the plasma pressure to the magnetic pressure.In particular, it has been shown that a monotonous increase of the plasma pressure in an open trap leads to the formation of a plasma filled region with a displaced magnetic field, the so-called diamagnetic "bubble".Further increase of plasma pressure results in an increase of the "bubble" radius.The plasma confinement time grows proportionally to the "bubble" radius and can significantly exceed the confinement time in a vacuum magnetic field.Furthermore, the diamagnetic confinement regime with intense off-axis injection of atomic beams at an angle to the magnetic field may additionally suppress losses of the background plasma, raise the temperature of electrons and the lifetime of fast ions arising from the trapping of injected atomic beams by the background plasma [5].In this paper we present a hybrid MHD-Particle-in-Cell (PIC) numerical model to study the formation of the diamagnetic mode in an axisymmetric open magnetic plasma trap with continuous injection of a plasma beam.The paper is organized as follows.In Section 1, we describe the main aspects of our implementation of the hybrid model, including parallelization strategy 3 .In Section 2, we present typical results from MHD-PIC simulations of the interaction of a plasma beam injected into a cylindrical trap with the background trap plasma.

Numerical Algorithms and Methods
In our hybrid model, we use a kinetic description for the positively charged ions and treat electrons as a massless charge neutralizing fluid.This approach is justified because the relevant time and space scales are determined by the ions.For simplicity in the description of our model we consider the case of hydrogen plasma.The evolution of the ion distribution function f ( r, v, t) is governed by the Vlasov equation while the motion of the electron fluid follows Here m i , m e and q i , q e = −q i are the ion and electron mass and charge, E and B are the electric and magnetic fields, n e is number density of the electrons.We take into account the resistive coupling between electrons and ions via R = νm e ( V i − V e ), with ν being the electron-ion collision frequency which does not depend on the plasma and magnetic field parameters4 , V e , V i = vf d v f d v denote the electron and ion mean velocities.Finally, we assume that the electron pressure p e is scalar, and that the quasineutrality condition n e = n i = f d v is fulfilled.Since we consider only low-frequency processes and do not take into account the displacement current, the total current density J = q i n i ( V i − V e ) is obtained from Amperes law and the magnetic field evolves according to Faradays law The electric field is calculated from the electron momentum equation ( 2) Finally, the equation for the electron temperature reads where Q e = J 2 σ is the heat generated in electrons, κ∇T e is the electronic heat flux, with κ being the thermal conductivity, σ = e 2 n e m e ν -electric conductivity and γ = 5 3 -the adiabatic index of an atomic gas.In our simulations lengths are given in units of c/ω pi , with ω pi = 4πn 0 e 2 /m i being the ion plasma frequency, and time in units of 1/ω iH , where ω iH = eB 0 /(m i c) is the ion cyclotron frequency with B 0 being the initial amplitude of the magnetic field in the center of the trap and n 0 is the initial background plasma density.
Our code employs state-of-the-art, widely used numerical algorithms and methods.In particular, the particle-in-cell (PIC) method, see e.g.[8], is used to solve the Vlasov equation for the ion components of the background plasma and of the injected beam 5 .The first-order weighting is applied to interpolate the fields to the ion positions and to obtain ion current and charge density from the positions of the ion-particles relative to the grid points.The Boris [2] pusher is used to propagate ion positions and velocities.The equations of the evolution of the electromagnetic fields and of the electron temperature are discretized using 2 nd -order regular finite difference stencils.More details about the implemented numerical schemes can be found in [7,12].

Initialization
At t = 0 the cold uniform hydrogen plasma is located inside a cylindrical chamber of radius R 0 and length L with a magnetic field 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 ϕ, and that its amplitude does not depend on the azimuthal angle.In this case at t = 0 the whole system has an axial symmetry.A beam of neutral hydrogen plasma with finite temperature is injected with constant velocity | v| into the chamber at the point r = 0, z = 0.The initial angular distribution of the velocity of the beam particles is randomized.In this paper we consider the reduced 2D 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] × [0, R 0 ], Fig. 1a.The distribution of the initial magnetic field along the trap axis is shown in Fig. 1b.The amplitude of the field in the center of the trap at the point r = 0, z = 0 is equal to B 0 , the fields B 1 at the points r = 0, z = ±L/2 are higher6 .

Parallelization
The computational effort in PIC simulations scales with the number of simulation particles.For large numerical boxes, especially in two-dimensional (2D) or three-dimensional (3D) simulations, this number can be as high as 10 9 ÷ 10 12 .This makes the large simulations unfeasible on a desktop computer.In order to use the parallel computing capabilities of modern supercomputers, an efficient parallelization of the numerical code is required.Our algorithms are local and allow parallelization via domain decomposition.In particular, we divide the simulation domain (a) Background magnetic field is created by two identical current coils, located at both ends of the chamber.Black lines are the magnetic field lines (b) The nonuniform magnetic field of a simple pair of current coils forms two magnetic mirrors between which a plasma can be trapped [4] Figure 1.Simulation setup and the distribution of the magnetic field along the trap axis into subdomains in z-direction.Each subdomain is assigned to a group of processor cores and the particles belonging to the subdomain are distributed among the cores of the group.The maximal number of groups is defined by the grid size and should not exceed N z /4, where N z is the number of grid cells in z-direction.At the initial stage, the background particles and the particles of the injected beams are distributed evenly between the cores of their group.While the particle positions and velocities can be updated on each subdomain independently, communication between neighboring domains is necessary.The particles crossing the domain boundaries are marked during the update of the positions, collected and then sent (all within a single MPI SEND call) to the respective neighboring domains.Additional communication routines are invoked after depositing the charge and current density on the grid to obtain the values at the boundaries.The same is done when calculating electromagnetic fields and electron temperature.

Load Balancing
From a computational point of view, numerical routines, handling particle-grid connection, are the most time-consuming.For this reason, the parallelization is only efficient if the number of particles on each core is kept similar during simulation.The large load imbalance, i.e. when some processes contain a much bigger number of simulation particles, will result in a longer duration of the particle propagation step for those processes and the other will be waiting.In this way, valuable computational resources are wasted.To ensure load balancing in our algorithm, we compute the average number of particles N k in one core per every k th group every few thousand time steps, increase on the basis of N k the number of cores in a group in dense regions (with larger N k ) to balance the number of particles per core, and redistribute the particles within the new group.More details about the particular implementation can be found in [7].

Simulation Results
In this Section, we present typical results from MHD-PIC simulations of the interaction of a proton beam injected into a cylindrical trap with the background trap plasma.The simulations were performed with a temporal resolution of ∆t = 10 −5 ω −1 iH and cell size of ∆z × ∆r = [0.05× 0.05](c/ω pi ) 2 .The size of the simulation domain was L × R 0 = [12 × 4](c/ω pi ) 2 with ∼ 2 × 10 5 background particles.The injection of a plasma beam normal to the z-axis takes place at a rate of 10 3 particles per time unit.The background magnetic field in the trap is created by a pair of current coils of radius 4.3 c/ω pi , placed at z = ±L/2 = ±6 c/ω pi .The amplitude of the current in the coils is such that the mirror ratio of the trap = 2, see Fig. 1b).
Note, that on the axis of the trap the magnetic field has only z-component, B z .All numerical parameters have been checked for convergence.As for the physical parameters, they are chosen to be close to the parameters of laboratory experiments on the CAT installation at the Budker Institute of Nuclear Physics of SB RAS (BINP SB RAS, Novosibirsk, Russia) [6].Namely, we consider an open magnetic trap with the length L = 60 cm and the radius R 0 = 15 cm.The number density of the background plasma n 0 = 4 • 10 13 cm −3 , the magnetic field strength in the center of the trap B 0 = 2 kG.The simulation was performed on 30 groups of cores of the computer facility of the Siberian Supercomputer Center of SB RAS (SSCC SB RAS, Novosibirsk, Russia).
Figure 2 shows the snapshots of the distribution of the magnetic pressure inside the trap at different times.The injection of plasma beam leads to the displacement of the magnetic field and to the formation and growth of the magnetic cavity, i.e. of the region with low amplitude of the field.The values of the magnetic field pressure inside the cavity are very small, < 5% of the initial value, and at the sharp boundaries of the cavity the amplitude of the magnetic field is in turn high.In order to investigate the temporal evolution of the cavity we plot in Fig. 3a the distribution of the amplitude of the magnetic field |B(r, t)| z=0 | in the central plane z = 0.The magnetic field in Fig. 3a is measured in units of B 0 .The size of the cavity in the radial direction grows over time, but this growth has an oscillatory character.The simulations with different mirror ratios R m have shown that the period of the oscillations is not constant but depends monotonously on the increasing radius of the cavity.Further investigation is needed to confirm this hypothesis and find the exact dependence of the oscillation period on the cavity radius if any, or to deny it.The decreasing rate of growth of the cavity radius with time implies that it is possible to reach a quasistationary regime in which the transverse size of the cavity will remain almost constant over time.Snapshots of the distribution of the injected ions displayed in Fig. 3b show the accumulation and capture of the plasma in the magnetic cavity region.

Conclusion
In this paper, we present a MHD-PIC hybrid numerical model to study the interaction of a beam of plasma injected into an axisymmetric magnetic trap with the background trap plasma.The background magnetic field is created by two identical current coils, located at both ends of the trap chamber.We assume that the field configuration and the evolution of the plasma do not depend on the azimuthal angle through the whole process and consider the reduced 2D cylindrical problem, with z-axis directed along the symmetry axis of the trap.We apply a kinetic description for the positively charged ions and treat electrons as a massless charge neutralizing fluid.The parallel numerical algorithm is based on the domain and particle decomposition.The model adequately describes nonlinear nonstationary evolution of the plasma and of the magnetic field and allows to follow this evolution for large simulation times in a wide range of the initial magnetic field and plasma parameters, as applied to the conditions of laboratory experiments at the BINP SB RAS.Our simulations show that the continuous injection of plasma beam into the trap leads to the formation of an extended region with the displaced magnetic field, the magnetic cavity.The value of the field pressure inside this cavity is less than 5% of the initial value.The space-temporal evolution of the cavity evidences that the growth of its size in the radial direction has an oscillatory character.Moreover, the simulations with different mirror ratios suggest that the period of the oscillations may monotonously depend on the increasing radius of the cavity.In addition, the decreasing growth rate of the cavity radius with time implies that it is possible to reach a quasistationary regime in which the transverse size of the cavity will remain constant or further change only insignificantly.

Figure 2 .Figure 3 .
Figure 2. Snapshots of the magnetic field pressure at t = 0, t = 24 ω −1 iH and t = 240 ω −1 iH shows the formation and growth of the magnetic cavity.Lengths are measured in units of c/ω pi