{'fraction_non_alphanumeric': 0.050198688992290186, 'fraction_numerical': 0.01389124749516014, 'mean_word_length': 4.37888198757764, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Alternating-Order Interpolation in a Charge-Conserving Scheme for Particle-In-Cell Simulations

We discuss the interpolation of the electric and magnetic fields within a charge-conserving Particle-In-Cell scheme. The choice of the interpolation procedure for the fields acting on a particle can be constrained by analyzing conservation of the energy and the particle generalized momentum. The better conservative properties are achieved, if the alternating-order form-factor is used for interpolation, which combines the lower-order and higher-order interpolation from integer and semi-integer points of a staggered grid. This approach allows us to significantly reduce noise in the charge conserving scheme and improves both the results quality and the computational efficiency.


Introduction: charge-conserving PIC schemes.
Here we discuss the conservative properties of the Particle-In-Cell (PIC) numerical schemes. Recall that usually the conservative schemes are employed to solve the system of conservation laws. For example, the continuity equation, ∂ρm ∂t + ∇ · (ρ m u) = 0, for the fluid mass density, ρ m , may be advanced through the time step, ∆t, using the conservative scheme: V i (ρ m ) n+1 i = V i (ρ m ) n i − ∆t j σ ij , with the computational domain split into a set of control volumes ("cells"), the mass density at a given time instant, t n = n∆t, averaged over the volume of the cell, i, and the mass flux through the ij face averaged over the time step: where u is the fluid velocity and the face area vector, S ij is directed from cell i, to cell j. Since σ ij = −σ ji , the total mass is conserved: i (ρ m ) n+1 This idea may be employed to assure the charge conservation within the PIC method, although the governing equations differ from conservation laws: where the index p enumerates particles (electrons, ions), w p is a momentum-tomass ratio, and other notations are usual. The particle charge density, ρ, obeys not only the continuity equation, ∂ρ ∂t + ∇ · J = 0, but also the Poisson equation. Therefore, the charge conservation property is formulated as the relationship between the charge density and the electric field, E, and/or the electric current, J.
In the present paper we discuss only the CHarge-Conserving PIC (ChCPIC) schemes. The need to introduce this family of numerical schemes is that to simulate fully relativistic motions of plasmas in strong rapidly varying electromagnetic fields (such as the laser field, for example), one may want to directly advance the electric field using the first Eqs. (3). The "charge conservation" property in this case allows one to fulfill the Poisson equation for this advanced field without solving the equation directly. The way to assure that the charge conservation was developed in [1, 2,3,4,5] (see also [6,7,8,9] and references therein). First, a staggered grid should be used as in [10] to ensure the finite-difference approximations for [∇ × B], [∇ × E] to be divergence-free, as we discuss briefly in Sec. 2. Second, the currents through the cell faces should be calculated in such way that their divergence balances the charge leakage from the cell. In Sec. 3, this is done using the virtual path integration to solve the time integral in Eq.(1). However, once the way to compute the particle current is modified, we should also modify accordingly the scheme for interpolating electric and magnetic fields acting on this particle, which is the goal of this paper.
Note an important distinction of the ChCPIC schemes from the general cloudin-cell scheme, which is explicitly pronounced once we follow the conservative schemes ideology. On one hand, a "cloud" within the framework of the cloudin-cell scheme can be thought of as some charge density distribution, ρ p (x, x p ), centered about the particle coordinate vector, x p . Following this ideology, the point-wise values of the shape-function (form-factor), ρ p (x, x p ), at the points x, where the electromagnetic fields are localized, should be used then as the interpolation coefficients to average the electromagnetic fields acting on the particle, and to calculate the Lorentz force. Contrarily, the conservative schemes approach assumes that the cloud charge density should be integrated over some control volume or its faces, to represent the contribution to the plasma density and the plasma electric current (see 1), thus making inapplicable the form-factor pointwise values. To constrain the field interpolation procedure, in Sec. 4 we discuss the accuracy of conservation, for the energy integral and for the particle generalized momentum. Their conservation in the ChCPIC scheme can be achieved, if the alternating-order form-factor is employed in the interpolation procedure. This alternating-order form-factor combines the interpolation of different order for different physical variables to interpolate. Equivalently, it integrates the "cloud" charge density over faces and edges of the control volume to interpolate facecentered and edge-centered electromagnetic field values.

Grid geometry and notations.
We use a 3D Cartesian grid in the domain 0 ≤ x ≤ L x , 0 ≤ y ≤ L y , 0 ≤ z ≤ L z , split for N x * N y * N z cells. The coordinates of the cell corners are (i ∆x, j ∆y, k ∆z), where i, j, k are integers and ∆x = L x /N x , ∆y = L y /N y , ∆z = L z /N z are the cell sizes. Then we introduce the normalized coordinates, x = x/∆x,ỹ = y/∆y,z = z/∆z, and time,t = t/∆t, and use them below with no tilde. Magnetic field, electric current, and particle momenta are defined at semi-integer time instants, t = n + 1/2, the electric field and the particle coordinates -at integer time instants, t = n. In the normalized coordinates, Eqs.(5) read: where V = ∆x∆y∆z and c x = c∆t/∆x, .... The grid functions are defined: the cell-centered charge density, ρ i+1/2,j+1/2,k+1/2 ; the electric field, E i+1/2,j+1/2,k , and the current density components, J i+1/2,j+1/2,k , defined at the centers of the faces, normal to the axis, x, y, z; and the magnetic field components, i,j,k+1/2 , defined at the midpoints of the edges directed along the axis, x, y, z. The subscript indexes denote coordinates of the point at which the grid function is defined.
The PIC scheme as taken from [11] with the suggested modifications is described in the Appendix. The algorithm involves interpolation for the fields acting on a particle: where α, β are the weights, their sums for a given particle should be equal to one: Herewith, we provide only the expressions for x component of vectors, whenever possible, denoting the generalization for the other components by '...'. The contribution to the current density from a charged particle can be expressed in terms of the particle position, x n p , and its velocity, u n+1/2 p , and/or x n+1 p : , , The advantage of the staggered grid is that the magnetic field divergence does not change and equals zero as long as it is initially equal to zero. Analogously, [∇ × B] term does not affect the electric field divergence. The Poisson equation, i,j+1/2,k+1/2 )/∆x + ..., is satisfied, if: 3. Charge density and charge conservation law.

