the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A deep insight into the ion foreshock with the help of test particle twodimensional simulations
Philippe Savoini
Bertrand Lembège
Twodimensional (2D) test particle simulations based on shock profiles issued from 2D full particleincell (PIC) simulations are used in order to analyze the formation processes of ions back streaming within the upstream region after interacting with a quasiperpendicular curved shock front. Two different types of simulations have been performed based on (i) a fully consistent expansion (FCE) model, which includes all selfconsistent shock profiles at different times, and (ii) a homothetic expansion (HE) model in which shock profiles are fixed at certain times and artificially expanded in space. The comparison of both configurations allows one to analyze the impact of the front nonstationarity on the backstreaming population. Moreover, the role of the space charge electric field E_{l} is analyzed by either including or canceling the E_{l} component in the simulations. A detailed comparison of these last two different configurations allows one to show that this E_{l} component plays a key role in the ion reflection process within the whole quasiperpendicular propagation range. Simulations provide evidence that the different fieldaligned beam (FAB) and gyrophase bunched (GPB) populations observed in situ are essentially formed by a E_{t}×B drift in the velocity space involving the convective electric field E_{t}. Simultaneously, the study emphasizes (i) the essential action of the magnetic field component on the GPB population (i.e., mirror reflection) and (ii) the leading role of the convective field E_{t} in the FAB energy gain. In addition, the electrostatic field component E_{l} is essential for reflecting ions at high θ_{Bn} angles and, in particular, at the edge of the ion foreshock around 70^{∘}. Moreover, the HE model shows that the rate BI_{%} of backstreaming ions is strongly dependent on the shock front profile, which varies because of the shock front nonstationarity. In particular, reflected ions appear to escape periodically from the shock front as bursts with an occurrence time period associated to the selfreformation of the shock front.
While upstream ions of the incoming solar wind interact with the curved terrestrial bow shock, a certain percentage is reinjected back into the solar wind and propagates along the interplanetary magnetic field (IMF); they form the socalled ion foreshock. This population has been extensively studied both with the help of experimental data (Tsurutani and Rodriguez, 1981; Paschmann et al., 1981; Bonifazi and Moreno, 1981a, b; Fuselier, 1995; Eastwood et al., 2005; Oka et al., 2005; Kucharek, 2008; Hartinger et al., 2013) and numerical simulations (BlancoCano et al., 2009; Lembège et al., 2004; Savoini et al., 2013; Kempf et al., 2015; Savoini and Lembège, 2015; Otsuka et al., 2018).
Even if we restrict ourselves to the quasiperpendicular region (i.e., for $\mathrm{45}{}^{\circ}\le {\mathit{\theta}}_{\mathrm{Bn}}\le \mathrm{90}{}^{\circ}$, where θ_{Bn} is the angle between the local shock normal and the IMF), different types of backstreaming ions are identified, and (a) some are characterized by a gyrotropic velocity distribution and form the fieldaligned ion beam population (hereafter FAB); conversely, (b) others exhibit a nongyrotropic velocity distribution and form the gyrophase bunched ion population (hereafter GPB). None of these populations yet have a wellestablished origin, and different mechanisms have been proposed for years (Möbius et al., 2001; Kucharek et al., 2004), including (i) scenarios based on the specular reflection (Sonnerup, 1969; Paschmann et al., 1980; Schwartz et al., 1983; Schwartz and Burgess, 1984; Gosling et al., 1982) with or without the conservation of the magnetic moment and (ii) scenarios which invoke the leakage of some magnetosheath ions producing a lowenergy FAB population (Edmiston et al., 1982; Tanaka et al., 1983; Thomsen et al., 1983). Nevertheless, the origin of FAB ions could be due to (iii) the diffusion of some reflected ions (called gyrating ions; these ions are reflected by the supercritical shock front but do not manage to escape into the upstream region and go into the downstream region after their initial gyration (Schwartz et al., 1983)). The diffusion can either be generated by upstream magnetic fluctuations (Giacalone et al., 1994) or more directly by the shock ramp itself (with a pitch angle scattering during the reflection process) (Kucharek et al., 2004; Bale et al., 2005). All scenarios have some drawbacks and are not able to clearly explain the origin of both populations. On the other hand, GPB are preferentially observed at some distance from the curved shock front (Thomsen et al., 1985; Fuselier et al., 1986a), and their synchronized nongyrotropic distribution is part of a lowfrequency monochromatic waves trapping (Mazelle et al., 2003; Hamza et al., 2006), or of beam plasma instabilities (Hoshino and Terasawa, 1985). As a conclusion, it is quite difficult to discriminate between these different scenarios which can be present simultaneously or separately in time.
Our previous papers (Savoini et al., 2013; Savoini and Lembège, 2015) were focused on the origin of these two populations. A largescale, twodimensional particleincell (PIC) simulation of a curved shock has been used, where full curvature and timeofflight effects for both electrons and ions are selfconsistently included. Our simulations have shown that both FAB and GPB populations and their typical associated pitch angle distributions observed experimentally (Fuselier et al., 1986b; Meziane, 2005) have been retrieved not far from the front (up to 2–3R_{E}, where R_{E} is the Earth's radius). Moreover, results have shown that these two populations can be generated directly by the macroscopic electric E and magnetic B fields present at the shock front itself. In other words, the differences observed between FAB and GPB populations are not the result of distinct reflection processes but are the consequence of the time history of ions interacting with the shock front. The FAB population loses their initial phase coherency by suffering several bounces along the front, which is in contrast with the GPB population which suffers mainly one bounce (i.e., mirror reflection process). This important result was not expected and greatly simplifies the question on each population origin (Savoini and Lembège, 2015)
Nevertheless, some further questions, which are difficult to investigate with full PIC simulations (because of the selfconsistency), still need to be answered in order to analyze several aspects of the reflection process. For this reason, we use complementary test particle simulations herein to clarify the respective impact of the shock curvature and the time variation in the macroscopic fields at the shock front on the backstreaming ion reflection process. The main questions presently addressed are summarized as follows:

Is the reflection process noncontinuous in time (bursttype reflection process) or not? In this case, how is it linked to the θ_{Bn} angle variation (i.e., space dependence) and/or to the shock profile variation (i.e., time dependence)?

What is the impact of the space charge electric field localized within the shock ramp on the reflection process?

What kind of reflection mechanisms can be identified in present simulations?
The paper is organized as follows. Section 2 briefly summarizes the conditions of the previous 2D PIC simulations (Savoini and Lembège, 2015) and of present particle test simulations. In Sects. 3 and 4, results of test particles are presented, and the ion reflection processes are investigated. The discussion and conclusions will be presented in Sects. 5 and 6, respectively.
The numerical conditions concerned in the present paper are similar to those described in Savoini et al. (2013) and Savoini and Lembège (2015). In short, we used a 2D, fully electromagnetic, relativistic particle code based on a standard finitesized particle technique (similar to Lembège and Savoini, 1992, 2002 for planar shocks).
2.1 Selfconsistent full PIC simulations
The code solves Maxwell and Poisson's equations in the Fourier space (socalled pseudospectral code), which allows one to separate the electric field contribution in two distinct parts, namely (i) a longitudinal or electrostatic component, hereafter denoted by subscript _{l} (built up by the space charge effects $\mathbf{\nabla}{\mathit{E}}_{\mathrm{l}}=\mathit{\rho}/{\mathit{\u03f5}}_{o}$), and (ii) a transverse or induced component, hereafter denoted by subscript _{t} (coming from the temporal variations in the magnetic field $\mathbf{\nabla}\times {\mathit{E}}_{\mathrm{t}}=\partial \mathit{B}/\partial t$). The longitudinal component is essentially built up within the shock front due to the different dynamics of ions and electrons, whereas the induced component is mainly generated by the propagating shock front itself (see Fig. 1, panel 2a) through the convective term ${\mathit{E}}_{\mathrm{t}}=\mathit{U}\times \mathit{B}$ (where U corresponds to the bulk shock front velocity since we are in the solar wind frame). In addition, the subscripts _{∥} and _{⟂} stand for parallel and perpendicular directions to the local magnetic field, respectively. In Fig. 1 and in following, the X–Y reference frame is the solar wind frame, with the third direction along Z pointing backward into the plot. Then, E_{t} has the direction of the increasing Y, and ∇B has the same direction as the present U vector.
In our configuration, the magnetostatic field is partially lying outside the simulation plane (see Savoini and Lembège, 2001 for more details). Then, the simulation is limited to the whole quasiperpendicular shock (i.e., for $\mathrm{45}{}^{\circ}\le {\mathit{\theta}}_{\mathrm{Bn}}\le \mathrm{90}{}^{\circ}$). We use a magnetic piston for which the geometry is adapted to initiate a shock front with a curvature radius large enough to be compared with the upstream ion Larmor gyroradius ${\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}_{\mathrm{ci}}$ (all normalized quantities are indicated with a tilde $\stackrel{\mathrm{\u0303}}{\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}}$); this is the same normalization which is used in the previous selfconsistent PIC simulations (Savoini and Lembège, 2001; Savoini et al., 2013; Savoini and Lembège, 2015). The curvature increases during the simulation. This configuration has two consequences; (i) first, as the time increases and the shock front expands, its velocity slightly decreases and so does the Alfvén Mach Number M_{A} from ≈5 to ≈3, where the velocity is measured at ${\mathit{\theta}}_{\mathrm{Bn}}=\mathrm{90}{}^{\circ}$ used as a reference angle; (ii) the timeofflight effects are selfconsistently included. Indeed, this ballistic process is observed when the upstream magnetic field lines are convected by the incoming solar wind. In present simulations (based on upstream rest frame), the curved shock front expands and scans different θ_{Bn} values. As a result, backstreaming particles, collected at a given upstream location, come from different parts of the curved shock front, depending on their respective velocity.
Initial plasma conditions are summarized as follows: light velocity $\stackrel{\mathrm{\u0303}}{c}=\mathrm{3}$ and temperature ratio between ion and electron population ${T}_{\mathrm{el}}/{T}_{\mathrm{io}}=\mathrm{1.58}$. A mass ratio ${m}_{i}/{m}_{e}=\mathrm{84}$ is used in order to save CPU time, and the Alfvén velocity is ${\stackrel{\mathrm{\u0303}}{v}}_{\mathrm{A}}=\mathrm{0.16}$. The number of cells in the simulation plane along each axis is NCX = NCY $=\mathrm{8192}\approx \mathrm{150}{\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}_{\mathrm{ci}}$, with the size of a grid cell ${\mathrm{\Delta}}_{x}=\mathrm{\Delta}y\approx \mathrm{1}{\stackrel{\mathrm{\u0303}}{\mathit{\rho}}}_{ce}$. The shock is in the supercritical regime. with a timeaveraged Alfvén Mach number M_{A}≈4 measured at ${\mathit{\theta}}_{\mathrm{Bn}}=\mathrm{90}{}^{\circ}$. In order to observe the early stage of the ion foreshock formation, the end time of the simulation is ${\stackrel{\mathrm{\u0303}}{t}}_{\mathrm{simul}}=\mathrm{5.4}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$ (where ${\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$ is the upstream ion gyroperiod), which is large enough to investigate the interaction of incoming ions with the shock front and the further formation of backstreaming ions.
2.2 Test particle simulations
In the present paper, we use all field components issued from the same previous PIC simulation as in Savoini and Lembège (2015); all components have been saved every $\mathrm{\Delta}\stackrel{\mathrm{\u0303}}{T}={\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}/\mathrm{20}$. Test particle simulations appear to be a straightforward way to evaluate the action of different field components on the ion dynamics. Indeed, the feedback effects of particles on electromagnetic fields are excluded in test particle simulations, and one can modify or cancel some field components independently from each other. This allows one to identify their specific actions on particles and on the resulting ion reflection processes.
Figure 1 plots an example of the two configurations used hereafter in this paper. Panels 2a–b show the fully consistent expansion model (hereafter named FCE), which corresponds to the results where test ions interact with the E and B fields issued from the selfconsistent simulation and where both spatial inhomogeneities and nonstationarities are fully included. If this configuration is easy to understand, the socalled particular approach named the homothetic expansion model (hereafter named HE) shown in panels 3a–b is complementary. In this case, particles interact with a propagating fixed front profile, i.e., alltime profile variations are excluded; only spatial inhomogeneities of the shock front profile chosen at the selected time are included, as detailed in Sect. 4.
In the two configurations of FCE and HE, we inject test particles distributed within 10 individual sampling boxes located along the curved shock front (Fig. 2). This procedure allows us to analyze the impact of the front curvature (local θ_{Bn}) on the formation of backstreaming ions. We follow a total of 1 million test particles. Then, each box has the same number of particles N=100 000, and they are initialized as a Maxwellian distribution with a thermal velocity v_{thi}, which is the same as in the selfconsistent simulation (Savoini and Lembège, 2015).
Let us point out that the use of finitesized sampling boxes at different initial θ_{Bn} angles only estimates the location where the particles hit the shock but does not provide an exact value for the local θ_{Bn} seen by the particles when these interact with the expanding shock front. Nevertheless, it is quite helpful when classifying the different types of particle interactions with the curved front. The sizes of all the identical boxes are chosen so that (i) along the curved shock front each box has an angular extension of $\approx \mathrm{4}{}^{\circ}$, which is small enough to scan the different orientations of θ_{Bn} and large enough for statistical constraints, and (ii) along the local shock normal, where each box has a length large enough (L_{size}=2000Δx) to ensure that most particles interact with the shock front during a noticeable time range (i.e., ${D}_{T}\approx \mathrm{3}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$).
Then, Sect. 3 will present the results obtained in the FCE model, which is the usual configuration representing the time evolution of test particles with the selfconsistent shock front profile. Section 4 will then introduce the more unusual HE expansion model (i.e., homothetic simulation approach).
To clarify the presentation, we will split this section into two parts in which we analyze, respectively, (i) the dynamics of the backstreaming ions (alias BI) in the different boxes and their main features, in terms of spatial and time evolution, and (ii) their behavior when some field components at the shock front are included/artificially excluded.
3.1 General features of the backstreaming ions
Figure 3 plots the spatial distribution of the backstreaming ions' density for the different sampling boxes (defined in Fig. 2) at the end time of the simulation for the FCE model. Different information can be summarized as follows:
 i.
The percentage of the backstreaming ions BI_{%} is obtained by computing the ratio of the backstreaming ions over the number of ions which have interacted with the shock front. This number increases when moving further into the foreshock (i.e., for decreasing θ_{Bn}) from BI_{%}=0.1 for N_{Box}=0 to BI_{%}≈14 for N_{Box}=9. This θ_{Bn} dependence is in agreement with previous experimental observations (Ipavich et al., 1981; Eastwood et al., 2005; Mazelle et al., 2005; Turner et al., 2014) and numerical simulations (Savoini and Lembège, 2015; Kempf et al., 2015).
 ii.
The upstream edge of the ion foreshock (mixed line in Fig. 3 for N_{Box}=10) is not parallel to the IMF but is the result of the timeofflight effect included in our simulations (Savoini et al., 2013). At the end of the simulation, this edge starts from the shock at the same critical angle, the socalled ${\mathit{\theta}}_{\mathrm{io},\mathrm{fore}}\approx \mathrm{70}{}^{\circ}$, as found in our previous selfconsistent simulations (Savoini et al., 2013; Savoini and Lembège, 2015).
 iii.
The backstreaming ion density is not uniform along the shock normal but exhibits different maxima, according to our observations. Not only is the spatial distribution not the same for all boxes but is even not uniform within a given same box, i.e., backstreaming ions do not escape uniformly away from or along the shock front. For instance, boxes N_{Box}=0–3 evidence two distinct spots near the shock front, as indicated by black arrows. As θ_{Bn} decreases (i.e., N_{Box}=5–9), the righthand spot disappears and the backstreaming population increasingly aligns along the upstream magnetic field B_{o}. Accordingly, the width of the reflection area (i.e., the angular extension of the ion foreshock defined very near the shock front) shrinks from ≈50ρ_{ci} (for N_{Box}=0) to ≈17ρ_{ci} (for N_{Box}=9).
These two distinct spots may be explained by the different time histories of the backstreaming ions within the shock front, as reported in Savoini and Lembège (2015). Actually, the interaction time strongly differs from one ion to another, depending on its gyrating feature when it hits the shock front for the first time. Short and long interactions over time can be defined, depending on whether the reflection process is respectively associated to a short or long displacement of the ion along the shock front before escaping upstream. If the individual trajectory of the reflected ions has been already evidenced in Savoini and Lembège (2015), present test particle simulations allow one to generalize the results via a statistical approach versus their initial angular locations (i.e., the N_{Box} number). Figure 4 plots the distribution of θ_{Bn} angles seen by the particles when these hit for the first time the shock front (hereafter named ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{hit}}$ in red) and when they finally exit the shock front to escape upstream (hereafter named ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ in blue). These statistical results are obtained by computing these angles for each particle. As a consequence, angle values are computed neither at the same time nor at the same location along the curved front (even if they are initially located in the same box). In other words, each particle sees different local shock front profiles in terms of the spatial inhomogeneity and time nonstationarity of the shock front.
First, let us note that the averaged value of ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{hit}}$ corresponds mainly to the initial location of the box, and therefore, the most important feature is not the angle values themselves but rather the difference between the averaged values of ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{hit}}$ and ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ when ions hit and leave the shock front, respectively. For this reason, we will use the angular range of the particles' interaction with the front as defined by ${\mathrm{\Delta}}_{\mathrm{int}}{\mathit{\theta}}_{\mathrm{Bn}}={\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}{\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{hit}}$. Obviously, ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{hit}}$ decreases as N_{Box} increases until it approaches the limit of the quasiperpendicular domain of propagation, i.e., ${\mathit{\theta}}_{\mathrm{Bn}}=\mathrm{45}{}^{\circ}$ for N_{Box}=9.
Second, the distribution functions of ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ strongly differ according to the respective box. For N_{Box}=0–2, two (blue) peaks occur, namely one for high ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ (for which Δ_{int}θ_{Bn}≈4–5^{∘}) and the other for low ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ (${\mathrm{\Delta}}_{\mathrm{int}}{\mathit{\theta}}_{\mathrm{Bn}}\approx \mathrm{15}{}^{\circ}$). In terms of the time trajectory, the presence of these two peaks suggests that some ions have spent different interaction times (subscript _{int}) within the shock front. Some escape after a short interaction time (i.e., small Δ_{int}θ_{Bn}), while others escape after a long interaction time (i.e., large Δ_{int}θ_{Bn}), where the terms short and long refer to a small and large drift along the shock front, as already analyzed in Savoini and Lembège (2015). In other words, a small drift refers to one bounce whereas a large drift refers to multiple bounces along the shock front.
Moreover, as N_{Box} increases (i.e., N_{Box}≥3), the lower ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ distribution (i.e., correspondingly the largest Δθ_{Bn}) decreases rapidly in amplitude and disappears from N_{Box}=6 (i.e., ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{hit}}\le \mathrm{56}{}^{\circ}$). Simultaneously, the other peak (i.e., correspondingly the smaller Δ_{int}θ_{Bn}) becomes dominant for all higher order boxes, meaning that fewer and fewer ions are associated to large drifts along the shock front.
Third, in order to complete information deduced from Fig. 4, Fig. 5 plots the number of reflected ions versus the time spent within the shock front. This interaction time $\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{T}}_{\mathrm{int}}$ is defined as the time difference between the time associated to ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ and to ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{hit}}$. Different main maxima of backstreaming ions' density are evidenced, namely f_{1}, f_{2}; a third maximum f_{3} can be also observed for boxes N_{Box}=0–4, but its amplitude is too weak to be relevant in this discussion. One important feature is that f_{1} and f_{2} appear in all boxes and are independent of the box number. More precisely, f_{1} appears about $\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{T}}_{\mathrm{int}}\approx \mathrm{0.25}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}\approx {\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}^{\mathrm{shock}}$ while f_{2} is observed at $\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{T}}_{\mathrm{int}}\approx \mathrm{1}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}\approx \mathrm{4}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}^{\mathrm{shock}}$, where ${\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}^{\mathrm{shock}}$ is the local gyroperiod estimated within the shock front (at the middle of the ramp). This indicates that the reflection process is not uniform in time but leads to the formation of ion bursts associated to the shock dynamics, even if the number of ions which spend several gyroperiods ${\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}^{\mathrm{shock}}$ (i.e., ≈4 bounces) within the shock front is rapidly negligible. In addition, for N_{Box}=0–2, f_{1} and f_{2} have a similar amplitude, which is not the case for N_{Box}=4–9. In fact, a close look at f_{1} and f_{2} shows that f_{2} does not decrease in magnitude but rather the amplitude of f_{1} drastically increases from 10 (N_{Box}=0) to 2500 (N_{Box}=9). Then, the f_{2} population is always present but becomes negligible for lower θ_{Bn} as compared with f_{1}; this explains why we do not observe two distinct spots for N_{Box}=4–9.
One helpful aspect of the test particle approach is including or excluding some electromagnetic field components in order to analyze their impact on the particles dynamics. Indeed, it is clear that some electric field component (i.e., E_{l}×B drift) and strong magnetic gradients drift (i.e., $\propto \mathbf{\nabla}B\times \mathit{B}$ drift) can be a prerequisite for a large drift along the front (i.e., ${\mathrm{\Delta}}_{\mathrm{int}}{\mathit{\theta}}_{\mathrm{Bn}}\approx \mathrm{15}{}^{\circ}$), whereas it could be unnecessary for the other case Δ_{int}θ_{Bn}≈4–5^{∘}. Unfortunately, the shock front magnetic gradient cannot be canceled without the shock itself; then, we will focus our study by including or excluding the electric components which will shed new light on the origin of backstreaming ions filling the foreshock.
3.2 Impact of electric field components
Savoini and Lembège (2015) have analyzed the impact of E×B drift velocity on the dynamics of backstreaming ions (Gurgiolo et al., 1983) and, more particularly, as a source of FAB and/or GPB populations. This study has shown that the origin of both populations can be easily explained in terms of E×B drift associated or not to a diffusion in the velocity space, but it was not able to explain the details of the reflection mechanism itself. Then, herein, we will focus on the role of the electrostatic field component ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ built up within the shock front (i.e., space charge effects). This longitudinal component, defined along the normal to the shock front, can be associated to the electrostatic potential wall responsible for some reflected ions. In the case of a constant shock profile in time with a planar geometry, this reflection does conserve the energy since the potential is the same before and after the reflection, and the total work of the electric force is canceled. Nevertheless, in more realistic conditions, this scenario is not valid anymore for ions which drift along the shock front and suffer both time and space electric and magnetic field variations. Then, in the following sections, we will preferentially use the field E_{l} rather than the potential Φ.
Figure 6a shows the percentage BI_{%} of backstreaming ions versus the box number where the black and red curves are defined for ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}\ne \mathrm{0}$ and ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}=\mathrm{0}$ respectively. The impact of the ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ field (i.e., the potential wall) on the reflection process is clearly apparent for all θ_{Bn} values (namely for each N_{Box} number). In particular, the percentage BI_{%} strongly decreases as ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}=\mathrm{0}$, which illustrates the dominant role of ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ field whatever the box is. This is especially true for the lower box number N_{Box}=0–3 (i.e., high θ_{Bn} approaches 90^{∘}) where very few backstreaming ions are observed. This point is not surprising if one remembers that this electrostatic field decelerates the incoming ions (i.e., accelerates ions in the present reference frame) and contributes to the reflection process. In other words, for N_{Box}=0 (i.e., largest θ_{Bn}), escaping ions have to be accelerated to higher parallel velocity, as reviewed in Burgess et al. (2012). Let us stress that Fig. 6a exhibits a clear change in the slope of BI_{%} increase at the box N_{Box}=2 centered around ${\mathit{\theta}}_{\mathrm{Bn}}=\mathrm{70}{}^{\circ}$. Herein, we will consider this value as the reference angle identifying the starting location of the ion foreshock edge attached to the shock front. This value is in reasonable agreement with the value (${\mathit{\theta}}_{\mathrm{io},\mathrm{fore}}\approx \mathrm{66}{}^{\circ}$) found approximately in the previous selfconsistent PIC simulations Savoini and Lembège (2015).
Another consequence is illustrated in Fig. 6b, which shows that the edge of the ion foreshock is shifted due to the lack of reflected ions, and it starts around ${\mathit{\theta}}_{\mathrm{Bn}}\approx \mathrm{55}{}^{\circ}$. Clearly, the contribution of the electric field is important for ions which populate the edge of the foreshock and need to escape at high θ_{Bn}.
Another way to observe the strong impact of ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ on the dynamics of reflected ions is illustrated in Fig. 7, which shows the two ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{hit}}$ and ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ distributions in the same format as that of Fig. 4. The number of backstreaming ions decreases drastically for all boxes and is zero for N_{Box}=0. Furthermore, the density is not uniform for all boxes and appears to be much more important for N_{Box}=5–9 than for N_{Box}=1–4. The ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ distribution is strongly modified and a comparison between Figs. 4 and 7 can be summarized as follows:

