Species

This page describes the different models that hPIC2 can use for a plasma species.

Full-orbit particles

The basic particle model of hPIC2 is a classical non-relativistic full-orbit Lagrangian particle of mass \(m\) and charge \(q\), with dynamics described by the classical Newton-Lorentz equation,

\[m \frac{d^2 \mathbf{x} }{dt^2} = q \mathbf{E} + q \mathbf{v} \times \mathbf{B}\]

where the electric field \(\mathbf{E}\) and the magnetic field \(\mathbf{B}\) are either calculated self-consistently by solving the Maxwell Equations, or provided externally by the user (eg. an imposed magnetic field).

The Newton-Lorentz equation is integrated by means of the Boris-Bunemann numerical method, with linearized tangent.

Simulation options in hPIC2 are specified in this section of the manual.

Fig. 1 Full-orbit dynamics of a single charged particle

Boltzmann electrons

On ion-transport time scales, the electron behavior can in first approximation be described simply considering a balance between electrostatic forces and pressure forces on an isothermal fluid: \(-k_B T_e \nabla n_e + e n_e\nabla \phi \approx 0\), with usual meaning of symbols as in Chen. Integrating the balance of forces leads to a relation between the electron particle density and the plasma potential in the form of equation

(1)\[\begin{aligned} n_e(\mathbf{x})=n_0 \exp( e \phi(\mathbf{x})/k_B T_e),\label{botlzmann.equation1} \end{aligned}\]

where \(n_0\) is the reference electron density corresponding to \(\phi=0\). Boltzmann electrons hold an advantage in terms of computational cost over the alternative approximations used in PIC simulations. While alternative methods capture the physical phenomena of electron motion to a higher degree of accuracy, the added simulation complexity makes it computationally expensive to run large timescale simulations. A description of how to use Boltzmann electrons in hPIC2 is provided here.

Time advancement schemes calculate unknown time-dependent variables at time \(t^{k+1} = t^k + \Delta t\) from known variables at time \(t^k\). Common time advancement algorithm in PIC codes calculates the ion density \(n_i^{k+1}\) using plasma potential \(\phi^k\). Subsequently, the plasma potential \(\phi^{k+1}\) is solved using the newly calculated ion density \(n_i^{k+1}\) and equation (1), i.e,;

(2)\[\begin{split}\begin{aligned} \epsilon_0 \nabla^2\phi^{k+1}(\mathbf{x})&=-\rho^{k+1}(\mathbf{x})\label{poisson.equation}\\ &=en_e^{k+1}(\mathbf{x})-en_i^{k+1}(\mathbf{x})\label{poisson.equation1}\\ &=en_0^{k+1} \exp(\phi^{k+1}(\mathbf{x})/T_e)-en_i^{k+1}(\mathbf{x})\label{poisson.equation2}. \end{aligned}\end{split}\]

Equation (2) can be solved using Newton-Raphson, or other methods, to calculate the plasma potential for the next iteration. Problems arise when the reference electron density \(n_0\) varies with time as is the case in the presence of a volumetric source/loss, or a boundary flux. A self-consistent numerical scheme to calculate \(n_0^{k+1}\) is required to maintain charge conservation. Breaking charge conservation leads to numerical oscillations and simulation divergence.

The adoption of Boltzmann electrons always require to enforce charge conservation through a dedicated scheme. Details of the charge conservation scheme are described in the paper Elias and Curreli [EC20], and are briefly described below. This is the "elias" option for the charge_conservation_scheme in hPIC2 input decks.

The charge conservation scheme is derived from the Ampere-Maxwell equation in differential form,

(3)\[\begin{aligned} \nabla \times \mathbf{B}&= \mu_0 \mathbf{J} + \epsilon_0 \mu_0 \frac{\partial \mathbf{E}}{\partial t}\label{max.equation1} \end{aligned}\]

As usual, local charge conservation is obtained by taking the divergence of equation (3) and calling the displacement current as \(\mathbf{J_D}=\epsilon_0 \frac{\partial \mathbf{E}}{\partial t}\)

(4)\[\begin{split}\begin{aligned} \nabla \cdot (\nabla \times \mathbf{B})&= \mu_0 \nabla \cdot \mathbf{J} +\mu_0 \nabla \cdot \left( \epsilon_0 \frac{\partial \mathbf{E}}{\partial t} \right) \label{max.equation2}\\ 0 &=\nabla \cdot \mathbf{J} + \nabla \cdot \mathbf{J}_D, \label{globalcharge.equation1} \end{aligned}\end{split}\]

where the conduction current \(\mathbf{J}=\mathbf{J}_i + \mathbf{J}_e\) is the sum of the contributions from the ion current \(\mathbf{J}_i\) and the electron current \(\mathbf{J}_e\). Equation (4) can equivalently be expressed as

