A full electromagnetic Particle-In-Cell code to model collisionless plasmas in magnetic traps

TL;DR

This paper presents a relativistic Particle-in-Cell code for simulating collisionless plasma behavior in magnetic traps under electron cyclotron resonance, demonstrating its application to plasma heating in a minimum-B trap configuration.

Contribution

The work introduces a specialized PIC code tailored for ECR plasma heating in magnetic traps, combining analytical and numerical methods for complex plasma physics problems.

Findings

Successful simulation of plasma heating in a minimum-B trap

Identification of electron energy groups during ECR heating

Validation of the PIC code for collisionless plasma modeling

Abstract

A lot of plasma physics problems are not amenable to exact solutions due to many reasons. It is worth mentioning among them, for example, nonlinearity of the motion equations, variable coefficients or non lineal conditions on known or unknown borders. To solve these problems, different types of approximations which are combinations of analytical and numerical simulation methods are put into practice. The problem of plasma behavior in numerous varieties of a minimum-B magnetic trap where the plasma is heated under electron cyclotron resonance (ECR) conditions is the subject of numerical simulation studies. At present, the ECR minimum-B trap forms the principal part of the multicharge ion sources. In this work, a scheme of the relativistic Particle-in-Cell (PIC) code elaborated for an ECR plasma heating study in minimum-B traps is presented. For a PIC numerical simulation, the code is applied to an ECR plasma confined in a minimum-B trap formed by two current coils generating a mirror magnetic configuration and a hexapole permanent magnetic bars to suppress the MHD instabilities. The plasma is maintained in a cylindrical chamber excited at $TE_{111}$ mode by $2.45$ $GHz$ microwave power. In the obtained magnetostatic field, the ECR conditions are fulfilled on a closed surface of ellipsoidal type. Initially, a Maxwellian homogeneous plasma from ionic temperature of $2$ $eV$ being during $81.62$ $ns$, that correspond to $200$ cycles of microwaves with an amplitude in the electric field of $1$ $kV/cm$ is heated. The electron population can be divided conditionally into a cold group of energies smaller than $0.2$ $keV$, a warm group whose energies are in a range of $3-10$ $keV$ and hot electrons whose energies are found higher than $10$ $keV$.

Equations

$$\frac{\partial f_{\alpha }}{\partial t}+\vec{v}\cdot {\nabla }_{\vec{r}}{f}_{\alpha }+\frac{q_{\alpha }}{m_{\alpha }}(\vec{E}+\vec{v}\times \vec{B})\cdot {\nabla }_{\vec{v}}\frac{f_{\alpha }}{\gamma }=0\text{mbox}\text{mbox}$$

$$\vec{\nabla }\times {\vec{E}}^{sc}=-\frac{\partial {\vec{B}}^{sc}}{\partial t},$$

$$\rho (\vec{r},t)=\sum _{\alpha }{q}_{\alpha }\int \gamma f_{\alpha }d\vec{v},$$

$$f_{\alpha }(\vec{r},\vec{v},t)=\sum _p{f}_p(\vec{r},\vec{v},t)$$

$$f_{p}(\vec{r},\vec{v},t)=N_{p}S_{\vec{r}}\big{(}\vec{r}-\vec{r}_{p}(t)\big{)}\delta\big{(}\vec{v}-\vec{v}_{p}(t)\big{)}$$

$$S_{\vec{r}}\big{(}\vec{r}-\vec{r}_{p}(t)\big{)}=\prod_{j=1}^{3}b_{0}\Bigg{(}\frac{x_{j}-x_{jp}(t)}{\Delta x_{j}}\Bigg{)}$$

