A Semi-implicit Particle-in-cell Expanding Box Model Code for Fully Kinetic Simulations of the Expanding Solar Wind Plasma

Fully kinetic simulations of space plasmas conducted with the Particle-In-Cell (PIC) approach have traditionally been limited in their capability of comparison against spacecraft observations by the computational cost of resolving both the ion and electron scales and the related constraints on domain size and temporal duration (Birdsall & Langdon 2004). Recently, however, these limitations have become less and less of a hindrance owing to a combination of factors. On one hand, supercomputers, largely thanks to multithreading and multicore technologies (Mallon et al. 2012; Intel 2018; Nvidia 2018), are approaching exascale, 10¹⁸ floating point operations per second (Top500, 2018). On the other hand, innovative algorithms (semi-implicit and implicit temporal discretizations; Brackbill & Forslund 1982; Langdon et al. 1983; Chen et al. 2011; Markidis & Lapenta 2011; Chen & Chacon 2015; Gonzalez-Herrero et al. 2018; grid refinement and grid remapping methods, Vay et al. 2004; Delzanno et al. 2008; Chacón et al. 2011; Fujimoto 2011; Innocenti et al. 2013, 2015; fluid-kinetic coupling strategies, Daldorff et al. 2014; Makwana et al. 2017) are pushing the boundaries of traditional PIC simulations.

Here we introduce the first fully kinetic Expanding Box Model semi-implicit code, obtained by introducing the Expanding Box Model (EBM; Velli et al. 1992; Grappin & Velli 1996; Liewer et al. 2001; Hellinger & Trávníček 2008; Hellinger et al. 2015) into the fully kinetic semi-implicit PIC code iPic3D (Markidis et al. 2010).

Our target environment is the expanding solar wind. This work is to be considered as a first step toward the self-consistent simulations of kinetic solar wind plasma evolution. In particular, the goal is to assess how the solar wind expansion drives the nonlinear evolution of coupled electromagnetic field fluctuations and the particle distribution functions of ions and electrons.

Electron and ion velocity distribution functions (VDFs) in the solar wind are far from thermodynamic equilibrium and exhibit distinct kinetic features (Marsch 2006, 2012). While some of these features may be of coronal origin (Pinfield et al. 1999; Esser & Edgar 2000; Chiuderi & Drago 2004), it is an established fact that propagation in the interplanetary space modifies particle distribution functions (Maksimovic et al. 2005), possibly as a result of wave/particle interaction through kinetic instabilities (Hellinger et al. 2006).

Our simulations will help to pinpoint the role of solar wind propagation and expansion in generating kinetic features in electron and ion VDFs. This work is particularly timely now after the launch of the Parker Solar Probe (PSP; Fox et al. 2016), which will explore the heliosphere from 1 au down to 9.86 solar radii and will help to clarify the "coronal origin" versus "propagation" conundrum.

PIC simulations have been used to study a variety of solar wind processes, including the constraining effect on the electron temperature of thermal anisotropy-induced instabilities (Camporeale & Burgess 2008; Hellinger et al. 2014), the role that scattering by high-frequency electromagnetic instabilities (e.g., whistler instability, electron cyclotron wave) has in producing suprathermal electrons (Saito & Gary 2007), and solar wind turbulence and its role in producing the observed spectra (Camporeale & Burgess 2011; Gary et al. 2012; Gary 2015; Pucci et al. 2017).

Solar wind expansion is almost never included in PIC simulations, with the exception of Camporeale & Burgess (2010), where a hot plasma expands radially in a cylindrical geometry, and Sironi & Narayan (2015), who studied plasma in an accretion flow (i.e., a compressing plasma). Ignoring plasma expansion removes from the simulations fundamental processes already partly uncovered by EBM hybrid simulations (Liewer et al. 2001; Hellinger & Trávníček 2008; Hellinger et al. 2015), where electrons, however, are treated as a massless isothermal fluid. One should also consider that adding expansion effects in MHD simulations allows us to better reproduce observed solar wind features, such as the local anisotropy of turbulent structures (Verdini & Grappin 2015, 2016), with respect to analogous simulations where expansion is not included. We can expect expansion in fully kinetic simulations to play an equally crucial role.

One of the reasons why solar wind propagation and expansion are not routinely included in PIC simulations may be the rapid increase of the domain size that comes with these processes if they are not incorporated into the evolution equations through a variable change, as done in the EBM. If one wants to simulate solar wind propagation and expansion with a standard PIC, one has to release the plasma parcel in a huge domain that already takes into account the distance ${\boldsymbol{U}}^{\prime} {T}_{\mathrm{sim}}$ that the parcel will travel owing to solar wind propagating with velocity ${\boldsymbol{U}}^{\prime}$ in the simulated time T_sim. At the beginning of the simulation, when the plasma has not propagated and expanded significantly yet, most of the cells will be empty of particles. However, they will nevertheless be included in code operations such as the field solver, thus unnecessarily increasing the communication burden of the system. Also, one has to consider the impact of such empty cells on a parallel code: with basic parallelization techniques such as domain decomposition, the cores hosting the empty cells would wait idly at the synchronization points unavoidable in PIC codes, thus wasting computational resources.

We circumvent these problems by using, in our code, the EBM, where plasma motion away from the Sun and nonradial expansion are introduced in the system via transformation of coordinates: this allows us to reduce the domain size and number of cells in our simulations from the "blue" grid to the smaller, stretching "red" grids in Figure 1. There, we depict respectively the domain and gridding of a "standard PIC" (blue) and "EBM-PIC" (red) simulation. The Sun is on the left of the figure.

Zoom In Zoom Out Reset image size

Using EBM-PIC already allows a significantly reduced computational cost for expanding solar wind plasma simulations. To further optimize our simulation efforts, the algorithm that we use in our code is the semi-implicit Implicit Moment Method (IMM; Brackbill & Forslund 1982; Lapenta et al. 2006). With a semi-implicit, rather than an explicit, PIC method, we can adjust the spatial and temporal resolution of the simulation to the spatial and temporal scales of the process we want to investigate, rather than having to resolve the smallest, fastest scales in the system out of stability concerns. This will allow us to target (moderately) large-scale solar wind processes, of the kind usually out of reach of explicit PIC codes, while retaining a kinetic description of the electrons.

The EBM IMM equations derived in this paper are implemented in the IMM, parallel, C++, Object-Oriented code iPic3D (Markidis et al. 2010). The EBM-IMM code is addressed as EB-iPic3D.

The structure of the paper is as follows: In Section 2 we describe the EBM and derive a formulation of Maxwell's and Newton's equations in an EBM system. In Section 3, the equations are discretized semi-implicitly following the blueprint of the IMM, and the equations for the EB-iPic3D code are obtained. In Section 4, we test EB-iPic3D against expected plasma behavior for unmagnetized and magnetized plasmas. Conclusions are drawn in Section 5, where we also illustrate our plans for the future.

The EBM, shown schematically in Figure 1, has been designed to describe the dynamics of a small sector of plasma (the simulation box) that is in spherical expansion owing to the advection, with speed ${\boldsymbol{U}}^{\prime}$, by the supersonic and super-Alfvénic average solar wind outflow. A change of frame of reference from a Sun-centered frame to a comoving frame is done, which allows us to follow the evolution of such a plasma sector in a semi-Lagrangian fashion by accounting for expansion effects via an explicit dependence on the average heliocentric distance R(t). The latter is taken along the x'-direction and can be considered as the Lagrangian variable of the background solar wind, evolving as R(t) = R₀ + U₀t. Both R₀ and U₀ are free parameters of the system, and R₀/U₀ = τ_e defines the expansion timescale. In this way, the EBM can describe in a locally plane geometry the dynamics occurring on timescales shorter than τ_e, while the effects of the expansion are taken into account as slow time variations of average quantities via a parametric dependence on R(t).

By assuming that spatial scales (i.e., the dimensions of the simulation box) are small with respect to R, the background solar wind velocity can be approximated as