\[\begin{aligned} \nabla \cdot (\mathbf{J}_e + \mathbf{J}_i + \mathbf{J}_D)&=0 \label{globalcharge.equation2} \end{aligned}\]

or using its integral form,

(5)\[\begin{aligned} \int_V \nabla \cdot (\mathbf{J}_e + \mathbf{J}_i + \mathbf{J}_D) dV&= 0 \label{displacemen.equation1} \end{aligned}\]

In the presence of volumetric source \(G\) and loss \(L\) terms, equation (5) becomes

(6)\[\begin{aligned} \int_V \nabla \cdot (\mathbf{J}_e + \mathbf{J}_i + \mathbf{J}_D) dV&= G-L \label{displacemen.equation2} \end{aligned}\]

The Boltzmann electron model described in equation (1) implicitly assumes the electron distribution is at a Maxwellian thermal equilibrium. For a Maxwellian thermal distribution, with a mean thermal electron velocity \(\mathbf{u_e}=\sqrt{\frac{8 K_b T_e}{\pi m_e}}\), the current density at the location \(\mathbf{x}\) can, as in Chen, be expressed as

(7)\[\begin{aligned} \mathbf{J}_e(\mathbf{x})=-e \boldsymbol{\Gamma}_e(\mathbf{x})=-e n_0 \mathbf{u}_e \exp(e\Phi(\mathbf{x})/T_e) \label{boundaryflux} \end{aligned}\]

By substituting Equation (7) into Equation (6) and solving for \(n_0\), immediately yields an expression for the reference Boltzmann electron density \(n_0\)

(8)\[\begin{aligned} n_0= \frac{\int_V \nabla \cdot (\mathbf{J}_i + \mathbf{J}_D) dV - G + L }{\int_V \nabla \cdot e \mathbf{u}_e \exp(e\Phi(\mathbf{x})/T_e) dV} \label{density_update} \end{aligned}\]

Equation (8) can be directly used to enforce global charge conservation in explicit PIC schemes with Boltzmann electrons. An example algorithm is discussed hereafter.

A simple explicit algorithm implementing Equation (8) for updating the Boltzmann density \(n_0\) from time step \(t^{k}\) to time step \(t^{k+1}\) is as follows.

1. Calculate ion density \(n_i^{k+1}\) using the plasma potential \(\phi^k\) at the previous time step, using the classical explicit PIC scheme;

2. Calculate reference Boltzmann electron density at \(n_0^{k+1}\) at time step \(t^{k+1}\) using equation (8) and boundary conditions for \(\phi^{k+1}\);

   \[ \begin{align}\begin{aligned}\begin{aligned} n_0^{k+1}= \frac{\int_V \nabla \cdot (\mathbf{J}_i^{k+1} + \mathbf{J}_D^{k}) dV - G^{k+1} + L^{k+1} }{\int_V \nabla \cdot e \mathbf{u_e} \exp(e\phi^{k+1}/T_e) dV} \label{density_update1}\\\end{aligned}\end{aligned}\end{align} \]

3. Solve the plasma potential \(\phi^{k+1}\) using ion density \(n_i^{k+1}\), boundary conditions for \(\phi^{k+1}\), the Poisson equation and reference Boltzmann electron reference density \(n_0^{k+1}\).

The algorithm can be equally applied to plasma domains of arbitrary dimensionality in 1D, 2D or 3D without any loss of accuracy. However, the conventional Courant–Friedrichs–Lewy (CFL) condition on the time step remains necessary to ensure accuracy on the particle pusher, and to resolve ion-timescale phenomena. In the next section we apply this algorithm to two cases, a steady-state plasma sheath and a radio-frequency plasma sheath.

Euler fluid

Many important physical quantities can be computed as moments of a distribution in velocity space. The number density \(n = n(\vec{x}, t)\) of a species described by the distribution \(f\) can be computed as

\[n = \int_{\mathbb{R}^3} f \, \mathrm{d} \vec{v};\]

the momentum density \(n m \vec{u} = n m \vec{u}(\vec{x}, t)\) is

\[n m \vec{u} = \int_{\mathbb{R}^3} m \vec{v} f \, \mathrm{d} \vec{v};\]

the stress tensor \(P_{ij} = P_{ij} (\vec{x}, t)\) is

\[P_{ij} = \int_{\mathbb{R}^3} m v_i v_j f \, \mathrm{d} \vec{v};\]

and the energy flux density \(\vec{Q} = \vec{Q}(\vec{x}, t)\) is

\[\vec{Q} = \int_{\mathbb{R}^3} \frac{1}{2} m v^2 \vec{v} f \, \mathrm{d} \vec{v}.\]

