Simulation of an inductively coupled plasma with a two-dimensional Darwin particle-in-cell code

Sydorenko, Dmytro; Kaganovich, Igor D.; Khrabrov, Alexander V.; Ethier, Stephane A.; Chen, Jin; Janhunen, Salomon

doi:10.1063/5.0241152

A two-dimensional particle-in-cell code for the simulation of low-frequency electromagnetic processes in laboratory plasmas has been developed. The code uses the Darwin method omitting the electromagnetic wave propagation. The Darwin method separates the electric field into solenoidal and irrotational parts. The solenoidal electric field is calculated with a new algorithm based on the equation for the electric field vorticity. The system of linear equations in the new algorithm is readily solved using a standard iterative method. The irrotational electric field is the electrostatic field calculated with the direct implicit algorithm. The code is verified by reproducing the two-stream instability, electron electromagnetic waves, and shear Alfvén waves. The code is applied to simulate an inductively coupled plasma with the driving current flowing around the plasma region. In this simulation, a ring of dense plasma forms at the initial stage but then the density becomes maximal in the center and decays monotonically toward the walls. The skin effect is in the transitional mode between local and non-local, and the electron velocity distribution function is non-Maxwellian.

I. INTRODUCTION

In laboratory plasma devices and plasma material processing reactors, the plasma is often dense (density $10^{17} m^{- 3}$ and up), relatively cold (electron temperature a few eV), and has dimensions of tens of centimeters. In this case, the plasma size is thousands of electron Debye lengths in one dimension. If the neutral pressure is a few mTorr, the electron mean free path may exceed the size of the plasma. Kinetic effects in such plasmas are important and can be studied numerically using particle-in-cell (PIC) codes. Explicit PIC algorithms are simple and robust, but they must resolve the electron Debye length and the plasma frequency.¹ Combination of the large number of cells in the numerical grid and the small time step makes the explicit algorithms numerically expensive, especially for simulations resolving more than one spatial dimension. Implicit PIC algorithms are free from these requirements and allow to reach the final state in simulation much sooner compared to explicit codes. On the other hand, fully implicit algorithms are complex and involve solutions of systems of nonlinear equations.^2,3 A simpler predictor–corrector approach is implemented in the direct implicit method.⁴ The direct implicit algorithm is not energy conserving, but effects of numerical heating or cooling may be minimized by proper selection of the cell size and time step.⁵

Electromagnetic simulations using explicit numerical schemes have a strong limitation on the simulation time step because they have to resolve waves propagating with the light speed.⁶ This limitation does not apply to fully implicit electromagnetic codes, which also have to solve a system of nonlinear equations.⁷ Another way to bypass this limitation is to omit the propagation of electromagnetic waves and use a method based on the Darwin approach where the inductive (solenoidal) electric field is neglected in the Ampère–Maxwell law.^8–12 The Darwin PIC algorithm described in Ref. 12 is based on the Hamiltonian approach and has good energy conservation properties.

The PIC code described below is developed based on the two-dimensional Darwin direct implicit PIC code DADIPIC of Gibbons and Hewett.¹³ The electrostatic part of the algorithm is similar to DADIPIC (though it is not an exact copy), while the electromagnetic part has significant differences. DADIPIC uses the Streamlined Darwin Field (SDF) Formulation to find solenoidal components of the electric field.¹⁴ An advantage of the SDF method is the simplicity of boundary conditions it uses. A drawback is that the method produces a system of linear equations that is difficult to solve using most standard iterative linear solvers. Previously, the SDF system of equations was solved using the dynamic alternating-directions implicit technique.^13,14 In this paper, the solenoidal electric field is obtained from an equation for the field vorticity. Nielson and Lewis suggested using such an equation but concluded that solving it is difficult.⁹ However, the system of linear equations resulting from this approach is reliably solved with the bi-conjugate gradient (KSPBCGSL) solver of the Portable Extensible Toolkit for Scientific Computation (PETSc)¹⁵ with a multigrid preconditioner (PCHYPRE) from the LLNL package HYPRE.¹⁶

The code is applied to simulate a side-coil inductively coupled plasma (ICP) system. Such devices are used for spectral chemical analysis, chemical synthesis, coating of materials, etc.^17–19 The side-coil configuration is chosen because the dominant solenoidal electric field components are the ones calculated with the new method. The simulated system has a rather high neutral density and antenna current. Such values are selected to test code stability and performance under extreme conditions rather than to investigate kinetic effects. Nevertheless, the simulation demonstrates the transition from the classical to anomalous skin effect as the plasma density grows and the skin layer size becomes comparable to the electron mean free path. The density is initially maximal near the walls surrounding the plasma but eventually, it acquires a maximum in the center. The electron velocity distribution function (EVDF) in the simulation is non-Maxwellian.

The paper is organized as follows. Section II contains a description of the electromagnetic part of the algorithm. Simulation parameters and numerical results of the side-coil ICP simulation are given in Sec. III. Conclusion is in Sec. IV. The electrostatic direct implicit algorithm is described in Appendix A. Calculation of moments of velocity distribution function used in the Darwin electromagnetic field equations is described in Appendix B. The finite difference form of the equation used to calculate solenoidal electric field components in the simulation plane is given in Appendix C. Appendix D explains how the parameters of the simulation described in Sec. III are selected. Appendix E describes simulations of the two-stream instability, electron electromagnetic waves, and shear Alfvén waves verifying the code.

II. MODEL DESCRIPTION

The code uses particle-in-cell formalism in Cartesian coordinates. Electrons and ions are represented as charged particles. Multiple ion species can be introduced. The code resolves two spatial coordinates (x and y) and three velocity components (⁠ $v_{x}$ ⁠, $v_{y}$ ⁠, and $v_{z}$ ⁠) for each charged particle. The simulation domain is a rectangle. For electrostatic simulations, the domain can be completely enclosed, periodic in one direction only, or periodic in two directions; rectangular metal or dielectric objects can be introduced inside the domain. For electromagnetic simulation, the domain is completely enclosed, without metal objects inside. Including metal objects inside the domain in electromagnetic simulation requires additional steps, which will be described in a separate publication.

The code includes models of electron–neutral collisions, charge-exchange ion–neutral collisions,²⁰ and electron emission from material surfaces due to electron and/or ion bombardment, similar to the models used in other codes developed by the authors.^21–23

The code is written in Fortran and parallelized with Message Passing Interface (MPI). Domain decomposition is used. Linear equation systems defining components of the electrostatic and electromagnetic fields are solved using the PETSc library. To ensure the efficient operation of parallel field solvers, the whole domain is divided into sub-domains (referred to as blocks) of the same size. The number of blocks equals the number of MPI processes. The code includes an algorithm ensuring approximately even the numerical load of each MPI process at the particle processing stage. Particles are advanced within sub-domains (referred to as clusters) combining several adjacent blocks. Particles belonging to a single cluster are processed by several MPI processes. The number of these processes changes dynamically, depending on the number of particles in the cluster, and is always no less than one.

The code uses four staggered structured uniform grids, similar to the Yee grids for electromagnetic simulation,⁶ shown in Fig. 1. The first grid has nodes with coordinates $x_{i} = i h$ ⁠, $y_{j} = j h$ ⁠, where i and j are integer numbers, $h$ is the cell size. Below, this grid is referred to as the centered grid. The second grid is shifted in the x-direction by a half-a-cell relative to the centered grid, it has nodes with coordinates $x_{i + 1 / 2} = (i + 1 / 2) h$ ⁠, $y_{j} = j h$ ⁠. The third grid is shifted in the y-direction and has nodes with coordinates $x_{i} = i h$ ⁠, $y_{j + 1 / 2} = (j + 1 / 2) h$ ⁠. The fourth grid is shifted in both x and y directions, coordinates of its nodes are $x_{i + 1 / 2} = (i + 1 / 2) h$ ⁠, $y_{j + 1 / 2} = (j + 1 / 2) h$ ⁠. In this paper, a node of a grid is specified by the pair of factors before the $h$ in the node coordinates x and y introduced above. One or both of these factors may be half-integer for a node of one of the shifted grids. For instance, node $(i, j)$ belongs to the centered grid and has coordinates $x_{i}$ and $y_{j}$ ⁠, node $(i + 1 / 2, j)$ belongs to the second grid and has coordinates $x_{i + 1 / 2}$ and $y_{j}$ ⁠, etc. Boundaries of the simulation domain and material objects inside the domain are aligned with the nodes of the centered grid. The electrostatic potential $Φ$ ⁠, charge density $ρ$ ⁠, vector components of electric field $E_{z}$ and electric current density $J_{z}$ are defined on the centered grid (solid purple circles in Fig. 1). Vector components of electric field $E_{x}$ ⁠, magnetic field $B_{y}$ ⁠, and electric current density $J_{x}$ are defined on the second grid (blue triangles in Fig. 1). Vector components of electric field $E_{y}$ ⁠, magnetic field $B_{x}$ ⁠, and electric current density $J_{y}$ are defined on the third grid (black diagonal crosses in Fig. 1). Vector component of magnetic field $B_{z}$ is defined on the fourth grid (red circled dots in Fig. 1).

FIG. 1.

View large Download slide

Markers show the positions of nodes of numerical grids on the coordinate plane x–y. The solid purple circles are nodes of the centered grid shown by the purple straight lines. The triangles and the diagonal crosses represent nodes of the grid shifted along the x and y directions, respectively. The red circled dots are nodes of the grid shifted along both the x and y directions.

Below, the electric field and current are often separated into their solenoidal (⁠ ${\vec{E}}_{sol}$ ⁠, ${\vec{J}}_{sol}$ ⁠) and irrotational (⁠ ${\vec{E}}_{irr}$ ⁠, ${\vec{J}}_{irr}$ ⁠) parts. This means that $\vec{E} = {\vec{E}}_{sol} + {\vec{E}}_{irr}$ where $\nabla \cdot {\vec{E}}_{sol} = 0$ and $\nabla \times {\vec{E}}_{irr} = 0$ ⁠. Similarly, $\vec{J} = {\vec{J}}_{sol} + {\vec{J}}_{irr}$ where $\nabla \cdot {\vec{J}}_{sol} = 0$ and $\nabla \times {\vec{J}}_{irr} = 0$ ⁠.

The numerical scheme is of a “leap-frog” type, with particle coordinates and velocities defined at times shifted by half of a time step $Δ t$ ⁠. In this paper, superscript n denotes values defined at time $t^{n} = n Δ t$ ⁠, superscript $n \pm 1 / 2$ denotes values defined at time $t^{n \pm 1 / 2} = t^{n} \pm Δ t / 2$ ⁠, and subscript p denotes values related to a single particle.

The direct implicit Darwin algorithm is well described in Ref. 13. A brief description of the electrostatic part of this algorithm is provided in Appendix A for convenience. The electromagnetic part is described below since it has significant differences from Ref. 13.

A. Electromagnetic equations in Darwin method

The Darwin method uses the following equations for solenoidal electric and self-consistent magnetic fields:

\nabla \times \vec{B} = μ_{0} {\vec{J}}_{sol},

(1)

\nabla^{2} {\vec{E}}_{sol} = μ_{0} \frac{\partial {\vec{J}}_{sol}}{\partial t},

(2)

where

μ_{0}

is the magnetic constant and

{\vec{J}}_{sol}

is the solenoidal part of the electric current. Equation follows from Ampère–Maxwell law where the solenoidal part of the displacement current

c^{- 2} \partial {\vec{E}}_{sol} / \partial t

is omitted and the relationship

\frac{\partial {\vec{E}}_{irr}}{\partial t} = - \frac{{\vec{J}}_{irr}}{ε_{0}},

(3)

where

ε_{0}

is the electric constant, is applied. Equation follows from Gauss's law and charge continuity. Equation is obtained by applying the curl operator to Faraday's law and using .

Ref. 13 mentions that finite time differencing cannot be applied in because it causes strong numerical instability. Instead, one can use Boltzmann's equation with the Bhatnagar-Gross-Krook (BGK) collision term

\frac{\partial (n m \vec{u})}{\partial t} = n q (\vec{E} + \vec{u} \times \vec{B}) - \nabla (m n ⟨ \vec{v} \vec{v} ⟩) - m n \vec{u} ν

(4)

and obtain the following equation for the time derivative of the total electric current:

μ_{0} \frac{\partial \vec{J}}{\partial t} = μ {\vec{E}}_{sol} + \vec{Q} .

(5)

The calculation of coefficients

μ

and

\vec{Q}

used in is described in Appendix B. Note that

\vec{Q}

accounts for particle scattering by neutrals, see Eqs. and .

B. Solenoidal electric field

The x and y components of (5) contain both the solenoidal and the irrotational parts. Separating these parts without knowing $E_{sol, x}$ and $E_{sol, y}$ is not straightforward. One way of doing this is the SDF method.^13,14 In the present paper, the following alternative approach is described.

One can apply the curl operator to and use

\nabla \times {\vec{J}}_{sol} = \nabla \times \vec{J}

to obtain

\nabla^{2} \nabla \times {\vec{E}}_{sol}^{n + 1} = μ_{0} \nabla \times \frac{\partial \vec{J}}{\partial t} .

(6)

In order to ensure that

\nabla \cdot {\vec{E}}_{sol} = 0

⁠, one can assume

{\vec{E}}_{sol}^{n + 1} = \nabla \times \vec{F},

(7)

where

\vec{F}

is some vector function. Substituting and into Eq. one obtains

\nabla^{2} (\nabla \times \nabla \times \vec{F}) = \nabla \times (μ \nabla \times \vec{F} + \vec{Q}) .

(8)

Note that Eq. is similar to Eq. of Ref. 9. In the selected 2D geometry, it is sufficient to know only the z component of

\vec{F}

⁠. Below, for the sake of brevity, this component is denoted F, with subscript z omitted. Then becomes

\begin{matrix} (\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}) (\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}) F - (\frac{\partial}{\partial x} μ \frac{\partial}{\partial x} + \frac{\partial}{\partial y} μ \frac{\partial}{\partial y}) F = \frac{\partial Q_{x}}{\partial y} - \frac{\partial Q_{y}}{\partial x} . \end{matrix}

(9)

Once the F values are known, the solenoidal electric field components in the simulation plane are found as

\begin{matrix} E_{sol, x}^{n + 1} = \frac{\partial F}{\partial y}, E_{sol, y}^{n + 1} = - \frac{\partial F}{\partial x} . \end{matrix}

