Computer simulation of diamagnetic regime in open magnetic trap

Preliminary results of numerical simulation of transition to regime of diamagnetic confinement in an axially-symmetric mirror trap are presented. Author’s two-dimensional modification of hybrid Particle-in-Cell model is used in the simulations. Formation of region with hot dense plasma and very small magnetic field is demonstrated. Evolution of plasma and magnetic field, as well as structure of the formed region are discussed and compared with simple analytical estimations.


Introduction
Present work is aimed at numerical simulation of newly proposed regime of plasma confinement in mirror trap, so-called diamagnetic confinement [1]. This method relates to confinement in axially-symmetric mirror trap of plasma with extremely high pressure which is equal to pressure of magnetic field of the trap. The magnetic field is fully displaced from the region occupied by the plasma because of plasma diamagnetism so that cluster of hot dense plasma with very weak magnetic field forms inside the trap (so-called diamagnetic "bubble"). The mirror ratio inside the "bubble" formally is extremely small so good particle and energy confinement is expected. The structure of the magnetic field in the diamagnetic trap has much in common with the structure of the field in elongated Field Reversed Configurations (FRCs) with very small field reversal [2]. So many of methods of MHD stabilization (in particular stabilization by conducting shell) verified in the experiments on FRCs can be used in the diamagnetic trap.
The attractive features of the diamagnetic trap (simple structure of magnetic system, stationarity, principal possibility of developing of compact reactor with high density of thermonuclear power) stimulate experimental investigation of the regime of the diamagnetic confinement. The diamagnetic-confinement-like regimes were demonstrated recently at the machine with neutral beam injection (NBI) C-2W [3]. The Compact Axisymmetric Toroid (CAT) device targeted on studying of formation of FRC and diamagnetic "bubble" in mirror trap with powerful NBI is being constructed in Budker Institute of Nuclear Physics [4]. Models describing equilibrium, evolution and stability of plasma with extremely high pressure are necessary for planning and adequate interpretation of results of experiments targeted to demonstration of diamagnetic confinement. At present one-and two-dimensional models describing equilibrium in the diamagnetic trap in the MHD approximation exist [1], [5]. The models predict the particle and the energy confinement time √ τ τ ⊥ , here τ is the particle lifetime in the gas-dynamic regime and τ ⊥ is the time of magnetic field diffusion across the margin of "bubble" [1]. The models also allow calculating the equilibrium magnetic field, as well as the distribution of the plasma density and temperature [5]. At the same time, some assumptions underlying the MHD models cannot be adequate to this task, in particular, the assumption of zero Larmor radius of ions. Therefore, hybrid MHD-kinetic models based on the Vlasov equations for ions and the MHD approximation for magnetized electrons should be developed. Even two-dimensional model without asymmetry on azimuthal angle can help to clarify such questions as prediction of the particles lifetime and the structure of the magnetic field generated by the plasma. The first is important for estimation of reactor perspectives of the diamagnetic trap, the second plays significant role in the investigation of stabilization of MHD instabilities by conducting shell. Also, the simulations can clarify some properties of the trajectories of the ions, for example, the rate of regular and stochastic orbits and the ions spin-up caused by preference losing of ions with positive sign of azimuthal momentum. These aspects of ions dynamic in FRCs are briefly discussed in [6]. Currently, there is a fairly large number of publications devoted to the problems of numerical modelling of nonlinear processes in high-pressure plasma. However, most of the articles are aimed either at the stability of the FRC, or at their formation by pulsed methods or a rotating magnetic field. Articles devoted to computer simulation of the configuration of high-beta plasma in the mirror traps with a powerful injection are virtually absent. It should be noted, in the paper [7] simulation of FRC formation in a mirror machine with powerful off-axis NBI is described.
The Particle-in-Cell method (PIC) was used for simulation of evolution of ion component but the electron currents were not taken into account. Formation of FRC in a mirror trap with offaxis NBI was simulated numerically with hybrid PIC code in paper [8]. Also, there are articles describing simulation of NBI in FRC, see for example [9]. Therefore, the development of a new numerical code for modeling the regime of diamagnetic plasma confinement in open magnetic systems is necessary.
Preliminary and test results of simulations of the regime of the diamagnetic confinement are presented in this article. Author's modification of hybrid PIC model with kinetic description for the ions and hydrodynamic description for the electrons is used in the simulations. The article is organized as follows. The formulation of the problem is discussed in the second section. Some analytical estimations of the injection power necessary for the transition regime of the diamagnetic confinement and the time required to reach steady-state are given in the third section. The numerical code is described in the fourth section. The results of the preliminary simulations are discussed in the fifth section.

