Three-dimensional hybrid numerical tool for collisionless plasma modeling

Recent laser-produced plasmas experiments open up new opportunities for the so-called laboratory astrophysics, allowing observation and studying a number of fundamental physical processes relevant to magnetized plasmas, such as thermo-magnetic instabilities leading to magnetic field generation, magnetic reconnection, collisionless shocks. In order to supplement those experiments with full-scale numerical simulations we develop a code AKA52 (Arbitrary-Kinetic-Algorithm) implementing a hybrid model that includes the dynamics of magnetic fields: advection by the ion flow and Hall effect, magnetic field generation by the Biermann battery effect and Weibel instability. The fully-parallelized high-performance hybrid algorithm includes Particle-in-Cell (PIC) formalism for ions and a 10-moment fluid model for electrons that are described by density, bulk velocity and the six-component pressure tensor evolution equation. Laser-plasma interaction is simulated by means of an ablation operator which imitates laser ionization and heating at critical density surface. As an example, we chose a problem of plasma expansion in the externally applied magnetic field perpendicular to the flow that is related to a number of recent laser-plasma experiments.


Introduction
Laboratory, space and astrophysical plasmas generally introduce multiple scales, both in space and time. The smallest scales are associated with light and fast electrons that require special treatment in full scale simulations. As most of the studied plasmas is magnetized, modeling of realistic the proton to electron mass ratio and the speed of light to the Alfvén speed ratio (both of the order of 10 3 ) require significant computational resources. Due to large scale separation both in time (between electron oscillation time and observation time) and in space (between Debye length and system size), various approximations have to be done to be able to run simulations for reasonable time. The prominent approach is a hybrid [1] that keeps the ions kinetic effects but removes fast electron gyration.
Hybrid codes, in which the heavier ions are treated as macro-particles following the PIC formalism and the electrons are assumed to be a massless neutralizing fluid, have been widely used in space physics over the past forty years. Such codes are used to model phenomena on ion inertial and gyro-radius scales which fall between longer scales obtained by magnetohydrodynamic simulations and shorter scales attainable by full particle simulations. From the first hybrid models applied to laboratory pinch experiments [2], it was subsequently developed for realistic space weather simulations [3]. Our motivation to develop the hybrid approach is coming from a request to model full scale laboratory astrophysics experiments that are developing in parallel with the high power laser facilities. Over the past few decades, a great experimental effort has been done in studying the expansion of laser-produced plasma in an external magnetic field. In the work we consider as an example a problem of plasma flow expansion in the externally applied magnetic field perpendicular to the flow. The choice of a such configuration is inspired by novel experiments in context of laboratory astrophysics [4; 5; 6], laser-driven inertial confinement fusion [7] and industrial applications [8; 9]. Such configurations are well know for being accompanied by rich physics, including many plasma instabilities such as lower hybrid drift instability [10] (LHDI), the magnetic Rayleigh-Taylor instability [11] (MRTI).
In the next section we discuss details of the implemented hybrid model. In Section 3 we overview numerical issues coming from interconnection between fluid and kinetic approaches. Section 4 is dedicated to performance analysis. In section 5 we present an example of modeling a real experiment performed in the Institute of Applied Physics RAS on PEARL laser facility [12].