$$\begin{eqnarray}&&{{\boldsymbol{U}}}^{{\prime} }={U}_{0}\left(1,\displaystyle \frac{{y}^{{\prime} }}{R},\displaystyle \frac{{z}^{{\prime} }}{R}\right).\end{eqnarray} \tag{ 1 }$$

As can be seen from Equation (1), the plasma parcel undergoes a stretching in the transverse directions y' and z' at a rate U₀/R because of the advection of the underlying solar wind velocity. Such an effect can be removed by the following change of variables that define the comoving frame:

$$\begin{eqnarray}\begin{array}{rcl}x & = & {x}^{{\prime} }-R(t),y={y}^{{\prime} }{R}_{0}/R(t),\\ z & = & {z}^{{\prime} }{R}_{0}/R(t),t={t}^{{\prime} }\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\boldsymbol{\nabla }}}^{{\prime} }=\left(\displaystyle \frac{\partial }{\partial x},\displaystyle \frac{{R}_{0}}{R}\displaystyle \frac{\partial }{\partial y},\displaystyle \frac{{R}_{0}}{R}\displaystyle \frac{\partial }{\partial z}\right)\equiv {\boldsymbol{\nabla }},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial {t}^{{\prime} }}=\displaystyle \frac{\partial }{\partial t}-{\boldsymbol{U}}\cdot {\boldsymbol{\nabla }},\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{U}}$ is the advection velocity written in terms of the new variables:

$$\begin{eqnarray}&&{\boldsymbol{U}}={U}_{0}\left(1,\displaystyle \frac{y}{{R}_{0}},\displaystyle \frac{z}{{R}_{0}}\right).\end{eqnarray} \tag{ 5 }$$

Similarly, the particle velocity ${\boldsymbol{v}}$ in the comoving frame is given by

$$\begin{eqnarray}&&\begin{array}{l}{v}_{x}={v}_{x}^{{\prime} }-{U}_{0},\quad {v}_{y}={v}_{y}^{{\prime} }-{U}_{0}\,y/{R}_{0},\quad {v}_{z}={v}_{z}^{{\prime} }-{U}_{0}\,z/{R}_{0}.\end{array}\end{eqnarray} \tag{ 6 }$$

The variable change of Equations (2) and (6) allows us to write Maxwell's and Newton's equation in the comoving frame. In the EB system (where ${U}_{0}/c\ll 1$), Maxwell's equations become

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\nabla }}\cdot {\boldsymbol{E}} & = & 4\pi \rho ,\\ {\rm{\nabla }}\cdot {\boldsymbol{B}} & = & 0,\\ {\rm{\nabla }}\times {\boldsymbol{E}} & = & -\displaystyle \frac{1}{c}\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}-\displaystyle \frac{1}{c}\overline{\overline{{\boldsymbol{P}}}}\cdot {\boldsymbol{B}},\\ {\rm{\nabla }}\times {\boldsymbol{B}} & = & \displaystyle \frac{4\pi }{c}{\boldsymbol{J}}+\displaystyle \frac{1}{c}\displaystyle \frac{\partial {\boldsymbol{E}}}{\partial t}+\displaystyle \frac{1}{c}\overline{\overline{{\boldsymbol{P}}}}\cdot {\boldsymbol{E}},\end{array}\end{eqnarray} \tag{ 7 }$$

where P_xx = 2U₀/R, P_yy = P_zz = U₀/R, and P_ij = 0 if $i\ne j$.

The continuity equation becomes

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-{\rm{\nabla }}\cdot {\boldsymbol{J}}-\rho {\rm{\nabla }}\cdot {\boldsymbol{U}}.\end{eqnarray} \tag{ 8 }$$

Newton's equations for particle motion read, in the EB system,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d{{\boldsymbol{x}}}_{p}}{{dt}} & = & {\overline{\overline{{\boldsymbol{P}}}}}_{x}\cdot {{\boldsymbol{v}}}_{p},\\ \displaystyle \frac{d{{\boldsymbol{v}}}_{p}}{{dt}} & = & \displaystyle \frac{{q}_{p}}{{m}_{p}}\left({{\boldsymbol{E}}}_{p}+\displaystyle \frac{{{\boldsymbol{v}}}_{p}\times {{\boldsymbol{B}}}_{p}}{c}\right)-{\overline{\overline{{\boldsymbol{P}}}}}_{v}\cdot {{\boldsymbol{v}}}_{p},\end{array}\end{eqnarray} \tag{ 9 }$$

where ${P}_{x,{xx}}=1$, ${P}_{x,{yy}}={P}_{x,{zz}}={R}_{0}/R$, and ${P}_{x,{ij}}=0$ if $i\ne j$. Similarly, ${P}_{v,{xx}}=0$, ${P}_{v,{yy}}={P}_{v,{zz}}={U}_{0}/R$, and ${P}_{v,{ij}}=0$ if $i\ne j$.

In this section, we derive the equations for the IMM in the Expanding Box (EB) system. The IMM equations in the nonexpanding system are derived in Brackbill & Forslund (1982), Ricci et al. (2002), and Lapenta et al. (2006).

For the sake of simplicity, this section is divided into four subsections. In Section 3.1 we derive the equation for the calculation of the electric field at the advanced time, ${{\boldsymbol{E}}}^{n+\theta }$, as a function of the particle moments at the advanced time, ${\rho }^{n+\theta }$ and ${{\boldsymbol{J}}}^{n+\theta }$. As customary in semi-implicit discretization, θ is the time-decentering parameter, with values $0.5\leqslant \theta \leqslant 1$. The time-advanced value of a dummy quantity ${\boldsymbol{\Xi }}$ is

$$\begin{eqnarray}&&{{\boldsymbol{\Xi }}}^{n+\theta }=\theta {{\boldsymbol{\Xi }}}^{n+1}+\left(1-\theta \right){{\boldsymbol{\Xi }}}^{n}.\end{eqnarray} \tag{ 10 }$$

In Section 3.2 we illustrate the core operation of the IMM method, the Taylor expansion of the moments at the advanced time, focusing on the specificity of the operation in the EB system. Then, we illustrate the semi-implicit approximation of particle motion in Section 3.3. Finally, in Section 3.4 we derive the equation used for the calculation of ${{\boldsymbol{E}}}^{n+\theta }$.

3.1. Derivation of the Function for the Time Update of the Electric Field as a Function of the Time-advanced Moments

The third and fourth lines of Equation (7) are discretized semi-implicitly as

$$\begin{eqnarray}&&{{\boldsymbol{B}}}^{n+1}={{\boldsymbol{B}}}^{n}-c{\rm{\Delta }}t{\rm{\nabla }}\times {{\boldsymbol{E}}}^{n+\theta }-{\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot {{\boldsymbol{B}}}^{n+\theta }\end{eqnarray} \tag{ 11 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{E}}}^{n+1} & = & {{\boldsymbol{E}}}^{n}+c{\rm{\Delta }}t\left({\rm{\nabla }}\times {{\boldsymbol{B}}}^{n+\theta }-\displaystyle \frac{4\pi }{c}{{\boldsymbol{J}}}^{n+\theta }\right)\\ & & -\ {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot {{\boldsymbol{E}}}^{n+\theta },\end{array}\end{eqnarray} \tag{ 12 }$$

respectively. We can see that the expansion introduces an extra time-advanced term on the right-hand side of each equation, $-{\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot {{\boldsymbol{B}}}^{n+\theta }$ and $-{\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot {{\boldsymbol{E}}}^{n+\theta }$, respectively.

The ${\rm{\nabla }}\times {{\boldsymbol{B}}}^{n+\theta }$ term in Equation (12) can be obtained by taking the curl of Equation (11) and using the vector identity

$$\begin{eqnarray}&&{\rm{\nabla }}\times \left({\rm{\nabla }}\times {{\boldsymbol{E}}}^{n+\theta }\right)={\rm{\nabla }}\left({\rm{\nabla }}\cdot {{\boldsymbol{E}}}^{n+\theta }\right)-{{\rm{\nabla }}}^{2}{{\boldsymbol{E}}}^{n+\theta }\end{eqnarray} \tag{ 13 }$$

and the time-discretized Gauss equation, the first line of Equation (7), at the advanced time, which reads

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {{\boldsymbol{E}}}^{n+\theta }=4\pi {\rho }^{n+\theta }.\end{eqnarray} \tag{ 14 }$$

We get

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\nabla }}\times {{\boldsymbol{B}}}^{n+\theta } & = & \displaystyle \frac{{\rm{\nabla }}\times {{\boldsymbol{B}}}^{n}}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}+\displaystyle \frac{c\theta {\rm{\Delta }}t}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}\\ & & \times \ \left({{\rm{\nabla }}}^{2}{{\boldsymbol{E}}}^{n+\theta }-{\rm{\nabla }}4\pi {\rho }^{n+\theta }\right),\end{array}\end{eqnarray} \tag{ 15 }$$