The boxes N_{Box}=1–2 evidence a total absence of reflected ions with a small range Δ_{int}θ_{Bn}, and only the ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ distribution around 60^{∘} persists. This result shows that ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ field plays a key role, more specifically on backstreaming ions suffering a one bounce reflection near the edge of the ion foreshock (i.e., Δ_{int}θ_{Bn}≈4–5^{∘}).

The ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ distribution in both cases (${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}\ne \mathrm{0}$ in Fig. 4 and ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}=\mathrm{0}$ in Fig. 7) is roughly similar in corresponding boxes for all N_{Box}≥6. This can be interpreted as either that the ions have been accelerated enough during their reflection at the shock front or that they need a lower parallel velocity to escape upstream. As a consequence, the E_{l} component is not mandatory anymore, and the mirror magnetic mechanism at the shock front can be invoked as the only reflection process.
In summary, the comparison between Figs. 4 and 7 evidences that ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ components are essential in the ion reflection for a high θ_{Bn} angle ($>\mathrm{56}{}^{\circ}$, i.e., N_{Box}=6–9) where ions need a strong acceleration process but play a less important role at lower θ_{Bn} angles. Conversely, for N_{Box}=6–9, the reflection process takes place with a very small Δ_{int}θ_{Bn}, with or without the electric field. In other words, the large shock drift invoked for ${\mathit{\theta}}_{\mathrm{Bn}}\le \mathrm{56}{}^{\circ}$ seems to be mainly supported by the convective electric field ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{t}}$ components present at the shock front.
Similarly, the one bounce reflection always occurs, even in the absence of ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ field (i.e., in the absence of the shock front potential wall), and can then be associated to a magnetic reflection which seems to be very efficient, especially at lower θ_{Bn}. This one bounce reflection (i.e., f_{1}) corresponds essentially to a short interaction time, as illustrated in Fig. 5 ($\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{T}}_{\mathrm{int}}\approx \mathrm{1}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}^{\mathrm{shock}}$). Then, the ion energy gain is essentially due to a Fermitype acceleration.
3.3 Impact of the shock front nonstationarity
Previous studies have largely evidenced that a quasiperpendicular shock front can be intrinsically nonstationary due to different mechanisms (for a review, see Lembège et al., 2004; Marcowith et al., 2016). Then, it is important to analyze the impact of such nonstationarity on the temporal ion foreshock dynamics. As a first step, we plot in Figure 8 the time evolution of the backstreaming ions' percentage BI_{%} as these leave the front and escape into the upstream region, where BI_{%} is the instantaneous rate computed during a short time range $\mathrm{\Delta}\stackrel{\mathrm{\u0303}}{T}={\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}/\mathrm{20}$. The time ${\stackrel{\mathrm{\u0303}}{T}}_{\mathrm{init}}=\mathrm{1248}$ is the initial time when test particles are launched into the timedependent simulation.
 i.