(10)

Equation (9) creates a system of linear equations for values of F defined in the nodes of the grid shifted relative to the centered grid in both the x and y directions. The finite difference form of these equations, the related boundary conditions, and the method of solving the system are described in Appendix C. Since the left-hand side (LHS) of (6) contains the curl of ${\vec{E}}_{sol}$ ⁠, below the method of calculation of $E_{sol; x, y}$ involving Eqs. (9) and (10) is referred to as the vorticity method.

The z-component of the time derivative of the predicted electric current is solenoidal. In combination with , it gives the following equation for the solenoidal electric field normal to the simulation plane:

(\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} - μ) E_{sol, z}^{n + 1} = Q_{z} .

(11)

The boundary condition for is

E_{sol, z} = 0

at the domain surface.

C. Self-consistent magnetic field

Prior to calculating the magnetic field, the predicted full electric current

{\vec{J}}_{pred}

must be found as described by in Appendix B. Equation requires the solenoidal part of the current. Since any vector normal to the simulation plane is solenoidal, one can use

J_{pred, z}

directly in the z-component of specifying

B_{x}

and

B_{y}

⁠. The situation with the

B_{z}

is less straightforward. The solenoidal current in the simulation plane required to find

B_{z}

can be obtained from

{\vec{J}}_{sol} = {\vec{J}}_{pred} - \nabla ϕ

where scalar function

ϕ

satisfies

\nabla^{2} ϕ = \nabla \cdot {\vec{J}}_{pred}

⁠. The corresponding equation

(\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}) ϕ = \frac{\partial J_{pred, x}}{\partial x} + \frac{\partial J_{pred, y}}{\partial y}

(12)

is solved with boundary condition

ϕ = 0

at the domain boundary.²⁴ The solenoidal current components are found as

\begin{matrix} J_{sol, x} = J_{pred, x} - \frac{\partial ϕ}{\partial x}, \\ J_{sol, y} = J_{pred, y} - \frac{\partial ϕ}{\partial y} . \end{matrix}

(13)

To calculate the magnetic field, a vector potential

\vec{A}

is introduced such that

\vec{B} = \nabla \times \vec{A},

(14)

and

\vec{A}

satisfies the Coulomb gauge,

\nabla \cdot \vec{A} = 0

⁠. With , Eq. gives

(\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}) A_{x} = - μ_{0} J_{sol, x},

(15)

(\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}) A_{y} = - μ_{0} J_{sol, y},

(16)

and

(\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}) A_{z} = - μ_{0} J_{pred, z} .

(17)

In this paper, Eqs. (15)–(17) are solved with boundary conditions $A_{x} = 0$ ⁠, $A_{y} = 0$ ⁠, and $A_{z} = 0$ at the domain boundary, respectively. The boundary condition for $A_{z}$ ensures that the magnetic field in the simulation plane is normal to the surface of the domain.

The magnetic field components are calculated as

\begin{matrix} B_{x}^{n + 1} = \frac{\partial A_{z}}{\partial y}, B_{y}^{n + 1} = - \frac{\partial A_{z}}{\partial x}, \end{matrix}

(18)

B_{z}^{n + 1} = \frac{\partial A_{y}}{\partial x} - \frac{\partial A_{x}}{\partial y} .

(19)

It is necessary to mention that the values of

A_{x}

are defined in the nodes of the grid shifted along the x-direction (same as

J_{x}

⁠),

A_{y}

are defined in the nodes of the grid shifted along the y-direction (same as

J_{y}

⁠), and

A_{z}

are defined in the nodes of the centered grid (same as

J_{z}

⁠).

An alternative to Eqs. , , and is to calculate the magnetic field normal to the simulation plane using the z-component of the curl of Eq. . This eliminates the need to calculate the solenoidal electric current:¹³

- (\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}) B_{z} = μ_{0} (\frac{\partial J_{pred, y}}{\partial x} - \frac{\partial J_{pred, x}}{\partial y}) .

(20)

Equation is solved with boundary condition

B_{z} = 0

at the simulation domain.

The method of calculation of $B_{z}$ has to be selected based on the peculiarities of a specific problem. One way is to choose the method, which satisfies best the energy conservation law, as described in Sec. III A.

D. External drivers

In an ICP simulation, the simulation domain is completely enclosed by grounded metal walls. An antenna driving the ICP is represented by a region where the electric current density is a given function of time, ${\vec{J}}^{ext} (t)$ ⁠. This external current must have zero divergences. Note that $\partial {\vec{J}}^{ext} / \partial t$ is also a known function of time. Such a region is considered transparent to the electromagnetic field. Below, it is referred to as the current patch. The electric current in the patch can be in the simulation plane, i.e., has components $J_{x}^{ext}$ and $J_{y}^{ext}$ ⁠, or be normal to the simulation plane, i.e., has component $J_{z}^{ext}$ only.

With the external currents, the right-hand sides (RHS) of Eqs. acquire an additional term

- μ_{0} J_{x, y, z}^{ext}

⁠, respectively. The RHS of acquires an additional term

μ_{0} (\frac{\partial J_{y}^{ext}}{\partial x} - \frac{\partial J_{x}^{ext}}{\partial y}) .

(21)

The RHS of Eq. acquires an additional term

μ_{0} (\frac{\partial}{\partial y} \frac{\partial J_{x}^{ext}}{\partial t} - \frac{\partial}{\partial x} \frac{\partial J_{y}^{ext}}{\partial t}) .

(22)

The RHS of Eq. acquires an additional term

μ_{0} \frac{\partial}{\partial t} J_{z}^{ext}

⁠.

E. Energy conservation

The code diagnostics includes an analysis of energy conservation. The total kinetic energy of all particles

W_{kin}

satisfies

\frac{\partial W_{kin}}{\partial t} = L_{wall} + L_{coll} + P_{field},

(23)

where

L_{wall}

and

L_{call}

represent energy losses due to collisions with walls and neutrals, respectively, and

P_{field}

is the power spent by the electric field to move particles. The code calculates the left-hand side (LHS) of ,

L_{wall}

⁠, and

L_{coll}

accounting for contributions from all particles and then finds the power spent by the electric field using

P_{field} = \frac{\partial W_{kin}}{\partial t} - L_{wall} - L_{coll} .

(24)

Since the electric field at the boundaries of the simulation domain is zero, the electromagnetic field energy flux through the boundary is also zero and the electromagnetic field energy

W_{field}

satisfies

\frac{\partial W_{field}}{\partial t} = - P_{field} - \int_{V} d^{3} r {\vec{J}}^{ext} \cdot \vec{E},

(25)

where

W_{field} = W_{E} + W_{B}

⁠,

W_{E} = \int_{V} d^{3} r \frac{ε_{0} {\vec{E}}^{2}}{2}

(26)

is the energy of the electric field,

W_{B} = \int_{V} d^{3} r \frac{{\vec{B}}^{2}}{2 μ_{0}}

(27)

is the energy of the magnetic field. The last term in the RHS of represents the power deposited in the system by the external electric current, the integration is performed over the volume of the system. Combining and gives the conservation law of the full energy

W_{full}

⁠,

\frac{\partial W_{full}}{\partial t} = L_{wall} + L_{coll} - \int_{V} d^{3} r {\vec{J}}^{ext} \cdot \vec{E},

(28)

where

W_{full} = W_{kin} + W_{field}

⁠. In view of , Eq. is always satisfied. Equations and , however, may have an imbalance between the LHS and RHS due to errors in the electromagnetic field and numerical heating or cooling included in the

P_{field}

obtained via . The rate of full energy change due to numerical effects can be calculated e.g., from as

{(\frac{\partial W}{\partial t})}_{num} = \frac{\partial W_{field}}{\partial t} + P_{field} + \int_{V} d^{3} r {\vec{J}}^{ext} \cdot \vec{E} .

(29)

Below this rate is referred to as the numerical energy rate.

III. SIMULATION OF AN INDUCTIVELY COUPLED PLASMA

A. System setup and numerical parameters

The simulation domain is a square surrounded by metal walls (purple square in Fig. 2). A relatively thin layer of dielectric with dielectric constant

ε_{r} = 1

is placed along the walls leaving an empty square area in the center of the system (shown by the green line in Fig. 2), this is the area where the plasma and the neutral gas exist. Inside the dielectric layer, there is a region with the driving electric current density

J_{x, y}^{ext}

bounded between two squares (shown by red lines in Fig. 2). The driving electric current flows in a closed loop around the plasma region. The electric current density components are defined as

\begin{matrix} J_{x}^{ext} (x, y, t) = I_{1} sin (ω t) \frac{\partial}{\partial y} g (x, y), \\ J_{y}^{ext} (x, y, t) = - I_{1} sin (ω t) \frac{\partial}{\partial x} g (x, y), \end{matrix}

(30)

where

I_{1}

is the amplitude of the total electric current flowing through the region's cross section per 1 m of length along the z-direction (in units of A/m),

ω

is the driving current frequency, and function g is zero outside the external boundary of the current region,

g = 1

inside the internal boundary of the current region, and inside the current region it is a solution of Poisson's equation

(\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}) g = 0,

(31)

solved with boundary conditions

g = 0

and

g = 1

at the external and internal boundaries of the driving current region, respectively. Note that the electric current with components defined by has zero divergence.

FIG. 2.

View large Download slide

Setup for simulation of an inductively coupled plasma with the driving current in the x–y plane. The cross marks the probe where the dependencies of amplitude and phase shift of $E_{sol, y}$ vs time shown in Fig. 5(b) are obtained.

The cell size of the numerical grid

h

and the time step

Δ t

are calculated from the following input parameters: the scale electron density

n_{e, 0}

⁠, the scale electron temperature

T_{e, 0}

(in units of eV), the ratio of the scale Debye length and the grid cell size

α

⁠, and the ratio of the maximal allowed particle velocity and the scale electron thermal velocity

β

⁠. The scale Debye length is calculated as

λ_{e, 0} = v_{t h, e, 0} / ω_{e, 0}

where

v_{t h, e, 0} = \sqrt{2 e T_{e, 0} / m_{e}}

is the scale electron thermal velocity, and

ω_{e, 0} = \sqrt{n_{e, 0} e^{2} / m_{e} ε_{0}}

is the scale electron plasma frequency. The grid cell size and the time step are calculated as

h = \frac{λ_{e, 0}}{α}, Δ t = \frac{h}{β v_{t h, e, 0}} .

(32)

In the simulation described in this section, $n_{e, 0} = 1.6 \times 10^{18} m^{- 3}$ ⁠, $T_{e, 0} = 4 eV$ ⁠, $α = 0.044194$ (the exact value is $\sqrt{2} / 32$ ⁠), $β = 9$ ⁠, $ω_{e, 0} = 7.136 \times 10^{10} s^{- 1}$ ⁠, $v_{t h, e, 0} = 1.186 \times 10^{6} m / s$ ⁠, $λ_{e, 0} = 1.662 \times 10^{- 5} m$ ⁠, $h = 3.761 \times 10^{- 4} m$ ⁠, and $Δ t = 3.523 \times 10^{- 11} s$ ⁠. The number of particles per cell for the scale density is $N_{ppc} = 8000$ ⁠. For a plasma with density and temperature as the scale values, the grid does not resolve the Debye length (⁠ $h / λ_{e, 0} = 22.63$ ⁠) and the time step barely resolves the period of electron plasma oscillations (⁠ $ω_{e, 0} Δ t = 2.51$ ⁠). The grid, however, resolves the skin layer length which for the plasma with the scale density is $c / ω_{e, 0} = 4.2 \times 10^{- 3} m$ ⁠. The reasoning behind the selection of $β$ and $N_{ppc}$ is given in Appendix D.

The centered grid has 266 nodes along the x and y directions. Correspondingly, the size of each side of the domain is $99.7 mm$ ⁠, see Fig. 2. The square plasma region in the center has $186 \times 186$ nodes, and the size of one side of the region is $69.6 mm$ ⁠. The external and internal square boundaries of the driving current region are $226 \times 226$ and $206 \times 206$ nodes, respectively. The sides of the external and internal driving current region boundaries are $84.6$ and $77.1 mm$ ⁠, respectively. The amplitude of the driving current is $I_{1} = 1000 A / m$ ⁠. The frequency of the driving current is $10 MHz$ ⁠. The plasma region is filled with uniform neutral gas Argon of density $10^{21} m^{- 3}$ and temperature $300 K$ ⁠, and the corresponding pressure is $31 mTorr$ ⁠. Elastic, inelastic, and ionization collisions between electrons and neutral atoms are accounted for. Charge-exchange collisions between ions and neutral atoms are accounted for. Particles colliding with the dielectric surrounding the plasma are considered absorbed by the dielectric and contribute to the surface charge. The ion mass is $40 AMU$ ⁠. The simulation starts with the uniform plasma of density $10^{16} m^{- 3}$ ⁠, electron temperature $4 eV$ ⁠, and ion temperature $0.03 eV$ ⁠. Initial electron and ion velocity distributions are Maxwellian, electromagnetic fields are zero, and particle accelerations (A10) used by the electrostatic algorithm are zeros. The simulated time interval is 30 000 ns. The simulation runs on a 96-core node of the Stellar cluster in PPPL and takes about 22 days of wall time to complete.

The simulation is initially attempted with the direct $B_{z}$ solver, which uses Eq. (20). It is found, however, that in this case there is a significant difference between the instant rate of the field energy change [LHS of Eq. (25)] and the corresponding balance of source and sink terms [RHS of Eq. (25)], compare curves 1 and 2 in Fig. 3(a). Then the simulation is attempted with $B_{z}$ solver using the vector potentials [Eqs. (15), (16), and (19)]. With this solver, the difference between the LHS and RHS of Eq. (25) is minimal, see virtually indistinguishable curves 1 and 2 in Fig. 3(b). The difference between the two solvers appears because the boundary condition $B_{z} = 0$ used by the direct solver results in overestimating the magnetic field and, correspondingly, the magnetic field energy, compare curves 1 and 2 in Fig. 3(c). At the same time, both simulations have very close values of the electric field energy, compare curves 1 and 2 in Fig. 3(d). In view of the energy conservation error introduced by the direct solver, the simulation described in Sec. III is carried out with $B_{z}$ calculated using the vector potentials.

FIG. 3.

View large Download slide