Form-factors.
To discretize the charge and current densities, one needs to specify the numerical representation for δ functions in Eqs. (7). We do this using the family of form-factor functions, f (l) (x, x p ), where the form-factor of a zero order is a cap-function: f (0) (x, x p ) = 1, and the higher-order form-factors are recursively defined: All form-factors: (1) are symmetric functions of x − x p ; (2) turn to zero at |x − x p | > (l + 1)/2; and (3): We are interested both in point values of the form-factor function, and in its integrals over the grid size. So, for a chosen form-factor, f (x, x p ) = f (l) (x, x p ), we introduce: By definition, ∆F i+1/2 (x p ) = f (l+1) (i + 1/2, x p ). The applicability of the formfactors for constructing the interpolation weights, which should satisfy Eq.(8), is ensured by the equation: From Eq. (12) we can obtain yet another identity to be used below: To prove Eq.(13) one can note that for x p = x ′ p both sides of the equation turn to zero and that on taking the derivative of this equation over x ′ p we obtain the earlier proven Eq.(12).

Conservative scheme for electric charge
A particle can be thought of as a cloud with the charge density, Following the conservative scheme idea, we define the contribution from this particle to the charge density grid function not as a point value of ρ p in the cell centers, but via Eqs.(1,11): Note that the integrated over the cell size form-factor value, ∆F i+1/2 (x n p ), ... is at the same time the point-wise value of the form-factor function of by unity higher order, . Assuming that within the time interval, (n, n + 1), the particle moves from the point x n p to the point x n+1 p along an arbitrary virtual path x v (t), we define the particle currents, ξ (x,y,z) , in Eq.(9) following Eqs.(1,11,14)): By definition, Eqs.(14-17) satisfy the charge conservation law as in Eq.(10). This can be verified by observing that the linear combination of Eqs.(15-17) as presented in Eq.(10), reduces to d Discuss a choice of the virtual path. For the straight path, (15)(16)(17) are piecewise polynomials of the order of 3l + 2. The integration was performed in [2], [5] only for the lowest order form-factor, l = 0. With a finite time step, especially for higher-order form-factors, the integration becomes sophisticated, because the form-factor is a piecewise smooth function. All ranges of smoothness should be accounted for separately while evaluating the integral analytically.
To overcome this difficulty, we observe that the integrals in Eqs. (15)(16)(17) can be solved, if the virtual path is composed of the edges of the rectangular box, such that the points are the opposite corners of this box and its edges are parallel to the coordinate axis. Upon calculating the integral in Eq.(15) as a sixth of a sum of the integrals along six possible virtual paths, we find: Here ∆F The scheme as in Eq.(18) is not new and was obtained from different considerations and in a different form in [6] (see also discussion in [8,9]).
Note a useful property of the interpolation coefficients. Summing up the contributions to the electric current at different faces and using Eqs. (12,13) we obtain the following exact relationship: and the analogous relationships for y and z components. These relationships for the interpolation coefficients for the electric current are similar to those for the electric and magnetic fields (see 8).