Numerical model
In the presented hybrid model, electrons are described by a 10-moment fluid model. Those 10 moments are density -n (equal to the total ion density by quasi-neutrality), bulk velocity -V e and the six-component pressure tensor -P e . To properly describe ion-scale instabilities, we keep kinetic description for ions. The hybrid model we use, solves the ion kinetic dynamics using the PIC method. The macro-particle with mass m and charge q obeys the equations: where x and v denote the particle position and velocity, E and B are the electric and magnetic fields interpolated to the particle location. The equation is solved using the Boris algorithm [13] that saves computation time compared to other schemes and is quite accurate. At each time step, the particle moments, namely the density n and the ion bulk velocity V i , are computed in each grid cell using a first order assignment function for each macro-particle. These ion moments are smoothed using a three-points stencil. Such smoothing helps to prevent the growth of smallscale electric fields in low-density areas and has limited consequences for the numerical diffusion processes.
The electromagnetic fields are treated in the low-frequency (Darwin) approximation, as we consider that the phase velocity of electromagnetic fluctuations is small compared to the speed of light. Neglecting the displacement current, we then write an electron Ohm's law : where J is the total current density being equal to the curl of magnetic field ∇ × B, ν is the hyper-viscosity. For the hyper-viscous term, we set ν = 10 −3 , that provides a dissipation process at sub-ion scales but does not affect studied physics. The reason to have such dissipation is coming from the limits of the fluid electron approximation that eliminates electron Landau damping, that in turn, should suppress whistler growth at short wavelengths. This is why we need to use an additional smoothing to reduce the amplitude of the short wavelength fluctuations that helps to prevent numerical instabilities. Electromagnetic fields are calculated on two staggered grid to preserve a time centering of the scheme.
The Ohm's law includes several effects which are important at scales ranging from global to kinetic one. Among them the advection of magnetic field by the ion fluid V i × B describing also the plasma diamagnetic effect, the Hall effect J × B describing sheath fields appearing when electrons are magnetized whereas ions are not, and the electron pressure divergence term ∇. P e . The last one can lead to magnetic field generation due to the Biermann battery effect in regions where the temperature and density gradients are not collinear. Also it contains Weibel type instabilities leading to current filamentation. The numerical issues associated with the inverse density dependence of the Hall and pressure divergence terms are discussed in Section 3.
For electrons description we use the six-component pressure tensor evolution equation: where the electron bulk velocity V e is expressed through the ion bulk velocity and the current density : V e = V i − J/n. Following the accepted terminology [14], the driver term [D] (first line of right hand side of Eq. 4) describes the transport of pressure at the electron fluid velocity, the cyclotron term [C] (second line) describes the rotation of pressure at the electron gyrofrequency and the divergence of electron heat flux (third line) is replaced by isotropization term [I].
As the model treats ion timescales the time step ∆t used in equations is small enough to resolve the gyro-motion of ions. To integrate the cyclotron term being on electron timescales we use an explicit scheme based on sub-cycling method. Defining µ = m i /m e as the ion-to-electron mass ratio, a time step ∆t/µ is appropriate for the electron magnetization treating. As we use a predictor-corrector scheme [1] to preserve a second order, the electron pressure tensor is defined at half time step, to obtain P n+1/2 we use µ sub-cycles with the following algorithm (for l ∈ [0, µ − 1]) In this equation we need to set the driver at even step, the term is extrapolated at the predictor phase and interpolated at the corrector one. The value of izotropization operator is defined using the last calculated pressure value.
In the code it is assumed that density and magnetic field are normalized to their characteristic values n 0 and B 0 respectively, lengths are normalized to the ion inertial length d 0 (calculated with n 0 ), times are normalized to the inverse of ion gyrofrequency Ω −1 0 (calculated with B 0 ) and velocities are normalized to the Alfvén velocity V 0 (calculated with B 0 and n 0 ). The normalization of the other quantities follows from these ones.
Laser effects are imitated by ablation operator, including heat operator and particle creation operator. Ablation operator works in localized area corresponding to the laser focal spot. Heat operator provides pressure increasing in the focal spot. Each time step we pump the driver part of the pressure evolution equation. Ions are accelerated by the arising electron pressure gradient. For electrons, we add a given constant to diagonal pressure tensor components imitating the temperature increase. Magnitude of the added constant delta is adjusted to obtain a desired temperature and speed of flow. The particle creation operator sustains constant target density, mimicking the reservoir of plasma provided by a solid-density target. Particles are added at a specified cold target temperature.
As the most challenging in the interconnection between kinetic and fluid models is a transition from discrete particle properties as position and velocity to continuous fluid moments, in the next section we discuss numerical issues arising while treating the electron pressure tensor evolution equation.

