High-order implicit particle-in-cell method for plasma simulations at solid densities

A high-order implicit multidimensional particle-in-cell (PIC) method is developed for simulating plasmas at solid densities. The space-time arrangement is based on Yee and a leapfrog algorithm for electromagnetic fields and particle advancement. The field solver algorithm completely eliminates numerical instabilities found in explicit PIC methods with relaxed time step and grid resolution. Moreover, this algorithm eliminates the numerical cooling found in the standard implicit PIC methods by using a pseudo-electric-field method. The particle pusher algorithm combines the standard Boris particle pusher with the Newton-Krylov iteration method. This algorithm increases the precision accuracy by several orders of magnitude when compared with the standard Boris particle pusher and also significantly decreases the iteration time when compared with the pure Newton-Krylov method. The code is tested with several benchmarks, including Weibel instability, and relativistic laser plasma interactions at both low and solid densities.

Figure 1

Spatial staggering in the LAPINS code. The charge density is located at cell centers. The electric fields and the current density are staggered upward to cell faces. The magnetic fields are placed at cell edges. Please note that in the original Petrov and Davis's algorithm, the electric fields, the current densities, and charge densities are located at the cell center, while the magnetic fields are placed at cell edges.

Figure 2

Time staggering and integration order in the LAPINS code. Note that $t^{\prime n-1/2}=t^{\prime n}-\delta t^{\prime }/2$, $t^{\prime n+1/2}=t^{\prime n}+\delta t^{\prime }/2$, and $t^{\prime n+1}=t^{\prime n}+\delta t^{\prime }$.

Figure 3

Red (three point) and blue (five point) shadows define the particles' weight functions. These two functions are all accessible in the code.

Figure 4

Example: Second-order and fourth-order difference operation ${\mathbf{\nabla }}^{\prime }\times$ in 1D geometry. Due to the Yee mesh staggered layout of variables, the central difference is computed exactly where needed.

Figure 5

The $\mathbf{E}\times \mathbf{B}$ drift motion for four cases: (a) ${\omega }_c^{\prime }\delta t^{\prime }=5$, (b) ${\omega }_c^{\prime }\delta t^{\prime }=10$, (c) ${\omega }_c^{\prime }\delta t^{\prime }=50$, and (d) ${\omega }_c^{\prime }\delta t^{\prime }=100$. The following parameters are the same for the simulation cases: ${\mathbf{u}}^{\prime }(t^{\prime }=0)=(0.1,0,0)$, ${\mathbf{B}}^{\prime }=(250,0,0)$, and ${\mathbf{E}}^{\prime }=(0,1,0)$. Here the (upper) red lines shows results from Newton-Krylov method, (lower) black lines are results from Boris's method, and the dashed lines are the theoretical predictions.

Figure 6

The first row shows energy distribution of electromagnetic fields as functions of the time and positions. The second row shows the value of $F_n^{\prime }$ (see text for explanation) as functions of time and positions. The first, second, and third columns correspond to the results carried out by LAPINS code with $d^{\prime }=0.0$, $d^{\prime }=0.5$, and $d^{\prime }=1$, respectively.

Figure 7

(a) The residue charge, i.e., integrating $F_n^{\prime }$ over the whole space, as a function of time. (b) The energy conservation as a function of time. Here (thick) blue line, (dashed) red line, and black line refer to results carried out by LAPINS code with $d^{\prime }=0.0$, $d^{\prime }=0.5$, and $d^{\prime }=1$, respectively.

Figure 8

The energy (transfer) as a function of time. Data are shown for plasma electron kinetic energy (dashed red), ion kinetic energy (black), and electromagnetic energy (thick blue). Here the electromagnetic energy is defined as $(1/2)\sum ({\mathbf{E}}^{\prime 2}+{\mathbf{B}}^{\prime 2})\delta x^{\prime }\delta y^{\prime }\delta z^{\prime }$. The straight line shows the analytical growth rate.

Figure 9

The evolution of plasma electron density (first column), beam density (second column), and plasma ion density (third column) as a function of time. The densities are in units of $1\times {10}^{21}{\mathrm{cm}}^{-3}$, and time is in units of 3.33 fs. For all simulations, the size of the simulation box is fixed as $4\mu \text{m}\times 4\mu \text{m}$ with $400\times 400$ cells and the time step is fixed the same. See text for the detailed explanation.

Figure 10

The energy density distributions of electromagnetic fields at $t=15$. Here the electromagnetic energy density is defined as $(1/2)({\mathbf{E}}^{\prime 2}+{\mathbf{B}}^{\prime 2})$. The different rows, (a), (b), (c), and (d), correspond to different simulation cases with background plasma densities $n_p=1\times {10}^{21}{\mathrm{cm}}^{-3}$, $n_p=1\times {10}^{22}{\mathrm{cm}}^{-3}$, $n_p=1\times {10}^{23}{\mathrm{cm}}^{-3}$, and $n_p=1\times {10}^{24}{\mathrm{cm}}^{-3}$, respectively. The first column is run by LAPINE code, the second column is run by high-order LAPINE code, and the third column is run by recently developed LAPINS code. For all simulations, the size of the simulation box is fixed as $5\mu \text{m}\times 5\mu \text{m}$ with $100\times 100$ cells and the time step is fixed the same. See text for the detailed explanation.

Figure 11

The energy (transfer) as a function of time. Data are shown for laser energy entering into the simulation box (black), electromagnetic energy (dashed red), and particle kinetic energy (thick blue).

Figure 12

The distribution of electron kinetic energy, electron density, and static magnetic field. The first row shows values at $t=35$ fs and the second row shows at $t=50$ fs. This simulation is produced by LAPINS code.

Figure 13

The distribution of electron kinetic energy, electron density, and static magnetic field. The first row shows values at $t=35$ fs and the second row shows at $t=50$ fs. This simulation is produced by high-order LAPINE code approximated with the layered-density method.