Panels (a) and (b) show time dependencies of the rate of change of field energy (curve 1, purple) and the RHS of Eq. (25) (curve 2, green) when $B_{z}$ is calculated directly (a) and via vector potentials (b). Panels (c) and (d) show time dependencies of the magnetic field energy (c) and electric field energy (d) when $B_{z}$ is calculated directly (curve 1, red) and via vector potentials (curve 2, black).

B. Time dependencies of global parameters

The plasma is sustained by the applied oscillating driving current. Oscillations of the driving current [see Fig. 4(a)] create a time-varying magnetic field $B_{z}$ [see Fig. 4(b)] and, correspondingly, a solenoidal electric field with components $E_{sol; x, y}$ [see $E_{sol, y}$ in Fig. 4(c)]. The other electromagnetic field components $B_{x, y}$ and $E_{sol, z}$ are relatively small and are not considered here. The electric field penetrating into the plasma creates an oscillating electric current $J_{x, y}$ [see $J_{y}$ in Fig. 4(d)]. The plasma electrons are energized and produce ionization, the average electron energy and ionization frequency make two oscillations per one period of the driving current, see Figs. 4(e) and 4(f).

FIG. 4.

View large Download slide

Driving electric current (a), magnetic field $B_{z}$ (b), solenoidal electric field $E_{sol, y}$ (c), plasma electric current density $J_{y}$ (d), average electron kinetic energy (e), and average ionization electron–neutral collision frequency (f) vs time. In (a), the positive values correspond to the current flowing in the counterclockwise direction. Curves shown in (b), (c), and (d) are obtained in a location (probe) shown by a diagonal cross in Figs. 15(b), 16(b), and 17(b), respectively. In (e) and (f), the averaging is performed over all electron particles. Note that the collision frequency in (f) is not averaged over time like the frequencies in Fig. 5(d).

The number of macroparticles grows with time and by the end of the simulation it exceeds 250 × 10⁶ for each charged species, see Fig. 5(a). The growth of the number of particles is initially exponential but it slows down after 10 000 ns. While the plasma density rapidly grows between $t = 5000$ and $10 000 ns$ ⁠, the amplitude of the solenoidal electric field near the driving current region decreases [see curve 1 in Fig. 5(b)] and the phase shift between the field and the driving electric current reduces from 90 ° to about 80 ° [see curve 2 in Fig. 5(b)]. In a real ICP device, such a change of the phase shift between the current and voltage in the antenna corresponds to the growth of irreversible power losses (Joule heating).

FIG. 5.

View large Download slide

Time dependencies of (a) number of electron macroparticles, (b) amplitude (curve 1) and phase shift (curve 2) of $E_{sol, y}$ in probe shown by the black diagonal cross in Fig. 2, (c) average kinetic electron energy, (d) ionization (curve 1) and elastic (curve 2) electron–neutral collisions frequencies, and (e) number of ionization collisions per time step (curve 1) and number of macroparticles escaped at the dielectric walls per time step (curve 2). Values in (d) and (e) are averaged over one driving current period (⁠ $100 ns$ ⁠). In (b) and (d), curves 1 and 2 correspond to the left and right vertical coordinate axes, respectively (note the purple and green arrows pointing left and right). In (b), the phase shift (curve 2) is calculated relative to the phase of the driving electric current $I_{1} sin (ω t)$ ⁠, i.e., it represents $E_{sol, y} \approx | E_{sol, y} | sin (ω t + ϕ_{E_{sol, y}})$ ⁠.

The average electron energy gradually decreases as more and more low-energy electrons are produced via ionization while the driving electromagnetic field penetrates only into a narrow boundary layer of the dense plasma, see Fig. 5(c). The average frequency of ionization collisions rapidly decreases during the period of exponential growth of the number of particles, between 3000 and 10 000 ns, see curve 1 in Fig. 5(d). After that, the average ionization collision frequency decreases much slower, in the end of the simulation it is about $4.6 \times 10^{4} s^{- 1}$ ⁠. The average frequency of elastic electron–neutral collisions is several orders of magnitude higher than the ionization frequency, compared curves 2 and 1 in Fig. 5(d). Although the ionization collision frequency becomes almost constant, the system is not in a quasi-stationary state yet. Final particle losses at the walls are about 63% of the total ionization rate, compared curves 1 and 2 in Fig. 5(e). If the slopes of these curves remain the same as they are at $t = 30 000 ns$ ⁠, the two curves will cross around $t = 56 000 ns$ and then the quasi-stationary state will be achieved.

Figure 6 represents various energy change rates averaged over the driving rf period. The rate of kinetic energy change [curve 1 in Fig. 6(a)] is positive and small compared to energy losses at walls and due to collisions with neutrals [curves 2 and 3 in Fig. 6(a), respectively]. The corresponding power spent by the electric field to compensate for the losses is shown by curve 1 in Fig. 6(b). It is very close to the power deposited in the system by the external rf electric current, compare curves 1 and 2 in Fig. 6(b). The difference between curves 2 and 1 in Fig. 6(b) is the RHS of (25). This difference is not equal to the actual rate of variation of the electromagnetic field energy [the LHS of (25)], compared curves 2 and 1 in Fig. 6(c). The corresponding numerical energy rate calculated using (29) is shown by curve 1 in Fig. 6(d).

FIG. 6.

View large Download slide

(a) Rate of change of kinetic energy of all particles (curve 1), rate of kinetic energy loss due to collisions with walls (curve 2), and rate of kinetic energy change due to collisions with neutrals (curve 3) vs time. (b) Power spent by the electric field to move particles (curve 1) and power deposited in the system by the external oscillating electric current (curve 2) vs time. (c) Rate of change of electromagnetic field energy (curve 1) and the difference between the external deposited power and the electric field power spent to move particles (curve 2) vs time. (d) Rate of full energy change due to the numerical scheme (curve 1) and the rate of change of the full energy (curve 2) vs time. All values are averaged over one rf period of 100 ns and are obtained for 1 m of length along the z-direction.

During the simulation, the numerical energy rate is between $- 25 %$ and $+ 15 %$ of the full energy change rate, compare curves 1 and 2 in Fig. 6(d). While these numbers look large, they are not of great concern because the numerical energy rate is only between $- 2 %$ and $+ 0.2 %$ of the external power deposited in the system, compare curve 1 in Fig. 6(d) with curve 2 in Fig. 6(b). Therefore, in this simulation, the violation of the energy conservation caused by numerical effects is insignificant.

The power deposited by the driving current is very small at the beginning of the simulation and grows until $t = 8800 ns$ ⁠, see Fig. 6(b). This matches the change of the phase shift between the driving current and the induced electric field shown by curve 2 in Fig. 5(b).

The rate of full energy change gradually decreases after $t = 7000 ns$ ⁠, see curve 2 in Fig. 6(d). If the slope of the curve for $t > 30 000 ns$ remains the same as at the end of the simulation, the rate of full energy change becomes zero around $t = 60 000 ns$ ⁠. This estimate of the time when the quasi-stationary state is achieved is close to the value of $56 000 ns$ obtained above from the extrapolation of the particle wall loss and ionization curves.

The analysis of particle sources and sinks and the analysis of the energy rates show that the quasi-stationary state will be achieved as follows. The growing plasma density will increase screening of the driving rf field, which will reduce the power deposited by the rf current in the plasma and reduce the ionization rate. At the same time, the losses of particles at the walls will grow. The quasi-stationary state will be achieved when the power deposited in the plasma is balanced by the wall losses and losses due to inelastic collisions with neutrals which includes ionization.

The simulation is not continued till the quasi-stationary state because of the high numerical cost of the late simulation stage caused by the large number of particles. In its turn, this is due to the large value of $N_{ppc}$ selected to minimize numerical effects during the initial stage of the simulation, as explained in Appendix D.

C. System evolution during the initial/transitional stage

The simulation begins with uniform low plasma density and electron temperature. In a collisional plasma, the classical skin depth is calculated as

δ = \frac{c}{ω_{e}} \frac{{(1 + ν_{e}^{2} / ω^{2})}^{1 / 4}}{cos [{tan}^{- 1} (ν_{e} / ω) / 2]},

(33)

where

ω_{e}

is the electron plasma frequency and

ω

is the electromagnetic wave frequency.²⁵ For the initial plasma parameters and the electron–neutral elastic collision frequency of

ν_{e} = 1.586 \times 10^{8} s^{- 1}

the classical skin depth is

δ = 106 mm

⁠. This length exceeds the size of the plasma, therefore, at the initial stage of simulation, the spatial variation of the solenoidal electric field is defined not by the skin effect but by the condition of asymmetry relative to the center of the system. Namely, the solenoidal electric field components

E_{sol; x, y}

are zero in the center and are approximately linear functions of the distance from the center, see Fig. 7.

FIG. 7.

View large Download slide

Solenoidal electric field at time $t = 49.7 ns$ ⁠. Panel (a) shows $E_{sol, y}$ vs the x-coordinate in the cross section with $y = 50 mm$ ⁠. Panel (b) shows the field vs the x- and y-coordinates. In (b), the color map shows the absolute value of the field, $| E_{sol, x y} | = {(E_{sol, x}^{2} + E_{sol, y}^{2})}^{1 / 2}$ ⁠, the arrows show the direction of the vector. The cross section shown in (a) is marked by the horizontal black solid line in (b).

The heating of the plasma and, correspondingly, the ionization rate are the strongest in the regions adjacent to the plasma boundary where the solenoidal electric field is the strongest. This causes the growth of plasma density within a region of the width of about $20 mm$ along the plasma boundary, see Fig. 8. Instantaneous 2D profiles of the electron density and temperature taken at the initial stage show a ring of dense and relatively warm plasma, see Figs. 9 and 10. The corresponding 2D profiles of the ionization rate also have the characteristic ring structure, see Fig. 11.

FIG. 8.

View large Download slide

Electron density in a cross section with $y = 50 mm$ vs the x-coordinate and time. The cross section is shown by the solid horizontal black line in Figs. 18(b) and 9(b).

FIG. 9.

View large Download slide

Electron density at time $t = 8094.4 ns$ ⁠: average density profile along the cross section with $y = 50 mm$ (a) and the density vs the x- and y-coordinates (b). The values in (a) are averaged over time interval $8000$ – $8095 ns$ ⁠. In (b), the purple and green squares mark the boundaries of the simulation domain and the plasma region, respectively. The cross section shown in (a) is marked by the horizontal black solid line in (b).

FIG. 10.

View large Download slide

(a) Profile of the electron temperature in the x-direction $T_{e, x}$ (red curve) and the y-direction $T_{e, y}$ (black curve) averaged over time interval $8000$ – $8095 ns$ in the cross section with $y = 50 mm$ ⁠. (b) The electron temperature in the y-direction $T_{e, y}$ vs x- and y-coordinates at time $t = 8094.4 ns$ ⁠. In (b), the purple and green squares mark the boundaries of the simulation domain and the plasma region, respectively. Rectangles 1 and 2 mark regions with particles used to calculate the electron velocity distribution function shown in Fig. 13. The horizontal solid black line shows the cross section plotted in (a).

FIG. 11.

View large Download slide

(a) Profile of the ionization rate averaged over time interval $7000$ – $8000 ns$ in a cross section with $y = 50 mm$ ⁠. (b) The ionization rate averaged over time interval $7900$ – $8000 ns$ vs the x- and y-coordinates. In (b), the purple and green squares mark the boundaries of the simulation domain and the plasma region, respectively. The horizontal solid black line shows the cross section plotted in (a).

The temperatures shown by red and black curves in Fig. 10(a) are related to the electron motion in the x-direction only or the y-direction only and are defined as $T_{e, x} = m_{e} (⟨ v_{e, x}^{2} ⟩ - {⟨ v_{e, x} ⟩}^{2})$ and $T_{e, y} = m_{e} (⟨ v_{e, y}^{2} ⟩ - {⟨ v_{e, y} ⟩}^{2})$ ⁠, respectively, here $⟨ ⟩$ means averaging over particles. The $T_{e, x}$ and $T_{e, y}$ temperatures show similar behavior, note that the red and black curves in Fig. 10(a) are virtually indistinguishable. Because of this, only the $T_{e, y}$ is represented by a 2D snapshot in Fig. 10(b).

In order to maintain plasma quasineutrality, the electrostatic potential at the initial stage has to acquire the ring-shaped profile as well. The potential is maximal near the plasma boundary where the plasma density is maximal and is much lower in the center, see Fig. 12. The difference between the potential in the maximum and the center can be as large as 10 V. Such a profile contains plasma electrons energized by the solenoidal electric field within the near-wall region of the dense plasma. As a result, not only the electron temperature in the center is about $1.7 eV$ lower than in the near-wall region, but also the highest electron energy in the center is about $10 eV$ smaller than in the near-wall region, compared EVDFs shown by curves 1 and 2 in Fig. 13. This results in the very low ionization rate in the center of the system during the initial stage of simulation.

FIG. 12.

View large Download slide

(a) Profile of the electrostatic potential averaged over time interval $8000$ – $8095 ns$ in the cross section with $y = 50 mm$ ⁠. (b) The potential vs the x- and y-coordinates at time $t = 8094.4 ns$ ⁠. In (b), the purple and green squares mark the boundaries of the simulation domain and the plasma region, respectively. The horizontal solid black line shows the cross section plotted in (a).

FIG. 13.

View large Download slide

Electron velocity distribution functions over the velocity along the x-direction (a) and y-direction (b) plotted vs the electron energy of motion along the corresponding direction, $w_{x} = m_{e} v_{e, x}^{2} / 2$ and $w_{y} = m_{e} v_{e, y}^{2} / 2$ ⁠. Negative values of the energy correspond to particles moving in the negative direction. Curves 1 (red) and 2 (black) are calculated using particles within regions 1 and 2 shown in Fig. 10(b) and averaged over time interval $8000$ – $8095 ns$ ⁠. Curves 3 and 4 show Maxwellian velocity distributions with temperatures $4.3$ and $2.6 eV$ ⁠, respectively.

The ions accumulated via ionization gradually accelerate away from the density maximum. This creates ion flows toward the plasma boundary outside the dense plasma ring and toward the center inside the dense plasma ring, see Fig. 14. The ion flow toward the center is converging, which increases the plasma density in the central region. During the first $15 000 ns$ this mechanism of density increase is much more important than the ionization.

FIG. 14.

View large Download slide

