Field initialization for relativistic species

As explained in PIC algorithms, if a net charge is present at the beginning of the simulation, the initial electromagnetic fields are computed. For static charge distributions, the solution of Poisson’s equation will be necessary to find the initial electrostatic field. If the initial charge has a non-zero initial speed, in general the electric and magnetic field should be computed solving the full set of Maxwell’s equations or equivalently the potentials equations. In some physical setups of interest, one or more relativistic species are injected in a plasma. In these cases, the computation of the initial electromagnetic fields can be reduced to the solution of a modified version of Poisson’s equation.

---

The relativistic Poisson’s equation

From the continuity equation ${\partial }_t\rho +\mathrm{\nabla }\cdot \mathbf{J}=0$ and Maxwell’s equations, it can be shown that the quantities $\mathrm{\nabla }\cdot \mathbf{B}$ and $\mathrm{\nabla }\cdot \mathbf{E}-\rho$ do not change over time:

$$\begin{array}{r}\begin{array}{rcl}{\partial }_t(\mathrm{\nabla }\cdot \mathbf{B})& =& 0,\\ {\partial }_t(\mathrm{\nabla }\cdot \mathbf{E}-\rho )& =& \mathrm{\nabla }\cdot {\partial }_t\mathbf{E}-{\partial }_t\rho =\mathrm{\nabla }\cdot (\mathrm{\nabla }\times \mathbf{B}-\mathbf{J})-{\partial }_t\rho =-(\mathrm{\nabla }\cdot \mathbf{J}+{\partial }_t\rho ).\end{array}\end{array}$$

Thus, if a simulation starts with $\rho \ne 0$, the electromagnetic fields must be properly initialized.

In the case of a static charge distribution, i.e. $\rho \ne 0$, $\mathbf{J}=0$, the initial electrostatic potential $\mathrm{\Phi }$ can be computed solving Poisson’s equation:

$${\mathrm{\nabla }}^2\mathrm{\Phi }=-\rho ,$$

and then be integrated to find the initial electric field: $\mathbf{E}=-\mathrm{\nabla }\mathrm{\Phi }$. The initial magnetic field $\mathbf{B}$ will be zero.

In general when the initial current $\mathbf{J}$ is not zero, the full set of fields equations should be solved to correctly initialize the electromagnetic fields.

However, if a species is already relativistic when the simulation starts, e.g. a relativistic electron bunch, its initial electromagnetic fields can be computed through a simplified procedure, described in [Vay2008], [Londrillo2014], [Massimo2016] and [Marocchino2018].

An important assumption of this calculation is that the species is highly relativistic, moving in the positive $x$ direction, with negligible momentum spread. Under this hypothesis, the transverse components of the species current density are neglected and the four-current quadrivector can be written as:

$$(\mathbf{J},\rho )=(\rho {\beta }_0,0,0,\rho ),$$

where ${\beta }_0$ is the initial mean velocity of the relativistic species. At least locally, the potentials $\mathbf{A}$, $\mathrm{\Phi }$ in the laboratory frame will be only function of $x-{\beta }_{0}t$, as they are propagating with the species at uniform relativistic velocity.

In the relativistic species rest frame $S^{\prime }$, the charge distribution is static and the electrostatic potential in that reference frame ${\mathrm{\Phi }}^{\prime }$ is related to the charge density in that reference frame ${\rho }^{\prime }$ through Poisson’s equation:

(68)$${\mathrm{\nabla }}^{\prime 2}{\mathrm{\Phi }}^{\prime }=-{\rho }^{\prime },$$

where the Laplacian operator is computed in the reference frame $S^{\prime }$:

$${\mathrm{\nabla }}^{\prime 2}={\partial }_{x^{\prime }}^2+{\partial }_{y^{\prime }}^2+{\partial }_{z^{\prime }}^2.$$

The vector potential in the species rest frame can be set to zero: ${\mathbf{A}}^{\prime }=0$. Through the above mentioned assumptions, it is possible to rewrite Eq. (68) only in terms of laboratory frame quantities.

Lorentz transformation of the four-vector $(\mathbf{J},\rho )$ yields ${\rho }^{\prime }=\rho /{\gamma }_0$, where ${\gamma }_0=1/\sqrt{1-{\beta }_0^2}$ is the average Lorentz factor of the relativistic species. Similarly, the potential ${\mathrm{\Phi }}^{\prime }$ can be rewritten in terms of the potential in the laboratory frame: ${\mathrm{\Phi }}^{\prime }=\mathrm{\Phi }/{\gamma }_0$. The Lorentz back-transformation of coordinates

$$x={\gamma }_0(x^{\prime }+{\beta }_0{t}^{\prime }),t={\gamma }_0(t^{\prime }+{\beta }_0{x}^{\prime }),y=y^{\prime },z=z^{\prime }$$

allows to transform the derivatives in Eq. (68) as

$${\partial }_{x^{\prime }}={\gamma }_0({\partial }_x+{\beta }_0{\partial }_t),{\partial }_{y^{\prime }}={\partial }_y,{\partial }_{z^{\prime }}={\partial }_z.$$

The partial derivative along the $x^{\prime }$ direction can be further simplified, through the hypothesis of temporary dependence of all quantities on $x-{\beta }_{0}t$, implying ${\partial }_t=-{\beta }_0{\partial }_x$:

$${\partial }_{x^{\prime }}=\frac{1}{{\gamma }_0}{\partial }_x.$$

Equation (68) can thus be rewritten as

(69)$$(\frac{1}{{\gamma }_0^2}{\partial }_x^2+{\mathrm{\nabla }}_{\perp }^2)\mathrm{\Phi }=-\rho ,$$

here informally referred to as the relativistic Poisson’s equation. In Smilei, as for Eq. (68), the solution of the relativistic Poisson’s equation is performed through the conjugate gradient method.

Once the potential $\mathrm{\Phi }$ is found, we can compute all the components of the electromagnetic field, using again the relations ${\partial }_t=-{\beta }_0{\partial }_x$, ${\mathrm{\Phi }}^{\prime }=-\mathrm{\Phi }/{\gamma }_0$ and the Lorentz back-transformation of the vector potential $\mathbf{A}$:

$$A_x={\gamma }_0(A_x^{\prime }+{\beta }_0{\mathrm{\Phi }}^{\prime })={\gamma }_0{\beta }_0{\mathrm{\Phi }}^{\prime }={\beta }_0\mathrm{\Phi },A_y=A_y^{\prime }=0,A_z=A_z^{\prime }=0.$$

From all these relations, the electromagnetic field can be computed as usual, through the definitions of potentials $\mathbf{E}=-\mathrm{\nabla }\mathrm{\Phi }-{\partial }_t\mathbf{A}$, $\mathbf{B}=-\mathrm{\nabla }\times \mathbf{A}$:

$$\begin{array}{r}\begin{array}{rcl}{E}_x& =& -{\partial }_x\mathrm{\Phi }-{\partial }_t{A}_x=-{\partial }_x\mathrm{\Phi }+{\beta }_0^2{\partial }_x\mathrm{\Phi }=-\frac{1}{{\gamma }_0^2}{\partial }_x\mathrm{\Phi },\\ E_y& =& -{\partial }_y\mathrm{\Phi }-{\partial }_t{A}_y=-{\partial }_y\mathrm{\Phi },\\ E_z& =& -{\partial }_z\mathrm{\Phi }-{\partial }_t{A}_z=-{\partial }_z\mathrm{\Phi },\\ \\ B_x& =& {\partial }_y{A}_z-{\partial }_z{A}_y=0,\\ B_y& =& {\partial }_z{A}_x-{\partial }_x{A}_z={\partial }_z{A}_x={\beta }_0{\partial }_z\mathrm{\Phi }=-{\beta }_0{E}_z,\\ B_z& =& {\partial }_x{A}_y-{\partial }_y{A}_x=-{\partial }_y{A}_x=-{\beta }_0{\partial }_y\mathrm{\Phi }={\beta }_0{E}_y,\end{array}\end{array}$$

or in more compact form: $\mathbf{E}=(-\frac{1}{{\gamma }_0^2}{\partial }_x\mathrm{\Phi },-{\partial }_y\mathrm{\Phi },-{\partial }_z\mathrm{\Phi })$, $\mathbf{B}={\beta }_0\hat{\mathbf{x}}\times \mathbf{E}$.

From the previous equations, it can be inferred that, in a 1D cartesian geometry, the fields computed through this procedure equal those obtained through the standard Poisson’s problem. This can also be inferred from the relativistic transformations of fields, which conserve the $x$ components of the electromagnetic fields for boosts in the $x$ direction.

---

Recommendations for relativistic species field initialization

In Smilei, each species can independently benefit from this field initialization procedure. Its field will be initialized when the species will start to move, in order not to interfere with the other species’ dynamics. The initialized fields will be superimposed to the electromagnetic fields already present in the simulation. To have physically meaningful results, we recommend to place a species which requires this method of field initialization far from other species, otherwise the latter could experience instantaneous unphysical forces by the relativistic species’ fields.

Remember that the transverse field of a moving charge with relativistic factor $\gamma$ is greater than the electrostatic transverse field of that charge, by a factor $\gamma$. This means that for highly relativistic particles, you will need to use a transversely large simulation window to let the field decrease enough to reduce border effects during its propagation. A complete absence of boundary effects in this case would be provided by perfectly absorbing boundary conditions, which are not implemented yet in the code. If the relativistic species propagates in a plasma, these border effects could be partially screened by the plasma.

A relativistic mean velocity in the $x$ or $-x$ direction and a negligible energy spread are assumed for each species in the hypotheses of this procedure. In particular it is assumed that all $p_x$ of the macro-particles of a species are of the same sign. The user must ensure these conditions when defining the species requiring field initialization in the namelist. In order to describe a non-monoenergetic species, it is necessary to divide it into into several sub-species for which the assumptions hold.

---

References

[Vay2008]

J.-L. Vay, Physics of Plasmas 15, 056701 (2008)

[Londrillo2014]

P. Londrillo, C. Gatti and M. Ferrario, Nucl. Instr. and Meth. A 740, 236-241 (2014)

[Massimo2016]

F. Massimo, A. Marocchino and A. R. Rossi, Nucl. Instr. and Meth. A 829, 378-382 (2016)

[Marocchino2018]

A. Marocchino, E. Chiadroni, M. Ferrario, F. Mira and A.R. Rossi, Nucl. Instr. and Meth. A (2018)