Results of Fig. 8 are obtained as the ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ field components are included (black curve) and artificially excluded (red curve). One concludes that the percentage BI_{%} strongly decreases as ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ components are excluded, and that the impact of ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ field is emphasized for lower N_{Box}. In other words, the backstreaming ions mainly appear for higher N_{Box}>5 (i.e., for lower θ_{Bn}), even in absence of ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ field.
 ii.
The different results may be classified into two groups, where (i) the first one concerns boxes N_{Box}=0–4 showing a slow increase (almost monotonic) in the reflection rate and (ii) a second group N_{Box}=5–9 evidences a steep increase followed by a flattop shape around BI_{%}≈1, even if it increases slightly with N_{Box}. At the end of the simulation, the strong decrease in BI_{%} observed for all boxes corresponds to the time when all ions of the different boxes have been swapped by the propagating shock front and then no more ions are back streaming.
For the first group N_{Box}=0–4, a delay is observed in the formation of backstreaming ions between the different boxes, although test particles are initially evenly distributed in the whole boxes. For N_{Box}=0, backstreaming ions appear around $\stackrel{\mathrm{\u0303}}{T}\approx \mathrm{5248}$, i.e., well after the initial release time ${\stackrel{\mathrm{\u0303}}{T}}_{\mathrm{init}}=\mathrm{1248}$ (i.e., $\mathrm{\Delta}\stackrel{\mathrm{\u0303}}{T}=\mathrm{5248}\mathrm{1248}$, which covers $\approx \mathrm{2.6}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$), whereas this time delay decreases to $\stackrel{\mathrm{\u0303}}{T}\approx \mathrm{2018}$ (i.e., $\mathrm{\Delta}\stackrel{\mathrm{\u0303}}{T}=\mathrm{2018}\mathrm{1248}=\mathrm{770}$, which covers $\approx \mathrm{0.5}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$) as N_{Box} increases. This illustrates the larger time delay of ions having interacted with the front to escape upstream at high θ_{Bn}. For N_{Box}=0–2, ions have to stay longer within the shock front to finally escape at lower ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ (Fig. 7), which is illustrated by the increase of BI_{%} as the time evolves (Fig. 8). The second group (N_{Box}=5–9) concerns boxes which are already at a lower angle θ_{Bn} with easier escaping conditions. In this case, ions are reflected continuously with some time variation in BI_{%}.
 iii.