which, inserted into Equation (12) and using Equation (10), gives

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{E}}}^{n+\theta }-{\left(\displaystyle \frac{c\theta {\rm{\Delta }}t}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}\right)}^{2}{{\rm{\nabla }}}^{2}{{\boldsymbol{E}}}^{n+\theta }\\ =\ \displaystyle \frac{{{\boldsymbol{E}}}^{n}}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}+\displaystyle \frac{c\theta {\rm{\Delta }}t}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}\\ \times \ \left(\displaystyle \frac{{\rm{\nabla }}\times {{\boldsymbol{B}}}^{n}}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}-\displaystyle \frac{4\pi }{c}{{\boldsymbol{J}}}^{n+\theta }\right)\\ -\ {\left(\displaystyle \frac{c\theta {\rm{\Delta }}t}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}\right)}^{2}4\pi {\rm{\nabla }}{\rho }^{n+\theta }.\end{array}\end{eqnarray} \tag{ 16 }$$

We can see that, to calculate the electric field at the advanced time, ${{\boldsymbol{E}}}^{n+\theta }$, we need particle moments at the advanced time, ${{\boldsymbol{J}}}^{n+\theta }$ and ${\rho }^{n+\theta }$. These quantities are defined as

$$\begin{eqnarray}&&{\rho }^{n+\theta }=\displaystyle \sum _{p}{q}_{p}S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n+\theta })\end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{J}}}^{n+\theta }=\displaystyle \sum _{p}{q}_{p}{\bar{{\boldsymbol{v}}}}_{p}S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n+\theta }),\end{eqnarray} \tag{ 18 }$$

where the index p marks a summation on all particles, q_p is the particle charge, ${{\boldsymbol{x}}}_{p}$ is the particle position, ${\boldsymbol{x}}$ is the grid point, $S({\boldsymbol{x}}-{{\boldsymbol{x}}}^{n+\theta })$ is the shape function of a particle normalized to the current cell volume at the advanced time n + θ, and ${\bar{{\boldsymbol{v}}}}_{p}$ is the intermediate particle velocity, i.e., the velocity at intermediate time, defined as

$$\begin{eqnarray}&&{\bar{{\boldsymbol{v}}}}_{p}=\displaystyle \frac{{{\boldsymbol{v}}}_{p}^{n+1}+{{\boldsymbol{v}}}_{p}^{n}}{2}.\end{eqnarray} \tag{ 19 }$$

Equations (17) and (18) can be rewritten as a sum of known quantities (field and particle variables at time n) plus a term dependent on ${{\boldsymbol{E}}}^{n+\theta }$, which can be moved to the left-hand side. This will be shown in Section 3.2.

3.2. Approximation of the Time-advanced Moments

The shape function centered around the particle position at the advanced time, $S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n+\theta })$, can be approximated around the shape function centered around the particle position at time tⁿ, $S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n})$. When applied to Equation (17), the approximation reads

$$\begin{eqnarray}\begin{array}{rcl}{\rho }^{n+\theta } & = & \displaystyle \sum _{p}{q}_{p}S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n})\\ & & +\ \theta {\rm{\Delta }}t{\left|\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \sum _{p}{q}_{p}S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n+\theta })\right)\right|}_{{\boldsymbol{x}}-{{\boldsymbol{x}}}^{n}}.\end{array}\end{eqnarray} \tag{ 20 }$$

The first term on the right-hand side is simply the definition of ρⁿ. The second term can be expanded as

$$\begin{eqnarray}\begin{array}{rcl}{\rho }^{n+\theta } & = & {\rho }^{n}-\theta {\rm{\Delta }}t{\rm{\nabla }}\cdot \left(\displaystyle \sum _{p}{q}_{p}{\bar{{\boldsymbol{v}}}}_{p}^{{\prime} }S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n+\theta })\right)\\ & & +\ \theta {\rm{\Delta }}t\displaystyle \frac{\partial {\boldsymbol{x}}}{\partial t}\cdot {\rm{\nabla }}\left(\displaystyle \sum _{p}{q}_{p}S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n+\theta })\right),\end{array}\end{eqnarray} \tag{ 21 }$$

taking into account that both particle positions (second term on the right-hand side in Equation (21)) and grid point positions (third term on the right-hand side) are dependent on time in the EBM and that ${\rm{\nabla }}S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n})\sim {\rm{\nabla }}S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n+\theta })$ in the approximation. Notice that ${\bar{{\boldsymbol{v}}}}_{p}^{{\prime} }$ is the particle intermediate velocity in the Sun-centered grid (refer to Section 2 for the notation in the expanded and nonexpanded frame), which can be mapped through

$$\begin{eqnarray}&&{\bar{{\boldsymbol{v}}}}_{p}^{{\prime} }={\bar{{\boldsymbol{v}}}}_{p}+{\boldsymbol{U}}\end{eqnarray} \tag{ 22 }$$

to the intermediate velocity ${\bar{{\boldsymbol{v}}}}_{p}$ in the comoving frame.

Hence, Equation (21) becomes

$$\begin{eqnarray}\begin{array}{rcl}{\rho }^{n+\theta } & = & {\rho }^{n}-\theta {\rm{\Delta }}t{\rm{\nabla }}\cdot {{\boldsymbol{J}}}^{n+\theta }-\theta {\rm{\Delta }}t{\rm{\nabla }}\cdot \left({\boldsymbol{U}}{\rho }^{n+\theta }\right)\\ & & +\ \theta {\rm{\Delta }}t{\boldsymbol{U}}\cdot {\rm{\nabla }}{\rho }^{n+\theta }\\ & = & {\rho }^{n}-\theta {\rm{\Delta }}t{\rm{\nabla }}\cdot {{\boldsymbol{J}}}^{n+\theta }-\theta {\rm{\Delta }}t{\rho }^{n+\theta }{\rm{\nabla }}\cdot {\boldsymbol{U}},\end{array}\end{eqnarray} \tag{ 23 }$$

where we recall that

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{x}}}{\partial t}={\boldsymbol{U}}\end{eqnarray} \tag{ 24 }$$

and the definition of ${{\boldsymbol{J}}}^{n+\theta }$ (Equation (18)). The reader may notice that Equation (24) corresponds to the semi-implicit discretization of the continuity equation in the EBM system, Equation (8).

Equation (18) can be approximated along similar lines to give

$$\begin{eqnarray}&&{{\boldsymbol{J}}}^{n+\theta }=\displaystyle \sum _{p}{q}_{p}{\bar{{\boldsymbol{v}}}}_{p}S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n})-\theta {\rm{\Delta }}t{\rm{\nabla }}\cdot \left(\displaystyle \sum _{p}{q}_{p}{\bar{{\boldsymbol{v}}}}_{p}{\bar{{\boldsymbol{v}}}}_{p}S({\boldsymbol{x}}-{{\boldsymbol{x}}}_{p}^{n})\right).\end{eqnarray} \tag{ 25 }$$

3.3. Semi-implicit Discretization of Particle Motion