(a) Profile of the ion flow velocity $V_{i, x}$ in the cross section with $y = 50 mm$ ⁠. (b) The ion flow velocity in the x–y plane vs the x- and y-coordinates. In (b), the color map shows the absolute value of the velocity, $| V_{i, x y} | = {(V_{i, x}^{2} + V_{i, y}^{2})}^{1 / 2}$ ⁠, the arrows show the direction of the vector. The horizontal solid black line shows the cross section plotted in (a). All values correspond to time $8094.4 ns$ ⁠.

As the density in the center grows, the magnitude of the electrostatic potential barrier gradually decreases, which allows more and more energetic electrons to enter the central region and produce ionization there. In addition, the high plasma density reduces the amplitude of the electric field in the plasma, see curve 1 in Fig. 5(b), which reduces the heating and ionization in the near-wall regions. The density becomes maximal in the center and decays monotonically toward the boundary of the plasma region at about $22 000 ns$ ⁠. After this time, the system behaves as described in Sec. III D.

D. Plasma state at the late stage of the simulation

In this section, the focus is on the plasma state at the late simulation stage (⁠ $t > 27 000 ns$ ⁠). Examples of instantaneous snapshots of the magnetic field $B_{z}$ ⁠, solenoidal electric field $E_{sol; x, y}$ ⁠, and solenoidal electric current density $J_{sol; x, y}$ are shown in Figs. 15, 16, and 17, respectively. The magnetic and electric fields in Figs. 15 and 16 correspond to different times when the magnetic and electric fields are at their maximum. The electric current in Fig. 17 corresponds to the time when the driving electric current is almost zero while the solenoidal electric field is maximal (the same time as in Fig. 16). This time is chosen to show rings of electric current with opposite direction which appear inside the plasma at certain times as discussed below. Note that the current density shown in Fig. 4(d) is the full unmodified electric current density $J_{y}$ while Fig. 17 represents the solenoidal electric current density obtained with Eqs. (12) and (13) and combined with the driving external current (30) and (31).

FIG. 15.

View large Download slide

Magnetic field $B_{z}$ at time $t = 28 024.8 ns$ ⁠. Panel (a) shows $B_{z}$ vs the x-coordinate in the cross section with $y = 50 mm$ ⁠. Panel (b) shows the field vs the x- and y-coordinates. The cross section shown in (a) is marked by the horizontal black solid line in (b). The diagonal cross marks the position of a probe where the time dependence shown in Fig. 4(b) is obtained. The white square marks the boundaries of the plasma region.

FIG. 16.

View large Download slide

Solenoidal electric field at time $t = 28049.6 ns$ ⁠. Panel (a) shows $E_{sol, y}$ vs the x-coordinate in the cross section with $y = 50 mm$ ⁠. Panel (b) shows the field vs the x- and y-coordinates. In (b), the color map shows the absolute value of the field, $| E_{sol, x y} | = {(E_{sol, x}^{2} + E_{sol, y}^{2})}^{1 / 2}$ ⁠, the arrows show the direction of the vector. The cross section shown in (a) is marked by the horizontal black solid line in (b). The diagonal cross marks the position of a probe where the time dependence shown in Fig. 4(c) is obtained. The white square marks the boundaries of the plasma region.

FIG. 17.

View large Download slide

Solenoidal electric current density in the simulation plane at time $t = 28 049.6 ns$ ⁠. Panel (a) shows $J_{sol, y}$ vs the x-coordinate in the cross section with $y = 50 mm$ ⁠. Panel (b) shows the current vs the x- and y-coordinates. In (b), the color map shows the absolute value of the current, $| J_{sol, x y} | = {(J_{sol, x}^{2} + J_{sol, y}^{2})}^{1 / 2}$ ⁠, the arrows show the direction of the vector. The cross section shown in (a) is marked by the horizontal black solid line in (b). The diagonal cross marks the position of a probe where the time dependencies shown in Fig. 4(d) are obtained. The white square marks the boundaries of the plasma region.

The plasma density is maximal in the center of the system and gradually decreases toward the dielectric walls, see Fig. 18. In the center, the density exceeds $2.5 \times 10^{18} m^{- 3}$ ⁠. The electron temperature is approximately uniform, in the central region it is slightly smaller than in regions near the plasma boundary (in the skin layer), see Fig. 19. Except for narrow near-wall regions, the temperature is between $3.3$ and $3.7 eV$ ⁠.

FIG. 18.

View large Download slide

Profile of the electron density averaged over time interval $28 000$ – $28 094 ns$ in a cross section with $y = 50 mm$ ⁠. (b) The density vs the x- and y-coordinates at time $t = 28 094.4 ns$ ⁠. In (b), the purple and green squares mark the boundaries of the simulation domain and the plasma region, respectively. The horizontal solid black line shows the cross section plotted in (a).

FIG. 19.

View large Download slide

(a) Profile of the electron temperature in the x-direction $T_{e, x}$ (red curve) and the y-direction $T_{e, y}$ (black curve) averaged over time interval $28 000$ – $28 095 ns$ in the cross section with $y = 50 mm$ ⁠. (b) The electron temperature in the y-direction $T_{e, y}$ vs x- and y-coordinates at time $t = 28 094.4 ns$ ⁠. In (b), the purple and green squares mark the boundaries of the simulation domain and the plasma region, respectively. Rectangles 1 and 2 mark regions with particles used to calculate the electron velocity distribution function shown in Fig. 23. The horizontal solid black line shows the cross section plotted in (a).

The average frequency of elastic electron–neutral collisions at $t = 28 000 ns$ is $ν_{e} = 1.249 \times 10^{8} s^{- 1}$ ⁠, see curve 2 in Fig. 5(d). Since the density and temperature are nonuniform, below the skin depth (33) is estimated using the electron density at $x = 20 mm$ where the electron temperature has a maximum. With $ν_{e}$ as above, $ω = 2 π \times 10^{7} s^{- 1}$ ⁠, and $ω_{e} = 6.176 \times 10^{10} s^{- 1}$ calculated for the density of $1.2 \times 10^{18} m^{- 3}$ at this point, see Fig. 18(a), the classical skin depth is $δ = 8.5 mm$ ⁠. For the electron temperature of $3.7 eV$ ⁠, the thermal velocity is $v_{t h} = 1.14 \times 10^{6} m / s$ and the effective electron mean free path is $v_{t h} / {(ν_{e}^{2} + ω^{2})}^{1 / 2} = 8.2 mm$ ⁠. Since $δ \approx v_{t h} / {(ν_{e}^{2} + ω^{2})}^{1 / 2}$ ⁠, the skin effect is in a transitional mode between the classical and anomalous modes.²⁶

The actual width [in the direction normal to the plasma region boundary] of the area in plasma where the electromagnetic field is significant is about $25 mm$ ⁠. In the center of the plasma, there is a region about $20 mm$ wide where the magnetic field amplitude is an order of magnitude less than the amplitude at the plasma boundary, see Figs. 15 and 20(a). A similar region exists for the solenoidal electric field, see Figs. 16 and 20(b). The solenoidal electric field there is small because the oscillating magnetic field is weak. The magnetic field of the driving electric current is screened by the electric current in the plasma. Interestingly, the plasma current canceling the magnetic field of the driving electric current at its maximum is generated earlier while the solenoidal electric field is the strongest.

FIG. 20.

View large Download slide

Magnetic field $B_{z}$ (a), solenoidal electric field $E_{sol, y}$ (b), and solenoidal electric current density $J_{sol, y}$ (c) in a cross section with $y = 50 mm$ vs the x-coordinate and time during one period of the driving current starting at $t = 28 000 ns$ ⁠. In (c), the values of the current density in the driving current region are multiplied by 0.1 to match the scale of the current density in the plasma region. The two vertical black lines show the boundaries of the plasma region.

The solenoidal electric current in the plasma flows along the plasma boundary forming a ring structure, see Fig. 17. The current induced near the plasma region boundary propagates toward the plasma center (i.e., the ring of current shrinks) with the speed about $5 \times 10^{5} m / s$ ⁠, see Fig. 20(c). Note that this speed is comparable to the electron thermal velocity used for the estimate of the skin depth above. The shrinking lasts for about half of the driving current period before the amplitude of the current in the ring reduces significantly due to collisions. While shrinking, the current travels about $25 mm$ in the radial direction, or about three effective electron mean free path lengths mentioned above. Before one current ring completely dissipates, another ring with the opposite current forms at the plasma edge. The two rings of plasma electric current with opposite directions co-exist when the driving electric current is minimal, see Fig. 20(c) for times 0, 50, and 100 ns. Such rings are shown in Fig. 17. When the driving electric current is maximal, there is only one ring of the plasma electric current directed against the driving current, see Fig. 20(c) for times 25 and 75 ns. The strong plasma current cancels most of the magnetic field of the driving current in the central region.

The ionization rate is maximal and approximately uniform in a wide region in the center of the plasma, see Fig. 21(b) and curve 1 in Fig. 21(a). The ionization collision frequency, however, is maximal in the skin layer and has a minimum in the center, see curve 2 in Fig. 21(a). This curve is obtained as a ratio of the ionization rate of curve 1 in Fig. 21(a) and the electron density curve in Fig. 18(a). The ionization frequency profile is qualitatively similar to the electron temperature profile but the relative difference between the near-wall peaks and the center is noticeably larger for the ionization frequency, compare curve 2 in Fig. 21(a) and the curves in Fig. 19(a). The reason for this is the non-Maxwellian electron velocity distribution function (EVDF).

FIG. 21.

View large Download slide

(a) Curve 1 (red) shows the profile of the ionization rate averaged over time interval $27 000$ – $28 000 ns$ in a cross section with $y = 50 mm$ ⁠, it corresponds to the left vertical coordinate axis (note the red arrow pointing left); curve 2 (black) is the ionization collision frequency, it corresponds to the right vertical coordinate axis (note the black arrow pointing right). (b) The ionization rate averaged over time interval $27 900$ – $28 000 ns$ vs the x- and y-coordinates. In (b), the purple and green squares mark the boundaries of the simulation domain and the plasma region, respectively. The horizontal solid black line shows the cross section plotted in (a).

The EVDF is depleted for energies above 16 eV, which is close to the ionization threshold for Argon, compare curves 1 and 2 with the Maxwellian EVDF curve 3 in Fig. 22. In the skin layer, the EVDF over velocity parallel to the direction of the solenoidal electric field [e.g., curve 1 (red) in Fig. 22(b)] has the number of particles with energy above the ionization potential [relative to the number of particles forming the EVDF maximum] larger than the same EVDF in the middle of the plasma region [curve 2 (black) in Fig. 22(b)]. The EVDFs over velocity normal to the plasma boundary are similar to each other (i.e., without relative difference in the number of electrons for all energies) in the skin layer and the center of the plasma region, compare curves 1 and 2 in Fig. 22(a).

FIG. 22.

View large Download slide

Electron velocity distribution functions over the velocity along the x-direction (a) and y-direction (b) plotted vs the electron energy of motion along the corresponding direction, $w_{x} = m_{e} v_{e, x}^{2} / 2$ and $w_{y} = m_{e} v_{e, y}^{2} / 2$ ⁠, respectively. Negative values of the energy correspond to particles moving in the negative direction. Curves 1 (red) and 2 (black) are calculated using particles within regions 1 and 2 shown in Fig. 19(b) and averaged over time interval $28 000$ – $28 095 ns$ ⁠. Curve 3 (green) is a Maxwellian velocity distribution for temperature 3.33 eV. In (b), curve 4 (blue) is the instantaneous velocity distribution in region 1 at time $28 064.5 ns$ ⁠, for clarity this curve is shifted downward by one order of magnitude.

In the center of the plasma, the EVDF is isotropic: curve 2 in Fig. 22(a) is very close to curve 2 in Fig. 22(b). In the skin layer, the EVDF is not isotropic: curve 1 in Figs. 22(a) and 22(b) are different, which is expected since the solenoidal electric field contributes mostly to the energy of electron motion in the direction parallel to the field. One may notice that the difference does not affect the values of temperatures $T_{e, x}$ and $T_{e, y}$ ⁠, which are very close to each other, compare the red and black curves in Fig. 19(a). The reason is that the EVDFs shown by curves 1 and 2 in Fig. 22 are averaged over one period of the driving electric current for the sake of improved statistics in depleted segments. Instantaneous EVDFs over $v_{y}$ in the skin layer in the region used to calculate curve 1 are asymmetric because of the induced electric current (i.e., electron flow) along the y-direction, an example of such an EVDF is curve 4 in Fig. 22(b). Temperatures shown in Fig. 19 are obtained from instantaneous EVDFs, the asymmetry is accounted for via the flow velocity in the definition of the temperature given above. The averaging of EVDF over the driving current period transformed the asymmetry of multiple EVDFs into the extended width of symmetric curve 1 in Fig. 22(b). The temperature of time-averaged EVDFs may be different from time-averaged temperatures obtained from instantaneous EVDFs.

The electrostatic potential is maximal in the center and decays monotonically toward the walls, see Fig. 23. In the center, the value of the potential averaged over the plasma period is about 19 V. The potential is exactly zero at the metal walls. Charge accumulation on the dielectric boundary results in the potential there being small, usually no more than $\pm 3 V$ ⁠, but non-zero. A sharp monotonic drop of the potential across a single-cell adjacent to the dielectric wall forms automatically and plays the role of a sheath. This potential difference between the center and the dielectric surface is sufficient to contain electrons with energy of motion in the direction normal to the wall about the threshold ionization energy (16 eV). For the average ionization frequency $5 \times 10^{4} s^{- 1}$ ⁠, the mean free path between ionization collisions for an electron with the energy $16 eV$ is about $47 m$ ⁠. This distance is much larger than the size of the plasma. Therefore, the tails of the EVDFs are depleted because energetic electrons escape to the walls, not because they all participate in ionization. On the other hand, for the frequency of elastic electron–neutral collisions of $1.24 \times 10^{8} s^{- 1}$ ⁠, the mean free path of an electron with thermal energy $3.5 eV$ is about $9 mm$ ⁠. This means that electrons participate in numerous collisions as they travel between the walls. It makes the EVDF more isotropic, and it also extends the time an energetic particle spends inside the plasma region before escaping at the wall, which increases the probability of an ionization collision. It is necessary to mention that the simulation is performed without the Coulomb collisions, which may be important for such a dense plasma and make the EVDF shape more Maxwellian.

FIG. 23.

View large Download slide

