The AMR Hybrid Particle-In-Cell formalism

From hybrid equations to AMR algorithms

The Hybrid formalism consists in modeling the plasma as a combination of constituants treated with a different physical models. This usually means that ions and electrons are treated differently. A rather complete description of the different ways a code can be "hybrid" is given in The Hybrid Multiscale Simulation Technology by A.S. Lipatov.
In astrophysical and space applications, the main application domains of PHARE, "hybrid" usually means that ions are considered at the kinetic level while electrons are considered as a fluid. This is the case for PHARE and this is what we mean by "hybrid" on this page.

The ion equations

The hybrid model consists in evolving in space \( r \) and time \(t\) the shape of the velocity distribution function \( f_p \) of each ion populations p under the influence of the electric \(E\) and magnetic field \(B\). This is done by solving the Vlasov equation for all ion populations when collisions are negligible.

\begin{equation} \frac{\partial f_p}{\partial t} + \mathbf{v}\cdot \frac{\partial f_p}{\partial \mathbf{r}} + \frac{q_p}{m_p}(\mathbf{\mathbf{E} + \mathbf{v}\times\mathbf{B}})\cdot \frac{\partial f_p}{\partial \mathbf{v}} = 0 \label{eq:vlasov} \end{equation}

Having the new distribution everywhere at \(t+\Delta t\), it is easy to calculate the ion moments as the sum of the moments of all populations. Namely, for the ion density \(n_i\) and bulk velocity \(\mathbf{u_i}\)

\begin{eqnarray} n_i(\mathbf{r},t) &= & \sum_p \int f_p(\mathbf{r}, \mathbf{v}, t) d\mathbf{v} \label{eq:density}\\ \mathbf{u}_i(\mathbf{r},t) &= & \frac{1}{n_i}\sum_p \int \mathbf{v} f_p(\mathbf{r}, \mathbf{v}, t) d\mathbf{v} \label{eq:bulk}\\ \end{eqnarray}

Electron momentum equation

What about the electrons? Remember? They are assumed to behave as a fluid. This is wrong of course in collisionless systems since nothing makes the density, velocity etc. of the electrons to depend on purely local physics as collisions would in a "real" fluid. But that's an approximation the hybrid formalism makes to simplify the physics compared to the fully kinetic system and that is already a much more realistic way of modeling the plasma and say, single fluid magnetohydrodynamics. Now there are subtleties. The electron momentum equation is:

\begin{equation} m_en_e \frac{d\mathbf{u_e}}{dt} = -\nabla\cdot \mathbf{P_e} - e n_e(\mathbf{E} + \mathbf{u_e}\times\mathbf{B}) \end{equation}

Electromagnetic field equations : Maxwell-Faraday and Ohm's law

"Treating electrons as a fluid", you probably think we solve that equation, in contrast to the Vlasov equation we used for ions. Well not really... But let's say we did. Now we would have to wonder where the magnetic field and electric field would come from. For the magnetic field, the answer is easy. We just use the Maxwell-Faraday equation:

\begin{equation} \frac{\partial \mathbf{B}}{\partial t} = -\nabla\times\mathbf{E} \label{eq:electronmomentum} \end{equation}

What about the electric field now? There is all the trick of Hybrid codes. We actually do not solve \ref{eq:electronmomentum} directly to get the new electron fluid momentum. Instead we make assumptions on the electron fluid, and use that momentum equation to calculate the electric field ! Thus, \ref{eq:electronmomentum} is re-written:

\begin{equation} \mathbf{E} = -\mathbf{u_e}\times\mathbf{B} - \frac{1}{en_e}\nabla\cdot \mathbf{P_e} +\frac{m_e}{e}\frac{d\mathbf{u_e}}{dt} \label{eq:ohmelectron} \end{equation}

Quasi-neutrality

