29 Particle-in-Cell
29.1 Phase Space Sampling
Let \(\hat{f}(\mathbf{z})\) be a general probability density function with normalization: \[ \int_\Omega \hat{f}(\mathbf{z})dV = 1 \]
Now sample the phase space with \(N\) markers according to \(\hat{f}(\mathbf{z})\) then the marker distribution, \(g\) is: \[ g(\mathbf{z}) = N\hat{f}(\mathbf{z}) \] and the phase space volume element occupied by a marker is \[ dV_i = \frac{1}{g(\mathbf{z}_i)} \] where \(\mathbf{z}_i\) is the phase space coordinate of the marker and \(dV_i\) the volume element occupied by it.
Define the weight of a marker as the number of physical particles in the corresponding phase space volume occupied by the marker: \[ w_i = f(\mathbf{z}_i)dV_i = \frac{f(\mathbf{z}_i)}{g(\mathbf{z}_i)} \]
This marker distribution can represent the total \(f\) or the \(\delta f\). (???)
We can now take moments of the distribution function in the phase space (such as density), and the Monte Carlo approximation to the integration is as follows: \[ M(A) = \int_\Omega A(\mathbf{z})f(\mathbf{z})dV = \sum_{i=1}^{N}A(\mathbf{z}_i)w_i + \mathcal{\frac{1}{\sqrt{N}}} \] where \(N\) is the number of markers in the phase space region \(\Omega\). This Monte Carlo approximation makes the particle simulation feasible for a number of problems, even in 3D, at the cost of introducing an error on the order of \(N^{-1/2}\).
29.2 Artificial Scalings
Common scales to characterize the plasma system:
Global scales of the system \(d_g\), characterized by the magnetopause standoff distance in the magnetosphere, some fraction of the solar radius in the corona, or system length for Tokamak.
Ion kinetic scales characterized by the ion inertial length \(d_i\) \[ d_i = \sqrt{\frac{m_i}{n_i q_i^2 \mu_0}} = \frac{m_i}{q_i}\sqrt{\frac{1}{\rho_i \mu_0}} \] and the ion gyroradius \[ r_i = \frac{m_i v_{\mathrm{th},i}}{q_i B} = \frac{m_i}{q_i}\frac{\sqrt{p_i/\rho_i}}{B} \]
Electron kinetic scales characterized by the electron skin depth \(d_e\) \[ d_e = \sqrt{\frac{m_e}{n_e q_e^2 \mu_0}} = \frac{m_e}{q_e}\sqrt{\frac{1}{\rho_e \mu_0}} \]
The smallest plasma scale characterized by the Debye length \(\lambda_D\) \[ \lambda_D = \sqrt{\frac{\epsilon_0 v_{\mathrm{th},e}^2}{q_e^2 n_e}} = \frac{m_e}{q_e}\frac{\epsilon_0 p_e}{\rho_e} \] where \(v_{\mathrm{th},e}=\sqrt{p_e / \rho_e}\) is the electron thermal velocity.
Typically, \(d_i\) are orders of magnitude smaller than \(d_g\), and \(d_e\) is even smaller. If we are interested in the interplay between the global plasma system and the kinetic processes, it is a natural question to ask how the behavior of the system depends on the ratio \(d_g∕d_i\). For example, in Ganymede’s magnetosphere \(d_g∕d_i \sim 10\) and the kinetic effects have non-negligible contributions to the global system. On the other hand, if \(d_g∕d_i\) is a very large number, then the kinetic effects will be mostly limited in a small region.1
For quasi-neutral plasmas, \(n_i = n_e \rightarrow \rho_i/\rho_e = m_i/m_e\), \(q_i = q_e\), \(p_i = p_e\). The ion scales are \(\sqrt{m_i/m_e}\) times larger than the corresponding electron scales. For a proton-electron plasma, \(\sqrt{d_i/d_e} \approx 43\). A standard trick to reduce the computational requirement is to artificially reduce the mass ratio to a smaller value (\(1836 \ge d_i/d_e \ge 25\)). In terms of reconnection, we have found only a relatively weak dependence of the process on the mass ratio. In practice almost all numerical studies, especially in 3-D, use a reduced ion-electron mass ratio.
Tóth et al. (2017) proposed another scaling methods to increase the ion and electron mass-to-charge ratios by a kinetic scaling factor \(f\) while keeping the MHD quantities, the mass densities \(\rho_e\) and \(\rho_i\), the pressures \(p_e\) and \(p_i\), the bulk velocities \(u_e\) and \(u_i\), the magnetic field B, and the various constants \(\epsilon_0, \mu_0\), and c unchanged. A nice property of this scaling method is that the characteristic speeds (bulk velocity, thermal velocity, Alfvén speed) are not modified; in other words, this kinetic scaling has no effect on ideal/resistive MHD. A direct implication of this scaling is the solution on the kinetic scales is similar for different values of f, but the spatial and temporal scales are proportional to f.
29.3 Numerical Heating and Thermalization
Both numerical heating and numerical thermalization are artifacts that can arise in PIC simulations, but they represent different phenomena. Note that “heating” is a general word: in fact, for implicit PIC schemes numerical cooling happens more often.
- Numerical Heating:
- Mechanism: This refers to the spurious increase in the overall kinetic energy of the system due to discretization errors. It often arises from the interaction of particles with the electric and magnetic fields that are defined on a grid. The finite grid size and time step can lead to inaccuracies in the force calculation, causing particles to gain energy artificially.
- Effect: Leads to an unphysical increase in the plasma temperature over time. This can affect the stability of the simulation and lead to inaccurate results.
- Example: Imagine a particle accelerating in an electric field. In a real plasma, the particle would gain energy smoothly. However, in a PIC simulation, the particle might “jump” over grid points, leading to a slightly higher energy gain than expected.
- Numerical Thermalization:
- Mechanism: This is a more subtle effect related to the artificial relaxation of the velocity distribution function (VDF) towards a Maxwellian distribution. It is caused by the “numerical collisions” between macroparticles due to their finite size and the grid. Even if the total energy is conserved, the way this energy is distributed among the particles can be artificially altered.
- Effect: Drives the system towards an artificial thermal equilibrium, even if the real plasma would not reach such a state. This can mask kinetic effects and lead to incorrect predictions of transport properties.
- Example: Consider a beam of electrons injected into a plasma. In reality, this beam might persist for a long time. However, in a PIC simulation, numerical collisions can cause the beam to spread out and lose its distinct identity much faster than expected.
Key Differences:
- Energy Conservation: Numerical heating violates energy conservation, while numerical thermalization can occur even when energy is conserved.
- Timescale: Numerical heating often occurs on faster timescales compared to numerical thermalization.
- Distribution Function: Numerical heating primarily affects the average energy, while numerical thermalization modifies the shape of the distribution function.
Mitigation Strategies:
- Increasing the number of macroparticles: This reduces the graininess of the simulation and the impact of numerical collisions. Refining the grid: A smaller grid size improves the accuracy of field calculations and reduces discretization errors.
- Employing higher-order shape functions: These provide a more accurate representation of the charge density and reduce numerical heating.
- Using advanced PIC algorithms: Some algorithms are designed to conserve energy or momentum more accurately, reducing numerical artifacts.
29.4 Computational Cost
The key feature that makes PIC more feasible than the brute-force nbody method is the reduction of computational cost. Given \(N\) macro-particles in the system, a brute force method has a complexity of \(\mathcal{O}(N^2)\), and a PIC method has \(\mathcal{O}(N\log N)\). A rough idea how this is possible can be obtained by ignoring the pair interactions and consider collective effects especially for those particles at a large distance.
A commonly used trick is to introduce a reduced ion-to-electron mass ratio. In a 3D simulation, the cost scales with \((m_i/m_e)^3\). However, note that in low mass ratio runs (e.g. \(m_i/m_e=25\)), the thermal speed of the electrons may be comparable to the Alfvén speed which causes higher wave decay through Landau damping on the electrons.
29.5 δf Method
Consider the growth rate of a single unstable wave in a 1D-1V Vlasov-Poisson system of plasma. The energy conservation leads to \[ \frac{\mathrm{d}}{dt}\Big(\frac{\epsilon_0}{2}E(x)^2\Big) + \int q v E(x)f(x,v)dv = 0 \]
For a single mode \(E(x,t) = A(t)\sin(kx-\omega t)\), the growth rate is \[ \gamma = \frac{1}{A}\frac{\mathrm{d}A}{\mathrm{d}t} = -\int qvEfdv = -\sum_j q_j w_j E(x_j) \] where \(w_j\) is the particle weight. All kinetic information lies in the marker particle position distribution.
The idea of \(\delta-f\) method is to decompose the total distribution function into two parts, \[ f = f_0 + \delta f,\quad |\delta f|\ll f_0 \] and let the particle weight be proportional to the perturbation part, \[ w\propto \delta f,\quad w_j = \frac{\delta f(x_j, v_j)}{f(x_j, v_j)} \]
If the contribution from \(f_0\) is known, for example \[ \int qvE f_0 dv = \int qvf_0(v) A \sin(kx-\omega t)dv = 0 \] then the growth rate can be written as \[ \gamma = -\int qvE\delta f dv = -\int qvE\frac{\delta f}{f}f dv = \sum_j q_j v_j w_j E(x_j) \]
Since now the particle weight is proportional to \(\delta f/f\), the system noise would be significantly reduced by \((\delta f/f)^2\); in other words, for the same noise level the number of macro-particles required would be reduced by \((\delta f/f)^2\), which is typically \(\mathcal{O}(10^6)\). Of course this requires the prior assumption \(|\delta f| \ll f\) to be valid.
Reconnection can lead to wave generation, which is a process of energy and momentum conversion from the kinetic region to the global scale. Foreshock waves are also triggered by the kinetic processes and have a global impact. Therefore, this statement requires extra consideration.↩︎