(a) Profile of the electrostatic potential averaged over time interval $28 000$ – $28 095 ns$ in the cross section with $y = 50 mm$ ⁠. (b) The potential vs the x- and y-coordinates at time $t = 28 095 ns$ ⁠. In (b), the purple and green squares mark the boundaries of the simulation domain and the plasma region, respectively. The horizontal solid black line shows the cross section plotted in (a).

IV. CONCLUSION

A two-dimensional particle-in-cell code in Cartesian geometry has been developed based on the direct implicit Darwin electromagnetic algorithm described in Ref. 13. Significant modifications of the original algorithm have been made. The code uses staggered grids convenient for electromagnetic simulation. The contribution of collisional scattering is included in the calculation of the solenoidal electric fields. The solenoidal electric field in the simulation plane is calculated with a new method using the equation for the vorticity of the solenoidal electric field.⁹ The linear system of equations in the vorticity method is solved using a standard iterative solver.

The self-consistent magnetic field is calculated using either vector potentials and the solenoidal part of the electric current obtained from a Poisson-like equation or the same direct method as Ref. 13. The difference between the methods is in the boundary conditions that can be formulated. Once one of these methods is selected for a particular simulation, the choice can be verified by checking the energy conservation. The code is verified by successful simulation of the two-stream instability, high-frequency plasma electromagnetic waves, and shear Alfvén waves. The verification simulations are described in Appendix E.

The code is applied to simulate an inductively coupled plasma system with the driving antenna wrapped around the plasma. Simulation resolves the plane normal to the axis of the system. In such a configuration, the components of the solenoidal electric field in the simulation plane are important. A simulation lasting 300 driving current periods is performed. In simulation, the plasma is sustained by the electromagnetic fields generated by the oscillating driving electric current in the antenna. Initially, a ring of dense plasma forms along the wall confining the plasma. By the end of the simulation, the plasma density increases two orders of magnitude compared to the initial density and acquires a maximum in the center of the system.

The electron mean free path between collisions with neutrals is small compared to the size of the plasma but is comparable to the classical skin depth for the achieved values of the plasma density. The skin effect is in a transitional mode between the classical (local) and anomalous (non-local) modes. On one hand, the electric current induced by the oscillating electromagnetic field near the plasma boundary is transported inside the plasma by particles moving normally to the boundary. On the other hand, the current quickly dissipates via collisions. The detailed investigation of the skin effect is out of the scope of the present paper.

The electron velocity distribution is non-Maxwellian, with the high-energy tails strongly depleted for energies above the threshold of confinement by the plasma potential. At the initial stage, the electrostatic potential has a non-monotonic profile: it is maximal near the wall confining the plasma (where the plasma density is maximal) and minimal in the center. This shape of the potential maintains plasma quasineutrality and prevents energetic electrons from reaching the center of the system. At the late stage of the simulation, the electrostatic potential is maximal in the center, and the corresponding threshold energy of electron confinement is slightly higher than the ionization threshold.

The simulation is in the implicit mode, the size of the numerical grid cell significantly exceeds the electron Debye length. Although the Debye length is not resolved, a drop of potential forms across a single-cell adjacent to the wall confining the plasma, which plays the role of the plasma sheath.

ACKNOWLEDGMENTS

This work at Princeton Plasma Physics Laboratory (PPPL) was supported by the US Department of Energy CRADA agreement between Applied Material Inc. and PPPL, Contract No. DE-AC02-09CH11466.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Dmytro Sydorenko: Investigation (lead); Software (lead); Visualization (lead); Writing – original draft (lead). Igor D. Kaganovich: Investigation (supporting); Project administration (lead); Supervision (lead). Alexander V. Khrabrov: Investigation (supporting); Software (supporting). Stephane A. Ethier: Software (supporting). Jin Chen: Software (supporting). Salomon Janhunen: Software (supporting).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A: ELECTROSTATIC DIRECT IMPLICIT ALGORITHM

In the beginning of a computational cycle advancing the system from time

t^{n}

to time

t^{n + 1}

⁠, one knows grid values of solenoidal electric field

{\vec{E}}_{sol}^{n}

⁠, irrotational electric field

{\vec{E}}_{irr}^{n}

⁠, magnetic field

{\vec{B}}^{n}

⁠, particle coordinates

{\vec{x}}_{p}^{n}

⁠, velocities

{\vec{v}}_{p}^{n - 1 / 2}

⁠, and accelerations

{\vec{a}}_{p}^{n - 1}

⁠. Particle advance is performed in two stages. First, particles are accelerated to “streaming” velocities and pushed to “streaming” positions:

{\tilde{\vec{v}}}_{p} = R^{n} \cdot {\vec{v}}_{p}^{n - 1 / 2} + \frac{Δ t}{2} (1 + R^{n}) \cdot (\frac{1}{2} {\vec{a}}_{p}^{n - 1} + \frac{q_{s}}{m_{s}} {\vec{E}}_{sol}^{n}),

(A1)

{\tilde{\vec{x}}}_{p} = {\vec{x}}_{p}^{n} + Δ t {\tilde{\vec{v}}}_{p},

(A2)

where

q_{s}

and

m_{s}

are the charge and mass of a particle, subscript s denotes particle species,

1

is a unitary matrix, matrix

R^{n}

describes rotation in the magnetic field

{\vec{B}}^{n}

and has coefficients

\begin{matrix} R_{11} = (1 - θ^{2} + 2 α^{2} B_{x}^{2}) / (1 + θ^{2}), \\ R_{12} = 2 (α^{2} B_{x} B_{y} + α B_{z}) / (1 + θ^{2}), \\ R_{13} = 2 (α^{2} B_{x} B_{z} - α B_{y}) / (1 + θ^{2}), \\ R_{21} = 2 (α^{2} B_{x} B_{y} - α B_{z}) / (1 + θ^{2}), \\ R_{22} = (1 - θ^{2} + 2 α^{2} B_{y}^{2}) / (1 + θ^{2}), \\ R_{23} = 2 (α^{2} B_{y} B_{z} + α B_{x}) / (1 + θ^{2}), \\ R_{31} = 2 (α^{2} B_{x} B_{z} + α B_{y}) / (1 + θ^{2}), \\ R_{32} = 2 (α^{2} B_{y} B_{z} - α B_{x}) / (1 + θ^{2}), \\ R_{33} = (1 - θ^{2} + 2 α^{2} B_{z}^{2}) / (1 + θ^{2}) \end{matrix}

(A3)

with

α = q Δ t / 2 m

and

θ^{2} = α^{2} (B_{x}^{2} + B_{y}^{2} + B_{z}^{2})

⁠, here superscripts denoting time in

B_{x, y, z}

and species subscripts in q and m are omitted for brevity.

Next, the “streaming” charge density

\tilde{ρ}

is calculated in the nodes of the centered grid using “streaming” particle coordinates. The following modified Poisson equation is solved to find the advanced electrostatic potential

Φ^{n + 1}

⁠:

\nabla \cdot [(ε_{r} 1 + X) \cdot \nabla Φ^{n + 1}] = - \frac{\tilde{ρ}}{ε_{0}},

(A4)

where

ε_{r}

is the relative dielectric constant, and

X

is the implicit susceptibility tensor defined as

X = \frac{Δ t^{2}}{4} \sum_{s} (1 + R^{n}) \frac{q_{s} {\tilde{ρ}}_{s}}{m_{s} ε_{0}},

(A5)

where the summation is performed over all charged species. Note that in plasma-filled areas

ε_{r} = 1

while inside the dielectric areas

\tilde{ρ} = 0

and

X = 0

⁠.

In the two-dimensional (2D) system, only four components of tensor

X

are required:

X_{11}

⁠,

X_{12}

⁠,

X_{21}

⁠, and

X_{22}

⁠. Like the charge density, these components are defined in nodes of the centered grid. Finite difference representation of Eq. , however, requires

X_{11}

and

X_{12}

in nodes of the grid shifted in the x-direction and

X_{21}

and

X_{22}

in nodes of the grid shifted in the y-direction. The code calculates the aforementioned tensor components in the centered grid nodes, tensor values in the shifted nodes are obtained by averaging of values in the two neighbor centered nodes. The electrostatic field is obtained as

{\vec{E}}_{irr}^{n + 1} = - \nabla Φ^{n + 1}

⁠. In finite differences, formulas for the electric field are

\begin{matrix} E_{irr, x; i + 1 / 2, j}^{n + 1} = (Φ_{i, j}^{n + 1} - Φ_{i + 1, j}^{n + 1}) / h, \\ E_{irr, y; i, j + 1 / 2}^{n + 1} = (Φ_{i, j}^{n + 1} - Φ_{i, j + 1}^{n + 1}) / h . \end{matrix}

(A6)

The modified Poisson's equation (A4) is solved with the values of the electrostatic potential given at metal walls enclosing the simulation domain (below referred to as the external walls) and at the surfaces of metal objects placed inside the simulation domain (below referred to as the inner objects). It is assumed that the metal is a perfect conductor, the potential is constant along the surface of metal walls and inner objects. Metal walls with different potential may be connected by segments with a linear potential profile along their surface, with the end values corresponding to the potentials of the metal walls. In electrostatic simulation, periodic boundary conditions can be applied at the opposite sides of the simulation domain.

At the second stage of particle advance, the electrostatic field is used to finish acceleration of particles to their final velocities:

δ {\vec{v}}_{p} = \frac{q_{s} Δ t}{4 m_{s}} (1 + R^{n}) {\vec{E}}_{irr}^{n + 1},

(A7)

{\vec{v}}_{p}^{n + 1 / 2} = {\tilde{\vec{v}}}_{p} + δ {\vec{v}}_{p} .

(A8)

In , the electrostatic field is taken in the “streaming” position

{\tilde{\vec{x}}}_{p}

of the particle. Then the particles are pushed to their final positions:

{\vec{x}}_{p}^{n + 1} = {\tilde{\vec{x}}}_{p} + Δ t δ {\vec{v}}_{p} .

(A9)

After this, the particle acceleration is updated as follows:

{\vec{a}}_{p}^{n} = \frac{1}{2} ({\vec{a}}_{p}^{n - 1} + \frac{q_{s}}{m_{s}} {\vec{E}}_{irr}^{n + 1}),

(A10)

where the electrostatic field is taken in the particle's final position

{\vec{x}}_{p}^{n + 1}

⁠. This concludes the electrostatic part of the algorithm.

APPENDIX B: COEFFICIENTS FOR THE DARWIN ALGORITHM

Below, various moments of the velocity distribution function are calculated, which requires to have particle velocities defined at the same time as the final coordinates . These velocities are obtained after the end of the electrostatic part as follows:

{\vec{v}}_{p}^{n + 1} = R^{n} \cdot {\vec{v}}_{p}^{n + 1 / 2} + \frac{q_{s} Δ t}{4 m_{s}} (1 + R^{n}) \cdot ({\vec{E}}_{irr}^{n + 1} + {\vec{E}}_{sol}^{n}) .

(B1)

Note that advances the velocity by half of the time step, which is why matrix

R^{n}

is calculated with

α = q Δ t / 4 m

in coefficients . Velocities are not used to advance particle coordinates. The full predicted electric current

{\vec{J}}_{pred}

required for calculation of the magnetic field is obtained as

{\vec{J}}_{pred, g} = \sum_{s} q_{s} \sum_{p} S ({\vec{x}}_{g} - {\vec{x}}_{p}^{n + 1}) {\vec{v}}_{p}^{n + 1},

(B2)

where subscript g denotes grid values, summation is over all particles (subscript p) of all charged species (subscript s),

{\vec{x}}_{g}

is the grid node coordinate vector, S is the function defining particle contribution to the density in surrounding nodes (corresponds to the bilinear interpolation in 2D).

Coefficients

μ

and

\vec{Q}

in Eq. for the time derivative of the total electric current are

\begin{matrix} μ = \sum_{s} \frac{ω_{p l, s}^{2}}{c^{2}}, \\ \vec{Q} = μ {\vec{E}}_{irr} + \vec{ζ} \times \vec{B} + \vec{K} + \vec{C}, \\ \vec{ζ} = \sum_{s} \frac{ω_{p l, s}^{2}}{c^{2}} {\vec{u}}_{s}, \\ \vec{K} = - μ_{0} \sum_{s} q_{s} \nabla (n_{s} ⟨ {\vec{v}}_{s} {\vec{v}}_{s} ⟩), \\ \vec{C} = - μ_{0} \sum_{s} q_{s} n_{s} {\vec{u}}_{s} ν_{s}, \end{matrix}

(B3)

where

ω_{p l, s}^{2}

is the plasma frequency,

{\vec{v}}_{s}

is the random velocity,

{\vec{u}}_{s}

=

⟨ {\vec{v}}_{s} ⟩

is the flow velocity,

n_{s}

is the number density, and

ν_{s}

is the momentum transfer collision frequency of charged species s, the summation is performed over all charged species. Field solver for

{\vec{E}}_{sol}

requires grid quantities

μ_{g}

and

{\vec{Q}}_{g}

⁠, which are obtained using direct accumulation of particle quantities in grid nodes:

\begin{matrix} μ_{g} = μ_{0} \sum_{s} \frac{q_{s}^{2}}{m_{s}} \sum_{p} S ({\vec{x}}_{g} - {\vec{x}}_{p}^{n + 1}), \\ {\vec{Q}}_{g} = {(μ {\vec{E}}_{irr}^{n + 1})}_{g} + {(\vec{ζ} \times {\vec{B}}^{n})}_{g} + {\vec{K}}_{g} + {\vec{C}}_{g}, \\ {(μ {\vec{E}}_{irr}^{n + 1})}_{g} = μ_{0} \sum_{s} \frac{q_{s}^{2}}{m_{s}} \sum_{p} S ({\vec{x}}_{g} - {\vec{x}}_{p}^{n + 1}) {\vec{E}}_{irr}^{n + 1} ({\vec{x}}_{p}^{n + 1}), \\ {(\vec{ζ} \times {\vec{B}}^{n})}_{g} = μ_{0} \sum_{s} \frac{q_{s}^{2}}{m_{s}} \sum_{p} S ({\vec{x}}_{g} - {\vec{x}}_{p}^{n + 1}) {\vec{v}}_{p}^{n + 1} \times {\vec{B}}^{n} ({\vec{x}}_{p}^{n + 1}), \\ {\vec{K}}_{g} = - μ_{0} \nabla \cdot \sum_{s} q_{s} \sum_{p} {\vec{v}}_{p}^{n + 1} {\vec{v}}_{p}^{n + 1} S ({\vec{x}}_{g} - {\vec{x}}_{p}^{n + 1}), \\ {\vec{C}}_{g} = μ_{0} Δ t^{- 1} \sum_{s} q_{s} \sum_{p_{c}} δ {\vec{v}}_{p_{c}} S ({\vec{x}}_{g} - {\vec{x}}_{p_{c}}^{n + 1}) . \end{matrix}