Formulation of the problem
To be able to compare the results of our and MHD simulations, we choose formulation of the problem close to that in MHD models [5]. Let's consider a mirror trap 0 < r < R 0 , −L/2 < z < L/2 with axially symmetrical vacuum magnetic field B = {B r , 0, B z }. Ions are continuously injecting into the trap from the point r = 0, z = 0. Under the assumption of quasineutrality of the plasma, the electron density inside the trap is equal to the sum of the ion densities of the background plasma and the injected ions. Here we do not consider specific methods of the matter and energy injection into the trap, for example, injection of a neutral beam, granules, etc. The probability of occurrence of ion-electron pair in the trap does not depend on the longitudinal coordinate (coordinate along the axis of the trap). The velocity distribution function of the electrons is Maxwellian with temperature T e , the injected ions can have arbitrary velocity distribution function. It should be noted, that in the MHD models the velocity distribution There are several typical times: ion cyclotron gyro-period 2π/Ω, here Ω 0 = eB 0 /m i c is ion cyclotron frequency (here B 0 is value of vacuum magnetic field in the center of the trap and m i is the mass of ions), period of bounce-oscillations of ions 2π/Ω b (which is of the order of the ratio of distance between mirrors to typical ion velocity) and times of ion-ion and ionelectron collisions, τ ii and τ d . With typical experimental parameters, there is a hierarchy of times . This hierarchy leads to some obvious difficulties in the simulation of formation of the "bubble"-like structure of the magnetic field. In the present model the ions are modelled by the PIC method, while the electrons are modelled as a liquid with density that coincides with the density of the ions. Typical times in this problem are the times of electron drag and L/v A , here L is the distance between mirrors and v A is the Alfvén velocity. The times are much greater than typical electron times, which are of the order of inverse electron cyclotron Ω −1 e , plasma ω −1 pe and bounce L/(T e /m e ) 1/2 frequencies. Thus using the hybrid model allows the simulation to be essentially accelerated. We assume m e → 0 in the fluid equations. This neglecting implies ω −1 pe → 0 and Ω −1 e → 0, so excitation of electrostatic and electroncyclotron waves is impossible in the simulation. These processes are insufficient for formation of the diamagnetic "bubble". This approach implies neglecting of all electron drifts, except of the electric drift. Strictly speaking, such approximation is correct if the electron temperature is much less than the mean energy of the ions.
3. Analytical estimations in the case of low plasma pressure Let us consider the evolution of ion population in a mirror trap with injected plasma (ion and electron pairs form inside the trap). For simplicity, we limit ourselves to the case when the injection power is small enough and the plasma pressure is small in comparison with the pressure of magnetic field.
The evolution of population of ions satisfies the kinetic equation Here B is magnetic field of the trap, S(r, v) is source of ions, term St ii describes ion-ion collisions and τ d is the time of electron drag. The boundary condition in the mirrors is ∂ v f i = 0 for trapped particles (here v is the component of velocity along magnetic field) and f i = 0 at ±v > 0 for passing particles (here signs "+" and "−" correspond the left and right mirrors). We neglect the density of passing ions further. Let's us neglect the ion-ion collisions now. After integration over speeds and space with the weights of 1 and v 2 here N and W are the total number of trapped particles in the trap and their transversal kinetic energy, J is the number of ions appearing in the trap per unit time, and m i w 2 ⊥ /2 is the mean transversal kinetic energy of appearing trapped ions. So plasma density grows linearly and plasma energy reaches a limit value τ d (m i w 2 ⊥ )J/3 in the time of the order of electron drag time if angle scattering of ions is neglected. So the system reaches a steady-state in a time of order τ d > τ ii . Also these calculations allow parameters of the diamagnetic "bubble" to be estimated through injection parameters: m i w 2 ⊥ J (B 2 0 /8π)La 2 /τ d , here B 0 is the vacuum magnet field in the center, L is the distance between mirrors and a is the radius of the "bubble". Accumulation of ions with low energy is a non-physical effect which is provoked by neglecting of ion-ion collisions. Actually, the angle scattering restricts the density of the trapped particles by value of the order of J/τ ii . Thus neglecting of the scattering does not prevent the formation of the "bubble"-like stationary structure, but it leads to some non-physical effects.