The equations for particle motion in the EB system, Equations (9), can be discretized semi-implicitly as follows:

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{x}}}_{p}^{n+1} & = & {{\boldsymbol{x}}}_{p}^{n}+{\overline{\overline{{\boldsymbol{P}}}}}_{x}^{n+\theta }\cdot {\bar{{\boldsymbol{v}}}}_{p}^{k}{\rm{\Delta }}t,\\ {{\boldsymbol{v}}}_{p}^{n+1} & = & {{\boldsymbol{v}}}_{p}^{n}+\displaystyle \frac{{q}_{p}{\rm{\Delta }}t}{{m}_{p}}\left({{\boldsymbol{E}}}_{p}^{n+\theta }\left({\bar{{\boldsymbol{x}}}}^{k}\right)+\displaystyle \frac{{\bar{{\boldsymbol{v}}}}_{p}^{k}\times {{\boldsymbol{B}}}_{p}^{n}\left({\bar{{\boldsymbol{x}}}}^{k}\right)}{c}\right)\\ & & -\ {\overline{\overline{{\boldsymbol{P}}}}}_{v}^{n+\theta }\cdot {\bar{{\boldsymbol{v}}}}_{p}^{k}{\rm{\Delta }}t,\end{array}\end{eqnarray} \tag{ 26 }$$

using a Crank–Nicolson predictor corrector scheme, where ${\bar{{\boldsymbol{x}}}}_{p}^{k}$ and ${\bar{{\boldsymbol{v}}}}_{p}^{k}$ are, respectively, the intermediate particle position and velocity at the kth iteration of the mover. Notice that ${{\boldsymbol{E}}}_{p}^{n+\theta }\left({\bar{{\boldsymbol{x}}}}^{k}\right)$ and ${{\boldsymbol{B}}}_{p}^{n}\left({\bar{{\boldsymbol{x}}}}^{k}\right)$ are the electric and magnetic fields calculated at the intermediate particle position.

After some straightforward manipulation, we obtain

$$\begin{eqnarray}&&\begin{array}{l}{\bar{{\boldsymbol{v}}}}_{p}+{\left(\displaystyle \frac{{q}_{p}{\rm{\Delta }}t}{2{m}_{p}c}\right)}^{2}{\bar{{\boldsymbol{v}}}}_{p}{\left|{{\boldsymbol{B}}}_{p}^{n}\right|}^{2}\\ +\ \displaystyle \frac{{\rm{\Delta }}t}{2}{\overline{\overline{{\boldsymbol{P}}}}}_{v}^{n+\theta }\cdot \left({\bar{{\boldsymbol{v}}}}_{p}+\displaystyle \frac{{q}_{p}{\rm{\Delta }}t}{2{m}_{p}c}\left({\bar{{\boldsymbol{v}}}}_{p}\times {{\boldsymbol{B}}}_{p}^{n}+\displaystyle \frac{{q}_{p}{\rm{\Delta }}t}{2{m}_{p}c}\left({\bar{{\boldsymbol{v}}}}_{p}\cdot {{\boldsymbol{B}}}_{p}^{n}\right)\cdot {{\boldsymbol{B}}}_{p}^{n}\right)\right)\\ =\ {\hat{{\boldsymbol{v}}}}_{p}+\displaystyle \frac{{q}_{p}{\rm{\Delta }}t}{2{m}_{p}c}\left({\hat{{\boldsymbol{v}}}}_{p}\times {{\boldsymbol{B}}}_{p}^{n}+\displaystyle \frac{{q}_{p}{\rm{\Delta }}t}{2{m}_{p}c}\left({\hat{{\boldsymbol{v}}}}_{p}\cdot {{\boldsymbol{B}}}_{p}^{n}\right)\cdot {{\boldsymbol{B}}}_{p}^{n}\right),\end{array}\end{eqnarray} \tag{ 27 }$$

where we have dropped heavy notation (e.g., ${{\boldsymbol{B}}}_{p}^{n}\left({\bar{{\boldsymbol{x}}}}^{k}\right)$ is now written as ${{\boldsymbol{B}}}_{p}^{n}$) and defined

$$\begin{eqnarray}&&{\hat{{\boldsymbol{v}}}}_{p}={{\boldsymbol{v}}}_{p}^{n}+\displaystyle \frac{{q}_{p}{\rm{\Delta }}t}{2{m}_{p}}{{\boldsymbol{E}}}_{p}^{n+\theta }.\end{eqnarray} \tag{ 28 }$$

Notice that the fourth and fifth terms on the left-hand side can be neglected, since they are proportional to ${U}_{0}/c\ll 1$. Hence, the formula for the intermediate velocity is obtained:

$$\begin{eqnarray}&&{\bar{{\boldsymbol{v}}}}_{p}=\displaystyle \frac{{\hat{{\boldsymbol{v}}}}_{p}+\tfrac{{q}_{p}{\rm{\Delta }}t}{2{m}_{p}c}\left({\hat{{\boldsymbol{v}}}}_{p}\times {{\boldsymbol{B}}}_{p}^{n}+\tfrac{{q}_{p}{\rm{\Delta }}t}{2{m}_{p}c}\left({\hat{{\boldsymbol{v}}}}_{p}\cdot {{\boldsymbol{B}}}_{p}^{n}\right)\cdot {{\boldsymbol{B}}}_{p}^{n}\right)}{1+\tfrac{{q}_{p}^{2}{\rm{\Delta }}{t}^{2}}{4{m}_{p}^{2}{c}^{2}}{\left|{B}_{p}^{n}\right|}^{2}+\tfrac{{\rm{\Delta }}t}{2}{\overline{\overline{{\boldsymbol{P}}}}}_{v}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}.\end{eqnarray} \tag{ 29 }$$

Taking advantage of Equations (28), Equation (29) can be conveniently separated as

$$\begin{eqnarray}&&{\bar{{\boldsymbol{v}}}}_{p}={\tilde{{\boldsymbol{v}}}}_{p}+{\beta }_{p}{\tilde{{\boldsymbol{E}}}}^{n+\theta },\end{eqnarray} \tag{ 30 }$$

with

$$\begin{eqnarray}&&{\beta }_{p}=\displaystyle \frac{{q}_{p}{\rm{\Delta }}t}{2{m}_{p}},\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{v}}}}_{p}={{\boldsymbol{R}}}_{p}\cdot {{\boldsymbol{v}}}_{p},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{E}}}}_{p}^{n+\theta }={{\boldsymbol{R}}}_{p}\cdot {{\boldsymbol{E}}}^{n+\theta },\end{eqnarray} \tag{ 33 }$$

and

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{p}\cdot =\displaystyle \frac{\overline{\overline{{\boldsymbol{I}}}}-\tfrac{{\beta }_{p}}{c}\overline{\overline{{\boldsymbol{I}}}}\times {{\boldsymbol{B}}}_{p}^{n}+\tfrac{{\beta }_{p}^{2}}{{c}^{2}}{{\boldsymbol{B}}}_{p}^{n}\cdot {{\boldsymbol{B}}}_{p}^{n}}{1+\tfrac{{\beta }_{p}^{2}}{{c}^{2}}{\left|{B}_{p}^{n}\right|}^{2}+\tfrac{{\rm{\Delta }}t}{2}{\overline{\overline{{\boldsymbol{P}}}}}_{v}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}\cdot \end{eqnarray} \tag{ 34 }$$

The ${\tilde{{\boldsymbol{v}}}}_{p}$ component of the intermediate velocity is calculated from known quantities, while the ${\tilde{{\boldsymbol{E}}}}_{p}^{n+\theta }$ part contains the unknown ${{\boldsymbol{E}}}^{n+\theta }$ term. Notice, in this regard, that ${\overline{\overline{{\boldsymbol{P}}}}}_{v}^{n+\theta }\cdot$ can be easily calculated from the x-position of the grid at the previous time step. Since we assume that the simulated plasma parcel (i.e., the simulated grid size) is small with respect to the distance from the Sun, ${\boldsymbol{R}}$, ${\overline{\overline{{\boldsymbol{P}}}}}_{v}^{n+\theta }\cdot$ has the same value for all grid points.

3.4. Derivation of the Function for the Time Update of the Electric Field as a Function of the Approximated Moments

After calculating ${\bar{{\boldsymbol{v}}}}_{p}$ in Section 3.3, we are now equipped with the necessary information to solve the equation for ${{\boldsymbol{E}}}^{n+\theta }$.

We can preliminarily see that ${{\boldsymbol{J}}}^{n+\theta }$ can be rewritten as