At this point, the equation for the electric field still has unknowns. The most obvious perhaps is \(n_e\) the electron particle density. This is where the hybrid formalism makes the assumption that at the scale we solve the equations, the plasma is quasineutral, and thus we can neglect the difference between \(n_i\) and \(n_e\) and have only one variable \(n\): the plasma density. Since we have the total ion density already, that's our \(n\). Quasineutrality enable us to get the electron bulk velocity from the known ion bulk velocity and the electric current density:

\begin{equation} \mathbf{u_e} = \mathbf{u_i} - \frac{\mathbf{j}}{en} \end{equation}

The total current density is obtained from the static Maxwell-Ampere equation, meaning we neglect the displacement current:

\begin{equation} \mu_0 \mathbf{j} = \nabla\times\mathbf{B} \end{equation}

Equation \ref{eq:ohmelectron} now reads:

\begin{equation} \mathbf{E} = -\mathbf{u_e}\times\mathbf{B} - \frac{1}{en}\nabla\cdot \mathbf{P_e} +\frac{m_e}{e}\frac{d\mathbf{u_e}}{dt} \label{eq:ohmelectron2} \end{equation}

Massless electrons

The next assumption usually made in Hybrid codes, that is also made in PHARE, is that the spatial and time scales at which we are interested in are much larger and longer that the scales at which the electron bulk inertia matters. The electrons being so light compare to even protons, that it is mostly ok to neglect the last term of \ref{eq:ohmelectron2}, which now reads:

\begin{equation} \mathbf{E} = -\mathbf{u_e}\times\mathbf{B} - \frac{1}{en}\nabla\cdot \mathbf{P_e} \label{eq:ohmelectron3} \end{equation}

Electron closure

Since we do not have an electron distribution function in hand, the last term of eq. \ref{eq:ohmelectron3} is not known a priori. Hybrid codes thus have to come with a so-called closure equation which role is to give us the pressure everywhere at time t, based on some assumption on the system. Usually, unless in very specific conditions, there is no rigorous way of getting such equation and most hybrid code assume a closure that is "reasonable" and above all "simple" to calculate. Perhaps the simplest and most used electron closure is the isothermal one. This simply say that the electron pressure \(P_e\) is given by the product of the density by some scalar constant that we call "the electron temperature".

\begin{equation} P_e= nT_e \label{eq:isothermal} \end{equation}

Dissipative terms

Solved as is, \ref{eq:ohmelectron3} and \ref{eq:isothermal} would result in current sheets to collapse at grid scale in the absence of an intrinsic dissipation scale in the system. Too ways are typically employed in Hybrid codes to include such a dissipation. Joule resistivity well known to be used already in MHD codes. It is a simple term \(\eta \mathbf{j}\) to add on the right hand side of equation \ref{eq:ohmelectron3}. This term adds diffusion of magnetic flux. However there is no scale at which this terms dominate over the electron ideal term \(-\mathbf{u_e}\times\mathbf{B}\), unless \(\eta\) is so large that ion scale structures are diffused away too.

Another term that can be employed is the so-called hyper-resistivity (sometimes called hyper-viscosity) that takes the form \(-\nu\nabla^2\mathbf{j}\) In contrast to classical resistivity, this terms (due to the second order derivative) comes with an intrinsic scale at which it is dominant over electron convection term and efficiently adds sub-ion scale dissipation.

PHARE include these two terms and eq. \ref{eq:ohmelectron3} that is solved in PHARE writes : \begin{equation} \mathbf{E} = -\mathbf{u_e}\times\mathbf{B} - \frac{1}{en}\nabla P_e +\eta\mathbf{j} - \nu\nabla^2\mathbf{j} \label{eq:ohmelectron4} \end{equation}

Numerical integration

The Particle-In-Cell formalism