Numerical issues
The classical laboratory astrophysics experiments [4; 5] are performed in high vacuum chambers or at least with background pressure being 5-6 order lower than the atmosphere one. Because of an uncertainty for zero density value in the Ohm's law (Eq. 3), more precisely for the Hall and pressure divergence terms, we cannot simulate an empty space with hybrid codes. We cannot simulate too low density neither as it potentially would increase these terms dramatically. For these reasons such simulations commonly use a background particles population that is uniformly distributed all over the simulation box. That makes the density level high enough to stabilize the integration scheme. For the electron pressure tensor integration (Eq. 4), we need to calculate the electron flow velocity that also depends on the inverse density value. As we are interested in comparison with fully kinetic models, we want to avoid using extra smoothing on pressure tensor components. In the work we do not smooth pressure tensor components while solving the evolution equation (4). However, to achieve numerical stability we smooth particle moments each time step and additionally current that is also used in the electron velocity expression. In practice, the most noisy components while integrating are electron velocity cross derivatives (∂ x V y , ∂ y V x , etc.) in driver part that are susceptible to the electron velocity gradients, that we want to keep smooth.
In the work we introduce an alternative method to avoid the use of high background density values while keeping the run stable. We use a multiplication factor in front of the terms in Ohm's law (Eq. 3) and in the part of electron flow velocity expression depending on inverse electron density. The factor is a function of the electron density that is zero when density is lower than a threshold and one otherwise. In between, it is a fifth order growing polynomial [15]. If the density is small at given points, it means we have an epsilon in front of the distribution function f (x, u, t), the same epsilon affects the higher order moments. For simulations in section 5 we use the high threshold 0.1 while the lowest is 10 −3 . The value 0.1 is a compromise for hybrid simulations, as the lower values could lead to the dramatic growth of the associated terms in the Ohm's law, as well as the electron velocity and, as a consequence, its gradients being important

Performance
The code AKA52 is written in C++11, for simulation tasks initialization we use a single script file written in python2. The problem is distributed among the nodes through MPI. To measure a strong scalability as a test problem we have chosen oscillations of the uniformly distributed demagnetized plasmas and used 128/256/512/1024 cores. The simulation box has 100×100×100 grid and we use 100 particles per cell number with total 10 8 particles in the box. We performed the experiments on the resources of Joint Supercomputer Center of the Russian Academy of Sciences (Intel Xeon E5-2697Av4, 128 GB RAM per node).
On Fig. 1 we can see that time per step (time counts both phases predictor and corrector one) decreases as 1/n, where n is a number of CPU cores. For the discussed in the next section simulations using an unrealistic ion-to-electron mass ratio µ = 10, the presented implementation spends 80% of the computational time for the particles pushing and gathering fluid moments, the rest for treating electron pressure tensor evolution and electromagnetic fields calculation including also MPI communication. Increasing ion-to-electron mass ratio up to realistic values µ = 2000, increases the part for electron pressure tensor integration (Eq. 4) up to 50 % of the total computational time.