$$\begin{eqnarray}&&{{\boldsymbol{J}}}^{n+\theta }={\hat{{\boldsymbol{J}}}}^{n}+\displaystyle \frac{1}{4\pi \theta {\rm{\Delta }}t}{\boldsymbol{\mu }}\cdot {{\boldsymbol{E}}}^{n+\theta },\end{eqnarray} \tag{ 35 }$$

where

$$\begin{eqnarray}&&{\hat{{\boldsymbol{J}}}}^{n}=\displaystyle \sum _{p}{q}_{p}{\tilde{{\boldsymbol{v}}}}_{p}S({\boldsymbol{x}}-{{\boldsymbol{x}}}^{n})-\theta {\rm{\Delta }}t{\rm{\nabla }}\cdot \displaystyle \sum _{p}{q}_{p}{\tilde{{\boldsymbol{v}}}}_{p}{\tilde{{\boldsymbol{v}}}}_{p}S({\boldsymbol{x}}-{{\boldsymbol{x}}}^{n}).\end{eqnarray} \tag{ 36 }$$

${\bar{{\boldsymbol{v}}}}_{p}{\bar{{\boldsymbol{v}}}}_{p}$ is approximated as ${\tilde{{\boldsymbol{v}}}}_{p}{\tilde{{\boldsymbol{v}}}}_{p}$ in the last term of Equation (36), neglecting the electric field contribution at this stage and removing the dependence of ${\hat{{\boldsymbol{J}}}}^{n}$ from the fields at time step n + θ.

The operator ${\boldsymbol{\mu }}$ is

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{\mu }}\cdot & = & \displaystyle \sum _{s}{{\boldsymbol{\mu }}}_{s}\cdot \\ {{\boldsymbol{\mu }}}_{s}\cdot & = & 4\pi \theta {\rm{\Delta }}t{\beta }_{s}{\rho }_{s}\displaystyle \sum _{p}{{\boldsymbol{R}}}_{p}\cdot ,\end{array}\end{eqnarray} \tag{ 37 }$$

where ρ_s is the density of the particle species s and we have taken advantage of the fact that the mass and charge are the same for particles belonging to the same species.

Similarly, we can identify in ${\rho }^{n+\theta }$ one component dependent on and another independent of ${{\boldsymbol{E}}}^{n+\theta }$. Defining the independent component as

$$\begin{eqnarray}&&{\hat{\rho }}^{n}={\rho }^{n}-{\rm{\Delta }}t\theta {\rm{\nabla }}\cdot {\hat{{\boldsymbol{J}}}}^{n}\end{eqnarray} \tag{ 38 }$$

and appropriately inserting the previous equations into Equation (16), we obtain an equation for the time update of the electric field:

$$\begin{eqnarray}&&\begin{array}{l}\left(\overline{\overline{{\boldsymbol{I}}}}+\displaystyle \frac{{\boldsymbol{\mu }}}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}\right)\cdot {{\boldsymbol{E}}}^{n+\theta }-{\left(\displaystyle \frac{c\theta {\rm{\Delta }}t}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}\right)}^{2}\\ \times \ \left({{\rm{\nabla }}}^{2}{{\boldsymbol{E}}}^{n+\theta }+\displaystyle \frac{{\rm{\nabla }}{\rm{\nabla }}\cdot \left({\boldsymbol{\mu }}\cdot {{\boldsymbol{E}}}^{n+\theta }\right)}{1+\theta {dt}{\rm{\nabla }}\cdot {\boldsymbol{U}}}\right)\\ =\ \displaystyle \frac{{{\boldsymbol{E}}}^{n}}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}+\displaystyle \frac{c\theta {\rm{\Delta }}t}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}\\ \times \ \left(\displaystyle \frac{{\rm{\nabla }}\times {{\boldsymbol{B}}}^{n}}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}-\displaystyle \frac{4\pi }{c}{\hat{{\boldsymbol{J}}}}^{n}\right)\\ -\ {\left(\displaystyle \frac{c\theta {\rm{\Delta }}t}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}\right)}^{2}4\pi {\rm{\nabla }}\displaystyle \frac{{\hat{\rho }}^{n}}{1+\theta {dt}{\rm{\nabla }}\cdot {\boldsymbol{U}}}.\end{array}\end{eqnarray} \tag{ 39 }$$

The temporal update of the magnetic field is obtained from Equation (11), conveniently rearranged after substituting ${{\boldsymbol{B}}}^{n+\theta }$ on the right-hand side following Equation (10):

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{B}}}^{n+1} & = & {{\boldsymbol{B}}}^{n}-\displaystyle \frac{{\rm{\Delta }}t}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot {{\boldsymbol{B}}}^{n}\\ & & -\ c{\rm{\Delta }}t\displaystyle \frac{{\rm{\nabla }}\times {{\boldsymbol{E}}}^{n+\theta }}{\overline{\overline{{\boldsymbol{I}}}}+\theta {\rm{\Delta }}t{\overline{\overline{{\boldsymbol{P}}}}}^{n+\theta }\cdot \overline{\overline{{\boldsymbol{I}}}}}.\end{array}\end{eqnarray} \tag{ 40 }$$

In this section, we test EB-iPic3D against expected plasma evolution in the EBM case. In particular, we test the evolution of the ion and electron thermal velocity and magnetic field (when present) in simulations without (Section 4.1) and with (Section 4.2) initial magnetic field, initialized in different directions. The expected evolution with R is recovered in all cases. Then, in Section 4.3, we focus on particle gyration in the case of radial and perpendicular magnetic field.

Before proceeding, it should be remarked that the parameters used in these tests (thermal and solar wind velocities, magnetic field intensity, distance from the Sun, etc.) are not necessarily realistic, since the aim at this stage is to test the implementation of the method rather than simulating solar wind plasma processes.

4.1. Nonmagnetized Plasmas

In this section, we test the evolution of the three components of the electron (e) and ion (i) thermal velocity in the direction parallel (${v}_{\mathrm{th},x}$) and perpendicular (${v}_{\mathrm{th},y}$, ${v}_{\mathrm{th},z}$) to the radial direction (x) in the case of nonmagnetized plasmas. 1D3V results are shown for the sake of simplicity at this stage.

Table 1 lists the parameters used for the different runs, with ${v}_{\mathrm{the},0}/c=0.01,0.1$ initial electron thermal velocity normalized to the speed of light c, R₀/d_e = 100, 1000 initial distance from the Sun normalized to the electron skin depth d_e, and U₀/c = 0.01, 0.1 plasma parcel velocity normalized to c. Notice, as remarked earlier, that these are convenient test values that do not reflect solar wind properties. Runs NoEB1 and NoEB2 are control runs, where the EBM is not used. For all simulations, the box length in the simulated x-direction is L_x/d_e = 10, 512 cells are used, the temporal resolution is ${\omega }_{\mathrm{pe}}{dt}=0.1$, with ω_pe the electron plasma frequency, T_i = T_e, with T the temperature, and 1024 particles per cell and a realistic mass ratio m_r = 1836 are used. The time-decentering parameter is θ = 0.5. A plasma composed by an electron and ion population with temperature ${T}_{s}={m}_{s}{v}_{{ths},0}^{2}$ and Maxwellian distribution function is loaded in the absence of background electric and magnetic field. The evolution of the thermal velocities is then followed, in the cases without (NoEB) and with (EB) expanding box, for a duration of ${\omega }_{\mathrm{pe}}t=1000$. Periodic boundary conditions are used for both fields and particles.

Table 1.  Parameters Used in the Simulations Analyzed in Section 4.1: Initial Electron Thermal Velocity Normalized to the Speed of Light, ${v}_{\mathrm{the},0}/c$, Initial Distance from the Sun Normalized to the Electron Skin Depth, R₀/d_e, Plasma Parcel Velocity Normalized to the Speed of Light, U₀/c

	${v}_{\mathrm{the},0}/c$	R0/de	U0/c
