A deep insight into the ion foreshock with the help of test particle two-dimensional simulations

. Two-dimensional (2D) test particle simulations based on shock proﬁles issued from 2D full particle-in-cell (PIC) simulations are used in order to analyze the formation processes of ions back streaming within the upstream region after interacting with a quasi-perpendicular curved shock front. Two different types of simulations have been performed based on (i) a fully consistent expansion (FCE) model, which includes all self-consistent shock proﬁles at different times, and (ii) a homothetic expansion (HE) model in which shock proﬁles are ﬁxed at certain times and artiﬁcially expanded in space. The comparison of both conﬁgurations allows one to analyze the impact of the front nonstationarity on the back-streaming population. Moreover, the role of the space charge electric ﬁeld E l is analyzed by either including or canceling the E l component in the simulations. A detailed comparison of these last two different conﬁgurations allows one to show that this E l component plays a key role in the ion reﬂection process within the whole quasiperpendicular propagation range. Simulations provide evidence that the different ﬁeld-aligned beam (FAB) and gyrophase bunched (GPB)

with or without the conservation of the magnetic moment and (ii) scenarios which invoke the leakage of some magnetosheath ions producing a low-energy 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) 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 low-frequency 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 large-scale, two-dimensional particle-in-cell (PIC) simulation of a curved shock has been used, where full curvature and time-of-flight effects for both electrons and ions are self-consistently 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  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 back-streaming ion reflection process. The main questions presently addressed are summarized as follows: 1. 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)?
2. What is the impact of the space charge electric field localized within the shock ramp on the reflection process?
3. 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  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.

Numerical simulation conditions
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 Savoini, 1992, 2002 for planar shocks).

Self-consistent full PIC simulations
The code solves Maxwell and Poisson's equations in the Fourier space (so-called 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 ∇E l = ρ/ o ), and (ii) a transverse or induced component, hereafter denoted by subscript t (coming from the temporal variations in the magnetic field ∇×E t = −∂B/∂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 E t = −U × 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. Panel (1) plots the simulation plane geometry; the B o magnetostatic field is mainly outside and directed downward from the plane. Panels (2a, b) illustrate the evolution of the magnetic field B tz in the fully consistent expansion (FCE) model (time dependent), respectively, at t init = 2.4 τ ci and t simul = 5.4 τ ci . In panel (2a), used as a reference, the vector velocity U = v shock has been superimposed to illustrate the shock front propagation (white thick arrows); the arrow length is not at the right scale. In addition, the projection of the B o magnetic field lines has been reported (oblique white thin lines). Panels (3a, b) illustrate one example of the curved magnetic field B tz in the homothetic expansion (HE) model (time independent), where the shock profile is fixed in time but expands in space via an expanding factor proportional to the shock front velocity v shock × t; this shock profile has been chosen at the time t init = 2.4 τ ci of the self-consistent simulation. From this time, the shock front dilates by a factor of 2.6 compared to its initial shape.
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 quasi-perpendicular shock (i.e., for 45 • ≤ θ Bn ≤ 90 • ). 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 ρ ci (all normalized quantities are indicated with a tilde˜); this is the same normalization which is used in the previous self-consistent 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 θ Bn = 90 • used as a reference angle; (ii) the time-of-flight effects are self-consistently 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, back-streaming 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 c = 3 and temperature ratio between ion and elec-tron population T el /T io = 1.58. A mass ratio m i /m e = 84 is used in order to save CPU time, and the Alfvén velocity is v A = 0.16. The number of cells in the simulation plane along each axis is NCX = NCY = 8192 ≈ 150 ρ ci , with the size of a grid cell x = y ≈ 1 ρ ce . The shock is in the supercritical regime. with a time-averaged Alfvén Mach number M A ≈ 4 measured at θ Bn = 90 • . In order to observe the early stage of the ion foreshock formation, the end time of the simulation is t simul = 5.4 τ ci (where τ ci is the upstream ion gyro-period), which is large enough to investigate the interaction of incoming ions with the shock front and the further formation of back-streaming ions.

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 T = τ ci /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 self-consistent 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., all-time 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 back-streaming 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 self-consistent simulation .
Let us point out that the use of finite-sized 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 Figure 2. Initial location of the 10 sampling boxes (labeled from zero to nine) which map the upstream ion foreshock region. All ions belonging to a given box are represented by the same color for statistical analysis only (Sects. 3 and 4). The θ Bn propagation angle, where each box is initially centered at time t = 0, is reported above the corresponding identification number of the box, but these colors will not be used anymore in this paper. 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 ≈ 4 • , 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 ≈ 3 τ 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 self-consistent shock front profile. Section 4 will then introduce the more unusual HE expansion model (i.e., homothetic simulation approach).