It is also useful to name some moments in the reference frame of the moving species. With \(\vec{w} = \vec{v} - \vec{u}\), let

\[p_{ij} = \int_{\mathbb{R}^3} m w_i w_j f \, \mathrm{d} \vec{v}\]

be the pressure tensor, and let

\[\vec{q} = \int_{\mathbb{R}^3} \frac{1}{2} m w^2 \vec{w} f \, \mathrm{d} \vec{v}\]

be the heat flux density. For convenience, let \(p = p_{ii}/3\) be the scalar pressure and decompose the pressure tensor as

\[p_{ij} = p \delta_{ij} + \pi_{ij},\]

where \(\pi_{ij}\) is the generalized viscosity tensor. Finally,

\[n m E = \int_{\mathbb{R}^3} H f \, \mathrm{d} \vec{v}\]

with the single-particle Hamiltonian \(H = \frac{1}{2} m v^2\) is the total energy density.

The Euler equations can be derived from the Boltzmann kinetic equation by computing moments as

\[\int_{\mathbb{R}^3} \psi \left[ \frac{\partial f}{\partial t} + \vec{v} \cdot \frac{\partial f}{\partial \vec{x}} + \frac{q}{m} \left( \vec{E} + \vec{v} \times \vec{B} \right) \cdot \frac{\partial f}{\partial \vec{v}} \right] \, \mathrm{d} \vec{v} = \int_{\mathbb{R}^3} \psi \mathcal{C} [f] \, \mathrm{d} \vec{v},\]

where \(\psi = \psi(\vec{v})\) is a polynomial. In particular, take \(\psi = m\), \(m \vec{v}\), and \(\frac{1}{2} m v^2\). This ultimately yields

\[ \begin{align}\begin{aligned}\frac{\partial}{\partial t} (nm) + \nabla \cdot (n m \vec{u}) = \int_{\mathbb{R}^3} m \mathcal{C}[f] \, \mathrm{d} \vec{v},\\\frac{\partial}{\partial t} (nmu_i) + \frac{\partial}{\partial x_j} P_{ij} - q n (\vec{E} + \vec{u} \times \vec{B})_i = \int_{\mathbb{R}^3} m \vec{v} \mathcal{C}[f] \, \mathrm{d} \vec{v},\\\frac{\partial}{\partial t} (nmE) + \nabla \cdot \left(nmE \vec{u} + \vec{q} + p \vec{u} + \pi_{ij} u_j \right) - q n \vec{u} \cdot \vec{E} = \int_{\mathbb{R}^3} H \mathcal{C}[f] \, \mathrm{d} \vec{v}.\end{aligned}\end{align} \]

These equations are closed by assuming that the heat flux density and generalized viscosity tensor are zero and relating the scalar pressure to the remaining fluid state variables through an equation of state (EOS), resulting in

\[ \begin{align}\begin{aligned}\frac{\partial}{\partial t} (nm) + \nabla \cdot (n m \vec{u}) = \int_{\mathbb{R}^3} m \mathcal{C}[f] \, \mathrm{d} \vec{v},\\\frac{\partial}{\partial t} (nmu_i) + \frac{\partial}{\partial x_j} \left( n m u_i u_j + p \delta_{ij} \right) = q n (\vec{E} + \vec{u} \times \vec{B})_i + \int_{\mathbb{R}^3} m \vec{v} \mathcal{C}[f] \, \mathrm{d} \vec{v},\\\frac{\partial}{\partial t} (nmE) + \nabla \cdot \left(nmE \vec{u} + p \vec{u} \right) = q n \vec{u} \cdot \vec{E} + \int_{\mathbb{R}^3} H \mathcal{C}[f] \, \mathrm{d} \vec{v}.\end{aligned}\end{align} \]

A common analytic EOS is the ideal gas law

\[p = n k T,\]

where \(k\) is the Boltzmann constant and \(T\) is the temperature, combined with the equipartition theorem for calorically perfect gases

\[n m E = \frac{1}{2} n m u^2 + \frac{1}{\gamma - 1} n k T,\]

which yields

\[p = (\gamma - 1) \left( n m E - \frac{1}{2} n m u^2 \right).\]

The specific heat ratio \(\gamma\) can be specified in hPIC2 using the gamma option.

Further description of the discretization is provided in Fluids. This section of the user manual described how a fluid may be used in hPIC2.

Uniform background

This model assumes that the species follows a Maxwellian distribution everywhere in space, so that the distribution is

\[f = n \sqrt{\frac{m}{2 \pi k T}} \exp \left( - \frac{m v^2}{2 k T} \right)\]

for a given number density \(n\) and temperature \(T\). The charge density is therefore simply \(\rho = q n\). A description of how to use uniform backgrounds in hPIC2 is provided here.