NoEB1	0.01	N/A	N/A
NoEB2	0.1	N/A	N/A
EB 1-1	0.01	1000	0.01
EB 1-2	0.01	100	0.01
EB 1-3	0.01	1000	0.1
EB 1-4	0.01	100	0.1
EB 2-1	0.1	1000	0.01
EB 2-2	0.1	100	0.01
EB 2-3	0.1	1000	0.1
EB 2-4	0.1	100	0.1

Download table as:  ASCIITypeset image

From the runs without EBM, we obtain a measure of the level of numerical cooling affecting the thermal velocity of electrons in the simulated direction, x. Numerical cooling is a well-known process in semi-implicit PIC codes (the analogous of numerical heating in explicit codes) and can be controlled by improving phase-space resolution. In the cases simulated, numerical cooling minimally reduces the electron thermal velocity in the simulated (radial) direction. We measure the level of numerical cooling by fitting the electron distribution functions obtained from the simulation with a Maxwellian distribution and examining the fitted thermal velocity values obtained in all directions, x, y, z, as time progresses. We see that radial thermal electron velocity (ions are not affected) at the end of the two non-EBM simulations is 99.75% and 99.15% of the respective initial values: the cooling is negligible in absolute value and is in any case much smaller than the cooling "legitimately" produced in the nonradial direction, y and z, by the expansion, when present.

In Figures 2 and 3 we show the variation of the electron and ion thermal velocity in the nonradial direction y, ${v}_{\mathrm{th},y}/c$, as a function of R/d_e for the EBM simulations listed in Table 1. As already done before, the thermal velocity values are obtained by fitting with a Maxwellian distribution the VDF obtained from the simulations at different times. Here, only ${v}_{\mathrm{th},y}$ is depicted; analogous results are obtained for the velocity in the other nonradial direction, ${v}_{\mathrm{th},z}$. The simulations with the same R₀/d_e and U₀/c but different initial thermal velocity ${v}_{\mathrm{the},0}/c$ are plotted in the same panel. They share the x-axis, while the y-axis is different (red curves map to the red, left axis, black curves to the right, black one). Each curve is fitted with a y = a/R curve (plotted with circular markers), and the fitted a parameters with confidence bars are reported in each case in the legend.

Zoom In Zoom Out Reset image size

As expected from the second line of Equation (9) in the absence of electric and magnetic field, the perpendicular velocity decreases in all cases as ${v}_{\mathrm{th}}={v}_{\mathrm{th},0}{R}_{0}/R$ (notice that the a parameter of the fits corresponds, within confidence bars, to ${v}_{\mathrm{th},0}{R}_{0}$ in the different cases), since there is no transfer of energy between waves and particles and ${P}_{v,{yy}}={P}_{v,{zz}}={U}_{0}/R$. ${P}_{v,{xx}}=0$. The velocity in the radial direction is not expected to change, and this is seen (within the margin defined by numerical cooling).

All simulations last ${\omega }_{\mathrm{pe}}t=1000;$ as expected from the initial parameters, a greater decrease in perpendicular velocities is registered in cases with lower R₀/d_e and higher U₀/c, depicted in panel (d) of Figures 2 and 3: at the end of the runs, the perpendicular velocity has dropped by almost 50% with respect to its initial value, as expected.

Zoom In Zoom Out Reset image size

Figure 4 shows the electron (panels (a), (b)) and ion (panels (c), (d)) v_x/c versus v_y/c velocity distribution at the beginning (panels (a), (c), ${\omega }_{\mathrm{pe}}t=0$, R/d_e = 100) and end (panel (b), (d), ${\omega }_{\mathrm{pe}}t=1000$, R/d_e = 200) of the simulation for the case with U₀/c = 0.1, R₀/d_e, v_the,0/c = 0.01. The black circle shows the Maxwellian fit used to calculate the thermal velocity in the different directions. The perpendicular cooling due to the expansion is quite visible at the end of the simulation for both electrons and ions.

Zoom In Zoom Out Reset image size

4.2. Magnetized Plasmas: Global Evolution

In this section, we analyze simulations of magnetized plasmas with initial magnetic field of equal magnitude and lying in the radial (B_x) or nonradial (B_y) direction.

We run two simulations initialized with the same initial parameters as EB 1-4 in Table 1, U₀/c = 0.1, R₀/d_e = 100, and v_the,0/c = 0.01. We add an initial magnetic field with magnitude ${B}_{0}/{\left(4\pi {m}_{e}{c}^{2}{n}_{0}\right)}^{1/2}=0.2$ in the radial, x-direction (run EB Bx) and in the perpendicular, y-direction (run EB By). ${\left(4\pi {m}_{e}{c}^{2}{n}_{0}\right)}^{1/2}$ is the normalization of the magnetic field in code units, where n₀ is the density at the beginning of the simulation. The density decreases with the heliospheric distance as $n={n}_{0}{R}_{0}^{2}/{R}^{2}$.

In Figures 5(a) and (b) we show the evolution of the B_x and B_y field components as a function of R/d_e for runs EB Bx and EB By. The solid lines are the measured values, while the superimposed lines with circular markers are fits whose dependence on R and coefficients are reported in the legend. One can see that B_x and B_y evolve (as expected from Equations (7)–(3)) as ${B}_{x}={B}_{x0}{R}_{0}^{2}/{R}^{2}$ and ${B}_{y}={B}_{y0}{R}_{0}/R$, respectively, since P_xx = 2U₀/R and P_yy = P_zz = U₀/R in Equations (7)–(3). Notice that $a=1999.67\sim {B}_{x0}{R}_{0}^{2}$ in the fit in panel (a) and $a=19.99\sim {B}_{y0}{R}_{0}$ in the fit in panel (b).

Zoom In Zoom Out Reset image size

In Figure 6 we report the evolution of the three components of the electron thermal velocity in simulations EB Bx (panel (a)) and EB By (panel (b)). We see that in the radial B case, panel (a), the x component of the thermal velocity (which is both parallel to the magnetic field and radial; blue) is not affected by the expansion, while the perpendicular components (perpendicular to both the magnetic field and the radial direction; red and black) decay with the usual ${v}_{\mathrm{the},\mathrm{perp}}={v}_{\mathrm{the}0,\mathrm{perp}}{R}_{0}/R$ dependence. The fit is superimposed with circular markers.

Zoom In Zoom Out Reset image size

In the nonradial B case, panel (b), the ${v}_{\mathrm{the},y}$ (red) velocity component is perpendicular to the radial direction but parallel to the magnetic field: it decays with the ${v}_{\mathrm{the},y}={v}_{\mathrm{the}0,y}{R}_{0}/R$ dependence expected for velocity in the nonradial direction (the fit is superimposed, red line with circular markers). The velocity in the directions perpendicular to the magnetic field, x (blue) and z (black), both evolve with a $1/\sqrt{R}$ dependence. This behavior is due to the conservation of magnetic moment: without magnetic field, ${v}_{\mathrm{the},x}$ would be unaltered, and ${v}_{\mathrm{the},z}$ would decay as R⁻¹. But the expansion is slow, compared to the cyclotron time, and the particle should spend only about 50% of the time moving in each direction. Then, both components end up decreasing at the "compromise" rate ${v}_{\mathrm{the},\mathrm{perp}}={v}_{\mathrm{the}0,\mathrm{perp}}\sqrt{B/{B}_{0}}\propto 1/\sqrt{R}$: the blue line with circular markers in Figure 6(b), ${v}_{\mathrm{the}0,\mathrm{perp}}\sqrt{B/{B}_{0}}$, shows the expected perpendicular velocity evolution from conservation of magnetic moment in the system. One can see that it superimposes excellently with the simulated values in the x- and z-direction.

4.3. Magnetized Plasmas: Particle Gyration

In Section 4.2, we focused on global characteristics of magnetized simulations, such as the evolution with R of the magnetic field and of the electron thermal velocities. Here we shift our attention to the motion of single electrons, captured by saving at high frequency the position and velocity of a randomly chosen subset of the simulated particles.

We start, for the sake of simplicity, from the radial B case EB Bx.