Interpolation procedure for the electric and magnetic fields.
To interpolate the charge density we had to integrate the form-factor over the cell volume. Here we show that to interpolate the electric and magnetic fields the form-factor should be integrated over the cell faces and edges correspondingly.

Energy integral
Consider the energy, E p , of the system consisting of the electric field and a single charged particle. Both the particle energy and the electric current contributing to the electric field energy change are additive by particles, hence, so is the error in the energy conservation and thus, we can calculate it for a single particle. For simplicity, we assume zero magnetic field because this field does not affect the particle energy and any change in the magnetic field energy is balanced with that for the electric field. At time instant, t = n, the energy, E n p , equals: the last term is to advance the particle energy through a half time step. The change in the particle energy with the use of Eq.(33), can be approximated as: with the error, O((∆t) 3 ). Neglecting this error, we derive the change in the total energy, using Eqs. (9,35): where: x face (x n+1 p ) , ..., "x face" stands for i, j + 1/2, k + 1/2. If the interpolation weights for the electric field match those for the current in the way as follows: then ∆ 1 E = 0 and the energy defect is small: The problem is that, in contrast with an ordinary PIC scheme, within the ChCPIC scheme one can hardly satisfy (19) exactly. Eq.(19) can be obtained as the trapezoidal estimate for the integrals in Eqs. (15)(16)(17), if the interpolation weights, α, are chosen as follows: The accuracy of Eq.(19) and the energy defect while using Eqs. (15)(16)(17)(18)20) are controlled by the choice of the form-factor order. Using the estimate: where we substituted and noting that for l = 0, 1, 2 correspondingly the form-factor, f (x − x p ), its first or second derivative are bounded, we can evaluate the accuracy of the energy conservation using the above estimate: is not fulfilled at all, the energy does not conserve in the ChCPIC scheme: ∆ 1 E ∼ O(∆t). We conclude that both the use of higher-order form-factor (l ≥ 1) and the interpolation following Eq.(20) are desirable.

Generalized momentum conservation
The governing equations (2-6) conserve the projection of the generalized particle momentum, m p w p + q p A/c, on the given direction g, if the electromagnetic field is constant along this direction, i.e. (g · ∇)A = 0 (see [12]). Indeed, and g · d dt (m p w p + q p A/c) = 0 as long as (g · ∇)A = 0. To verify the generalized momentum conservation within the ChCPIC scheme, the latter should be formulated in terms of the vector potential. The grid functions A i+1/2,j+1/2,k are introduced at the same points as E i+1/2,j+1/2,k and the time derivative of the vector potential grid function may be expressed in terms of the electric field: in accordance with the first of Eqs. (6). The electric field and the vector potential at x = x p are also linked via this equation, hence, the vector potential in a particle location, A(t, x p ) should be interpolated with the same weights, α, as we use for interpolating the electric field: On the other hand, the magnetic field acting on the particle can be interpolated via the grid function of the magnetic field, In turn, the latter grid function can be expressed in terms of that for the vector potential, as the discretization of the second of Eqs.(6), on the staggered grid: i+1/2,j,k−1/2 , ..., so that: Now, we introduce the generalized momentum, Evaluating the difference P n+1/2 − P n−1/2 , we see that the transformation as in Eq.(21), which allows the generalized momentum conservation, is possible, if the differential equation, B = ∇ × A, is exactly fulfilled with the interpolated values of the magnetic field, and vector potential, i.e.
where ∇ p = ( 1 ∆x ∂ ∂xp , ..). Applying the operator, ∇ p × to Eq.(23) and comparing the result with Eq.(24) we find that Eq.(25) is fulfilled if the following set of equations holds: etc. These equations as well as Eq.(20) dictate the following choice for the magnetic field weights:

Electrostatic limit and a link to the finite element approach
To discuss the relationship with the important concepts of the momentum conservation and self-force (see detail in [11]) we now consider the electrostatic limit in which the fields acting on a given macroparticle can be assumed steady-state. Within this approximation, the electric field may be expressed in terms of the scalar electric potential, Φ: The grid function for the scalar potential is cell-centered, Φ i+1/2,j+1/2,k+1/2 . The relationships for the face-centered electric field components are as follows: The Poisson equation should be fulfilled in the following form: With this observations, the contribution from the electric field to the energy integral may be transformed using Eqs.(28,29) as follows: The key point is that in the last sum the charge density can be expressed in terms of the contributions from the given particles, then resulting in the field energy transformation to the sum over all charged particles: where the static potential at the particle location is introduced as follows: If the energy is conserved, then the energy integral may play the role of the Hamiltonian function, so that the electric force acting on a particle can be expressed in terms of partial derivatives of the Hamiltonian function over particle coordinates. As the result, the electric field acting on the particle should obey Eq. (27): But this is achieved with the choice of Eq.(20) for the weight coefficients, to interpolate the electric field acting on a particle! Indeed, Here, the conclusion of Langdon is reproduced, that the "energy-conserving" plasma simulations models require the contribution from the potential electric field to the force acting on a particle to be a gradient of the grid potential (see [13]).
From here yet another way to arrive at the suggested scheme for interpolating the electric and magnetic force may be recalled, specifically, the finite elements. The nodal amplitudes for fields and the form-factors are equivalent to local polynomial representation for the field using finite elements. The introduction of incomplete polynomial face (div-conforming) and edge (curl-conforming) elements is generally attributed to Nédélec [14]. For the "energy-conserving" scheme presented by Lewis in [15] (see also [16]) the use of such finite elements allows maintaining some of the Maxwell equations as the exact relationships between the finite elements for the fields. In the application to arbitrary curvilinear grids made in [5] this approach resulted in the same form-factors as in Eqs.(20,26), the form-factor order being l = 1. Nevertheless, the new derivation of Eqs.(20,26) for newer ChCPIC deserves the attention, moreover that the formulation of the conservation laws for such schemes is different from the formulation of [15]. We also emphasized that for the form-factor order of l = 1 there is no point in applying the interpolation scheme of Eqs.(20,26), as long as the energy (nor the generalized momentum) is not conserved anyway with this form-factor.