Another interesting feature is the presence of different modulations which are superimposed on a timeaveraged reflection ion rate (blue dashed line), especially for the boxes N_{Box}=1–5 and the boxes N_{Box}=7–8. Nearly all boxes exhibit these modulations which represent about 40 % of the averaged BI_{%}, although the amplitude of these modulations varies versus time and N_{Box}. These evidence a nonstationary ion escaping rate. These modulations almost disappear in the case of ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}=\mathrm{0}$ (red line). This illustrates that the presence of the ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ field component is a key ingredient in the formation of these modulations since this electric field component is also involved in the shock front selfreformation, as described in previous works (Lembège and Savoini, 1992; Scholer et al., 2003; Matsukiyo and Scholer, 2006).
This result confirms the importance of the electrostatic field component at the shock front in the reflection processes of the backstreaming ions and, most importantly, on the ion foreshock nonstationarity behavior as described in a previous paper (Savoini and Lembège, 2015). Nevertheless, it is quite difficult to establish a onetoone correspondence between these modulations and the nonstationarity of the shock front because, during the sampling time interval $\mathrm{\Delta}\stackrel{\mathrm{\u0303}}{T}={\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}/\mathrm{20}$, the front nonstationarity and the timeofflight effects have mixed ions coming from either different times and/or different ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ regions (even if they are in the same box). So this selfconsistent approach is not totally adapted to resolve this question. A complementary approach is necessary, based on simplified simulations with fixed shock front profiles in expansion (nonstationary effects are excluded). This motivates the homothetic expansion model (HE) described in the next section.
4.1 Descriptions of the HE model
Let remember that all simulations are made in the solar wind reference frame (i.e., the curved shock expands into the upstream region where the solar wind is at rest). As a consequence, if one follows test particles within this configuration, we have to mimic this behavior. In order to proceed, we apply a homothetic transformation (homogeneous dilatation in all directions) with an expansion factor deduced from the shock front velocity determined at selected times, as illustrated in panels 3a–b of Fig. 1. Special attention has been taken in the determination of this homothetic factor $\mathit{\lambda}={v}_{\mathrm{shock}}\times t$, since the shock front velocity v_{shock} at a given time must fit with the corresponding value issued from the PIC simulations (Savoini and Lembège, 2015). With this information, we are able to expand the shock front through a cubic interpolation as it propagates with v_{shock} in an expanding simulation plane (i.e., the grid cells stay constant ${\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}}}_{x}={\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}}}_{y}=cte$ but the number of the grid cells increases accordingly). In other words, all points of electromagnetic fields at the shock front follow the relation $\mathit{O}\mathit{M}\u27fc({v}_{\mathrm{shock}}\times t)\mathit{O}\mathit{M}$, where v_{shock} is the value of the shock velocity, as computed from our selfconsistent 2D PIC simulations at the selected time, and OM is the vector between the initial location of the shock front (i.e., the point O) and any point of the field array (i.e., the point M). At this stage, we have to point out that the velocity v_{shock} remains artificially constant during the whole simulation, which is not the case for the FCE model where v_{shock} slightly decreases. Then, the same procedure is repeated for another selected time, so that one can analyze the impact of different shock front inhomogeneities and curvature on ion dynamics; let us note that timeofflight effects are always included. Each front profile is analyzed within a similar simulation time range $\approx \mathrm{3}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$. In summary, similar simulations are performed for 174 different times in order to simulate all the different shock profiles provided by the 2D PIC selfconsistent simulation from ${\stackrel{\mathrm{\u0303}}{t}}_{\mathrm{init}}=\mathrm{1.2}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$ to ${\stackrel{\mathrm{\u0303}}{t}}_{\mathrm{simul}}=\mathrm{5.4}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$.
4.2 General features of the backstreaming ions
Figure 9 has been achieved in the same format as Fig. 8 by performing 100 independent simulations (i.e., we take only the first 100 simulations so that all test particles hit the propagating shock front). During this range, the shock front is forced to expand (see Sect. 4.1). For example, the chosen time $\stackrel{\mathrm{\u0303}}{T}=\mathrm{4456}\approx \mathrm{4.2}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$ on the abscissa axis corresponds to one particular shock front profile (i.e., including all electromagnetic field components issued from the previous selfconsistent PIC simulations) from which we have measured the instantaneous shock front velocity and that we follow during the time range covering $\approx \mathrm{3}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$.
The purpose is to determine (i) whether some shock profiles are more appropriate than others for the formation of backstreaming ions and (ii) if, yes, whether better reflection takes place for some particular angular range of θ_{Bn}.
The comparison of Figs. 8 and 9 provides the following information. The maximum BI_{%} value is much higher for the HE model than for the previous FCE model for each corresponding box. For example, for N_{Box}=9, BI${}_{\mathit{\%}}^{\mathrm{max}}\approx \mathrm{1.2}$ in the FCE simulations as compared with BI${}_{\mathit{\%}}^{\mathrm{max}}\approx \mathrm{15}$ in the HE configuration. In fact, one has to remember that for the FCE model the BI${}_{\mathit{\%}}^{\mathrm{max}}$ value represents an instantaneous reflection rate versus the shock front evolution, whereas this rate is a timeintegrated value for the HE model. Indeed, in this model, the shock front profile stays the same during the whole simulation, and then, if this profile allows the reflection of some incident ions, they will be reflected continuously, leading to a high BI_{%} number.
Obviously, the main information is not the BI_{%} value itself but rather its evolution versus time and for different shock profiles. Other main results issued from Fig. 9 may be summarized as follows:

The N_{Box}=0 box evidences almost no reflection for the majority of the shock front profiles, which indicates that nonstationarity effects present in the FCE configuration (i.e., Fig. 3) are needed for feeding backstreaming ions along the edge of the foreshock.

Boxes N_{Box}=1–8 clearly show some strong modulations in the percentage BI_{%} versus the shock profile of concern, which corresponds to a quasiperiodic bursty emission of backstreaming ions. The BI_{%} rate periodically reaches a maximum value followed by a minimum around zero. The corresponding time period $\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{T}}_{\mathrm{max}}$ is about $\approx \mathrm{460}=\mathrm{0.5}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$ (between two successive maxima). The temporal width of each maximum is about $\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{T}}_{\mathrm{range}}\approx \mathrm{256}=\mathrm{0.25}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}$. These modulations mean that conditions for the formation of backstreaming particles are not continuous but correspond to some specific shock front profiles. In addition, these modulations appear synchronized in time for the different boxes 1–7, which implies that the local reflection conditions are not strongly dependent on the θ_{Bn} angle but rather depend on the shock profile at certain times.

In contrast, the boxes N_{Box}=8–9 also evidence the same kind of modulations but with greatly reduced amplitude; these are even nonexistent between $\stackrel{\mathrm{\u0303}}{T}\approx \mathrm{3456}$ and 4456, which indicates a low sensitivity to the shock front profile when approaching ${\mathit{\theta}}_{\mathrm{Bn}}=\mathrm{45}{}^{\circ}$.

Similarly, the maximum values BI${}_{\mathit{\%}}^{\mathrm{max}}$ (black curve) change drastically with box numbers, namely from small amplitudes, BI_{%}≤6 for N_{Box}=1–2, to very high values, BI_{%}≈30 for N_{Box}=3–5, before decreasing again for N_{Box}=6–9. These variations may be understood by taking into account the reflection processes present at these different θ_{Bn} and, more specifically, in regards to the electric field components. For boxes N_{Box}=0–2, reflection is almost impossible without the ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ component (i.e., electric potential wall). As θ_{Bn} decreases (for N_{Box}=3–5), the reflection becomes easier and both magnetic and electric fields contribute to the percentage of reflected ions. Finally, for the last boxes N_{Box}=6–9, the peak amplitude decreases but the contribution of ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ becomes less important in the reflection process, as evidenced by comparing both black (${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}\ne \mathrm{0}$) and red (${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}=\mathrm{0}$) curves. Instead, another process, essentially driven by the magnetic field (mirror reflection) contributes more since the peak amplitude of the red curve increases progressively as N_{Box} increases from 6 to 9.

Only the last peak around $\stackrel{\mathrm{\u0303}}{T}\approx \mathrm{5456}$ has a different behavior in comparison with others. In particular, we observe that, for these different times or profiles, the presence of the ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ leads to higher BI_{%}. This behavior can be understood because v_{shock} is lower for these times, and the ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ component is necessary for decelerating ions and reflecting them. Without this component, only the magnetic reflection is present, and BI_{%} has the same amplitude as the previous maximum.