Figure 7 shows the trajectory in the perpendicular yz-plane (panel (a)) and in the 3D space (panel (b)) of an electron from run EB Bx. The electron, chosen randomly among the electron population, has initial coordinates in phase space v_x/c = 0.0073, v_y/c = −0.0248, v_z/c = −0.0034 (initial perpendicular velocity v_perp/c = 0.0250), x/d_e = 1.0647, y/d_e = 0.4975, and z/d_e = 0.4997.

Zoom In Zoom Out Reset image size

Figure 7(a) (blue line) shows the gyration motion in the perpendicular plane, fitted with a red circle with center $y/{d}_{e}=0.5147$, z_c/d_e = 0.3757 (red circle) and radius r/d_e = 0.1260 (the expected initial gyroradius). The black circle marks the position of the electron at the beginning of the run, R/d_e = 100, ${\omega }_{\mathrm{pe}}t=0$. One can notice that the gyroradius does not change for the duration of the simulation.

The electron trajectory in the 3D space (panel (b)) shows gyration motion with fixed gyroradius and increasing gyroperiod (the increase in the gyroperiod is more visible plotting the y or z trace as a function of time, not shown here).

During the simulated time, considering that the magnetic field (and consequently the gyrofrequency) decreases as (R/R₀)⁻² (see Figure 5(a), radial B field) and that the perpendicular thermal velocities, ${v}_{\mathrm{th},y}$ and ${v}_{\mathrm{th},z}$, decrease as (R/R₀)⁻¹ (see Figure 6(a), radial B field), one might expect the gyroradius to increase linearly with R.