Numerical results: the fully consistent expansion model (FCE)
To clarify the presentation, we will split this section into two parts in which we analyze, respectively, (i) the dynamics of the back-streaming 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. Figure 3 plots the spatial distribution of the back-streaming 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:

General features of the back-streaming ions
Ann. Geophys., 38, 1217-1235, 2020 https://doi.org/10.5194/angeo-38-1217-2020 Figure 3. FCE configuration: spatial distribution of the back-streaming particle density within the simulation X-Y plane at the end of the simulation time ( t end ≈ 5.4 τ ci , where τ ci is the upstream ion cyclotron period). All boxes are plotted from N Box = 0 to 9, defined initially at θ Bn ≈ 90 • and ≈ 45 • , respectively; the bottom panel labeled N Box = 10 shows aggregate boxes in which we have reported the edge of the ion foreshock (dashed line) and the angle θ Bn = 45 • (dotted line) for reference. Considering the small number of ions involved in the reflection process, we have used a Gaussian interpolation which gives the relative density weight of each ion. Then, the color coding (vertical bar) gives only an indication of the relative density amplitude. The location of the curved shock front is defined at the middle of the front ramp (thick black line) at the last time t end . Moreover, we have reported the space-integrated percentage value BI % of back-streaming ions within each corresponding box. In order to exclude the gyrating ions present near the front from the back-streaming population, we have eliminated ions found within a small area ≈ 2-3ρ ci upstream of the shock front. For this reason, a very thin white area is visible along the curved front where no particles are present. For N Box = 0-3, the arrows point to the two spots (see the text for explanations). The edge of the foreshock is represented by a mixed line at θ Bn = 70 • . The dashed line (defined for θ Bn = 45 • ) is indicated for reference. P. Savoini and B. Lembège: Two-dimensional simulations of ion foreshock acceleration processes i. The percentage of the back-streaming ions BI % is obtained by computing the ratio of the back-streaming 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 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 time-of-flight 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 so-called θ io,fore ≈ 70 • , as found in our previous self-consistent simulations (Savoini et al., 2013;Savoini and Lembège, 2015).
iii. The back-streaming 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., back-streaming 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 right-hand spot disappears and the back-streaming 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 back-streaming 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 θ hit Bn in red) and when they finally exit the shock front to escape upstream (hereafter named θ exit Bn in blue). These statistical results are obtained by computing these angles for each particle. As a consequence, angle val-ues 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 θ hit Bn 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 θ hit Bn and θ exit Bn 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 int θ Bn = θ exit Bn − θ hit Bn . Obviously, θ hit Bn decreases as N Box increases until it approaches the limit of the quasi-perpendicular domain of propagation, i.e., θ Bn = 45 • for N Box = 9.
Second, the distribution functions of θ exit Bn strongly differ according to the respective box. For N Box = 0-2, two (blue) peaks occur, namely one for high θ exit Bn (for which int θ Bn ≈ 4-5 • ) and the other for low θ exit Bn ( int θ Bn ≈ 15 • ). 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 θ exit Bn distribution (i.e., correspondingly the largest θ Bn ) decreases rapidly in amplitude and disappears from N Box = 6 (i.e., θ hit Bn ≤ 56 • ). 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 T int is defined as the time difference between the time associated to θ exit Bn and to θ hit Bn . Different main maxima of back-streaming 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, , where τ shock ci 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 τ shock ci (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., ∝ −∇B ×B drift) can be a prerequisite for a large drift along the front (i.e., int θ Bn ≈ 15 • ), 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 back-streaming ions filling the foreshock. Savoini and Lembège (2015) have analyzed the impact of E × B drift velocity on the dynamics of back-streaming 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 E 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 Figure 5. FCE configuration: plots of the ion distribution function (for each box N Box = 0 − 9) versus the interaction time range T int spent by each particle within the shock front. As shown, this interaction time range is not continuous, but it evidences distinct bursts of reflected ions (hereafter named f 1 and f 2 ), respectively, defined at T int ≈ 250 ≈ 0.25 τ ci and T int ≈ 950 ≈ 1 τ ci , where τ ci is the upstream cyclotronic period. A third burst f 3 of reflected ions can be identified around T int ≈ 1500 ≈ 1.5 τ ci for N Box = 0 − 2 but can be neglected as compared with the others. 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 back-streaming ions versus the box number where the black and red curves are defined for E l = 0 and E l = 0 respectively. The impact of the E 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 E l = 0, which illustrates the dominant role of E 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 back-streaming 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 ), es- caping 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 θ Bn = 70 • . 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 (θ io,fore ≈ 66 • ) found approximately in the previous self-consistent PIC simulations Savoini and Lembège (2015).