$$b_0(\xi )=\{\begin{array}{ll}1& \text{mbox}if\text{arrowvert}\xi \text{arrowvert}<1/2\\ 0& \text{mbox}otherwise\end{array}$$

$$\frac{d{\vec{r}}_p}{dt}={\vec{v}}_p,$$

$${\vec{E}}_p=\sum _g{\vec{E}}_{g}W({\vec{r}}_g-{\vec{r}}_p)d\vec{r}$$

$${\vec{B}}_p=\sum _g{\vec{B}}_{g}W({\vec{r}}_g-{\vec{r}}_p)d\vec{r}$$

$$W\big{(}\vec{r}_{g}-\vec{r}_{p}\big{)}=\prod_{j=1}^{3}b_{1}\Bigg{(}\frac{x_{jg}-x_{jp}(t)}{\Delta x_{j}}\Bigg{)}$$

$$b_1(\xi )=\{\begin{array}{ll}1-\text{arrowvert}\xi \text{arrowvert}& \text{mbox}if\text{arrowvert}\xi \text{arrowvert}<1\\ 0& \text{mbox}otherwise\end{array}$$

$$\text{centering}\vec{\nabla }\times {\vec{H}}^{sc\ast }={\vec{J}}^{\ast }+j\omega {\epsilon }_0\overline{\overline{S}}{\vec{E}}^{sc\ast }\text{@add@centering}$$

$$\text{centering}\vec{\nabla }\times {\vec{E}}^{sc\ast }=-j\omega {\mu }_0\overline{\overline{S}}{\vec{H}}^{sc\ast }\text{@add@centering}$$

$$\overline{\overline{S}}=\left[\begin{array}{ccc}{s}_x^{-1}{s}_y{s}_z& 0& 0\\ 0& s_x{s}_y^{-1}{s}_z& 0\\ 0& 0& s_x{s}_y{s}_z^{-1}\end{array}\right],$$

$$s_x=k_x+\frac{{\sigma }_x}{j\omega {\epsilon }_0};s_y=k_y+\frac{{\sigma }_y}{j\omega {\epsilon }_0};s_z=k_z+\frac{{\sigma }_z}{j\omega {\epsilon }_0}$$

$$D_x^{sc}={\epsilon }_0\frac{s_z}{s_x}{E}_x^{sc},D_y={\epsilon }_0\frac{s_x}{s_y}{E}_y^{sc},D_z^{sc}={\epsilon }_0\frac{s_y}{s_z}{E}_z^{sc}$$

$$E_x^{n+1}(i+1/2,j,k)=\frac{\left[1-\frac{{\sigma }_z(k)\Delta t}{2{\epsilon }_0{k}_z(k)}\right]}{\left[1+\frac{{\sigma }_z(k)\Delta t}{2{\epsilon }_0{k}_z(k)}\right]}{E}_x^n(i+1/2,j,k)$$

$$+\frac{k_x(i+1/2)}{{\epsilon }_0{k}_z(k)}\frac{\left[1-\frac{{\sigma }_x(i+1/2)\Delta t}{2{\epsilon }_0{k}_x(i+1/2)}\right]}{\left[1+\frac{{\sigma }_z(k)\Delta t}{2{\epsilon }_0{k}_z(k)}\right]}{D}_x^{n+1}(i+1/2,j,k)$$

$$-\frac{k_x(i+1/2)}{{\epsilon }_0{k}_z(k)}\frac{\left[1-\frac{{\sigma }_x(i+1/2)\Delta t}{2{\epsilon }_0{k}_x(i+1/2)}\right]}{\left[1+\frac{{\sigma }_z(k)\Delta t}{2{\epsilon }_0{k}_z(k)}\right]}{D}_x^n(i+1/2,j,k)$$

$$D_x^{n+1}(i+1/2,j,k)=\frac{\left[1-\frac{{\sigma }_y(j)\Delta t}{2{\epsilon }_0{k}_y(i)}\right]}{\left[1+\frac{{\sigma }_y(j)\Delta t}{2{\epsilon }_0{k}_y(j)}\right]}{D}_x^n(i+1/2,j,k)+\frac{\Delta t}{k_y(j)\left[1+\frac{{\sigma }_y(j)\Delta t}{2{\epsilon }_0{k}_y(j)}\right]}\cdot$$

$$[\frac{H_z^{n+1/2}(i+1/2,j+1/2,k)-H_z^{n+1/2}(i+1/2,j-1/2,k)}{\Delta y}$$

$$-\frac{H_y^{n+1/2}(i+1/2,j,k+1/2)-H_y^{n+1/2}(i+1/2,j,k-1/2)}{\Delta z}]$$

$$-J^{n+1/2}(i+1/2,j,k)$$

$$B_x^{sc}={\mu }_0\frac{s_z}{s_x}{H}_x^{sc},B_y^{sc}={\mu }_0\frac{s_x}{s_y}{H}_y^{sc},B_z^{sc}={\mu }_0\frac{s_y}{s_z}{H}_z^{sc},$$

Full text

There are different numerical methods to model plasmas. Depending of both temperature and concentration, these can be classified in three main groups: fluid models, kinetic models and hybrid models [1]. The fluid models are the most simple way to describe the plasma from macroscopic quantities, which are used for the study of highly collisional plasmas where the mean free path is much smaller than size of plasma ($l_{mfp}\ll L$). The kinetic models are the most fundamental way to describe plasmas through the distribution function in phase-space for each particle specie; which are used for the study of weakly collisional ($l_{mfp}\sim L$) or collisionless plasmas ($l_{mfp}\gg L$) from the solution of the Boltzmann or Vlasov equation, respectively [2]. For kinetic simulations there are different method to solve the Boltzmann or Vlasov equation, being the Particle-In-Cell (PIC) codes one the most popular [3, 4, 5] The hybrid model combine both the fluid and kinetic models, treating some components of the system as a fluid, and others kinetically; which are used for the study of plasmas, may use the PIC method for the kinetic treatment of some species, while other species (that are Maxwellian) are simulated with a fluid model.

In this work, a scheme of the relativistic Particle-in-Cell (PIC) code elaborated for an ECR plasma heating study in minimum-B traps is presented. For a PIC numerical simulation, the code is applied to an ECR plasma confined in a minimum-B trap formed by two current coils generating a mirror magnetic configuration and a hexapole permanent magnetic bars to suppress the MHD instabilities. The plasma is maintained in a cylindrical chamber excited at $TE_{111}$ mode by $2.45$ $GHz$ microwave power. In the obtained magnetostatic field, the ECR conditions are fulfilled on a closed surface of ellipsoidal type. Initially, a Maxwellian homogeneous plasma from ionic temperature of $2$ $eV$ being during $81.62$ $ns$, that correspond to $200$ cycles of microwaves with an amplitude in the electric field of $1$ $kV/cm$ is heated. The electron population can be divided conditionally into a cold group of energies smaller than $0.2$ $keV$, a warm group whose energies are in a range of $3-10$ $keV$ and hot electrons whose energies are found higher than $10$ $keV$.

K****eywords Electron Cyclotron Resonance (ECR), minimum-B magnetic trap, Particle-In-Cell (PIC) method.

1 Introduction

Actually, an electron cyclotron resonance ion source (ECRIS) is widely used for multicharge ions production because in the ECR conditions the plasma electrons can be energized by microwaves without any trouble [6]. In these systems, a minimum-B magnetic trap formed by the superposition of an axial magnetic field produced by two or more axisymmetric D.C. current coils and a radial magnetic field formed by a cusp multipole system is harnessed for plasma confining. Impossibility to apply the analytical methods for the ECRIS plasma study and practical inaccessibility of those plasmas for experimental diagnostics tools compel investigators to elaborate numerical methods so that they could obtain information about the behavior of both the plasma and individual particles in ECRIS [7, 8, 9].

In this work, a numerical scheme developed as a relativistic full electromagnetic PIC code for plasma heating by microwaves under the ECR conditions in minimum-B traps is offered. In the PIC method chosen for our numerical simulations due to uncollisionality of the ECRIS plasma, the evolution of the plasma components is described by using superparticles (SP) [10, 11, 12]. The simulations are carried out in three stages:

Simulation of the stationary microwave field which is excited in the chamber through a port made in its cylindrical surface [13]. 2. 2.

Fill the chamber with a homogeneous plasma ball in the instant $t=t_{0}$ and a Maxwellian distribution for electrons and ions velocities in such way that $\vec{J}(\vec{r},t=t_{0})=0$ in order to not violate the continuity equation in the instant $t=t_{0}$. 3. 3.

Numerical study of the SP self-consistent dynamics.

The chamber of $4.54$ $cm$ in diameter and $10$ $cm$ long is excited at the electric transverse $TE_{111}$ mode by $2.45$ $GHz$ microwaves of $1$ $kV/cm$ in amplitude. The initial plasma is assumed maxwellian of a electronic temperature of $5.44\times 10^{-4}$ $eV$ and a density of $2\times 10^{10}$ $cm^{-3}$.

The system can be described in the framework of the Vlasov-Maxwell equation and solved by using the particle-in-cell (PIC) method [14, 15, 16].

This paper is organized as follows:

Section 2 describes the physical model and an electromagnetic PIC approach to the simulation of the plasma-microwave field interaction. In Section 3, the results concerning the plasma space distribution and the energy spectrum of electrons are depicted.

2 Physical system and numerical model

2.1 Physical scheme model

A physical model of the system under study is shown in Fig.1. The ECR plasma in the non-magnetic metal chamber (1) is confined by the minimum-B magnetic field formed by two current coils (7) and a cusp hexapole (6). The microwaves generated by a magnetron (5) are injected into a chamber (1) through waveguides (4) and (2).

2.2 Numerical model

2.2.1 Electromagnetic particle-in-cell method

The self-consistent simulation of the plasma-microwave interaction in our system is described in the framework of the Vlasov-Maxwell equation system:

[TABLE]

where $\gamma=[1-(v_{p}/c)^{2}]^{-1/2}$ is the relativistic factor, $f_{\alpha}(\vec{r},\vec{v},t)$ is a six-dimentional phase space distribution for a specie $\alpha$, electron or ion in our case, whose charge and mass are $q_{\alpha}$ and $m_{\alpha}$, respectively.

$\vec{E}$ and $\vec{B}$ are the total electric and magnetic fields, respectively. In our system, $\vec{E}=\vec{E}^{sc}$ correspond to the self-consistent electric field, which is the superposition of the microwave field $E^{hf}$ and the plasma self-generated electric field $E^{sg}$.

$\vec{B}=\vec{B}^{s}+\vec{B}^{sc}$, where $\vec{B}^{s}$ is the magnetostatic field of the minimum-B magnetic trap is calculated through a procedure described in [17]. $\vec{B}^{sc}$ is the self-consistent magnetic field produced due to the plasma particle motions. which has a similar meaning of the electric field component previously described. The evolution of the $\vec{E}^{sc}$ and $\vec{B}^{sc}$ are given by the following Maxwell equations:

[TABLE]

The charge density and current density are determined by the expressions:

[TABLE]

respectively.

In PIC simulations, the distribution functions of each specie is given by the superposition of SPs which are the clouds of real particles [14]:

[TABLE]

In our simulation the state of Np-particle ensemble at time t is specified by the exact one-SP distribution function:

[TABLE]

Here $N_{p}$ is the number of real particles that form a SP; \delta\big{(}\vec{v}-\vec{v}_{p}(t)\big{)} is the Dirac delta function for the velocities distribution, being $\vec{v}_{p}$ the velocity distribution. The spatial shape factor is given by,

[TABLE]

where the $j$-index refer to the $x$, $y$ and $z$ coordinates, $\vec{r}_{p}$ is the SP position; $\Delta x_{j}$ is the spatial step in $j$ direction, and

[TABLE]

corresponds to the first b-spline.

The SP dynamics determined by the relativistic Newton-Lorentz equation is solved numerically through the Boris leapfrog procedure [18]:

[TABLE]

where $\gamma=[1-(v_{p}/c)^{2}]^{-1/2}$ is the relativistic factor, $q_{s}=N_{p}q_{\alpha}$ and $m_{s}=N_{p}m_{\alpha}$ are the charge and mass of a SP. Such equations are very similar to those describing the motion of the real particles but the fields $\vec{E}_{p}$ and $\vec{B}_{p}$ are calculated as an average over the SP to preserve the moments of the Vlasov equation. These fields can be calculated from the values calculated in gridpoints, $\vec{E}_{g}$ and $\vec{B}_{g}$, as

[TABLE]

and

[TABLE]

respectively. The interpolation function $W$ depends on the choosen spatial shape factor. For the present case

[TABLE]

where

[TABLE]

corresponds to the second b-spline.

In our numerical model, it is supposed that the considered chamber is made of a perfect electric conductor (PEC) and in order to avoid nonphysical reflections at the microwave port, the Uniaxial Perfectly Matched Layer (UPML) method is used (See Fig.2).

The Maxwell rotational equations (2) at UPML conditions are expressed in a phasor-like form [18]:

[TABLE]

and

[TABLE]

where $\overline{\overline{S}}$ is a diagonal tensor defined by

[TABLE]

where:

[TABLE]

In our simulations $k_{x}$, $k_{y}$ y $k_{z}$ are chosen equal to the unity and the electrical conductivities $\sigma_{x},\sigma_{y},\sigma_{z}$ are chosen equal to zero in all space except the microwave port where $\sigma_{x}=(x/d)^{3}\sigma_{x}^{max}$, $d=10\Delta x$ is the width of the UPML with the spatial step $\Delta x$ in $x$ direction, $x$ is the distance from the border towards the interior of the UPML and $\sigma_{x}^{max}=0.8(m+1)/(\Delta x\sqrt{\mu_{0}/\varepsilon_{0}})$ with $m=3$, because 10 cells are used in the UPML zone.

Starting from the law of Ampere-Maxwell (13), defining the constitutive relations as:

[TABLE]

and applying a Fourier inverse transform, a set of equivalent equations in the time domain are obtained. In a finite difference form for the $x$ component, it can be written:

[TABLE]

and

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Where $n$ denote the index for the time $t_{n}=n\Delta t$. For simplicity, in the previous equations the superscript $"sc"$ has been omitted.

Similarly, from the law of Faraday (14) and the constitutive relations:

[TABLE]

the equations for the evolution of the rectangular components of the magnetic field are obtained.

In Fig.3, the computational cycle steps are shown:

calculation of the current densities in the mesh nodes beginning with the SP positions and velocities data; 2. 2.

computation of the total fields in the mesh nodes and their velocities data by using the charge conservation method [15]; 3. 3.

calculation of the total electric and magnetic fields in the mesh nodes, 4. 4.

calculation of the new positions and velocities of the SPs through integration of their motion equations.

The electric and magnetic fields calculated in the mesh nodes and then interpolated to the SP positions make it possible to determine the forces acting on particles that predetermines the SP positions and their impulses on the next time step (See Fig.3).

3 Results

The process of propagation of $290$ $kW$ microwave power along a $TE_{10}$ rectangular waveguide and its penetration into the chamber through a port made in its cylindrical surface is numerically simulated. In the stationary regime, a microwave field of $1$ $kV/cm$ tension is set. The chamber quality factor is of $50$.

The simulations are fulfilled on a rectangular $3D$ mesh at the spatial steps $\Delta x$ $=$ $\Delta y$ $=$ $0.07$ $cm$, $\Delta z$ $=$ $0.2$ $cm$ at a time step $\Delta t$ $=$ $2.07$ $ps$, that are chosen in accordance with the Courant stability condition. The initial plasma configuration presents a sphere of $2.27$ $cm$ in radius spaced in the center of the chamber.The energy distribution function is maxwellian at a maximum of $2$ $meV$ (See Figs.4 and 5).

The total number of $1.5\times 10^{12}$ particles is divided into $3\times 10^{6}$ superparticles that reciprocally interact through their electric and magnetic fields. In the course of simulations, the plasma particles are steeply redistributed in the chamber volume, but their total number remains fixed.

The calculated vector image of the magnetostatic field in the transverse plane is shown in Fig.6, and in the longitudinal plane in Fig.7. In Fig.6, we can see that the magnetic field generated by the hexapolar system, has the typical configuration of a minimum-B magnetic trap which contributes to the radial confinement of the plasma. On the other hand, in Fig.7, one can see the magnetistatic field lines in a longitudinal plane.

In Fig.8 and Fig.9, we show the results concerning the energy spatial distributions for electrons after 200 cycles of microwaves, where one can see that the core plasma density exceeds by far the density of the corona region. In general, the plasma is accumulated in the trap axis vicinity and symmetrically to the center of the trap.

The data shown in Fig.10 and Fig.11 evidence that the electron population can be divided into 3 groups: the cold electrons with energies less than $0.2$ $keV$, the electrons with intermediate energy between $0.2$ and $3$ $keV$, the warm electrons whose energies range from $3$ to $10$ $keV$ and the hot electrons with energies higher than $10$ $keV$. An estimation shows that the intermedite and warm electrons constitute 90.3% with reference to the total electron population, the cold electrons contribution is 9.5% and the hot electron fraction is 0.2%.

4 Conclusions

A numerical scheme for simulating the behavior of collisionless plasmas confined in the ECR minimum-B magnetic traps which include the simulation of chamber excitation is elaborated. This scheme implemented for determining the characteristics of ECR plasma heating in a minimum-B magnetic trap shows that the electron component is formed by three groups: cold, warm and hot. Appearing these groups can be attributed to different modes of the microwave electron interaction. This point will be elucidated through some more long-run simulations. The hot group, which is of interest to ECRIS is localized predominantly in the trap center.

Acknowledgments

The work is supported by the Universidad Industrial de Santander under the program of mobility VIE and for the GridUIS-2 testbed of SC3UIS. One of the authors (Alex Estupiñán) would like to thank the Universidad Autónoma de Bucaramanga.