Finally, Fig. 9 confirms the key role of E_{l} field in backstreaming ion formation, except when approaching ${\mathit{\theta}}_{\mathrm{Bn}}=\mathrm{45}{}^{\circ}$ (N_{Box}=8–9) while another reflection process is also at work (in absence of E_{l}). This represents an indirect way of stressing the noticeable impact of the magnetic field in this angular range. This magnetic reflection process is more evidenced at lower θ_{Bn} since ions need less parallel velocity to be reflected back into the upstream region. This statement can be quantified more precisely, as explained in Sect. 5.
This result is an illustration of the impact of the electrostatic component at the shock front. As is well known, this component works to decelerate incoming ions (i.e., accelerate in our solar wind frame) and to accelerate electrons (i.e., decelerate in our solar wind frame) to the downstream region (Savoini and Lembège, 1994; Bale et al., 2005). As a consequence, this electrostatic component is an essential ingredient in the formation of backstreaming ions, especially at higher θ_{Bn}.
A previous paper (Savoini and Lembège, 2015) demonstrated that all reflected ions suffered the same E×B drift in the velocity space, which can account for the pitch angle distributions observed at the shock front. In fact, the key point is the time spent by particles within the front shock, which finally decides whether ions will escape to form the FAB (with a pitch angle $\mathit{\alpha}\approx \mathrm{0}{}^{\circ}$) or GPB populations (with a pitch angle $\mathit{\alpha}\ne \mathrm{0}{}^{\circ}$), where α is the angle between the velocity vector and the magnetic field. In other words, the FAB population may be associated to a large drift along the shock front (and/or long interaction time) during which particles see a timevarying shock front and lose their phase coherency; this case corresponds to a large angular range Δ_{int}θ_{Bn} mentioned in Sect. 3.1. In contrast, the ions of the GPB population have a shorter interaction time with the shock front associated to a small angular range Δ_{int}θ_{Bn}. Present test particle simulations allow one to have a deeper insight to the spatial origin of the observed FAB and GPB populations. Then, we have to analyze the ion velocity distribution more carefully.
Figure 10 plots the local perpendicular velocity distribution functions $f({\mathit{v}}_{\u27c2\mathrm{1}},{\mathit{v}}_{\u27c2\mathrm{2}})$ in both FCE and HE approaches (where v_{⟂1} and v_{⟂2} refer the ion perpendicular velocity components defined with respect to the local magnetic field). All plots are obtained at the end of the simulations and take into account the different populations observed in Figs. 4 and 7 (i.e., both distinct peaks of ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ angle for lower N_{Box} are included); results issued from FCE (Fig. 10, left panel) and HE (Fig. 10, right panel) configurations are considered. For the HE configuration, we chose an initial time $\stackrel{\mathrm{\u0303}}{T}=\mathrm{4848}$ (see Fig. 9 for reference) which corresponds to a maximum of BI_{%} in order to have enough reflected ions in the velocity space. Results from the three different boxes, ${N}_{\mathrm{Box}}=\mathrm{1},\phantom{\rule{0.25em}{0ex}}\mathrm{4}$ and 8, are represented in order to give an overview of the whole ion foreshock components.
Results of the FCE configuration (left panels) can be analyzed in Figs. 4, 7 and 10. Plots of the E_{l}≠0 case (Fig. 10) show that the N_{Box}=1 has a low number of upstream reflected ions, which leads to poor statistics and a noisy $f({\mathit{v}}_{\u27c2\mathrm{1}},{\mathit{v}}_{\u27c2\mathrm{2}})$ distribution. Nevertheless, it evidences approximately a distribution with a maximum slightly noncentered at ${\mathit{v}}_{\u27c2}=\mathrm{0}$. Then, this distribution can be viewed as a mixing of GPB and FAB populations, even if the GPB population with a pitch angle different zero is the largest one. When moving further into the foreshock region (i.e., lower θ_{Bn} angle with N_{Box}=4), the number of reflected ions drastically increases and we observe more clearly the characteristic partial ring of the GPB population (as in Savoini and Lembège, 2015). In addition, the center of the ring is also partially filled in because of partial diffusion due to particles having a large Δ_{int}θ_{Bn} range, corresponding to the ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ peak around 60^{∘} in Fig. 4, and/or by the intrinsic time fluctuations of the front which tends to blur out the velocity distribution both in perpendicular and parallel directions. At last, in agreement with the associated small ${\mathit{\theta}}_{\mathrm{Bn}}^{\mathrm{exit}}$ of Fig. 4, N_{Box}=8 (in Fig. 10) also evidences a nonMaxwellianlike distribution ($\mathit{\alpha}\ne \mathrm{0}{}^{\circ}$). When we look at the E_{l}=0 case, we observe roughly the same behavior for all boxes, even if the decrease in the backstreaming ions number in box N_{Box}=8 makes the comparison difficult.
A further analysis requires a similar approach with the HE configuration in which we follow a succession of independent expanding shock profiles in order to exclude the impact of the time fluctuations on the velocity distribution $f({\mathit{v}}_{\parallel},{\mathit{v}}_{\u27c2})$; then, no ion diffusion associated to these fluctuations is allowed. Results of the HE configuration (right panel, Fig. 10) show reflected ions for N_{Box}=1 when ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}\ne \mathrm{0}$, but once again, no reflected ions can be seen when ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}=\mathrm{0}$. This evidences the importance of the electrostatic potential wall in order to reflect upstream ions for high θ_{Bn}. On the other hand, if the amplitude of the perpendicular velocity is roughly the same for the two different configurations (FCE and HE), the HE case shows a very wellformed ring, in contrast with the FCE case, which exhibits a diffuse velocity space. This illustrates that field time variations (FCE configuration) are much more efficient at diffusing particles than the fields spatial variations (HE configuration). Similarly, N_{Box}=4 exhibits a clear ring, which is a feature of the GPB population. The center of the ring is not partially filled in since time velocity diffusion is excluded. These results demonstrate that the formation of the FABlike population is also mainly due to ion velocity diffusion related to the time fluctuations of the shock front which can have different origins, as described in previous works (Kucharek et al., 2004; Bale et al., 2005). Finally, for N_{Box}=8, the velocity distribution drastically changes from a ring (${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}\ne \mathrm{0}$) to a localized bump (${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}=\mathrm{0}$) roughly similar to the FCE case. It is clear that the number of reflected ions decreases drastically as ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ field components are artificially suppressed. But, more important is that the formation of a nongyrotropic distribution does not depend strongly on the ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ component and is mainly controlled by the convective electric field through the E_{t}×B drift in the velocity space, as already described in Savoini and Lembège (2015).
Let us remember that each distribution results from a combination of all particles originating from one given box (no matter where they end up spatially). This differs from the more common strategy based on measurements of local ion distributions, as performed in Savoini and Lembège (2015), but which did not establish, at that time, which part of the curved shock front the FAB and GPB ions are issued from. A further analysis is needed and is left for later work.
Present test particle simulations reinforce the scenario described in Savoini and Lembège (2015) and have allowed us to investigate the formation of the ion foreshock more deeply. In summary, the impact of the three quantities or effects has been identified as follows: (i) the electric field (separately the ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ and the ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{t}}$ components), (ii) the magnetic field and (iii) the shock front nonstationarity. The synoptic of Fig. 11 summarizes the importance of each impact versus the shock front curvature from the edge of the ion foreshock (${\mathit{\theta}}_{\mathrm{Bn}}\approx \mathrm{70}{}^{\circ}$) to ${\mathit{\theta}}_{\mathrm{Bn}}=\mathrm{45}{}^{\circ}$. In this sketch, the colors vary from strong (full color) to weak (white) intensities as a function of their respective influence on the ion dynamics. These different effects are the following:

Impact of the E_{l} field on the ion reflection process. As is well known, the builtup potential wall at the shock front (i.e., the electric field E_{l} along the shock normal n) is mainly responsible for the deceleration (i.e., acceleration in our reference frame) of the incoming upstream ions by the shock front. The E_{l} component has essentially two distinct impacts where, (i) without this electric component, no reflected ions are observed for ${\mathit{\theta}}_{\mathrm{Bn}}>\mathrm{62}{}^{\circ}$, whereas in presence of this electric component reflected ions (having suffered even one bounce) can be observed at the edge of the ion foreshock ($\approx {\mathit{\theta}}_{\mathrm{Bn}}\le \mathrm{70}{}^{\circ}$); (ii) at lower angles (${\mathit{\theta}}_{\mathrm{Bn}}\le \mathrm{50}{}^{\circ}$) many ions are reflected without the help of the ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ component and can be associated to a magnetic mirror reflection. Then, in Fig. 11, ${\stackrel{\mathrm{\u0303}}{E}}_{\mathrm{l}}$ is only reported as being strong around the edge of the ion foreshock to emphasize its mandatory action for high θ_{Bn} angles.

Impact of the E_{t} field on the ion reflection process. Figure 10 evidences that the convective electric component E_{t} is always present in our simulation (we are in the solar wind reference frame and then E_{t}≠0 within the curved propagating shock front). Our previous work (Savoini and Lembège, 2015) was only able to show that the E×B drift scenario in the velocity space foreseen by Gurgiolo et al. (1983) was at the origin of two distinct GPB (i.e., one bounce) and FAB (i.e., multibounces) populations only separated from the particle time history within the shock front.
But this scenario was not able to distinguish the relative importance between the two electric field components of E_{l} and E_{t}, respectively. This question has been clarified in the present paper since FAB and GPB populations formed by the E×B drift in the velocity space are evidenced with or without the E_{l} component (Fig. 10). Then, the E_{l} field seems to be less dominant than the E_{t} field. Finally, Fig. 11 illustrates the E×B drift impact by a dark color that is almost uniform within the whole quasiperpendicular region.

Impact of the B field on the ion reflection process. The magnetic field is important for several reasons: (i) its increase at the shock front builds up the E_{l} component (space charge effect), which reflects back incoming lowenergy ions, and (ii) more importantly, it is also directly responsible for the reflection of some ions (i.e., through the magnetic mirror reflection) and for the drift along the shock front of the multibounce ions (i.e., FAB population). Then, this population gains energy as ions propagate in the E_{t} direction along the shock front. Nevertheless, as θ_{Bn} decreases from 90 to 45^{∘}, the ion reflection becomes easier since the parallel guiding center velocity needed to overcome the shock front velocity decreases (Paschmann et al., 1980). This behavior is clearly illustrated herein by the increase in the percentage of reflected ions BI_{%} as θ_{Bn} decreases (i.e., N_{Box} increases). This behavior persists even in absence of E_{l}, where BI_{%} is only reduced by a factor of 2.5, as illustrated by Fig. 6. Then, if the magnetic field is important in the whole quasiperpendicular region, we emphasize, in Fig. 11, its stronger impact near ${\mathit{\theta}}_{\mathrm{Bn}}\approx \mathrm{45}{}^{\circ}$ where it is mainly responsible for the reflection (i.e., magnetic mirror) and acceleration (i.e., Fermitype) of ions (Webb et al., 1983).