Plasma expansion in magnetic field
In this section we analyze a simulation of plasma expansion in the externally applied magnetic field perpendicular to the flow. The chosen simulation parameters are close to ones from the experiment on PEARL laser facility, the results of which have been discussed in context of three-dimensional magnetohydrodynamic simulations of the similar experiment on TITAN and ELFIE laser facilities [12]. The comparison of parameters in the experiment and in our model is given in the Table. 1. In our model we use the characteristic magnetic field B 0 = 100 Tesla as in laser plasma experiments the Biermann magnetic field generation sets its level on the Megagauss one [16], being stronger than externally applied field. The external magnetic field then equals to 0.135B 0 and is at X direction. We use a near critical density value for n 0 . Also for simplicity we use ions mass and charge, both equal to 1. The diameter of the focal spot D spot , where we apply the ablation operator, is 20d 0 . The operator works during all the simulation imitating the target cooling.
Firstly, we look at the classical simulation case using background population. In the model initially we set the dense target density to 1 and background population density to 0.2. The simulation domain is 3D rectangular box with sizes L X = 100, L Y = 100, L Z = 125. We use a 125×125×150 grid corresponding to a mesh size equal to 0.8d 0 in all directions. The time-step is 5 × 10 −3 in order to satisfy the CFL conditions for the fastest whistler modes. We use periodic boundary conditions in X and Y directions and perfect conducting boundaries in Z direction. We use 50 particles per cell. For the cyclotron term integration in the pressure tensor evolution equation we use a reduced ion-to-electron mass ratio µ = 10 while keeping the relaxation parameter τ on electron timescales (τ = 0.01).
In experiments [17; 18; 12] analysed with the three-dimensional MHD simulations, plasma evolution follows three distinct phases. In the beginning plasma expands as a hemisphere. Therefore, the plasma initially expands as a diamagnetic cavity, canceling the magnetic field in the plume interior. The second phase is observed when magnetic field pressure starts to be comparable with plasma thermal pressure, the plasma plume expansion decelerates. Also at this phase small-scale structures on the plume front start to appear that are supposed to be seeds for further evolution. One of the seeds leads to the formation of a narrow directed slab which travels across the magnetic field. As was concluded [12] morphology of the slab is well reproducible and does not depend on small variations of the laser-target interaction from facility to facility. While  the small scale filaments in perpendicular to slab direction depend on the external magnetic field strength.
Our model well captures the diamagnetic cavity, Fig. 2 (a,b) shows B x in two cross sections. We see that magnetic field is pressed out the plume interior and compressed on the front. The general picture in perpendicular to external field cross section is asymmetric, plasma is drifting in Y direction. While the interior magnetic field is quite chaotic, we clearly recognize the Biermann battery magnetic field in the focal spot region (d) where pressure gradient is in the target plane and density gradient is at Z direction, non-collinear to each other. Also Weibel type instabilities are growing on the plume front (c) forming filaments. Both effects were also presented in 2D fully kinetic simulations [19] where pressure tensor is calculated from kinetic electrons.
In contrast to the recent 3D magnetohydrodynamics simulations [12] we do not see growth of small scale filament on the plume surface, one of the reason could be a dense background population. To circumvent this discrepancy we use the previously introduced technique of density factor multiplication. For such runs we use background density level 10 −6 that is comparable with experimental conditions. The main difference is an origin of the shock wave on the electron pressure. Fig. 3 shows pressure tensor component for cases with dense background (a) and without (b). Instead of an uniformly compressed cavity we see a thin sheath on the front. We still do not observe a formation of the narrow plasma slab that could be a consequence of the simulation domain size which doesn't allow to resolve well the focal spot diameter. For future work we plan to use an extended box and greater spot diameter that requires significant increase of the computational resources and is out of scope of this paper. The presented simulations have been performed using an unrealistic ion-to-electron mass ratio for the pressure tensor evolution equation integration. Normally, increasing the ratio also requires to decrease the spatial resolution, as the electron spatial scales are √ µ times smaller than the proton one. Because of the quadratic dependence of the time step on the cell size, it requires also to decrease significantly the time step. As the features we would capture were on ion scales, such a choice allowed us to use a reduced spatial and time resolution that saved computational resources. For the future work the choice of the ratio requires further analysis depending on the investigated problem and the observed gradients.

Conclusion
We have developed a 3D fully-parallelized high-performance hybrid code for collisionless plasmas modeling. We tested the performance on oscillations of the uniformly distributed plasmas and observed a quite good strong scalability increasing CPU cores number. Also we have presented the results of three-dimensional hybrid simulations of experiments designed to study how the external magnetic field applied perpendicular to plasma flow affects its dynamics. The observed structure consists of a plasma plume with a magnetic cavity. A shock envelope surrounds the magnetic cavity. To observe fine structures on the pressure tensor components, we introduce new technique that helps to reduce the density of background population.