There are two ways to solve equations \ref{eq:vlasov} for ion populations. It can be solved calculating eulerian derivatives, i.e. discretizing velocity and spatial dimensions and solving the equation at those fixed location. This is called a "Vlasoc Hybrid code". It is generally complex and require lots of computational resources. The other way consists in adopting a Lagrangian standpoint. That is, cutting the initial distribution function in N weighted bins and follow the dynamic of those chunks in phase space. The little pieces of distributions are called "macro-particles". Decomposing the distribution function \(f_p(\mathbf{r},\mathbf{v}, t)\) of the population \(p\) into the contribution of \(N_p\) macro-particles in the following way is the base of the so-called "Particle-in-Cell" (PIC) method

\begin{equation} f_p(\mathbf{r}, \mathbf{v}, t) = \sum_m^{N_p} w_m S(\mathbf{r} - \mathbf{r}_m(t))\delta(\mathbf{v}-\mathbf{v}_m(t)) \end{equation}

where \(\mathbf{r}_m\) and \(\mathbf{v}_m\) are the position and velocity of the \(m_{th}\) macro-particle. \(w_m\) represents the weight of that macro-particle, i.e. how much it counts in the evaluation of \(f_p\). \(\delta\) is the Dirac function, which says that a macro-particle represent a specific velocity in the distribution. In contrast, the function \(S\) is a finite support function representing the "shape" of the macro-particle in the spatial dimension. This function tells us how far from the macro-particle a local distribution sees its influence. In PHARE we use b-spline functions to model \(S\). PHARE uses b-splines of the first, second and third order. The higher the order the further a macro-particle influences the distribution, but the longer it takes to compute it.