1D scheme and self-force
To illustrate the entire scheme, consider one-dimensional (1D) motion of a single particle in the given steady-state fields depending on x-coordinate only.
Assume zero transverse components of the electric field as well as a zero longitudinal component of the magnetic field. Eqs.(20,26) in 1D read: where the subscript ⊥ denotes the tranverse components of the vector potential and magneitc field and n x is the vector along x-coordinate. In the limit of an infinitely small timestep and in negecting for a while the particle charge effect on the fields, we find that three integrals fully describe the mactroparticle motion: where, again, the continuous potentials at the particle position are expressed in terms of cell-centered grid function values. Within this approximation, however, the motion description is perfectly accurate, as long as the conservation laws are fulfilled which fully control the particle motion. However, as noticed in [11], [13], the energy conservaion property of the numerical scheme results in the side effect in a form of a self-force. The conserved energy includes the contribution from the particle's own field energy, which is weakly dependent on the particle location relative to the grid nodes. As a result the weak force appears tending to pull the particle from cell centers towards faces. Now one can calculate the self-force effect on the 1D macroparticle motion. The current through x-faces, − ∆X ∆t qp V F i (x n+1 p ) − F i (x n p ) , can be integrated over time, giving the electric field: E 2 . Now we apply the weight coefficients as in Eq.(20) to obtain the self-force acting on the macroparticle in the form of the gradient of the particle's own field energy: With this accounting, the energy conservation law provides the distorted result: E + q p Φ(x p ) + Φ (own) (x p ) = const. Upon calculating the transverse magnetic field and vector potential induced by the particle slow motion, one can find that the transverse particle momentum is distorted accordingly, so that the ratio, which is the transverse velocity, is not modified with the particle's own field. So, in 1D non-relativistic motion the self-force does not affect the transverse velocity.
To evaluate the self-force effect on the longitudinal velocity, note that the dimensional factor in the formula for q p Φ (own) (x p ) as related to temperature equals ∆ 2 x/(Nr 2 D ), where r D is the Debye length and N is the number of macroparticles per cell. Usually (but not in the numerical tests discussed below) this multiplier is small, which reduces the effect of the self-force comparing to the macroparticle interaction with the regular electromagnetic field and with the neighboring particles. The dimensionless factor was evaluated in [11] for the lowest-order form-factor and the discontinuous character of the self-force was noticed for this case, resulting in potential problems. However, for the choice of a higher order form-factor (which is called "form-factor of the order l = 1 1 2 " below), the field energy is not only a smooth function of x p , but its variance is very small: The variance from minimum to maximum is 1/64. We see that the self-force is rather small and purely potential and its effect was not revealed in the numerical tests discussed below.

Alternating-order form-factor
Eqs. (14,20,26) may be interpreted in terms of the alternating-order formfactor, ϕ(x, x p )ϕ(y, y p )ϕ(z, z p ), with the value of this function at the grid point, x, y, z, giving the interpolation weight for the grid function defined at this point. Although the form-factor is a continuous function, only semi-integer and integers x, y, z matter, to which the grid functions are assigned. Thus, we define we alternate the form-factor order. As long as the alternating-order form-factor combines the values of functions, f (l) and f (l+1) , we would characterize it by the semi-integer order, l a = l + 1/2 Finally, we provide the values of the alternating-order form-factor for l a = 1 1 2 : and for l a = 2 1 2 , which really ensures high accuracy and good energy conservation:

Numerical tests
The tests we present combine some features of the tests for ChCPIC as developed in [8,9]. The tests are highly demanding: the mesh size is as large as the Debye length, the time step is as large as a half of ∆x/c, the macroparticle density is as small as two per cell in a three-dimensional case, thus increasing the macroparticle charge. Increased effects of particle-particle and particle-grid interactions allow for producing a noticeable numerical error within a short test runs.
First, we study the noise in the plasma with comparatively high level of thermal fluctuations. We use a three-dimensional 64*64*64 rectangular grid. The average number of electron particles per grid cell is N = 2. The electron thermal speed is: v T e = T e /m e = 0.05c, where T e is the electron temperature. At the initial time instant 2 * 64 3 electron particles are randomly distributed over the computational domain, the averaged plasma frequency being equal to ω pe = 1.  Figure 2: Evolution of the total energy normalized per the initial total energy. Blue curves demonstrate perfect energy conservation as achieved with the alternating-order form-factor of the order l a = 1 1 2 . The pink curves show much worse result obtained with higher (but not alternating) order l = 2. The left panel is for the Maxwellian plasma simulated through two plasma wave periods. The right panel is for two-stream instability followed till ω pe t = 200 (8000 iterations). equal number of immovable ions are put to the same locations as the electrons, so that the electromagnetic field is zero initially. The mesh size is ∆x = ∆y = ∆z = v T e /ω pe . With the time step, ω pe ∆t = 0.025, the plasma evolution has been simulated through two plasma wave periods, nω pe ∆t ≈ 4π, n = 503. Fortran 90/95 code is compiled and run with the double precision.
In Fig.2 (left panel) we present the evolution of the total energy normalized per the initial total energy. The blue curve demonstrates an almost perfect energy conservation as achieved with the alternating-order form-factor of the order l a = 1 1 2 . The pink curve shows, for comparison, the result obtained with the "conventional" way to interpolate the fields in the particle location, as used, for example in [9] as well as in older works, including the TRISTAN code [3]. Within this approach the cell-centered form-factor provides the interpolation weight, which is applied to the electric field vector as averaged over the faces of the given cell, which results in the following expression for α; and the analogous principle is applied to construct β-coefficients, with Eq.(18) being used to interpolate the electric currents, With this approach higher order l = 2 is used to interpolate both the charge density and the electric and magnetic field.
We clearly see that the increase in the interpolation order for the fields, com-paring to the alternating-order form-factor results in the degrade of an accuracy, as long as the energy defect becomes two orders of magnitude higher. This numerical result entirely agrees with our theoretical claims that within ChCPIC scheme the energy conservation degrades unless the alternating-order form-factor is used. Note also that the worse quality of the results with the uniform order l = 2 form-factor is accompanied by the degrade in efficiency. The stencil for l = 2 form-factor is wider than that for l a = 1 1 2 , therefore, more search/algebra operations are involved in this case. Thus, the efficiency drops from ≈ 2.1 · 10 5 particles advanced per second of CPU time per processor, for l a = 1 1 2 case, down to ≈ 1.5 · 10 5 particles/s/processor, for l = 2 case.
In the second test taken from [9] the two-stream instability is studied. The initial distribution is modified as follows. First, the magnetic field is applied along the direction of x-axis, such that the cyclotron frequency of electrons in this field equals one: ω Be = q e B/(m e c) = 1 = ω pe . To set the unstable distribution function with two streams, for 50% of the electron particles the directed stream velocity, v = 2v T e = 0.1c, along x-direction is added to the thermal random velocity. For the other half of the electron particles the negative of the first stream velocity, −v = −2v T e = −0.1c, is added. The instability evolution is traced through the time interval, 0 ≤ ω pe t ≤ 200. In the right panel of Fig.2, the check of the total energy conservation, again, demonstrates the advantage of the alternating-order form-factor (herewith the contribution from the initial magnetic field is excluded from the total energy integral).
In Fig.3 we present two test results obtained with the increased particle charge, so high that N = 1 particle per cell simulates the plasma of unity plasma frequency ω pe = 1, i.e. the charge is two times higher compared to the previous tests. Left panel presents the Langmuir oscillations. The initial distribution is simular to that for studying the two-stream instability, but the number of particles (of two times higher charge) is by a factor of 0.5 less. The unidirectional stream velocity, v = 2v T e = 0.1c, is added to all electron particles. The left panel of Fig.3 shows the total energy of electrons (upper curve) and the electric field energy, as functions of ω pe t. The test parameters (such as the low number of particles, resulting in comparatively strong electron-ion correlations) are intentionally chosen to make the errors noticeable, so we can observe gradual attenuation of the oscillations, as the result of electron-ion interactions. However, the defect in the total energy is as small as ∼ 10 −5 , thanks to using the alternating order from-factor. This defect is negligible compared to the changes in the partial energy integrals. In the right panel, we show the trajectories of particles gyrating in the magnetic field. The magnetic field is increased ten times: ω Be = 10, and it is applied along   z-axis. Two particles are located initially in the center of the computational domain, with oppositely directed initial momenta, w = ±0.5c, the particle charges have opposite signs (electron and positron particles). The particle trajectories are shown in the right panel of Fig.3 for the time interval, 0 ≤ ω Be t ≤ 400 with the timestep, ω Be ∆t = 0.25. One can observe a small gradual decrease in the particle kinetic energy. It reveals itself in the curve "thickness", with its inner radius, which is the Larmor radius at the final particle energy, being somewhat less than the outer radius, which is the Larmor radius at the initial particle energy. Our comparison with the field interpolation scheme presented in Eq.(31) does not aim to criticize [3,9], moreover that Eq.(31) in TRISTAN is used in a combination with low order l = 1 form-factor for the charge density, so that reducing this order following Eq.(20). would be impractical and still would not allow us to reach the energy conservation. However, with the higher order form-factor, as we see, the development of the interpolation scheme should go along the direction different from that presented in Eq.(31).

Conclusion
The pulsed electric field in a focus of a high-field laser may cause a charge separation even in the target of a solid-state density. To simulate the laser-plasma interaction in that strong fields, the ChCPIC scheme could be the best choice (see [6]). However, these schemes are believed to be too noisy and the energy nonconservation, as we see, can potentially be a source for such noise.