Basic equations
The original system of equations includes the Vlasov equations for the ion component of the plasma and the incoming ions, the equations of magnetic hydrodynamics for the electrons and Maxwell equations [11]: Here, f i is the ion distribution function, B and E are the magnetic and the electric fields, v e is the electron velocity, v i = f i (r, v, t)vdv/n i is the average velocity of ions, p e and T e are the pressure and temperature of the electrons. The collision force R i = −R e takes into account the exchange of pulses between the electron and ion components of the plasma and the beam, It is assumed here that the thermal conductivity coefficient k and the collision frequency ν = 1/τ d are independent of the plasma and magnetic field parameters (ν = const and k = const).
The heat generated in the electrons is Q e = j 2 /σ, the electron heat flux is q e = −k∇T e . In the calculations, the adiabatic index γ = 5/3 is used. To solve the Vlasov kinetic equation, the author's modification of PIC is used [12]. To solve the MHD and Maxwell equations, finitedifference schemes of the first order of accuracy are used. To solve the heat balance equation for the electrons the up-wind scheme is used. In hybrid models, the condition of plasma quasineutrality is assumed to be n e = n i (the scales under consideration are larger than the Debye radius), and the displacement currents are also neglected, only low-frequency processes are considered. The electric field is determined from the equation of motion of the electrons under the assumption m e = 0.
A detailed description of the used numerical methods and algorithms for the problem solution is presented in the report [13].

Numerical results
Results of verification of the numerical code are described in this section. All quantities are normalized on the ion cyclotron frequency Ω 0 = eB 0 /(m i c), the Alfvén velocity v A = B 0 / √ 4πm i n 0 and the ion dispersion size c/ 4πn 0 e 2 /m i , here B 0 is the magnitude of vacuum magnetic field in the center and n 0 is the initial concentration of the background ions. We choose following normalizing parameters: B 0 = 2kG and n 0 = 10 12 cm −3 , so that Ω 0 ≈ 1.9 · 10 7 s −1 , v A ≈ 4.4 · 10 8 cm/s and the ion dispersion size is 22.8cm. The electron drag time is chosen nonphysically large to decrease time necessary for reaching stationary state. The ion-ion collisions are neglected in this simulation. The vacuum magnetic field is created by the two coaxial coils with equal radii and current. The mirror ratio is r mirror = 2. The distribution function of the injected ions is proportional to  Formed configuration of the magnetic field consists of inner space with approximately zero magnitude of the field and outer region with weakly modified field (see figure 2). These areas are separated by a transition region where the magnetic field changes sharply. The width of the transition region is of the order of the average Larmor radius of the beam ions. The radius of the transition region, defined as the radius at which the magnitude of the magnetic field is equal to one percent of B 0 , monotonically grows and reach a stationary value in a time of the order of τ d (see figure 2).  The changing of the magnetic field leads to the generation of the azimuthal component of the electric field. This vortex field decelerates the beam ions so that resulting changing in the kinetic energy of the ions is equal to the changing in the magnetic field energy. The magnitude of the vortex electric field decreases with time and is equal to zero when a stationary state is reached. The azimuthal momentum of the beam ions is close to zero because it was born on the axis. So the density of the beam ions inside the "bubble" decreases approximately inversely with the distance to the axis (see figure 2). The beam ions are reflected from the transition region, so their density is very small outside of the "bubble". The background plasma is extruded from the region occupied by the beam ions (see figure 2). There are two extrusion mechanisms. Firstly, an increase in the beam ion density leads to the formation of an electron density gradient due to quasi-neutrality. The gradient of electron pressure leads to arising potential electric field (see third equation in (3)), which ejects cold background ions along magnetic field lines. This effect is analogical to ejecting of cold ions from an open trap by ambipolar potential. On the other hand, an increase in the pressure of the beam ions leads to a change in the magnetic field and a generation of the azimuthal component of the electric field localized near the edge of the "bubble" (see figure3). This field removes the background ions and electrons in the radial direction due to the electric drift.

Conclusion
Preliminary results of simulations of transition to the regime of diamagnetic plasma confinement in open magnetic trap are presented. Formation of quasi-stationary structure of magnetic field with region of extremely small magnetic field is demonstrated. Typical times of processes and structure of magnetic field are consistent with qualitative conception given in [1] and analytical estimations. The numerical model is applicable for simulation with parameters corresponding to real experimental conditions.