$$ S_1(x) = \left\{ \begin{array}{ll} 1-|x| & \mbox{si } |x|\le 1 \\ 0& \mbox{sinon.} \end{array} \right. $$

$$ S_2(x) = \left\{ \begin{array}{ll} 3/4-x^2 & \mbox{si } |x|\le 1/2 \\ 1/2(3/2-|x|)^2 & \mbox{si } 1/2 \lt |x| \le 3/2 \\ 0 & \mbox{si } |x| \gt 3/2 \end{array} \right. $$

$$ S_3(x) = \left\{ \begin{array}{ll} 1/2|x|^3 -x^2 + 2/3& \mbox{si } |x|\le 1 \\ 4/3 (1 - 1/2|x|)^2 & \mbox{si } 1 \lt |x| \le 2 \\ 0 & \mbox{si } |x| \gt 2 \end{array} \right. $$


In 3D the shape function is just the product of the 1D functions \(S(\mathbf{\mathbf{r}_m}) = S(x-x_m)S(y-y_m)S(z-z_m)\)

Spatial discretization

In PHARE, equations are written independentely of the concrete discretization so it can be changed if needed. The discretization that is concretely implemented is the one commonly known as the Yee Lattice.

Time integration

The time evolution of Hybrid equations in PHARE is an asbtract object, which allows to isolate any concrete implementation. The current implementation is a predictor-predictor-corrector scheme, also used in Kunz et al. 2014. Knowing the magnetic field \(\mathbf{B}^n\) and electric field \(\mathbf{E}^n\) at time \(t=n\), we compute a first prediction of the magnetic field at \(t=n+1\), called \(\mathbf{B}_{p1}^{n+1}\).

\begin{equation} \mathbf{B}^{n+1}_{p1} = \mathbf{B}^{n} - \Delta t \nabla \times \mathbf{E}^n \label{eq:faraday_disc} \end{equation}

Then, having now \(\mathbf{B}_{p1}^{n+1}\) and the ion moments at \(t=n\), we compute equation \ref{eq:ohmelectron4} to get a first prediction of the electric field \(\mathbf{E}_{p1}^{n+1}\):

\begin{equation} \mathbf{E}^{n+1}_{p1} = -\mathbf{u_e}^n\times\mathbf{B}_{p1}^{n+1} - \frac{1}{en^n}\nabla P_e^n + \eta \mathbf{j}^{n+1} - \nu\nabla^2\mathbf{j}^{n+1} \label{eq:ohm_disc} \end{equation}

Now, we compute the electromagnetic field at \(t = n+1/2\) as the midpoint of the particle push that will take place just after.

\begin{eqnarray} \mathbf{E}^{n+1/2}_{p1} &=& \frac{1}{2}(\mathbf{E}^n + \mathbf{E}^{n+1}_{p1}) \\ \mathbf{B}^{n+1/2}_{p1} &=& \frac{1}{2}(\mathbf{B}^n + \mathbf{B}^{n+1}_{p1}) \end{eqnarray}

At this point we have the field required to accelerate and push the particles from \(t=n\) to \(t=n+1\), as a first prediction of their position and velocity.

\begin{equation} \mathbf{r}^n \to \mathbf{r}^{n+1}_p\\ \end{equation} \begin{equation} \mathbf{v}^n \to \mathbf{v}^{n+1}_p\\ \end{equation}

From these new positions and velocities, the density \(n\) and the ion bulk velocity \(\mathbf{u_i}\) are computed on the mesh using the b-spline function at chosen order.

Now we have \(n^{n+1}\) and \(\mathbf{u_i}^{n+1}\), and also \(\mathbf{E}^{n+1/2}\) and \(\mathbf{B}^{n+1/2}\), we can start again from \(t=n\) to get a better predictions of things at \(t=n+1\). Let's go first with Faraday's law, that can now be written in a time centered fashion:

\begin{equation} \mathbf{B}^{n+1}_{p2} = \mathbf{B}^{n} - \Delta t \nabla \times \mathbf{E}^{n+1/2} \label{eq:faraday_p2} \end{equation}

Then again Ohm's law, but this time using \(\mathbf{B}_{p2}^{n+1}\) and moments at \(t=n+1\)

\begin{equation} \mathbf{E}^{n+1}_{p2} = -\mathbf{u_e}^{n+1}\times\mathbf{B}_{p2}^{n+1} - \frac{1}{en^{n+1}}\nabla P_e^{n+1} + \eta \mathbf{j}^{n+1} - \nu\nabla^2\mathbf{j}^{n+1} \label{eq:ohm_p2} \end{equation}

From which we can again take the time centered average :
\begin{eqnarray} \mathbf{E}^{n+1/2}_{p1} &=& \frac{1}{2}(\mathbf{E}^n + \mathbf{E}^{n+1}_{p2}) \\ \mathbf{B}^{n+1/2}_{p1} &=& \frac{1}{2}(\mathbf{B}^n + \mathbf{B}^{n+1}_{p2}) \end{eqnarray}

And push again the particles, to \(t=n+1\) one last time: \begin{equation} \mathbf{r}^n \to \mathbf{r}^{n+1}\\ \end{equation} \begin{equation} \mathbf{v}^n \to \mathbf{v}^{n+1}\\ \end{equation}
And update the total density and ion bulk velocity at \(t=n+1\).

At this point since we just got a better time centered electric field, we compute Faraday's law one last time to get the final version of the magnetic field

\begin{equation} \mathbf{B}^{n+1} = \mathbf{B}^{n} - \Delta t \nabla \times \mathbf{E}^{n+1/2} \label{eq:faraday_c} \end{equation}

And using the just obtained magnetic field, together with newly calculated moments, we compute Ohm's law one last time to get our final version of the electric field:

\begin{equation} \mathbf{E}^{n+1} = -\mathbf{u_e}^{n+1}\times\mathbf{B}^{n+1} - \frac{1}{en^{n+1}}\nabla P_e^{n+1} + \eta \mathbf{j}^{n+1} - \nu\nabla^2\mathbf{j}^{n+1} \label{eq:ohm_c} \end{equation}

Adaptive Mesh Refinement

Spatial refinement ratio and refined time stepping

Refinement ratio 2 and recursive time stepping

Field refinement

Linear field interpolation...

Particle splitting

Exact and approximative splitting

Coarsening and synchronization points

coarsening algorithm and which quantity is coarsenned

Refinement criteria

tagging for refinement...