(B4)

Term

{\vec{C}}_{g}

in represents variation of the electric current due to collisions. In this term, summation over

p_{c}

is summation over particles collided at this time step,

δ {\vec{v}}_{p_{c}} = {\vec{v}}_{p_{c}}^{a} - {\vec{v}}_{p_{c}}^{b}

is the particle velocity change in the collision,

{\vec{v}}_{p_{c}}^{b}

and

{\vec{v}}_{p_{c}}^{a}

are the particle velocity vectors before and after the collision. Collision procedures for electrons and ions are involved after the electrostatic part of the algorithm is completed. Contributions to

{\vec{C}}_{g}

are accounted for every time a collision occurs.

APPENDIX C: FINITE DIFFERENCE FORM OF EQUATION FOR F

The finite difference form of Eq. involves 13 nodes of the F grid positioned according to the stencil shown by red circles in Fig. 24. The equation with the central node of the stencil in node

(i + 1 / 2, j + 1 / 2)

of the F grid is

\begin{matrix} F_{i + \frac{1}{2}, j + 2 \frac{1}{2}} + 2 F_{i - \frac{1}{2}, j + 1 \frac{1}{2}} + (- 8 - h^{2} μ_{i + \frac{1}{2}, j + 1}) F_{i + \frac{1}{2}, j + 1 \frac{1}{2}} + 2 F_{i + 1 \frac{1}{2}, j + 1 \frac{1}{2}} + F_{i - 1 \frac{1}{2}, j + \frac{1}{2}} + (- 8 - h^{2} μ_{i, j + 1 / 2}) F_{i - \frac{1}{2}, j + \frac{1}{2}} + [20 + h^{2} (μ_{i + \frac{1}{2}, j + 1} + \\ μ_{i, j + \frac{1}{2}} + μ_{i + 1, j + \frac{1}{2}} + μ_{i + \frac{1}{2}, j})] F_{i + \frac{1}{2}, j + \frac{1}{2}} + (- 8 - h^{2} μ_{i + 1, j + \frac{1}{2}}) F_{i + 1 \frac{1}{2}, j + \frac{1}{2}} + F_{i + 2 \frac{1}{2}, j + \frac{1}{2}} + 2 F_{i - \frac{1}{2}, j - \frac{1}{2}} + (- 8 - h^{2} μ_{i + \frac{1}{2}, j}) F_{i + \frac{1}{2}, j - \frac{1}{2}} + 2 F_{i + 1 \frac{1}{2}, j - \frac{1}{2}} + F_{i + \frac{1}{2}, j - 1 \frac{1}{2}} = h^{3} (Q_{x; i + \frac{1}{2}, j + 1} - Q_{x; i + \frac{1}{2}, j} + Q_{y; i, j + \frac{1}{2}} - Q_{y; i + 1, j + \frac{1}{2}}) . \end{matrix}

(C1)

FIG. 24.

View large Download slide

Markers show the positions of nodes of numerical grids on the coordinate plane x–y. The solid red circles belong to the grid where F is defined. Nodes of the grid for $Q_{x}$ and $Q_{y}$ used in the right-hand side of (C1) are shown by the blue triangles and black diagonal crosses, respectively. Note that the four coefficients $μ$ in the left-hand side of (C1) are defined in the same nodes as the $Q_{x}$ and $Q_{y}$ ⁠. The solid purple lines go through nodes of the centered grid.

Note that since the F grid nodes are shifted along both directions, they are not on the domain boundary aligned with the nodes of the centered grid. The field solver includes F grid nodes, which are half-a-cell outside the domain. The boundary conditions for (C1) require that $F = 0$ in the nodes half-a-cell outside the domain and in the nodes, which are inside the domain and half-a-cell from the domain boundary. With (10), this condition ensures that both $E_{sol, x}$ and $E_{sol, y}$ are zero at the domain boundary, similar to Ref. 13. The system of equations (C1) is solved using the KSPBCGSL solver of the PETSc library¹⁵ with the multigrid preconditioner from the LLNL package HYPRE.¹⁶

APPENDIX D: SELECTION OF NUMERICAL PARAMETERS

The direct implicit algorithm may introduce numerical heating or cooling.^5,27 However, for given grid spatial resolution and time step, there are optimal values of electron density and temperature for which the numerical heating and cooling are negligible. The optimal set of numerical and plasma parameters is not universal though and depends on details of the simulation setup and the numerical scheme. For example, Ref. 27 considering 2D electrostatic simulation of a double-periodic collisionless plasma with density $n_{e, 0}$ and initial temperature $T_{e, 0}$ suggests using $ω_{e, 0} Δ t λ_{e, 0} / h \approx 0.1$ ⁠. For comparison, Ref. 5 considering 1D simulation of a plasma slab expanding into a vacuum half-space concludes that the optimum is when $ω_{e, 0} Δ t λ_{e, 0} / h \approx 0.3$ ⁠. If in a simulation both the density and temperature of electrons are nonuniform, for a given numerical scheme it is necessary to know, first, how strong the numerical heating or cooling may be when and where the plasma parameters differ from the optimal ones. Second, what can be done to mitigate this effect.

To understand the effect of the direct implicit algorithm on plasmas with different parameters, a set of electrostatic simulations is performed as follows. The electromagnetic part of the algorithm is turned off. The simulation domain is a square filled with plasma only. Periodic boundary conditions are applied at all domain boundaries. The numerical grid is created with the same scale values of $n_{e, 0}$ ⁠, $T_{e, 0}$ ⁠, and $α$ as in Sec. III, therefore, the grid resolution $h$ is the same. The centered grid has 66 nodes along the x and y directions. There are no collisions with neutrals and no external fields. Initial plasma density and temperature are uniform. The ion mass is 40 AMU, the initial ion temperature is $0.03 eV$ ⁠, initial electron and ion velocity distributions are Maxwellian. In the simulation set, the initial plasma density $n_{e}$ is either $1.6 \times 10^{18} m^{- 3}$ or $1.6 \times 10^{16} m^{- 3}$ ⁠. The initial electron temperature $T_{e}$ is 3, 4, or 5 eV. Parameter $β$ defining the time step (32) is 8, 9, or 10. A simulation lasts for $50 ns$ ⁠. The first set of simulations is performed with the number of particles per cell for the scale density $N_{ppc} = 1000$ ⁠. Note that this value is 8 times lower than the value used in Sec. III.

In each simulation of the set, the average kinetic electron energy is, approximately, a linear function of time. The slope of this function represents the rate of average kinetic electron energy change with time. The positive or negative rate corresponds to the numerical heating or cooling, respectively. The rates from the first set are assembled in Table I. The numerical cooling is dominant if $β = 8$ ⁠. The numerical heating is the dominant process for $β = 10$ ⁠. For all values of $β$ ⁠, the rates obtained for the small plasma density are either much larger than the rates obtained for the large density and the same initial temperature, or change the sign. It is necessary to point out that the low density simulations suffer from poor statistics—there are only 10 particles in a cell for $n_{e} = 1.6 \times 10^{16} m^{- 3}$ ⁠.

TABLE I.

The topmost two lines and the leftmost column contain simulation parameters in the set of electrostatic simulations with $N_{ppc} = 1000$ ⁠. The bottom 3 lines, columns 2–7, contain the rate of average kinetic electron energy change (in units of $eV / μ s$ ⁠) for each simulation of the set. The table should be read as follows: for example, the rate value $+ 1.511 eV / μ s$ in the bottom line of column 2 corresponds to simulation with $n_{e} = 1.6 \times 10^{18} m^{- 3}$ ⁠, $T_{e} = 3 eV$ ⁠, $β = 10$ ⁠, etc.

$n_{e} (m^{- 3})$	$1.6 \times 10^{18}$			$1.6 \times 10^{16}$
$T_{e}$ (eV)	3.0	4.0	5.0	3.0	4.0	5.0
$β = 8$	+0.519	−0.261	−1.324	−0.260	−3.697	−8.956
$β = 9$	+1.036	+0.492	−0.415	+3.270	−0.491	−4.067
$β = 10$	+1.511	+1.186	+0.486	+4.816	+2.064	−0.837

$n_{e} (m^{- 3})$	$1.6 \times 10^{18}$			$1.6 \times 10^{16}$
$T_{e}$ (eV)	3.0	4.0	5.0	3.0	4.0	5.0
$β = 8$	+0.519	−0.261	−1.324	−0.260	−3.697	−8.956
$β = 9$	+1.036	+0.492	−0.415	+3.270	−0.491	−4.067
$β = 10$	+1.511	+1.186	+0.486	+4.816	+2.064	−0.837

Simulations with $β = 8$ and $β = 9$ have minimal (in absolute value) energy change rates at $T_{e} = 4 eV$ ⁠, see Table I. If the expected range of electron temperatures is between 3 and 5 eV, one can choose either $β = 8$ or $β = 9$ for the actual simulation if the plasma density is large. The value used in Sec. III is $β = 9$ ⁠. Note, however, that with $N_{ppc} = 1000$ as in the first set, the simulation is affected significantly by numerical cooling in areas with $T_{e} > T_{e, 0}$ or heating in areas with $T_{e} < T_{e, 0}$ at the transitional stage when the plasma density is low.

Experimenting with the code found that the rate of the numerical heating or cooling reduces for larger values of $N_{ppc}$ ⁠. Therefore, another simulation set is performed in the same way as the set described above, but with $N_{ppc} = 8000$ ⁠. The results of the second set are assembled in Table II. The energy change rates in the second set are significantly lower than those in the first set. Moreover, rates in the second set obtained for the low density are often smaller than the rates from the first set obtained for the large density (with the same initial electron temperature). The value $N_{ppc} = 8000$ used in the second set and selected for Sec. III is the largest value, which does not push the numerical cost too far for the computational hardware available to the authors.

TABLE II.

Same as Table I, for the set of electrostatic simulations with $N_{ppc} = 8000$ ⁠.

$n_{e} (m^{- 3})$	$1.6 \times 10^{18}$			$1.6 \times 10^{16}$
$T_{e}$ (eV)	3.0	4.0	5.0	3.0	4.0	5.0
$β = 8$	+0.058	−0.044	−0.191	+0.095	−0.499	−1.116
$β = 9$	+0.130	+0.054	−0.071	+0.482	−0.018	−0.587
$β = 10$	+0.187	+0.145	+0.049	+0.749	+0.285	−0.107

$n_{e} (m^{- 3})$	$1.6 \times 10^{18}$			$1.6 \times 10^{16}$
$T_{e}$ (eV)	3.0	4.0	5.0	3.0	4.0	5.0
$β = 8$	+0.058	−0.044	−0.191	+0.095	−0.499	−1.116
$β = 9$	+0.130	+0.054	−0.071	+0.482	−0.018	−0.587
$β = 10$	+0.187	+0.145	+0.049	+0.749	+0.285	−0.107

APPENDIX E: VERIFICATION OF THE CODE

Simulations of known physical phenomena described in this appendix are performed to verify the code. In all simulations, collisions with neutrals are turned off. Electromagnetic simulations of electron and shear Alfvén waves are performed with the direct solver for $B_{z}$ using Eq. (20). The solver was selected after the energy conservation law analysis similar to the one described in Sec. III A.

1. Two-stream instability

Simulation of the two-stream instability of a cold electron beam in a cold plasma is performed as follows. The electromagnetic module is turned off, the simulation is purely electrostatic. The grid and the time step are calculated as described in Sec. III A with $n_{e, 0} = 10^{16} m^{- 3}$ ⁠, $T_{e, 0} = 1 eV$ ⁠, $α = 4$ ⁠, and $β = 9$ ⁠. The resulting grid cell size is $h = 2.628 \times 10^{- 5} m$ and the time step is $Δ t = 4.924 \times 10^{- 12} s$ ⁠. The number of particles per cell is $N_{ppc} = 9000$ ⁠. The initial plasma density equals to the scale density. The ion mass is 40 AMU, the initial energy of ambient electrons and ions is zero. The electron beam initial energy is $10 eV$ and the beam-to-plasma density ratio is $10^{- 3}$ ⁠. The beam electrons have initial zero velocities in the y and z directions and the velocity in the x-direction $v_{b} = 1.875 529 \times 10^{6} m / s$ ⁠. The rectangular domain has periodic boundary conditions on all sides. The periodicity length along the y-direction is 10 cells. No initial perturbation is introduced. The wavelength of the instability along the x-direction is defined by the periodicity length of the domain in this direction. Five simulations are performed where the periodicity length along the x-direction is 72, 80, 88, 96, and 104 cells. The selected parameters ensure that only one wavelength fits within the system.

In each simulation, after a short transitional period the electrostatic potential shows exponentially growing oscillations with a constant growth rate, see Fig. 25(a). The exponential growth is followed by saturation and oscillation of the amplitude, as expected for the nonlinear stage of the instability.²⁸ The frequency and growth rate of the instability are defined by the dispersion equation²⁹

1 - \frac{ω_{p e}^{2}}{ω^{2}} - \frac{ω_{b e}^{2}}{{(ω - k v_{b})}^{2}} = 0,

(E1)

where the frequency

ω

is a complex number, k is the wave number,

ω_{p e}

and

ω_{b e}

are the electron plasma frequencies of the ambient and beam electrons, respectively. The frequencies and growth rates of the instability obtained in the five simulations mentioned above have a good quantitative match with the linear theory, compare black markers and red curves in Figs. 25(b) and 25(c). This test verifies the electrostatic part of the code.

FIG. 25.

View large Download slide

Simulation of the two-stream instability. (a) Electrostatic potential vs time for $k v_{b} / ω_{e, 0} = 0.993$ ⁠, the red straight line represents exponential growth with the corresponding theoretical growth rate. Panels (b) and (c) show the real (b) and imaginary (c) parts of the frequency vs the wavenumber. In (b) and (c), the red curves are obtained from the dispersion equation (E1) and the black markers are obtained in simulations. The markers, from left to right, correspond to simulations with the periodicity length along the x-direction of 104, 96, 88, 80, and 72 cells.

2. Electromagnetic electron R-wave