In the simulation (comoving) frame this does not happen, but when the coordinates of the plane perpendicular to the radial direction are transformed back from comoving (O) to the Sun-centered (O') system, ${{\boldsymbol{x}}}_{\mathrm{perp}}^{{\prime} }=R/{R}_{0}{{\boldsymbol{x}}}_{\mathrm{perp}}$, the gyroradius in the Sun-centered system increases, as expected, as R/R₀. This is illustrated in Figures 7(c) and (d), which are otherwise analogous to panels (a) and (b). The red line in Figures 7(c) and (d) depicts the motion of the center of gyration in the Sun-centered system. If the center of gyration of the particle is not in the center of the comoving frame, the expansion of the perpendicular direction in the Sun-centered system translates into a drift of the gyrocenter, which moves along the ${y}_{c}^{{\prime} }\,={y}_{c}R(t)/{R}_{0}$, ${z}_{c}^{{\prime} }={z}_{c}R(t)/{R}_{0}$ line (Figures 7(c) and (d)).

In Figure 8 we show the trajectory in space (panels (a), (b)) and in velocity space (panels (c), (d)) of an electron with the same initial position and velocity as the electron depicted in Figure 7 (v_x/c = 0.0073, v_y/c = −0.0248, v_z/c = −0.0034, x/d_e = 1.0647, y/d_e = 0.4975, z/d_e = 0.4997), but this time belonging to run EB By. In panels (a) and (c) we depict the plane perpendicular to the magnetic field (xz plane), and in panels (b) and (d) the 3D space and velocity space, respectively. The traces recorded from the simulated particle are in blue. In red, we show the particle path resulting from the numerical integration of the equations of electron motion in the EBM, Equation (9). We have used as initial conditions for the integrated particle the position and velocity of the simulated particle, and we have assumed that the only field component different from zero is B_y, which varies with R as B_y = B_y0R₀/R (as expected from the equation for the magnetic field evolution and as observed in the simulation).

Zoom In Zoom Out Reset image size

One can see that the trajectories recorded from the simulations and integrated numerically differ significantly in space (we observe a drift of the gyrocenter in the z-direction), not so much in velocity space.

The reason for the discrepancy is clarified in Figure 9(a), analogous to Figure 8(b), with the difference that, instead of plotting the integrated solution for particle motion in a B_y = B_y0R₀/R field, we plot in black the integrated solution for a B_y = B_y0R₀/R field and an electric field in the radial direction of the form depicted in Figure 9(b) (blue line). Electric field oscillations in the simulated direction (x, in this case) are to be expected in a PIC simulation and do not constitute a numerical artifact, but a natural consequence of physical density oscillations through Gauss's law. The shape of the electric field introduced in the integration has been selected to match the E_x field sampled by the particle, ${E}_{x,p}$, recorded from the simulation (red line in panel (b)). The radial electric field experienced by each electron determines whether the electron deviates from the integrated B_y-only solution in the $+z$- or $-z$-direction.

Zoom In Zoom Out Reset image size

We observe that electric field oscillations in the simulated direction x develop also in the case of radial magnetic field. There, however, they do not impact the observed particle trajectory strongly, as evident from Figure 7 and its discussion. One possible reason for this may be that the electric field impact on particle trajectory is mediated by its perturbing the velocity in the radial direction, v_x. This, in turn, can affect the other velocity components through the ${\boldsymbol{v}}\times {\boldsymbol{B}}$ term, which is different from 0 in the EB ${{\rm{B}}}_{{\rm{y}}}$ but not in the EB ${{\rm{B}}}_{{\rm{x}}}$ case.

In this paper, we have introduced a new tool for fully kinetic simulations of the expanding solar wind. The solar wind constitutes a challenge for kinetic PIC approaches for several reasons: as in most plasma environments, the characteristic temporal and spatial scales that one either may want to (e.g., electron or ion skin depth, d_e or d_i) or is forced (to resolve, e.g., Debye length λ_D) in a fully kinetic simulation are spread over several orders of magnitude. This is seen in Figure 10, which shows the characteristic lengths and times in solar wind plasma at 1 au. There, together with the characteristic scales, we have plotted a realistic but very ambitious (and unattainable) system size and process duration of the order of 1 au and 3 days, respectively.

Zoom In Zoom Out Reset image size

Two particular challenges, however, are characteristic of the solar wind. First, the system size is larger, in terms of electron or ion skin depths, than most of the other plasma environments, such as planetary magnetospheres, which are already a challenging feat for kinetic simulations. Second, fundamental solar wind processes, i.e., the propagation away from the Sun and radial expansion of solar wind parcels, increase the size of the domain that one has to simulate in order to capture the fundamental dynamics. This introduces new characteristic lengths in the system. Returning to Figure 10, one sees that the characteristic electron (blue) and ion (red) length scales (depicted as vertical lines in the length vs. timescale plot) are three, rather than the usual two, per species, i.e., skin depths (solid lines) and gyroradii (dotted lines). The third length scale is the distance that a plasma parcel (moving at a characteristic velocity V, here V = 800 km s⁻¹, the average fast wind speed) travels in a characteristic timescale. Taking as characteristic timescales the electron and ion gyroperiods, we obtain the thick dashed blue and red lines, respectively, for electron and ions. While for electrons this new characteristic length scale is shorter than the others (in particular, V/Ω_ce < r_e, since ${v}_{\mathrm{th},e}\sim 2000\,\mathrm{km}\,{{\rm{s}}}^{-1}\gg V$), this is not the case for ions, for which $V\gg {v}_{\mathrm{th},i}\sim 35\,\mathrm{km}\,{{\rm{s}}}^{-1}$. This means that the window that one needs to simulate to capture, together, the relevant electron and ion scale processes becomes larger, due to V/Ω_ce pushing for higher resolution and V/Ω_ci for larger domain sizes, thus increasing computational cost. The EBM formulation shrinks the window again with the use of the comoving frame, which hides the V/Ω_ce, V/Ω_ci length scales in the variable change used to define the comoving frame from the Sun-centered frame. Also, one has to notice, as already commented on in Section 1, that global kinetic solar wind simulations, going up to fractions of L, become significantly cheaper with EBM rather than without, since the gridding to be simulated covers only the moving parcel, rather than the entire distance L, hence the need for using, together, a semi-implicit algorithm and the EBM method in solar wind kinetic simulations.

In Figure 10, we show as a black rectangle the scales that we address in the simulations shown in this paper. The left and bottom sides of the rectangle are the spatial and temporal resolutions, respectively (${dx}=0.02{d}_{e}$, ${\omega }_{\mathrm{pe}}{dt}=0.1$); the right and top sides are the domain size in space and time, respectively (L = 10d_e, Tω_pe = 1000). These simulations are easy feats in a 1D3V configuration, and in fact they run on a laptop. In blue and red we show the same rectangles for typical production runs done with the iPic3D code we started from in developing EB-iPic3D. The blue and red rectangles correspond to simulations addressing processes at the electron and ion scale, respectively. Here, temporal and spatial resolutions are fractions of the respective inverse plasma frequency and skin depth, and domain sizes are hundreds of skin depths and thousands of the inverse plasma frequency. Simulations encompassing these kinds of scales can be done on a small cluster in 1D. In 2D, they require high-performance computing (HPC) resources of the order of tens or hundreds of thousands of core-hours. 3D analogous simulations would require one to shrink box dimensions and duration and reduce spatial and temporal resolution.

The EB evolution does not alter significantly the performance of the iPic3D code in terms of execution time with respect to the non-EB version. Hence, we expect to be able to reach the same resolutions, box sizes, and simulation duration with EB-iPic3D.

Separate considerations should be made regarding the typical EB parameters, U₀ and R₀, and the relation between the expansion time, ${\tau }_{\exp }=R/{U}_{0}$, and the other characteristic timescales. The EBM is valid if the characteristic process of the system is fast compared to the characteristic time of the expansion, i.e., if ${\omega }_{\mathrm{ch}}\gt 1/{\tau }_{\exp }$. This inequality is very easy to satisfy, e.g., for electron (ion) scale processes with ω_ch ∼ ω_pe (ω_pi), R/d_e (R/d_i) of the order of tens or hundreds and U₀ a small fraction of the speed of light. It is in fact too easy to satisfy with realistic parameters, namely, U₀/c ∼ 2 × 10⁻³: one has to wait "too long" to see the effects of the expansion in a realistic simulation.

For this reason, in the work presented here and in future work, we choose to compress the temporal scales by resorting to using U₀/c and R/d_e (R/d_i) values higher and lower, respectively, than in reality. This artificially increases the $1/{\tau }_{\exp }$ simulated and makes it closer to the characteristic frequency of the processes of interest. This operation is somehow analogous to one widely used in kinetic PIC simulations, i.e., decreasing the simulated ion to electron mass ratio with respect to the real one to compress ion and electrons scales.

To sum up, we address the particular challenges of kinetic solar wind simulations with a two-pronged approach.

On the one hand, we use the EBM, recapped in Section 2, to address solar wind propagation and expansion. The EBM allows us to follow the evolution of a solar wind plasma parcel as it moves away from the Sun and expands in the nonradial directions. A fixed velocity in the radial direction is assumed (but this restriction can be removed and an accelerating EB can be introduced instead, as done in Tenerani & Velli 2017 in the MHD context), and the perpendicular directions expand at a rate proportional to the distance from the Sun.

The EBM is implemented here for the first time in a semi-implicit PIC code, iPic3D. The derivation of the EBM IMM equations is recapped in Section 3. We use a semi-implicit PIC approach to be able to tailor the spatial and temporal resolution of the simulation to the specific process of interest, rather than having to resolve the fastest and smallest scales in the system out of stability concerns. This will allow us to increase the domain size and temporal duration of our simulations with respect to what is feasible with an explicit PIC code. This will limit the severity of the "domain size" challenge.

In Section 4, we test the EBM-IMM code EB-iPic3D against expected plasma evolution in EBM case. Our testing activities focus on the integration of the EBM into the preexisting IMM code rather than on testing the IMM code itself. This is because the IMM algorithm has been thoroughly validated since its early formulations more than 30 years ago (Brackbill & Forslund 1982), including by comparison against explicit codes (e.g., Ricci et al. 2004).

In Sections 4.1 and 4.2 we compare the evolution of our simulations against the expected thermal velocity (Sections 4.1 and 4.2) and magnetic field (Section 4.2) evolution with the distance from the Sun, R. The "expected evolution" is simply obtained from the EBM equations governing the velocity and magnetic field temporal update, as seen in the second line of Equation (9) and third line of Equation (7), discretized semi-implicitly in the second line of Equation (26) and (11). In both sections, electrons and ions are initially distributed following a Maxwellian VDF, which is altered during the course of the simulation. In Section 4.1 the plasma is not magnetized. A radial and nonradial magnetic field is then introduced in Section 4.2.

In Section 4.3 we study how solar wind expansion affects electron gyration around a decaying radial and nonradial magnetic field. We observe that, in the radial case, the gyroradius increases as expected with the heliospheric distance if the particle trajectory in space and velocity space is transformed from the comoving to the Sun-centered frame.

In the case of nonradial magnetic field, we observe a significant deviation of the electron trajectory in space from the integrated solution obtained assuming that the only field acting on the particle is a nonradial magnetic field exhibiting the dependence on the heliocentric distance expected from the magnetic field equation and effectively observed in the simulation. We explain this deviation in terms of the radial (x) velocity oscillations provoked by radial density oscillations, which give rise to radial electric field oscillations through Gauss's law. Radial velocity oscillations are able to significantly affect particle motion in nonradial directions in the nonradial magnetic field case because they enter the force equation through the ${\boldsymbol{v}}\times {\boldsymbol{B}}$ term, which is nonzero with nonradial magnetic field. The deviation of electron trajectory with respect to the integrated, B_y-only solution can be interpreted physically as the electrons "jumping" from a magnetic field line to another, i.e., diffusing across magnetic field lines. This process appears to be easier, in a 1D3V EBM geometry, in the case of nonradial rather than radial magnetic field for the reasons just mentioned. This can be taken as a small example of a wider range of processes captured by self-consistent PIC simulations, but not immediately evident from analytic calculations.

As a final remark, we would like to remind the reader that, since the aim of this paper is to test EB-iPic3D rather than to simulate reality, we have kept the problems artificially simple by numerically eliminating instabilities that should occur in the simulations as a result of the expansion-driven velocity and magnetic field evolution. The simulation of those same instabilities in expanding plasmas will be addressed in future work.

A temperature anisotropy such as the one registered at R/d_e = 200 in Section 4.1, unmagnetized plasmas, should give rise to a Weibel instability. The Weibel instability occurs in the presence of a thermal anisotropy and has a wavenumber in the "cold" direction (Fried 1959; Weibel 1959; Innocenti et al. 2011). In our simulations, the Weibel instability is suppressed by using a 1D3V reduced geometry: only instabilities with wavenumber in the simulated x-direction can develop, and cooling due to expansion affects the perpendicular y- and z-directions. Hence, the Weibel instability is prevented from occurring.

In Section 4.2, on magnetized plasma, we have suppressed the onset of the electron firehose instability, which develops in magnetized plasmas when the parallel temperature exceeds the perpendicular temperature, as in our case with EBM. The electron firehose instability is mostly parallel to the magnetic field and can therefore develop also in our reduced geometry in the radial B case. However, it has a wavenumber of the order ${{kd}}_{i}\gt 1$ (Gary 2005), which corresponds to relatively long wavelengths in our case, where the mass ratio is realistic (and therefore ${d}_{i}/{d}_{e}=\sqrt{1836}$) and the simulation box is L_x = 10d_e long. Larger simulation boxes allow the firehose instability to develop. With this relatively small box, the instability is artificially suppressed.

In future work, we intend to harness the full potential of the EBM-IMM EB-iPic3D code by addressing the large-scale, longtime nonlinear evolution of thermal-anisotropy-driven kinetic instabilities. Specifically, we will be able to simulate the dynamics of the core-strahl-halo electron distribution function and associated heat-flux-driven Alfvén instabilities constraining the electron core drift. Our code will be used to simulate the stability, for typical solar wind conditions, of core-strahl electron distributions to kinetic Alfvén and magnetosonic modes and whistler modes, and, at intermediate scales, the possible coupling of proton-distribution-function-driven instabilities, thus providing an understanding of kinetic-scale turbulence drivers and dissipation in the expanding solar wind.

M.E.I.'s research is supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, administered by Universities Space Research Association under contract with NASA. M.V. and A.T. acknowledge support from the PSP Observatory Scientist grant NASA SPP NNX15AF34G. M.E.I. is an FWO postdoctoral fellow currently on leave.