Impact of electric field components
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 θ Bn ≈ 55 • . 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 E l on the dynamics of reflected ions is illustrated in Fig. 7, which shows the two θ hit Bn and θ exit Bn distributions in the same format as that of Fig. 4. The number of back-streaming 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 θ exit Bn distribution is strongly modified and a comparison between Figs. 4 and 7 can be summarized as follows: 1. The boxes N Box = 1-2 evidence a total absence of reflected ions with a small range int θ Bn , and only the θ exit Bn distribution around 60 • persists. This result shows that E 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 θ exit
Bn distribution in both cases ( E l = 0 in Fig. 4 and E l = 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 E l components are essential in the ion reflection for a high θ Bn angle (> 56 • , 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 θ Bn ≤ 56 • seems to be mainly supported by the convective electric field E t components present at the shock front.
Similarly, the one bounce reflection always occurs, even in the absence of E 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 ( T int ≈ 1 τ shock ci ). Then, the ion energy gain is essentially due to a Fermi-type acceleration.

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 back-streaming 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 T = τ ci /20. The time T init = 1248 is the initial time when test particles are launched into the time-dependent simulation.
i. Results of Fig. 8 are obtained as the E l field components are included (black curve) and artificially excluded (red curve). One concludes that the percentage BI % strongly decreases as E l components are excluded, and that the impact of E l field is emphasized for lower N Box . In other words, the back-streaming ions mainly appear for higher N Box > 5 (i.e., for lower θ Bn ), even in absence of E 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 flat-top 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 back-streaming ions between the different boxes, although test particles are initially evenly distributed in the whole boxes. For N Box = 0, backstreaming ions appear around T ≈ 5248, i.e., well after the initial release time T init = 1248 (i.e., T = 5248 − 1248, which covers ≈ 2.6 τ ci ), whereas this time delay decreases to T ≈ 2018 (i.e., T = 2018 − 1248 = 770, which covers ≈ 0.5 τ 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 θ exit Bn (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 E l = 0 (red line). This illustrates that the presence of the E 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 self-reformation, 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 back-streaming ions and, most importantly, on the ion foreshock nonstationarity behavior as described in a previous paper . Nevertheless, it is quite difficult to establish a one-to-one correspondence between these modulations and the nonstationarity of the shock front because, during the sampling time interval T = τ ci /20, the front nonstationarity and the time-of-flight effects have mixed ions coming from either different times and/or different θ exit Bn regions (even if they are in the same box). So this self-consistent 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.

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 λ = v shock × t, since the shock front velocity v shock at a given time must fit with the corresponding value issued from the PIC simu-lations . 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 x = 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 OM −→ (v shock × t)OM, 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 time-of-flight effects are always included. Each front profile is analyzed within a similar simulation time range ≈ 3 τ 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 t init = 1.2 τ ci to t simul = 5.4 τ ci . 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 T = 4456 ≈ 4.2 τ 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 ≈ 3 τ ci .

General features of the back-streaming ions
The purpose is to determine (i) whether some shock profiles are more appropriate than others for the formation of back-streaming 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 max % ≈ 1.2 in the FCE simulations as compared with BI max % ≈ 15 in the HE configuration. In fact, one has to remember that for the FCE model the BI max % value represents an instantaneous reflection rate versus the shock front evolution, whereas this rate is a time-integrated 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: 1. 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 back-streaming ions along the edge of the foreshock.
2. Boxes N Box = 1-8 clearly show some strong modulations in the percentage BI % versus the shock profile of concern, which corresponds to a quasi-periodic bursty emission of back-streaming ions. The BI % rate periodically reaches a maximum value followed by a minimum around zero. The corresponding time period T max is about ≈ 460 = 0.5 τ ci (between two successive max-ima). The temporal width of each maximum is about T range ≈ 256 = 0.25 τ ci . These modulations mean that conditions for the formation of back-streaming 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.
3. 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 T ≈ 3456 and 4456, which indicates a low sensitivity to the shock front profile when approaching θ Bn = 45 • .
4. Similarly, the maximum values BI 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 E 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 E l becomes less important in the reflection process, as evidenced by comparing both black ( E l = 0) and red ( E l = 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. 5. Only the last peak around T ≈ 5456 has a different behavior in comparison with others. In particular, we observe that, for these different times or profiles, the presence of the E l leads to higher BI % . This behavior can be understood because v shock is lower for these times, and the E 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.
6. Finally, Fig. 9 confirms the key role of E l field in back-streaming ion formation, except when approaching θ Bn = 45 • (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 Figure 9. HE configuration: percentage BI % of back-streaming ions measured at the end of each simulation, where each time corresponds to a given fixed shock front. For each shock profile in homothetic expansion, the simulation covers 3 τ ci , allowing one to obtain a welldeveloped ion foreshock, and BI % represents the ratio of the back-streaming ions over the total number of upstream ions which are released at the beginning of the simulation within a given box. As in Fig. 7, black and red lines correspond to the case when the E l field is included and artificially excluded, respectively. The concerned shock profiles are chosen only at late times of the full PIC simulation (from T = 3456 to 5474), where the curvature radius of the shock front is relatively large ( R > 70 ρ ci ).
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 back-streaming ions, especially at higher θ Bn .

Discussions
A previous paper  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 α ≈ 0 • ) or GPB populations (with a pitch angle α = 0 • ), 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 time-varying 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 (v ⊥1 , v ⊥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 θ exit Bn 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 T = 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 Box = 1, 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 (v ⊥1 , v ⊥2 ) distribution. Nevertheless, it evidences approximately a distribution with a maximum slightly noncentered at v ⊥ = 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 θ exit Bn 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 θ exit Bn of Fig. 4, N Box = 8 (in Fig. 10) also evidences a non-Maxwellian-like distribution (α = 0 • ). 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 (v , v ⊥ ); 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 E l = 0, but once again, no reflected ions can be seen when E l = 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 well-formed 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 FAB-like 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 Bale et al., 2005). Finally, for N Box = 8, the velocity distribution drastically changes from a ring ( E l = 0) to a localized bump ( E l = 0) roughly similar to the FCE case. It is clear that the number of reflected ions decreases drastically as E l field components are artificially suppressed. But, more important is that the formation of a nongyrotropic distribution does not depend strongly on the E 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.

Conclusions
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 E l and the E 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 (θ Bn ≈ 70 • ) to θ Bn = 45 • . 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: 1. Impact of the E l field on the ion reflection process. As is well known, the built-up 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 Figure 10. Local perpendicular ion velocity space (v ⊥1 , v ⊥2 ) of all back-streaming ions computed at the end time of the simulations (i.e., after 3τ ci for all simulations) for boxes N Box = 1, 4 and 8 in the FCE approach when the E l is included (left panel); the case E l = 0 is not plotted since the percentage BI % is very weak (see Fig. 8). The right panel shows similar results for the same boxes in the HE configuration (corresponding to T = 4848) in both cases where E l is included and artificially excluded; statistical results where BI % is too weak are not shown. Red and blue hold for maximum and minimum density value in the velocity space.
essentially two distinct impacts where, (i) without this electric component, no reflected ions are observed for θ Bn > 62 • , whereas in presence of this electric component reflected ions (having suffered even one bounce) can be observed at the edge of the ion foreshock (≈ θ Bn ≤ 70 • ); (ii) at lower angles (θ Bn ≤ 50 • ) many ions are reflected without the help of the E l component and can be associated to a magnetic mirror reflection. Then, in Fig. 11, E l is only reported as being strong around the edge of the ion foreshock to emphasize its mandatory action for high θ Bn angles.
2. 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  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 quasi-perpendicular region.
3. 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 Figure 11. Sketch of the ion foreshock in the quasi-perpendicular shock region, illustrating the angular areas along the curved front where the four main identified processes contributing to and/or impacting on the formation of back-streaming ions apply (namely the longitudinal electric field E l , the magnetic field, the shock front nonstationary and the convective electric field E t ). Each process is illustrated by different thick band along the curved front which is shifted, one from the other, in order to avoid overwhelming the sketch. One color (red, blue, green and black, respectively) is associated to each process. The varying intensity of the color indicates where the process is strong (full color) or weak (white). This allows one to identify, at a glance, the angular areas where the different processes are complementary or accumulating. The upstream interplanetary magnetic field (IMF) is reported in gray (dotted and dashed lines), as is the shock front itself. Typical directions θ Bn = 90 and 70 • (blue and red dashed lines) defined between the local shock normal n and the IMF indicate the location where the electron and ion foreshock edge initiates from the curved shock front, respectively. The electron foreshock edge is indicated as a reference.
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 quasi-perpendicular region, we emphasize, in Fig. 11, its stronger impact near θ Bn ≈ 45 • where it is mainly responsible for the reflection (i.e., magnetic mirror) and acceleration (i.e., Fermi-type) of ions (Webb et al., 1983).
4. 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 back-streaming 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 1 τ shock ci (where τ shock ci 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 back-streaming 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 back-streaming 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 Ann. Geophys., 38, 1217Geophys., 38, -1235Geophys., 38, , 2020 https://doi.org/10.5194/angeo-38-1217-2020 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 time-offlight 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. Data availability. 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.
Author contributions. 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.
Competing interests. The authors declare that they have no conflict of interest.