Dispersion equation for the high-frequency electromagnetic waves propagating parallel to the magnetic field is³⁰

\frac{c^{2} k^{2}}{ω^{2}} = 1 - \frac{ω_{p e}^{2}}{ω (ω - ω_{c})},

(E2)

where c is the light speed, k is the wave number,

ω

is the frequency,

ω_{p e}

and

ω_{c}

are the electron plasma frequency and cyclotron frequency, respectively. The wave described by has right-hand circular polarization. The wave is transverse, the ratio of electric and magnetic field amplitudes equals the wave phase speed,

v_{p h} = ω / k

⁠. The dispersion curve of such a wave in a plasma with density

10^{15} m^{- 3}

in magnetic field

5 \times 10^{- 4} T

is shown by the red curve in Fig. 26.

FIG. 26.

View large Download slide

Dispersion of the high-frequency electromagnetic plasma R-wave. The cross marks parameters of the wave excited in the simulation are shown in Figs. 27–29.

In simulation described below, an electron R-wave is excited with a given frequency and wave number. The objective of the simulation is to demonstrate that the wave propagates along the ambient magnetic field with the correct phase speed, the wave is right-hand polarized, and the ratio of electric and magnetic field amplitudes equals the phase speed.

The grid size and the time step are calculated using (32) with $n_{e, 0} = 10^{15} m^{- 3}$ ⁠, $T_{e, 0} = 1 eV$ ⁠, $α = 0.05$ ⁠, and $β = 9$ ⁠. The grid cell size is $h = 6.649 \times 10^{- 3} m$ and the time step is $Δ t = 1.246 \times 10^{- 9} s$ ⁠, the number of particles per cell is $N_{ppc} = 500$ ⁠, the ion mass is 40 AMU. Initially, the plasma is uniform, the density equals the scale density, the electron temperature is $1 eV$ ⁠, the ion temperature is $0.03 eV$ ⁠. The simulation domain has 1441 and 241 cells along the x and y directions, respectively, the corresponding length is $L_{x} = 9.58 m$ and width is $L_{y} = 1.6 m$ ⁠. The domain boundaries are grounded metal walls. The ambient magnetic field of $5 \times 10^{- 4} T$ is in the positive x-direction.

The driving antenna is the rectangular contour with the current in the x–y plane, similar to the one used in Sec. III. The driving frequency $ω = 4.397 \times 10^{7} s^{- 1}$ equals to one half of the electron cyclotron frequency, $ω = ω_{c} / 2$ ⁠. The corresponding wave number calculated with (E2) is $k = 5.952 m^{- 1}$ ⁠, the wavelength is $λ = 1.055 m$ ⁠, and the phase speed is $v_{p h} = 7.387 \times 10^{6} m / s$ ⁠. The wave parameters are marked by the black cross in Fig. 26. The size of the antenna along the x-direction equals $0.532 m$ and is close to $λ / 2$ ⁠, the size along the y-direction is $1.037 m$ ⁠, the width of the electric current channel is $0.03324 m$ (5 cells). The center of the antenna is placed on the horizontal symmetry axis $y = L_{y} / 2$ at a distance of $L_{x} / 3$ from the left border, see Fig. 27. The amplitude of the electric current in the antenna is $3 A$ per 1 m of length along the z-direction. The simulation time is $711 ns$ ⁠.

FIG. 27.

View large Download slide

Simulation of the high-frequency electromagnetic plasma R-wave. 2D profiles of $E_{sol, x}$ (a), $E_{sol, y}$ (b), and $E_{sol, z}$ (c) at $t = 699 ns$ ⁠. The black rectangle represents the region with the driving external electric current in the simulation plane. The triangle marks the location of the probe where time dependencies shown in Fig. 29 are obtained.

The electric fields considered below are the solenoidal electric fields. The oscillating current in the antenna excites waves propagating leftward and rightward away from the antenna. The wave propagating leftward reflects off the left border, but the reflected wave does not enter the region rightward from the antenna during the simulated time. Since the system is bounded, the wave is not purely transverse and has the electric field component $E_{sol, x}$ parallel to the ambient magnetic field, see Fig. 27(a). On the horizontal symmetry axis $y = L_{y} / 2$ ⁠, the transverse (relative to the ambient magnetic field) electric field components $E_{sol, y}$ and $E_{sol, z}$ are maximal, see Figs. 27(b) and 27(c), while the parallel electric field is zero. The electromagnetic fields for comparison with the one-dimensional theory are taken on this axis.

Evolution of the solenoidal electric field $E_{sol, y}$ in space and time along the symmetry axis is shown in Fig. 28. In this figure, the colored bands represent waves propagating away from the antenna. The slope of these bands is constant, which means that the waves propagate with constant phase speed. The phase speed estimated from the figure is $7.513 \times 10^{6} m / s$ ⁠, which is very close to the theoretical value given above. Amplitudes of $E_{sol, y}$ and $E_{sol, z}$ are close to each other and phase shift between them is 90 ° corresponding to circular polarization, see Fig. 29(a). The direction of rotation is that of the right-hand polarized wave (same as the direction of rotation of electrons).³⁰ The transverse magnetic field rotates in the same direction as the electric field, amplitudes of $B_{y}$ and $B_{z}$ are very close to each other (circular polarization), and the ratio of electric and magnetic field amplitudes is close to the theoretical phase speed, compare Figs. 29(a) and 29(b). Thus, the code successfully reproduced a high-frequency electromagnetic plasma wave. This verifies electron dynamics in electromagnetic simulation.

FIG. 28.

View large Download slide

Simulation of the high-frequency electromagnetic plasma R-wave. Solenoidal electric field $E_{sol, y}$ vs the x-coordinate and time along the cross section of the system shown in Fig. 27 with $y = 0.798 m$ ⁠.

FIG. 29.

View large Download slide

Simulation of the high-frequency electromagnetic plasma R-wave. Time dependencies of the solenoidal electric field (a) and wave magnetic field (b) in the probe marked by the triangle in Fig. 27. In (a), curves 1 and 2 correspond to $E_{sol, y}$ and $E_{sol, z}$ ⁠, respectively. In (b), curves 1 and 2 correspond to $B_{y}$ and $B_{z}$ ⁠, respectively, multiplied by the theoretical wave phase velocity $v_{p h} = 7.387 \times 10^{6} m / s$ ⁠.

3. Shear Alfvén waves

Dispersion equation for shear Alfvén waves propagating parallel to the magnetic field is³⁰

\frac{c^{2} k^{2}}{ω^{2}} = 1 + \frac{Ω_{p}^{2}}{Ω_{c} (Ω_{c} \pm ω)},

(E3)

where

Ω_{p}

and

Ω_{c}

are the ion plasma frequency and cyclotron frequency, respectively. In the denominator of the last term in , the plus and minus signs correspond to the right- and left-hand polarized waves (or simply R- and L-waves), respectively. An example of dispersion curves of these waves for a plasma with density

10^{17} m^{- 3}

in magnetic field

10^{- 2} T

with ion mass

0.05 AMU

is shown in Fig. 30. The waves selected for simulation have the same wavelength but different frequencies, see the black crosses in Fig. 30. Note that the ion mass is reduced to save the numerical cost of simulations.

FIG. 30.

View large Download slide

Dispersion of the shear Alfvén R-wave (top curve) and L-wave (bottom curve). The straight cross marks parameters of the R-wave are shown in Figs. 31–33. The diagonal cross marks parameters of the L-wave shown in Figs. 34–37.

Simulations of both the R- and L-waves are carried out in the same system. The grid cell size and the time step are calculated using (32) with $n_{e, 0} = 10^{17} m^{- 3}$ ⁠, $T_{e, 0} = 1 eV$ ⁠, $α = 0.025$ ⁠, and $β = 9$ ⁠, resulting in $h = 1.33 \times 10^{- 3} m$ and $Δ t = 2.491 \times 10^{- 10} s$ ⁠. The size of the simulation domain is 4000 cells in the x-direction and 1200 cells in the y-direction corresponding to $L_{x} = 5.319 m$ and $L_{y} = 1.596 m$ ⁠. Due to the large number of cells, the number of particles per cell is reduced to $N_{ppc} = 100$ ⁠. As mentioned above, the ion mass is reduced to $0.05 AMU$ and the ambient magnetic field along the x-direction is $10^{- 2} T$ ⁠. Initially the plasma is uniform, the electron temperature is $1 eV$ ⁠, the ion temperature is $0.1 eV$ ⁠.

The wavelength selected is $λ = 1.064 m$ ⁠. For the R-wave, the frequency is $ω = 1.4971 Ω_{c} = 2.889 \times 10^{7} s^{- 1}$ ⁠, the corresponding wave period is $217.49 ns$ ⁠, the phase velocity is $v_{p h} = 4.892 \times 10^{6} m / s$ ⁠. The driving antenna is the rectangular contour with the current in the x–y plane. The size of the contour along the x-direction equals $0.532 m$ ⁠, which is $λ / 2$ ⁠, the size along the y-direction is $1.065 m$ ⁠, the width of the electric current channel is $0.0545 m$ (41 cells), the left edge of the contour is placed at a distance of $3 λ / 2$ from the left boundary (1200 cells), the center of the contour is on the horizontal symmetry line $y = L_{y} / 2$ ⁠, see Fig. 31. The driving current amplitude is $50 A$ per 1 m along the z-direction, the initial phase is $- π / 2$ ⁠. Since the domain length along the direction of wave propagation $L_{x}$ is only $5 λ$ ⁠, the distance between the right edge of antenna and the right border of the domain is $3 λ$ ⁠, which is why the reflections start distorting the field in the area on the right of the antenna after 3 wave periods. To avoid this, the simulation stops at $t = 660 ns$ ⁠.

FIG. 31.

View large Download slide

Simulation of the shear Alfvén R-wave. 2D profiles of $E_{sol, x}$ (a), $E_{sol, y}$ (b), and $E_{sol, z}$ (c) at $t = 658 ns$ ⁠. The black rectangle represents the region with the driving external electric current in the simulation plane. The triangle marks the location of the probe where time dependencies shown in Fig. 33 are obtained.

Similar to the electron wave simulation in Appendix B, the driving current in the antenna excites waves propagating leftward and rightward from the antenna along the x-direction. The wave has not only the transverse electric field components $E_{sol, y}$ and $E_{sol, z}$ ⁠, but also the component $E_{sol, x}$ parallel to the ambient field, see Fig. 31. On the horizontal symmetry axis, the transverse fields have maxima while the parallel field is zero. The phase speed of propagating waves derived from the evolution of $E_{sol, y}$ along the horizontal symmetry axis shown in Fig. 32 is $5.59 \times 10^{6} m / s$ ⁠, which is 14% larger than the theoretical value. The wave polarization is elliptical rather than circular since the amplitude of $E_{sol, z}$ exceeds the amplitude of $E_{sol, y}$ ⁠, see Fig. 33(a). The wave, however, is right-hand polarized (note the same phase shift between curves 1 in 2 in Figs. 33 and 29), and the ratio of electric and magnetic field amplitudes is close to the phase speed [compare Figs. 33(a) and 33(b)].

FIG. 32.

View large Download slide

Simulation of the shear Alfvén R-wave. Solenoidal electric field $E_{sol, y}$ vs the x-coordinate and time along the cross section of the system shown in Fig. 31 with $y = 0.798 m$ ⁠.

FIG. 33.

View large Download slide

Simulation of the shear Alfvén R-wave. Time dependencies of the solenoidal electric field (a) and wave magnetic field (b) in the probe marked by the triangle in Fig. 31. In (a), curves 1 and 2 correspond to $E_{sol, y}$ and $E_{sol, z}$ ⁠, respectively. In (b), curves 1 and 2 correspond to $B_{y}$ and $E_{z}$ ⁠, respectively, multiplied by the theoretical wave phase velocity $v_{p h} = 4.891 \times 10^{6} m / s$ ⁠.

For the selected wavelength, the frequency of the L-wave is $ω = 0.5995 Ω_{c} = 1.157 \times 10^{7} s^{- 1}$ ⁠, the corresponding wave period is $543.13 ns$ ⁠, the phase velocity is $v_{p h} = 1.959 \times 10^{6} m / s$ ⁠. Multiple external driving currents are required to excite the L-wave. First, there is the contour with the current in the x–y plane, with the same dimensions and same position as the antenna used for the simulation of the shear Alfvén R-wave above, see Figs. 34(a) and 34(b). Second, there are two rectangular regions (patches) where the external oscillating electric current flows along the z-direction: $I_{z}^{ext} (t) = I_{z} sin (ω t + φ_{z})$ ⁠. These patches have dimensions $0.0545 m$ along the x-direction and $1.065 m$ along the y-direction, see the two narrow black rectangles in Fig. 34(c), and are superimposed over the vertical (aligned along y) segments of the contour with the current in the x–y plane. The phases of the currents are adjusted so that the total external electric current vector in areas where the currents are superimposed rotates in the same direction as the field vector in the left-hand polarized wave (the direction of rotation of a positive ion in the ambient magnetic field), see Fig. 35. The amplitude of the current in the x–y plane is $100 A$ ⁠, the amplitude of the currents along the z directions is $106.4 A$ ⁠, these values ensure that the amplitudes of the external electric current density along the y and z directions are the same. The simulation lasts for $1651 ns$ ⁠.

FIG. 34.

View large Download slide

Simulation of the shear Alfvén L-wave. 2D profiles of $E_{sol, x}$ (a), $E_{sol, y}$ (b), and $E_{sol, z}$ (c) at $t = 1644 ns$ ⁠. In (a) and (b), the black rectangle represents the region with the driving external electric current in the simulation plane. In (c), the two black rectangles show the regions with the driving external electric current normal to the simulation plane. The triangle marks the location of the probe where time dependencies shown in Fig. 37 are obtained.

FIG. 35.

View large Download slide

Simulation of the shear Alfvén L-wave. The driving electric current in the simulation plane (curve 1) and normal to the simulation plane (curve 2) vs time. Curve 1 corresponds to the current in the y-direction in the right vertical segment of the contour shown in Figs. 34(a) and 34(b). Curve 2 shows the current in the z-direction in the rightmost external current area shown in Fig. 34(c).