Impact of the shock front nonstationarity. Present simulations show that the reflection process is not continuous in both time and space but strongly depends on the local shock front profile met by incoming ions at their hitting time. This behavior is difficult to identify in experimental measurements since the particles coming from different shock locations and at different times are mixed; in contrast, this can be easily evidenced in our HE test particles configuration. This configuration evidences that particular shock profiles are more suitable for the formation of backstreaming ions than other ones. Indeed, we observe modulations of the BI_{%} percentages in Fig. 9 which are much more pronounced than in our FCE configuration (Fig. 8). These modulations are so strong that BI_{%} drops to zero periodically, which means that for certain shock front profiles no ion can escape into the upstream region. This behavior is observed for all N_{Box} at the same time (i.e., same shock profile), which implies that the ion reflection does not depend on the location along the shock front but essentially on the global profile of the shock at a given time. In the present simulations, we can identify four distinct and noticeable bursts (i.e., maxima BI_{%} values) with an average cyclic occurrence period of $\mathrm{1}{\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}^{\mathrm{shock}}$ (where ${\stackrel{\mathrm{\u0303}}{\mathit{\tau}}}_{\mathrm{ci}}^{\mathrm{shock}}$ is defined at the shock ramp). Surprisingly, N_{Box}=0 (Fig. 9) does not show the same bursts as the others, which suggests that the time variations in the shock front (and associated particle diffusion, as suggested by Kucharek et al., 2004) are mandatory for obtaining backstreaming ions around the edge of the ion foreshock (i.e., high θ_{Bn}). This point will require a further investigation.
In summary, present results show that the formation of the ion foreshock is not a continuous process but must be considered as time dependent, which leads to bursty emissions of backstreaming ions. Three different contributions have been evidenced, namely (i) the E_{l} component in the ion global reflection process, in particular for high θ_{Bn}, (ii) the magnetic field B essentially observed when E_{l}=0 for lower θ_{Bn}, such as the magnetic mirror reflection, and (iii) finally, the E_{t}×B drift in the velocity space mainly sustained by the convective electric field which is necessary to generate both FAB and GPB populations as described in Savoini and Lembège (2015).
Unfortunately, the impact of the shock front nonstationarity on the ion foreshock is difficult to analyze (see, for example, Fig. 7) for two different reasons, namely (i) the timeofflight effects mix reflected ions coming from different shock profiles, and (ii) even if some shock profiles are more efficient than others at reflecting ions, their respective impacts disappear rapidly since they are being blurred out by the impact of less efficient profiles on particles as time evolves.
Test particle data have been deposited in the archive of the LPP laboratory at https://ao.lpp.polytechnique.fr/index.php/s/YtN3Bk4LttjpdD7 (Savoini, 2020). A “Readme” file has been added, which includes the creators, title, and date of last access as well as the description of the HDF5 data structure produced and a Python program that can be used to read the data directly and display the particle trajectories.
PS and BL contributed to the design and implementation of the research, the analysis of the results and the writing of the paper. The data were produced by PS.
The authors declare that they have no conflict of interest.
Numerical simulations were performed at the TGCC highperformance computer center located in BruyèresleChâtel (near Paris), and we thank the staff for their support (DARI project; grant no. A0050400295). ISSI (Bern, Switzerland) is thanked for supporting the collaboration network for “Resolving the Microphysics of Collisionless Shock Waves”. The authors wish to thank Yann PfauKempf and the anonymous referee for their helpful comments.
This research has been supported by Labex PLAS@PAR (grant no. ANR11IDEX000402) and by the Centre National d'Études Spatiales (CNES; grant nos. APRWEEXP/10010105 and APRZETPE0010/010105).
This paper was edited by Minna Palmroth and reviewed by Yann PfauKempf and one anonymous referee.
Bale, S. D., Balikhin, M. A., Horbury, T. S., Krasnoselskikh, V. V., Kucharek, H., Möbius, E., Walker, S. N., Balogh, A., Burgess, D., Lembège, B., Lucek, E. A., Scholer, M., Schwartz, S. J., and Thomsen, M. F.: QuasiPerpendicular Shock Structure and Processes, Space Sci. Rev., 118, 161–203, https://doi.org/10.1007/s1121400538270, 2005. a, b, c
BlancoCano, X., Omidi, N., and Russell, C. T.: Global Hybrid Simulations: Foreshock Waves and Cavitons under Radial Interplanetary Magnetic Field Geometry, J. Geophys. Res., 114, 01216, https://doi.org/10.1029/2008JA013406, 2009. a
Bonifazi, C. and Moreno, G.: Reflected and Diffuse Ions Backstreaming from the Earth's Bow Shock: 2. Origin, J. Geophys. Res., 86, 4405–4413, 1981a. a
Bonifazi, C. and Moreno, G.: Reflected and Diffuse Ions Backstreaming from the Earth's Bow Shock. I Basic Properties, J. Geophys. Res., 86, 4397–4404, 1981b. a
Burgess, D., Möbius, E., and Scholer, M.: Ion Acceleration at the Earth's Bow Shock, Space Sci. Rev., 173, 5–47, https://doi.org/10.1007/s1121401299015, 2012. a
Eastwood, J., Lucek, E., Mazelle, C., and Meziane, K.: The Foreshock, Space Sci. Rev., 118, 41–94, https://doi.org/10.1007/s1121400538243, 2005. a, b
Edmiston, J. P., Kennel, C. F., and Eichler, D.: Escape of Heated Ions Upstream of QuasiParallel Shocks, Geophys. Res. Lett., 9, 531–534, https://doi.org/10.1029/GL009i005p00531, 1982. a
Fuselier, S. A.: Ion Distributions in the Earth's Foreshock Upstream from: The Bow Shock, Adv. Space Res., 15, 43–52, 1995. a
Fuselier, S. A., Thomsen, M., Gary, S., and Bame, S.: The Phase Relationship between GyrophaseBunched Ions and MHDlike Waves, Geophys. Res. Lett., 13, 60–63, 1986a. a
Fuselier, S. A., Thomsen, M. F., Gosling, J. T., Bame, S. J., and Russell, C. T.: Gyrating and Intermediate Ion Distributions Upstream from the Earth's Bow Shock, J. Geophys. Res., 91, 91–99, https://doi.org/10.1029/JA091iA01p00091, 1986b. a
Giacalone, J., Jokipii, J. R., and Kota, J.: Ion Injection and Acceleration at QuasiPerpendicular Shocks, J. Geophys. Res., 99, 19351–19358, 1994. a
Gosling, J. T., Thomsenn, M. F., Bame, S. J., Feldman, W. C., Pashmann, G., and Sckopke, N.: Evidence for Specularly Reflected Ions Upstream from the QuasiParallel Bow Shock, Geophys. Res. Lett., 9, 1333–1336, 1982. a
Gurgiolo, C., Parks, G. K., and Mauk, B. H.: Upstream Gyrophase Bunched Ions: A Mechanism for Creation at the Bow Shook and the Growth of Velocity Space Structure through Gyrophase Mixing, J. Geophys. Res., 88, 9093–9100, https://doi.org/10.1029/JA088iA11p09093, 1983. a, b
Hamza, A. M., Meziane, K., and Mazelle, C.: Oblique Propagation and Nonlinear Wave Particle Processes, J. Geophys. Res., 111, 04104, https://doi.org/10.1029/2005JA011410, 2006. a
Hartinger, M. D., Turner, D. L., Plaschke, F., Angelopoulos, V., and Singer, H.: The Role of Transient Ion Foreshock Phenomena in Driving Pc5 ULF Wave Activity, J. Geophys. Res., 118, 299–312, https://doi.org/10.1029/2012JA018349, 2013. a
Hoshino, M. and Terasawa, T.: Numerical Study of the Upstream Wave Excitation Mechanism. INonlinear Phase Bunching of Beam Ions, J. Geophys. Res., 90, 57–64, 1985. a
Ipavich, F. M., Galvin, A. B., Gloeckler, G., Scholer, M., and Hovestadt, D.: A Statistical Survey of Ions Observed Upstream of the Earth's Bow Shock – Energy Spectra, Composition, and Spatial Variation, (Upstream Wave and Particle Workshop, 86, 4337–4342, https://doi.org/10.1029/JA086iA06p04337, 1981. a
Kempf, Y., Pokhotelov, D., Gutynska, O., Wilson, III, L. B., Walsh, B. M., von Alfthan, S., Hannuksela, O., Sibeck, D. G., and Palmroth, M.: Ion Distributions in the Earth's Foreshock: HybridVlasov Simulation and THEMIS Observations, J. Geophys. Res., 120, 3684–3701, https://doi.org/10.1002/2014JA020519, 2015. a, b
Kucharek, H.: On the Physics of Collisionless Shocks: Cluster Investigations and Simulations, J. Atmos. Sol.Terr. Phys., 70, 316–324, https://doi.org/10.1016/j.jastp.2007.08.052, 2008. a
Kucharek, H., Möbius, E., Scholer, M., Mouikis, C., Kistler, L. M., Horbury, T., Balogh, A., Réme, H., and Bosqued, J. M.: On the origin of fieldaligned beams at the quasiperpendicular bow shock: multispacecraft observations by Cluster, Ann. Geophys., 22, 2301–2308, https://doi.org/10.5194/angeo2223012004, 2004. a, b, c, d
Lembège, B. and Savoini, P.: Nonstationarity of a Twodimensional Quasiperpendicular Supercritical Collisionless Shock by SelfReformation, Phys. Fluids, 4, 3533–3548, 1992. a, b
Lembège, B. and Savoini, P.: Formation of Reflected Electron Bursts by the Nonstationarity and Nonuniformity of a Collisionless Shock Front, J. Geophys. Res., 107, 1037, https://doi.org/10.1029/2001JA900128, 2002. a
Lembège, B., Giacalone, J., Scholer, M., Hada, T., Hoshino, M., Krasnosel'skikh, V., Kucharek, H., Savoini, P., and Terasawa, T.: Selected Problems in CollisionlessShock Physics, Space Sci. Rev., 110, 161–226, https://doi.org/10.1023/B:SPAC.0000023372.12232.b7, 2004. a, b
Marcowith, A., Bret, A., Bykov, A., Dieckman, M. E., Drury, L. O., Lembège, B., Lemoine, M., Morlino, G., Murphy, G., Pelletier, G., Plotnikov, I., Reville, B., Riquelme, M., Sironi, L., and Novo, A. S.: The Microphysics of Collisionless Shock Waves, Rep. Prog. Phys., 79, 046901, https://doi.org/10.1088/00344885/79/4/046901, 2016. a
Matsukiyo, S. and Scholer, M.: On Reformation of QuasiPerpendicular Collisionless Shocks, Adv. Space Res., 38, 57–63, https://doi.org/10.1016/j.asr.2004.08.012, 2006. a
Mazelle, C., Meziane, K., Lequéau, D., Wilber, M., Eastwood, J. P., Rème, H., Sauvaud, J. A., Bosqued, J. M., Dandouras, I., McCarthy, M., Kistler, L. M., Klecker, B., Korth, A., BavassanoCattaneo, M. B., Pallocchia, G., Lundin, R., and Balogh, A.: Production of Gyrating Ions from Nonlinear WaveParticle Interaction Upstream from the Earth's Bow Shock: A Case Study from ClusterCIS, Planet. Space Sci., 51, 785–795, https://doi.org/10.1016/S00320633(03)001077, 2003. a
Mazelle, C., Meziane, K., and Wilber, M.: Fieldaligned and gyrating ion beams in a planetary foreshock, AIP Conf. Proc., 781, 89–94, https://doi.org/10.1063/1.2032680, 2005. a
Meziane, K.: AIP Conference Proceedings, in: THE PHYSICS OF COLLISIONLESS SHOCKS: 4th Annual IGPP International Astrophysics Conference, AIP, 116–122, https://doi.org/10.1063/1.2032683, 2005. a
Möbius, E., Kucharek, H., Mouikis, C., Georgescu, E., Kistler, L. M., Popecki, M. A., Scholer, M., Bosqued, J. M., Rème, H., Carlson, C. W., Klecker, B., Korth, A., Parks, G. K., Sauvaud, J. C., Balsiger, H., BavassanoCattaneo, M.B., Dandouras, I., DiLellis, A. M., Eliasson, L., Formisano, V., Horbury, T., Lennartsson, W., Lundin, R., McCarthy, M., McFadden, J. P., and Paschmann, G.: Observations of the spatial and temporal structure of fieldaligned beam and gyrating ring distributions at the quasiperpendicular bow shock with Cluster CIS, Ann. Geophys., 19, 1411–1420, https://doi.org/10.5194/angeo1914112001, 2001. a
Oka, M., Terasawa, T., Saito, Y., and Mukai, T.: FieldAligned Beam Observations at the QuasiPerpendicular Bow Shock: Generation and Shock Angle Dependence, J. Geophys. Res., 110, 05101, https://doi.org/10.1029/2004JA010688, 2005. a
Otsuka, F., Matsukiyo, S., Kis, A., Nakanishi, K., and Hada, T.: Effect of Upstream ULF Waves on the Energetic Ion Diffusion at the Earth's Foreshock. I. Theory and Simulation, Astrophys. J., 853, 117, https://doi.org/10.3847/15384357/aaa23f, 2018. a
Paschmann, G., Sckopke, N., and Asbridge, J.: Energization of Solar Wind Ions by Reflection from the Earth's Bow Shock, J. Geophys. Res., 85, 4689–4693, 1980. a, b
Paschmann, G., Sckopke, N., Papamastorakis, I., Asbridge, J. R., Bame, S. J., and Gosling, J. T.: Characteristics of Reflected and Diffuse Ions Upstream from the Earth's Bow Shock, J. Geophys. Res., 86, 4355–4364, https://doi.org/10.1029/JA086iA06p04355, 1981. a
Savoini, P.: Data for self consistent test particles, available at: https://ao.lpp.polytechnique.fr/index.php/s/YtN3Bk4LttjpdD7, last access: 13 November 2020. a
Savoini, P. and Lembège, B.: Electron Dynamics in Two and One Dimensional Oblique Supercritical Collisionless Magnetosonic Shocks, J. Geophys. Res., 99, 6609–6635, 1994. a
Savoini, P. and Lembège, B.: TwoDimensional Simulations of a Curved Shock: SelfConsistent Formation of the Electron Foreshock, J. Geophys. Res., 106, 12975, https://doi.org/10.1029/2001JA900007, 2001. a, b
Savoini, P. and Lembège, B.: Production of Nongyrotropic and Gyrotropic Backstreaming Ion Distributions in the QuasiPerpendicular Ion Foreshock Region, J. Geophys. Res., 120, 7154–7171, https://doi.org/10.1002/2015JA021018, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x
Savoini, P., Lembège, B., and Stienlet, J.: On the Origin of the QuasiPerpendicular Ion Foreshock: FullParticle Simulations, J. Geophys. Res., 118, 1132–1145, https://doi.org/10.1002/jgra.50158, 2013. a, b, c, d, e, f
Scholer, M., Shinohara, I., and Matsukiyo, S.: QuasiPerpendicular Shocks: Length Scale of the CrossShock Potential, Shock Reformation, and Implication for Shock Surfing, J. Geophys. Res., 108, 1014, https://doi.org/10.1029/2002JA009515, 2003. a
Schwartz, S. J. and Burgess, D.: On the Theoretical/Observational Comparison of FieldAligned Ion Beams in the Earth's Foreshock, J. Geophys. Res., 89, 2381–2384, https://doi.org/10.1029/JA089iA04p02381, 1984. a
Schwartz, S. J., Thomsen, M. F., and Gosling, J. T.: Ions Upstream of the Earth's Bow Shock – A Theoretical Comparison of Alternative Source Populations, J. Geophys. Res., 88, 2039–2047, https://doi.org/10.1029/JA088iA03p02039, 1983. a, b
Sonnerup, B. U. Ö.: Acceleration of Particles Reflected at a Shock Front, J. Geophys. Res., 74, 1301–1304, https://doi.org/10.1029/JA074i005p01301, 1969. a
Tanaka, M., Goodrich, C. C., Winske, D., and Papadopoulos, K.: A Source of the Backstreaming Ion Beams in the Foreshock Region, J. Geophys. Res., 88, 3046–3054, https://doi.org/10.1029/JA088iA04p03046, 1983. a
Thomsen, M., Gosling, J., and Bame, S.: Gyrating Ions and LargeAmplitude Monochromatic MHD Waves Upstream of the Earth's Bow Shock, Jo. Geophys. Res., 90, 267–273, 1985. a
Thomsen, M. F., Gosling, J. T., and Schwartz, S. J.: Observational Evidence on the Origin of Ions Upstream of the Earth's Bow Shock, J. Geophys. Res., 88, 7843–7852, https://doi.org/10.1029/JA088iA10p07843, 1983. a
Tsurutani, B. T. and Rodriguez, P.: Upstream Waves and Particles – An Overview of ISEE Results, J. Geophys. Res., 86, 4317–4324, 1981. a
Turner, D., Angelopoulos, V., Wilson, L., Hietala, H., Omidi, N., and Masters, A.: Particle Acceleration during Interactions between Transient Ion Foreshock Phenomena and Earth's Bow Shock, Geophys. Res. Abstr., EGU20142276, EGU General Assembly 2014, Vienna, Austria, 2014. a
Webb, G., Axford, W., and Terasawa, T.: On the Drift Mechanism for Energetic Charged Particles at Shocks, Astrophys. J., 270, 537–553, 1983. a
 Abstract
 Introduction
 Numerical simulation conditions
 Numerical results: the fully consistent expansion model (FCE)
 Numerical results: the homothetic expansion (HE) model
 Discussions
 Conclusions
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Numerical simulation conditions
 Numerical results: the fully consistent expansion model (FCE)
 Numerical results: the homothetic expansion (HE) model
 Discussions
 Conclusions
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References