The wave excited in this simulation has the electric field component $E_{sol, x}$ parallel to the ambient magnetic field, see Fig. 34(a). The fields for comparison with the linear theory are taken along the horizontal symmetry axis where $E_{sol, x}$ is zero while $E_{sol, y}$ and $E_{sol, z}$ have maximum. During the first period of the driving rf current, two types of waves are excited. There is a wave with phase speed of about $5.7 \times 10^{6} m / s$ identified as the shear Alfvén R-wave. This wave is seen in Fig. 36 for $x > 3 m$ and $t < 650 ns$ ⁠. There is also a much slower wave, which gradually becomes dominant, see Fig. 36 for $2 m < x < 4 m$ ⁠. The phase speed of this wave is $2.04 \times 10^{6} m / s$ ⁠, which is only 5% larger than the theoretical value for the expected L-wave. The transverse solenoidal electric field shown in Fig. 37(a) is not sinusoidal, probably because it includes the field of the fast R-wave reflected from the left and right domain boundaries. The amplitude of $E_{sol, y}$ is about 30% larger than the amplitude of $E_{sol, z}$ ⁠, so the wave polarization is elliptical rather than circular, compare curves 1 and 2 in Fig. 37(a). Despite these differences, the phase shift between $E_{sol, y}$ and $E_{sol, z}$ corresponds to the left-hand polarized wave, see Fig. 37(a). Also, the ratio of amplitudes of the transverse electric and magnetic fields is close to the phase velocity of the L-wave, compare Figs. 37(a) and 37(b).

FIG. 36.

View large Download slide

Simulation of the shear Alfvén L-wave. Solenoidal electric field $E_{sol, y}$ vs the x-coordinate and time along the cross section of the system shown in Fig. 34 with $y = 0.798 m$ ⁠.

FIG. 37.

View large Download slide

Simulation of the shear Alfvén L-wave. Time dependencies of the solenoidal electric field (a) and wave magnetic field (b) in the probe marked by the triangle in Fig. 34. In (a), curves 1 and 2 correspond to $E_{sol, y}$ and $E_{sol, z}$ ⁠, respectively. In (b), curves 1 and 2 correspond to $B_{y}$ and $B_{z}$ ⁠, respectively, multiplied by the theoretical wave phase velocity $v_{p h} = 1.959 \times 10^{6} m / s$ ⁠.

Differences between the waves obtained in simulations and theoretical shear Alfvén R- and L-waves, such as the elliptical polarization rather than the circular one and the slightly higher phase speed, are minor and can be attributed to the presence of boundaries in the system, relatively low number of particles per cell (due to high numerical cost) increasing the electrostatic noise, thermal effects. With these reservations, one can conclude that the code successfully reproduced left- and right-hand polarized shear Alfvén waves. This verifies ion dynamics in electromagnetic simulation.

REFERENCES

1.

C. K.

Birdsall

and

A. B.

Langdon

,

Plasma Physics via Computer Simulations

(

IOP Publishing

,

Bristol/Philadelphia

,

1991

).

Google Scholar

Crossref

2.

G.

Chen

,

L.

Chacon

, and

D. C.

Barnes

, “

An energy- and charge-conserving, implicit, electrostatic particle-in-cell algorithm

,”

J. Comput. Phys.

230

,

7018

–

7036

(

2011

).

https://doi.org/10.1016/j.jcp.2011.05.031

Google Scholar

Crossref

3.

J. R.

Angus

,

W.

Farmer

,

A.

Friedman

,

V.

Geyko

,

D.

Ghosh

,

D.

Grote

,

D.

Larson

, and

A.

Link

, “

An implicit particle code with exact energy and charge conservation for studies of dense plasmas in axisymmetric geometries

,”

J. Comput. Phys.

519

,

113427

(

2024

).

https://doi.org/10.1016/j.jcp.2024.113427

Google Scholar

Crossref

4.

A. B.

Langdon

,

B. I.

Cohen

, and

A.

Friedman

, “

Direct implicit large time-step particle simulation of plasmas

,”

J. Comput. Phys.

51

,

107

–

138

(

1983

).

https://doi.org/10.1016/0021-9991(83)90083-9

Google Scholar

Crossref

5.

B. I.

Cohen

,

A. B.

Langdon

, and

D. W.

Hewett

, “

Performance and optimization of direct implicit particle simulations

,”

J. Comput. Phys.

81

,

151

–

168

(

1989

).

https://doi.org/10.1016/0021-9991(89)90068-5

Google Scholar

Crossref

6.

K.

Yee

, “

Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media

,”

IEEE Trans. Antennas Propag.

14

,

302

–

307

(

1966

).

https://doi.org/10.1109/TAP.1966.1138693

Google Scholar

Crossref

7.

S.

Mattei

,

K.

Nishida

,

M.

Onai

,

J.

Lettry

,

M. Q.

Tran

, and

A.

Hatayama

, “

A fully-implicit particle-in-cell monte carlo collision code for the simulation of inductively coupled plasmas

,”

J. Comput. Phys.

350

,

891

–

906

(

2017

).

https://doi.org/10.1016/j.jcp.2017.09.015

Google Scholar

Crossref

8.

C. G.

Darwin

, “

The dynamical motions of charged particles

,”

London, Edinburgh, Dublin Philos. Mag. J. Sci.

39

,

537

–

551

(

1920

).

https://doi.org/10.1080/14786440508636066

Google Scholar

Crossref

9.

C. W.

Nielson

and

H. R.

Lewis

, “

Particle-code models in the nonradiative limit

,”

Methods Comput. Phys.: Adv. Res. Appl.

16

,

367

–

388

(

1976

).

https://doi.org/10.1016/B978-0-12-460816-0.50015-4

Google Scholar

Crossref

10.

W.

Lee

,

R. C.

Davidson

,

E. A.

Startsev

, and

H.

Qin

, “

The electromagnetic Darwin model for intense charged particle beams

,”

Nucl. Instrum. Methods Phys. Res. A

544

,

353

–

359

(

2005

).

https://doi.org/10.1016/j.nima.2005.01.233

Google Scholar

Crossref

11.

D.

Eremin

,

T.

Hemke

,

R. P.

Brinkmann

, and

T.

Mussenbrock

, “

Simulations of electromagnetic effects in high-frequency capacitively coupled discharges using the Darwin approximation

,”

J. Phys. D: Appl. Phys.

46

,

084017

(

2013

).

https://doi.org/10.1088/0022-3727/46/8/084017

Google Scholar

Crossref

12.

D. C.

Barnes

, “

Time-explicit Darwin PIC algorithm

,”

J. Comput. Phys.

462

,

111151

(

2022

).

https://doi.org/10.1016/j.jcp.2022.111151

Google Scholar

Crossref

13.

M. R.

Gibbons

and

D. W.

Hewett

, “

The Darwin direct implicit particle-in-cell (DADIPIC) method for simulation of low frequency plasma phenomena

,”

J. Comput. Phys.

120

,

231

–

247

(

1995

).

https://doi.org/10.1006/jcph.1995.1160

Google Scholar

Crossref

14.

D. W.

Hewett

,

D. J.

Larson

, and

S.

Doss

, “

Solution of simultaneous partial differential equations using dynamic ADI: Solution of the streamlined Darwin field equations

,”

J. Comput. Phys.

101

,

11

–

24

(

1992

).

https://doi.org/10.1016/0021-9991(92)90039-2

Google Scholar

Crossref

15.

S.

Balay

,

S.

Abhyankar

,

M. F.

Adams

,

S.

Benson

,

J.

Brown

,

P.

Brune

,

K.

Buschelman

,

E. M.

Constantinescu

,

L.

Dalcin

,

A.

Dener

,

V.

Eijkhout

,

J.

Faibussowitsch

,

W. D.

Gropp

,

V.

Hapla

,

T.

Isaac

,

P.

Jolivet

,

D.

Karpeev

,

D.

Kaushik

,

M. G.

Knepley

,

F.

Kong

,

S.

Kruger

,

D. A.

May

,

L. C.

McInnes

,

R. T.

Mills

,

L.

Mitchell

,

T.

Munson

,

J. E.

Roman

,

K.

Rupp

,

P.

Sanan

,

J.

Sarich

,

B. F.

Smith

,

S.

Zampini

,

H.

Zhang

,

H.

Zhang

, and

J.

Zhang

, see https://petsc.org/ for “

PETSc Web page

” (

2025

).

Google Scholar

16.

hypre

, see https://llnl.gov/casc/hypre, https://github.com/hypre-space/hypre for “

hypre: High performance preconditioners

” (2025).

17.

J.

Amorim

,

H. S.

Maciel

, and

J. P.

Sudano

, “

High-density plasma mode of an inductively coupled radio frequency discharge

,”

J. Vac. Sci. Technol. B

9

,

362

–

365

(

1991

).

https://doi.org/10.1116/1.585576

Google Scholar

Crossref

18.

J.

Hopwood

, “

Review of inductively coupled plasmas for material processing

,”

Plasma Sources Sci. Technol.

1

,

109

–

116

(

1992

).

https://doi.org/10.1088/0963-0252/1/2/006

Google Scholar

Crossref

19.

T.

Okumura

,

“Inductively coupled plasma sources and applications

,”

Phys. Res. Int.

2010

,

164249

.

https://doi.org/10.1155/2010/164249

Crossref

20.

S. A.

Maiorov

, “

Ion drift in a gas in an external electric field

,”

Plasma Phys. Rep.

35

,

802

–

812

(

2009

).

https://doi.org/10.1134/S1063780X09090098

Google Scholar

Crossref

21.

D.

Sydorenko

, “

Particle-in-cell simulations of electron dynamics in low pressure discharges with magnetic fields

,” Ph.D. thesis (

University of Saskatchewan

,

2006

).

Google Scholar

22.

D.

Sydorenko

and

A.

Khrabrov

, see https://github.com/PrincetonUniversity/EDIPIC for “

Edipic-1d Web page

” (

2021

).

Google Scholar

23.

D.

Sydorenko

,

A.

Khrabrov

,

W.

Villaphana

,

A.

Powis

, and

S.

Ethier

, see https://github.com/PrincetonUniversity/EDIPIC-2D for “

Edipic-2d Web page

” (

2021

).

Google Scholar

24.

E. A.

Startsev

,

R. C.

Davidson

, and

H.

Qin

, “

Numerical studies of the electromagnetic Weibel instability in intense charged particle beams with large temperature anisotropy using the nonlinear best Darwin delta-f code

,” in

2007 IEEE Particle Accelerator Conference (PAC)

(

IEEE

,

2007

), pp.

4297

–

4299

.

Google Scholar

Crossref

25.

V. I.

Kolobov

and

D. J.

Economou

, “

The anomalous skin effect in gas discharge plasmas

,”

Plasma Sources Sci. Technol.

6

,

R1

–

R17

(

1997

).

https://doi.org/10.1088/0963-0252/6/2/003

Google Scholar

Crossref

26.

V. A.

Godyak

,

B. M.

Alexandrovich

, and

V. I.

Kolobov

, “

Lorentz force effects on the electron energy distribution in inductively coupled plasmas

,”

Phys. Rev. E

64

,

026406

(

2001

).

https://doi.org/10.1103/PhysRevE.64.026406

Google Scholar

Crossref

27.

H.

Sun

,

S.

Banerjee

,

S.

Sharma

,

A. T.

Powis

,

A. V.

Khrabrov

,

D.

Sydorenko

,

J.

Chen

, and

I. D.

Kaganovich

, “

Direct implicit and explicit energy-conserving particle-in-cell methods for modeling of capacitively coupled plasma devices

,”

Phys. Plasmas

30

,

103509

(

2023

).

https://doi.org/10.1063/5.0160853

Google Scholar

Crossref

28.

N. G.

Matsiborko

,

IN.

Onishchenko

,

V. D.

Shapiro

, and

V. I.

Shevchenko

, “

On non-linear theory of instability of a mono-energetic electron beam in plasma

,”

Plasma Phys.

14

,

591

–

600

(

1972

).

https://doi.org/10.1088/0032-1028/14/6/003

Google Scholar

Crossref

29.

A. B.

Mikhailovskii

,

Theory of Plasma Instabilities

(

Consultants Bureau

,

New York

,

1974

), Vol.

1

.

Google Scholar

Crossref

30.

R. J.

Goldston

and

P. H.

Rutherford

,

Introduction to Plasma Physics

(

Institute of Physics Publishing

,

Bristol/Philadelphia

,

1995

).

Google Scholar

Crossref

© 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC) license (https://creativecommons.org/licenses/by-nc/4.0/).

This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License .

Simulation of an inductively coupled plasma with a two-dimensional Darwin particle-in-cell code

I. INTRODUCTION

II. MODEL DESCRIPTION

A. Electromagnetic equations in Darwin method

B. Solenoidal electric field

C. Self-consistent magnetic field

D. External drivers

E. Energy conservation

III. SIMULATION OF AN INDUCTIVELY COUPLED PLASMA

A. System setup and numerical parameters

B. Time dependencies of global parameters

C. System evolution during the initial/transitional stage

D. Plasma state at the late stage of the simulation

IV. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: ELECTROSTATIC DIRECT IMPLICIT ALGORITHM

APPENDIX B: COEFFICIENTS FOR THE DARWIN ALGORITHM

APPENDIX C: FINITE DIFFERENCE FORM OF EQUATION FOR F

APPENDIX D: SELECTION OF NUMERICAL PARAMETERS

APPENDIX E: VERIFICATION OF THE CODE

1. Two-stream instability

2. Electromagnetic electron R-wave

3. Shear Alfvén waves

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Simulation of an inductively coupled plasma with a two-dimensional Darwin particle-in-cell code Open Access

I. INTRODUCTION

II. MODEL DESCRIPTION

A. Electromagnetic equations in Darwin method

B. Solenoidal electric field

C. Self-consistent magnetic field

D. External drivers

E. Energy conservation

III. SIMULATION OF AN INDUCTIVELY COUPLED PLASMA

A. System setup and numerical parameters

B. Time dependencies of global parameters

C. System evolution during the initial/transitional stage

D. Plasma state at the late stage of the simulation

IV. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: ELECTROSTATIC DIRECT IMPLICIT ALGORITHM

APPENDIX B: COEFFICIENTS FOR THE DARWIN ALGORITHM

APPENDIX C: FINITE DIFFERENCE FORM OF EQUATION FOR F

APPENDIX D: SELECTION OF NUMERICAL PARAMETERS

APPENDIX E: VERIFICATION OF THE CODE

1. Two-stream instability

2. Electromagnetic electron R-wave

3. Shear Alfvén waves

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Simulation of an inductively coupled plasma with a two-dimensional Darwin particle-in-cell code