Evidence of the Nonstationarity of the Terrestrial Bow Shock from Multi- Spacecraft Observations: Methodology, Results and Quantitative Comparison with PIC Simulations

. The nonstationarity of the terrestrial bow shock is analyzed in detail from in situ magnetic field measurements issued from the FGM experiment on board of Cluster mission. Attention is focused on statistical analysis of quasiperpendicular 10 supercritical shock crossings. The present analysis stresses for the first time the importance of a careful and accurate methodology in the data processing which can be a source of confusion / misunderstanding if not treated properly. The analysis performed using 96 shock front crossings shows evidence of a strong variability of the microstructures of the shock front (foot and ramp) which are analyzed in deep details. Main results are : (i) most statistics clearly evidence that the ramp thickness is very narrow and can be as low as a few c/  pe (electron inertia length), (ii) the width is narrower when the angle θ Bn (between 15 the shock normal and the upstream magnetic field) approaches 90°, (iii) the foot thickness strongly varies but its variation has an upper limit provided by theoretical estimates given in previous studies (e.g., Schwartz et al., 1983; Gosling and Thomsen, 1985; Gosling and Robson, 1985); (iv) the presence of foot and overshoot, as shown in all front profiles confirms the importance of dissipative effects. Present results indicate that these features can be signatures of the shock front self-reformation among a few mechanisms of nonstationarity identified from numerical simulation/theoretical works. 20 A comparison 2D PIC simulation for a perpendicular supercritical shock (used as reference), has been performed and it shows that: (a) the ramp thickness varies only slightly in time over a large fraction of the reformation cycle and reaches a lower bound value of the order of a few electron inertial length, (ii) in contrast, the foot width selected shock crossings, the  Bn values extend from almost 90° to 75° but mostly above 84° (80%), the M A values of the resulting shocks are equally distributed from 2 to 6.5, the ratio v sh / v thi between 2 and 30 with a peak within the range 12-24, and  i extends from very low values to 0.6 but with 67% 400 values less than 0.1. the shock is strictly perpendicular ( θ Bn = 90°) where the upstream magnetostatic field B 0 is along z axis, and propagates within the (x-y) simulation plane; this perpendicular shock is used as a reference case. All quantities are dimensionless, are indicated by a tilde ‘ ~” and are normalized as follows : the spatial

Third, an improved approach has been supported by the quadri-satellite Cluster mission which not only allowed to analyze in great details the fine structures within the shock front itself, but also clearly evidenced that the terrestrial bow shock can be strongly nonstationary (Horbury et al., 2001(Horbury et al., , 2002Moullard et al., 2006;Lobzin et al. 2007;Lefebvre et al., 2009;Meziane et al., 2014Meziane et al., , 2015. More generally, the nonstationarity of quasi-perpendicular shocks in a supercritical regime has been also evidenced in plasma laboratory experiments several decades ago (Morse et al., 1972), 70 and largely analyzed in numerical simulations (Biskamp and Welter, 1972;Forslund et al., 1984;Lembège and Dawson, 1987;Winske and Quest, 1988;Lembège and Savoini, 1992;Lembège et al., 2008Lembège et al., , 2009Nishimura et al., 2003;Schöler and Matsukiyo, 2004;Chapman et al., 2005). Particle-in-cell (PIC) and hybrid simulations clearly evidenced different processes as sources of the shock front nonstationarity which are summarized as follows (Lembège and al., 2004): a) the self-reformation driven by the accumulation of the reflected ions at a foot distance from the ramp. During their 75 coherent gyromotion, the reflected ions accumulate at roughly the same location and build up a foot whose amplitude largely increases in time until its upstream edge suffers an important steepening and becomes a new ramp. New incoming ions are reflected at this new ramp and the same process continues cyclically. This process has been shown to be quite robust since it was evidenced in hybrid and PIC simulations as well as in 1D (Biskamp and Welter, 1972;Lembège and Dawson, 1987;Lembège and Savoini, 1992;Hada et al., 2003;Schöler et al., 2003;Schöler and Matsukiyo, 2004), in 2D (Lembège and 80 Savoini, 1992), and 3D PIC simulations (Shinohara et al., 2007). The self-reformation persists for oblique propagation and disappears (i) when the oblique direction of the shock front normal is below a critical angle for which the density of newly reflected ions is too weak to feed the ions accumulation (Lembège and Savoini, 1992), and (ii) for shocks with higher βi (the ratio between thermal ion pressure and magnetic pressure) (see Sect. 5.1) as shown by Schmitz et al. (2002) and Hada et al. (2003), or equivalently as the ratio vsh/vthi, where vsh is the shock velocity in the solar wind rest frame and vthi is the upstream 85 ion thermal velocity (solar wind protons), is relatively weak (a few units) (Schöler et. al., 2003). b) a similar self-reformation is also observed in 1D PIC simulation of an oblique propagating shock (in a configuration not exactly perpendicular) but has a different origin since it is mainly driven by a microinstability (MTSI or Modified Two Stream Instability) excited by the relative drift between incoming electrons/incoming ions within the foot region. This instability develops quickly enough to thermalize the ions and generates a pressure increase which contributes to the build-up 90 of a new ramp (Matsukyio and Schöler 2003;Schöler and Matsukiyo, 2004). In addition to low βi, high/realistic mass ratio and a slightly oblique propagation are necessary conditions (small but finite component of wave vector parallel to the ambient magnetic field) for the MTSI to be excited. c) the so-called gradient catastrophe process (Krasnoselskikh et al., 2002) which results only from an imbalance between nonlinear and dispersive effects (Galeev et al., 1999), is expected as the Alfven Mach number MA is above a critical threshold 95 (or equivalently the shock front propagation is below a critical value). When this condition is satisfied, large amplitude nonlinear whistler waves are emitted from the ramp and are responsible for the front nonstationarity. This process also requires https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. 4 a slightly oblique propagation necessary for the emission of these nonlinear whistler waves. However, this theoretical process has very strong constraints as it is based on 1D model only and excludes any dissipative effect. Indeed, reflected ions (which play a key role in supercritical Mach regime) are neglected and their impact on the shock dynamics is totally absent. 100 d) the emission of large amplitude whistler waves within the front which propagate obliquely both to the shock normal direction and to the static magnetic field; this mechanism has initially been evidenced for a perpendicular shock (Hellinger et al., 2007;Lembège et al., 2009), and later for a quasi-perpendicular shock by Yuan et al. (2009). It is basically a 2D mechanism, which has been retrieved both in 2D-hybrid and 2D-PIC simulations and has been observed in 3D PIC simulation (Shinohara et al., 2007). One expects that the presence of these waves within the front enlarges the width of the ramp (since these propagate 105 almost at the same velocity of the shock front), but no detailed measurement has been reported yet. One proposed mechanism could be the ion Weibel instability triggered by the reflected gyrating ions (Burgess et al., 2016). Up to now, the mechanism responsible for this wave emission has not been clearly identified yet. This paper represents an extensive and more detailed analysis of previous works (Mazelle et al., 2007(Mazelle et al., , 2008a(Mazelle et al., , 2008b(Mazelle et al., , 2010a(Mazelle et al., , 2010b(Mazelle et al., , 2011(Mazelle et al., , 2012. Its aim is motivated by the following questions: (i) what important cautions are necessary when 110 measuring carefully the spatial widths of the fine structures (foot and ramp) of the shock front from experimental data? (ii) what are the sources of confusion or misunderstanding when comparing the present results (based on statistics of shock crossings) with interpretations of experimental data made in previous analysis mainly based on one or a few shock crossings only ?, (iii) among the different mechanisms proposed as sources of the shock front nonstationary is it possible to extract some features from the shock front microstructures (observed in experimental data) allowing to clearly identify the dominant 115 nonstationary processes?
The paper is organized as follows. Section 2 presents the motivation, data selection and methodology. A detailed analysis of individual shock crossings is presented in Sect. 3. The results of a statistical analysis and the comparison with PIC numerical simulations results are reported in Sect. 4. Section 5 is devoted to the comparison with previous studies (both experimental and numerical) while Sect. 6 draws some conclusions. 120 2 Motivation, data selection and methodology

Experimental data
The constraint on time resolution necessitates the sole use of the magnetic field data for the present detailed analysis. Particle data will be used only to provide the necessary "contextual" plasma parameters, since their time resolution is commonly much lower than the magnetic field data. The latter are provided by the fluxgate magnetometer (FGM) (Balogh et al., 2001). Since 125 the thicknesses of the shock front micro-structures must be determined accurately, the use of the best time resolution is vital, i.e., 22 Hz to 66 Hz depending on the available telemetry. But, lower time-resolution data (5 Hz and 1 Hz) are also used for comparison and to confirm the temporal extension of shock front micro-structures (mainly the ramp and the magnetic foot; the https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. 5 overshoot/undershoot are not analyzed in detail herein) when the level of microturbulence is relatively high. This is also very important for validating the determination of the shock normal vector (see Sect. 2.2). 130 Using the appropriate data from both Cluster and ACE satellites, the major shock parameters have been derived: the Alfvenic and magnetosonic Mach numbers MA and MMS, the angle Bn between the upstream interplanetary magnetic field and the local shock normal, and the ion i. We use herein Cluster particle data used from the Cluster Ion Spectrometer (CIS) experiment which is extensively described in Rème et al. (2001). It consists of two sensors: (i) a Hot Ion Analyzer (HIA) measuring particles within the energy range 0.005-26 keV, and (ii) a time-of-flight mass spectrometer (CODIF), combining a 135 top-hat analyzer with a time-of-flight section to measure the major species H+, He+, He++ and O+ covering an energy range extending from 0.02 to 38 keV/q. Both instruments allow to measure full 3D distributions with a typical angular resolution of 22.5°×22.5° within one satellite spin period (4 s). Two operational modes are possible: in "normal" telemetry modes, one distribution is transmitted every two or three spins (depending on the time of measurements), whereas in "burst" telemetry modes the distribution is transmitted at every single spin. With the instrumental characteristics mentioned above, both the solar 140 wind plasma and the energetic particles are detected. The measurements of the HIA analyzer are performed with high geometry factor (HIA-G) suitable for upstream ion measurements as well as with low geometry factor (HIA-g) suitable for the solar wind plasma measurements which use specific narrower angular bins (5.6°x5.6°). We are primarily interested in the moments of the velocity distribution obtained from both 'solar wind' modes and 'magnetospheric' modes of the HIA instrument. The best determination of the upstream solar wind parameters (density, velocity, and temperature) is obtained for a 'solar wind' 145 mode where HIA-g data are available. For a 'magnetospheric' mode, the solar wind velocity is still correctly determined though less accurately than in a 'solar wind' mode; however, the density is often underestimated, and the temperature is overestimated by a large factor when compared to other data sets (for details, see Rème et al., (2001)). The density underestimate is due to a saturation effect for the solar wind beam while the temperature overestimate is due to the use of wide angular bins of the High-G part of the HIA instrument. In this case, we use ACE satellite data considering the proper time-shift between ACE and the 150 Cluster constellation. show the profiles of the magnetic field magnitude measured by FGM at the spin time resolution (4s-average). The profiles seem very similar at first glance and do not allow any precise determination. In contrast, the bottom panels show the same shock crossing with high time resolution data. Many differences are then clearly revealed, which illustrate the necessity to use 155 the best possible time resolution. The methodology used here needs a relatively small separation distance between all four satellites in term of ion inertial length. This was the case only during the spring seasons of 2001 and 2002 for which all the inter-distances lied between 100 to 600 km typically, which is no more the case for more recent satellite configurations (Escoubet et al., 2015). Moreover, the Earth's bow shock is often crossed many times by the four satellites in this configuration as the quartet remains nearby the shock front itself while the shock is globally moving inward and outward. 165 https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.

Methodology
We detail for the first time a new methodology to carefully analyse experimental data for the study of super-critical quasiperpendicular shocks. The step-by-step analysis method we have developed is illustrated on Figure 2 and is conceptually straightforward. First, for each satellite, the temporal extension of each micro-structure of the shock front is determined in the time series of the magnetic field. For each of them, beginning and end times are first estimated by visual inspection. For many 170 previous studies in the literature, the process usually stops at this level leaving the possibility of strong bias. Here, this first step is made only to initialize later automatic algorithms which will be used eventually to get spatial error bars (both for the outer and inner edge of each micro-structure). For each beginning or end time of a micro-structure, we get a first estimate of an error bar on it by defining two associated times: one time (e.g. ta) when the spacecraft has already likely entered the shock front or left the previously estimated substructure; a second time (e.g. tb) is defined as when the spacecraft is inside the present 175 substructure for sure. The order of these two times obviously depends whether there is an inbound or outbound bow shock crossing. This may appear arbitrary and occasionally difficult to estimate but the experience shows that the risk to miss the real beginning/end time inside the defined interval is low.
The determination of upstream parameters entails a steady asymptotic upstream interplanetary magnetic field (IMF) B0 associated with each shock crossing by the four spacecraft. For highly perpendicular shocks, it is however often possible to 180 define an unperturbed upstream field on a time scale around a few minutes. For each satellite, we must define an upstream time interval displaying a magnetic field quite steady both in components and magnitude. For each spacecraft s/ci, where s/ci holds for spacecraft "i" among the Cluster quartet (where i = 1, 4), an average upstream field B0i is computed. Its magnitude is shown in the lower green horizontal line in Figure 2a for one example of a single spacecraft shock crossing. Again, an initial estimated time interval is obtained from the visual inspection of the time series. The time duration cannot be too long since the 185 usual stability time scale of the IMF due to the intrinsic solar wind variability (turbulence and directional discontinuities) is rarely more than a few tens of minutes (e.g., Burlaga, 1971) which we can take as an upper limit Δtu. In order to be physically reliable, Δtu cannot be too short. One reasonable lower limit Δtl is the time interval during which a solar wind parcel travels a distance of one ion inertial length (about 100 km usually) that is about one minute typically, considering the mean spacecraft velocities are all consistent within a reasonable range. In this purpose, the six relative angles ij measured between upstream 190 fields B0i and B0j for s/ci and s/cj with (i,j = 1 to 4) are computed. The upstream vector B0i measured by satellite "i" must be as possible close to the vector B0j measured by another satellite "j". If all angles ij are less than 5° (arbitrary value chosen herein from empirical considerations), then the procedure is pursued, otherwise, the upstream time intervals (where these ij values are calculated) are automatically modified by slightly shifting and reducing them until the above constraint on ij is satisfied.
The shock event is rejected from the present analysis when the adjustment fails or not satisfactory. 195 When possible, the use of a larger time interval improves the statistics. Using the initial estimate of Δt and constraining it to be between Δtl and Δtu, the finest upstream interval corresponds to the lowest values of the standard deviations  for both the magnitude and components of the field. Though this can be easily verified by visual inspection, it allows to ensure the https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. 8 absence of any solar wind disturbance. The following step deals with the determination of a unique unperturbed ambient IMF B0 for every shock crossing by the four-spacecraft constellation. First, it is necessary to compare the individual upstream fields 200 B0i. This determination of B0 is of crucial importance for the present study since it impacts both the value of the shock Alfven Mach number through the field magnitude used and the determination of the angle Bn0 between B0 and the global shock normal n0. This is obviously the first possible source of error in the determination of both the upstream Alfven velocity, the Alfven and magnetosonic Mach numbers and the angle Bn0.
For each satellite, the standard deviations  obtained above for the magnetic field components and magnitude allow the 205 determination of the time when the spacecraft penetrates inside the magnetic foot. We determine, for each component and magnitude, the time when it is above 3- threshold level and close to the initially estimated time of the entry in the foot obtained by visual inspection (ta blue dashed line in Figure 2a). Then, we select the 'external' limit of this entry (t1 in Figure   2a) among these four times as the closest to the external estimate ta and the 'internal limit' (t2 in Figure 2a) as the closest to the estimated internal estimate (tb blue dashed line in Figure 2a). The derived times t1 and t2 then define the error bar of the 210 foot entry in the time domain. For the present case, the visual estimate tb nearly coincides with the time t2 automatically determined. Moreover, the possible presence of solar wind or foreshock transients as well as detached whistler wave packet must be avoided. The later occurs quite often even in the nearly perpendicular configuration and can easily leads to a wrong determination. Also, solar wind transients are always possible while foreshock transients may only be an issue for more 'oblique' configuration which is not addressed in the present study. Nevertheless, it is necessary to make a cautious post visual 215 inspection.
Similarly, the initial estimated times of the other substructures are labelled with letters and the final times determined from the present methodology are labelled with numbers. The times tc and td are the initial 'visual' estimates for the exit from the foot and entry in the ramp respectively. They will be only used for initializing the automated algorithm described below.
A major caveat in many previous attempts to compute the ramp thickness in the literature is to include in the ramp the part of 220 the magnetic field profile up to the maximum of the first overshoot. To avoid that and to well define at least the first overshoot we need to determine a downstream (asymptotic) level which would be required for instance for checking the Rankine-Hugoniot conditions. For that purpose, it is necessary to define a downstream interval in the magnetosheath far enough from the consecutive overshoots/undershoots displaying a quasi-steady magnetic field. This choice is highly important and not always straightforward in the case of very close consecutive crossings (multiple crossings) where only an average over the 225 region of the sequence of overshoots / undershoots will be estimated. When a "clean" time interval downstream from this overshoots / undershoots interval is available, the average values obtained for B magnitude and components can be compared to the values inside the overshoots/undershoots interval. If these values are in good agreement, the downstream interval is appropriate for the determination of the asymptotic value. When a "clean" time interval downstream from this overshoots/undershoots interval is available, the average values obtained for B magnitude and components can be compared 9 to the values inside the overshoots/undershoots interval. If these values are in good agreement, the downstream interval is appropriate for the determination of the asymptotic value. It is then used to give an initial estimate (te, tf) of the final time  interval (t5, t6) for the exit from the ramp/entry in the overshoot by using +/-5% of the computed asymptotic value (which is arbitrary but appears to be a good empirical compromise). It should be mentioned that a 1- level can also be determined here but it will usually be large in this region (local microturbulence in the magnetosheath) and not appropriate for the present purpose. The two times (te, tf) are again only used for initializing the next automated step.
From the initial time estimates (tc, td) and (te, tf), the determination of the final times for the entry in and exit from the 245 ramp (t3, t4) and (t5, t6) respectively, is made with an automatic algorithm. In order to characterize a steep gradient, we perform a linear fit of the data points between tc and tf. The time range for the regression can vary on both sides. The process is initialized by the largest possible time interval with t3=tc and t6= tf. Then, the entry time can vary between tc and td, while the exit time can vary between te and tf. The choice of a linear fit mainly excludes any "pollution" of the ramp region by a part of the extended foot (around the upstream edge of ramp) and by downstream fluctuations around the overshoot (around the 250 downstream edge of the ramp). The resulting time interval is associated to the best linear regression coefficient and the bounding times (tc and tf) explored during the automatic procedure are kept to estimate the error bars. In other words, we impose t3 = tc and t6 = tf . The automatic procedure only defines t4 and t5 (internal limits of the ramp). Then, the crossing of the downstream asymptotic value closer to the first estimate defines the exit from the overshoot (nearly equivalent to the middle of the gradient between the first overshoot and the first undershoot, not shown in detail here). The times t7 and t8 for the exit 255 from the first overshoot are defined in the same way as the entry (Figure 2a). The time tg shown by a vertical blue dashed line is -as for the other microstructures-simply the initial rough 'visual' estimate and is just shown for reference for a later visual check. For every shock crossing by the Cluster quartet, we make the same procedure for each satellite. The satellite that is associated with the steeper ramp from this analysis is used to arbitrarily define what is called the 'reference satellite' for a shock crossing event. Then, we determine the middle time tramp,i of the ramp crossing time interval for each individual satellite 'i' and 260 the 'reference time' for a shock crossing is simply this middle time tref for the arbitrarily chosen 'reference satellite'. We then determine the shock velocity and normal from the four different middle times tramp,i and the four spacecraft positions in the GSE frame by using the classical 'timing-method' (e.g., Schwartz, 1998and Horbury et al., 2001. This method works if the separation between the spacecraft is not too large (otherwise the hypothesis of a planar structure breaks down) and provided that the four spacecraft are not coplanar (otherwise the method fails) and that the shape of the tetrahedron is not very 265 elongated (which is not the case for the present study). This method provides a 'global' normal vector n0 valid only at the length scale of the spacecraft tetrahedron (i.e. the largest inter-distance between two satellites) and at the temporal scale defined by the time interval between the first crossing and the last crossing among the four satellites. The error bars on the entry/exit times are used to determine an error bar on the normal vector. It must be mentioned that each local normal vector ni for the crossing by s/ci can differ from the global n0 because of the small-scale turbulence of the front (shock front ripples) but cannot 270 be determined experimentally as reliably as n0 by any classical single-spacecraft method. Moreover, as it will be discussed later, we intend here to compare the values determined experimentally with simulations results for which the normal vector n0 and thus the Bn0 angle are 'globally' defined for a shock front both averaged spatially perpendicular to the normal and in time (see Sect. 4.4). This will also conversely justify the averaging along the shock front made in the simulations (that will be https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. described in Sect. 4.4). The local effect due to the variation of the local normal direction ni with respect to B0 (local curvature 275 effect due to the shock rippling at a lower scale) is not analysed in the present work.
For all the shock front micro-structures, the above analysis gives an 'apparent' temporal extension along each satellite path and allows to compare them between the four satellites. In every individual spacecraft frame, we compute an individual shock velocity vector. Since the spacecraft velocities in the Earth frame are typically around a few km/s as also found for many slowly moving shocks in the same frame, it is essential to also take this small spacecraft velocities into account. The major 280 aim of the present analysis is to reconstruct a local spatial profile along the 'global' shock normal n0 as shown in Figure 2b.
For this purpose, the angle between each spacecraft 'i' trajectory and n0 has been used. Then, it is possible to measure the local physical thicknesses along n0 of the microstructures of the shock front and compare them between the different satellites. It must be mentioned that the shock velocity determined by the 4-spacecraft timing-method here is an averaged velocity within the time interval separating the first and the last shock crossings. Then, it is supposed to be constant during this time interval 285 and homogeneous over the spatial scale of the tetrahedron. One cannot preclude that the shock velocity can vary inside the considered time interval and that the shock motion is subject to some acceleration inward and outward. However, this possible effect can be assumed as weak if the temporal profiles of the shock crossing are quite similar in time series taken at low temporal resolution (such as on Figure 1a). On the contrary, when a sudden shock front acceleration/deceleration (due to a solar wind fluctuation) occurs between the crossings by two satellites, there is a clear signature in the time series (see in Figure  290 3 derived from Figure 1 in Horbury et al., (2002)). Panels 3b and 3c are typical cases of observations which are eliminated from our present analysis. We use the same shock velocity in the Earth's frame for the 4 crossings to derive the spatial widths along the global normal n0 for any case like in panel 3a.
Then, we can compare the obtained thicknesses of the shock microstructures with physical length values predicted either by theory or from numerical simulations results. It must be mentioned that using only time series can lead to wrong comparative 295 classification of the observed crossings. A large angle between the spacecraft velocity vector and the shock front normal can produce a long crossing duration for a real very small physical thickness and also will produce a small shock speed. In summary, a particular caution is necessary to avoid conclusions deduced exclusively from the 'temporal' profiles, especially in the selection process which is often based in the literature solely on the steepness of the temporal 'gradient'.
In order to have a qualitative check of the determination of the global normal n0, the times series of the component Bn of 300 the instantaneous magnetic field along this normal is systematically examined. First, this normal component is supposed to be conserved on both sides of the shock front. Second, for a nearly perpendicular shock, this component is expected to vanish (exactly zero for 90° shock) upstream from the shock front and inside the ramp, since the vector of the magnetic field is supposed to be strictly along the shock front. Figure 4 displays the temporal evolution of Bn for a shock front crossing on 03/31/2001 around 1738 UT for one of the four satellites (magenta lower curve) in comparison with the total magnitude B 305 (black upper curve). For this shock, the Bn angle is about 86°. If one takes into account the fluctuations observed in and around the ramp, the conservation in the front is well respected and the upstream value is also very low as expected for a nearly https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. perpendicular case. The downstream evolution in the overshoot is also consistent with expected oscillations around the average downstream value (e.g., Kennel and Sagdeev, 1967).

One typical case
A typical example of the spatial profiles is shown in Figure 5. It is obtained for one supercritical quasi-perpendicular shock 320 crossing by Cluster on 03/31/2001 at 1738 UT (Bn=86±2°, MA=4.1, i=0.05) when the s/c separation was about a few ion inertial lengths c/pi (herein c/pi = 80-120 km, typically). The profiles are ordered by the sequence of crossings starting from satellites C4 to C3 (top to bottom in Figures 5a and 5b). The relative locations of the four s/c defined at the 'reference time' given by satellite C4 (see Sect. 2.2) are represented in the plane containing both the XGSE direction and the normal n0 ( Figure   5c) and in the plane perpendicular to the normal n0 ( Figure 5d). 325 In a first approach, the differences between the 4 s/c spatial profiles in panel 5b are quite obvious in particular for the magnetic foot and the overshoot widths; however, these are less evident for the ramp. We emphasize that all upstream conditions remain unchanged for all four satellites. As a consequence, this evidences the strong variability of the quasi-perpen-https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. dicular shock front itself. We stress that the ramp widths can be very narrow and can reach values as low as a few electron inertial lengths with the thinnest one reaching 5 c/pe for satellite C4. This result strongly disagrees with the typical scale length of the ramp commonly considered as around c/pi (e.g., Russell and Greenstadt, 1979;Scudder et al., 1986;Bale et al., 2005) although in a very good agreement with some "exceptional" ISEE results Russell, 1996, 1998) as well as 335 theoretical studies predicting the existence of small scales in quasi-perpendicular shocks as a signature of the front nonstationarity or results issued from PIC numerical simulations (Lembège and Dawson, 1987;Lembège and Savoini, 1992;Lembège et al., 1999Lembège et al., , 2013aLembège et al., , 2013bSchöler et al. 2003, Marcowith et al, 2016 andreferences therein).
This nonstationarity can be also evidenced in the content of Figure 5. The satellites C3 and C4 locations are almost aligned with the direction n0 ( Figure 5c) and are close one to each other (over a distance less than c/ pi) in the plane perpendicular to 340 n0 ( Figure 5d). If the shock front was stationary, after propagating from C4 to C3, the profiles would be expected nearly identical while they appear significantly different. For instance, the spatial ramp thickness for C3 is about twice the one for https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. C4, and the shape and thickness of the overshoot are also very different. Moreover, C2 and C1 locations are not far to be along the n0 direction ( Figure 5c) but they well depart one from each other in the perpendicular front (i.e. along the shock front) by 350 several c/pi (Figure 5d). The differences mentioned here can be related to the shock front nonuniformity due to its rippling; the impact of this rippling on the front variations will be analysed in a future work.

Impact of the shock velocity
A careful determination of the shock velocity is necessary in order to discuss the physical spatial thicknesses of the shock front substructures. This also means that it is not possible to consider only times series to discuss the relative properties of different 355 shock crossings. Moreover, as already stated previously, the time resolution of the data needs to be considered before any https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.
conclusion. This is illustrated in Figure 6. The left part shows the time series of the magnetic field components and magnitude for two shock front crossings and for two different time resolution each (panels a1 and b1): the highest available resolution (15 or 60 ms depending on the mode) and the spin resolution of 4s. These two shocks have been chosen as they are characterized with similar parameters Bn, MA, and i. The first shock 'a' (panel a1) is a 'slow' shock with a velocity of about 11 km/s while 360 shock 'b' (panels b1) is a 'fast' shock with a velocity of about 78 km/s. When considering the lowest time resolution (averaged data on 4s), the two profiles in the upper plots of panels 'a1' and 'b1' appear very. If only data with this time resolution were considered, one should have concluded at first glance that the two shock crossings are very similar and display very similar properties with relatively smooth time gradients of a similar shape. But, the highest time resolution on the lower panels reveals 365 Figure 6: Comparison of a slow moving shock (panels (a1) and (a2), where vsh= 11 km/s, MA= 3.8, Bn = 89°, and I =0.04) and a rapid shock (panels (b1) and (b2), where vsh = 78 km/sec, MA = 3.5, Bn = 89°, and i = 0.045) defined for low (dt=4 s) and high resolution rates (dt=15 or 60 ms). Panels (a1) and (b1) provide the temporal profiles of the total magnetic field amplitude (defined respectively for C1 and C3, while panels (a2) and (b2) provide correspondingly the spatial profiles of the measured magnetic field (the distance is normalized with respect to 370 the ion inertial length defined upstream). The spatial origin corresponds to the location at the middle time of the ramp crossing interval for each satellite.
the large difference between shock profiles when changing the data sampling: shock 'b' appears steeper than shock 'a'.
Moreover, for shock 'a', the gradient remains nearly unchanged with only superimposed higher frequency fluctuations which 375 are not observed upstream from the sharp front displayed by shock 'b'. The full shock properties are revealed only if the spatial profiles, as shown on panels 'a2' and 'b2', are considered. The general shapes of the shock front substructures are very similar https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. for the two shocks, but the steeper spatial ramp is obtained for shock 'a' while shock 'b' displays a ramp spatial thickness more than twice that of shock 'a' (right hand panels (a2) and (b2) respectively). Moreover, it also reveals that the structure of the shock is better sampled up to much higher frequencies for a slow shock velocity. In the spacecraft frame, any instrument like 380 the magnetometer uses a constant sampling frequency and thus obviously gets much less measurement points inside any shock front substructure when this one is propagating faster. This example illustrates quite well the advantages of the methodology developed in the present study to correctly address the spatial scales of a shock front and the deep caution to keep in mind when analysing only the temporal shocks profiles.

Selection criteria
An analysis similar to the one described in Sect. 3 has been performed for a list of shocks after a severe selection procedure.
Several criteria have been used: (i) limitation to shocks normal direction Bn approaching 90° as much as possible (more than 75°) in order to avoid more complex mechanisms which take place at oblique shocks (Lembège and Savoini, 1992;Lembège et al, 2004;Marcowith et al. 2016), (ii) existence of well-defined upstream and downstream ranges for each of the 4 s/c, (iii) 390 stability of the upstream averaged magnetic field from one s/c to another, (iv) validity of the normal determination by checking the weak variability of the normal field component Bn around the ramp and low value of Bn/|B| upstream for Bn close to 90°, (v) availability of all data necessary to compute the shock parameters and (vi) a maximum spacecraft separation distance limited to a few c/pi between the four spacecraft in order to consider the crossing of a same "individual" shock front by all Cluster spacecraft. From an initial ensemble of 455 shock crossings (including essentially data for 2001 and 2002 for which 395 spacecraft separation distances are all small enough), only 24 shock events have been selected after the validation of the aforementioned criteria. This implies a detailed analysis of 96 individual shock crossings. Figure 7 displays the distributions of four main parameters Bn, MA, i and vsh/vthi (see Sect. 1). For the selected shock crossings, the Bn values extend from almost 90° to 75° but mostly above 84° (80%), the MA values of the resulting shocks are equally distributed from 2 to 6.5, the ratio vsh/vthi between 2 and 30 with a peak within the range 12-24, and i extends from very low values to 0.6 but with 67% Figure 7: Distribution of the 96 (selected) shock crossings versus (a) the angle Bn (defined between the shock normal direction and the upstream magnetic field), (b) the Alfven Mach number (MA= Vsw,n / vA, where Vsw,n is the solar wind bulk velocity component defined along the shock normal, and vA is the upstream Alfven velocity respectively), (c) the ratio i (between the ion thermal pressure and the magnetic 405 pressure), and (d) the ratio vsh /vthi where the vsh is the velocity of the shock front defined in the rest plasma frame (to be compared with present simulation results) and vthi is the upstream solar wind proton velocity.

Ramp thickness
The obtained ramp thickness distribution shown on Figure 8a appears Gaussian, indicating that most of the shocks are 410 associated with a narrow ramp. Interestingly, the distribution reveals a cut-off mark at c/pi. Figure 8b   Although Figure 8c shows no clear dependence of the ramp width with Bn , it however indicates that the narrowest ramps 425 are associated with nearly perpendicular geometries as pointed out by the red arrow. Figure 8c also shows that the ramp thickness variation is not solely controlled by the shock-geometry. It can vary a lot for a single value of Bn except maybe https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.
when close to 90. No less important, the ramp thickness obtained from the analysis of a same crossing with each single s/c appear inconsistent (up to a factor 7) although the upstream conditions remain similar. This illustrates quite well the signature of an intrinsic shock front nonstationarity. Figure 8d does not reveal any simple relation of the shock ramp thickness with MA. 430 It may only illustrate that the maximum observed thickness reduces as the Mach number increases, although the scatter plot and the statistics do not allow to make any conclusive result.
The statistics confirms the important result evidenced in the individual case study of Sect. 3, that the magnetic field ramp of a supercritical quasi-perpendicular shock often reaches a few c/pe. The time variability of scales observed for the same upstream conditions is also consistent with a signature of the front nonstationarity due to the self-reformation fed by the 435 accumulation of reflected ions in agreement with PIC and hybrid simulations results (Lembège and Dawson, 1987;Biskamp et Welter, 1972;Lembège et al., 2004;Lembège and Savoini, 1992;Hellinger et al., 2002). Other mechanisms raising from the modified two stream instability (triggered by the relative drift between incoming electrons and reflected ions within the foot region of the front) for oblique shocks  may also be responsible for the nonstationarity of the shock front as discussed in Sect. 5. 440

Magnetic foot
The left panel of Figure 9 displays the distribution of the magnetic foot thickness for all the 96 shock crossings by any individual spacecraft versus the foot thickness (normalized with the upstream ion gyroradius ) whereas the right panel shows the maximal value among the four satellites for each selected shock crossed by the quartet (labelled with the numbers on the horizontal axis) versus the foot thickness normalized with . 445 First, the distribution (Figure 9a) clearly shows that the foot thickness is always found less than the upstream gyroradius.
Second, most of the foot thicknesses (81%) are found to be less than 0.4 , and the median value is 0.2 , . Most of the observed shocks are characterized by a foot with a thickness (Lfoot < 0.1 , ) indicating that the foot width rarely reaches its maximal value. The literature provides different theoretical determinations of the foot thickness. Gosling and Robson (1985) used the results of Schwartz et al. (1983) to show that for arbitrary angle , the upstream extent of the foot is: 450 (2) and is the angle between the shock normal and the upstream solar wind velocity vector. 455 Another previous determination was given by Livesey et al. (1984) as: For a strictly perpendicular shock ( = 90°) and a normal incident upstream solar wind ( = 0°), the foot width then reaches the distance for a reflected ion at the time of its turn-around defined by dxn/dt = 0, where xn is the coordinate of the ion along the normal, given by Woods (1969;: 460 https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.
The values of d obtained from Eq. (1) -Eq.
(2) are shown by red points on Figure 9b for each shock event. For each event, the blue value with error bar is the largest experimental foot thickness for one among the four satellites. The horizontal blue dashed line shows the theoretical value given by Eq. (4) and clearly stresses that the maximal experimental foot thickness is usually below the 90°-critical value and the value obtained from the experimental derivation of and (considering 465 the error bars). Thus, for each specific shock crossing and associated parameters, the above two theoretical values which are supposed to be fixed by the two aforementioned angles (and thus not varying with time) can be generally considered as an upper bound of the observed "instantaneous" foot thickness.  Gosling and Robson, 1985). For reference, the value Lfoot / ci,us = 0.68 defined for Bn = 90° (Woods, 1971) is indicated by the blue horizontal dashed line in panel (b). 475

Comparison with numerical simulation results.
To compare more precisely with numerical simulation results, two-dimensional PIC simulations (where both ions and electrons are treated as individual macroparticles) have been performed like those of previous works (Lembège and Savoini, 1992;Lembège et al., 2009), where the planar shock is initialized by a magnetic piston (applied current pulse). Therefore, the shock 480 geometry is always defined in the upstream frame: the shock propagates along the x direction. Periodic conditions are used along the direction of the planar shock front (y-axis). Herein, the shock is strictly perpendicular (θBn = 90°) where the upstream magnetostatic field B0 is along z axis, and propagates within the (x-y) simulation plane; this perpendicular shock is used as a reference case. All quantities are dimensionless, are indicated by a tilde ' ~" and are normalized as follows : the spatial https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. coordinate is = x/Δ; velocity ̃ = v/ωpeΔ; time ̃ = ωpet, electric field ̃ = eE/meωpe 2 Δ; magnetic field ̃ = eB/mewpe 2 Δ. The 485 parameters e, me , Δ, and ωpe are respectively the electric charge, the electron mass, the numerical grid size and the electron plasma frequency. These definitions are identical to those used in previous 1-D PIC (Lembège and Dawson, 1987), and 2-D PIC simulations (Lembège and Savoini, 1992). The plasma conditions and shock regime used herein are similar to those used in Lembège et al. (2009). All basic parameters are summarized as follows: the plasma simulation box has 6144 x 256 grids with a spatial resolution Δ=Δx=Δy=Δz=1/60 ( / pi ) = 1/3 ( / pe) (where / pi and / pe are the upstream ion and 490 electron inertial lengths respectively), which is high enough to involve all microstructures of the shock front. Initially, the number of particles per cell is 4 for each species. Velocity of light = 3, and mass ratio of proton and electron Mi/me = 400. To achieve reasonable run times and simulation domains, a ratio of ̃pe/̃ce = 2 had been used as in Hada et al. (2003) and Schöler and Matsukiyo (2004). The electron/ion temperature ratio is Te/Ti = 1.58; upstream ions and electrons are isotropic so that thi= = thi,x,y,z and the = the,x,y,z respectively. The ambient magnetic field is 0= 1.5. The shock has an averaged 495 Alfvenic Mach number, MA = shock/ A = 5.06 where the upstream Alfven velocity A is equal to 0.075 and the shock velocity shock= 0.38. The ratio of upstream plasma thermal pressure to magnetic field pressure is taken as i = 0.10 for protons and e = 0.16 for electrons. For these initial conditions, all other upstream plasma parameters are detailed in for both electrons and protons.
In present plasma conditions, cycles of the shock front self-reformation are evidenced like those of previous works 500 (Lembège and Dawson, 1987;Savoini, 1992, Lembège et al., 2009). An enlarged view of the time stackplot of the (y-averaged) main magnetic field component is shown over one cycle in Figure 10a. Note that the local direction of the front normal versus the upstream magnetostatic field is varying along a front ripple. Then, considering an average along the shock front (y-axis) is equivalent to averaging over these local normal directions and agrees with the method used in experimental data for determining the normal direction of the shock front crossed by the Cluster tetrahedral configuration for 505 small satellite inter-distance (i.e. average over local normal directions given individually by each satellite).
As shown in previous works, the self-reformation is mainly characterized by (i) a strong cyclic magnetic field amplitude variation at the overshoot with a time period ref (= 1523) of the order of 1/3 ci (where ci = 5027 is the upstream ion cyclotron period) so that several self-reformations can take place within one upstream ion cyclotron period; (ii) a temporal anticorrelation of the B field amplitude measured at the foot and at the overshoot respectively, and (iii) a strong time variability of the shock 510 front width. These features persist quite well even when using a realistic mass ratio as shown by al. (2013a, 2013b); in particular, the self-reformation cyclic period does not depend almost on the mass ratio. Moreover, the features of the other self-reformation process due to the MTSI (for slightly oblique shocks, see Sect. 1) also persists for realistic mass ratio as shown by Schöler and al. (2003). Herein, the spatial ramp width is identified within the shock front by a linear fit defined in the steepest part of the magnetic field gradient (Figure 10b). However, the procedure to define this width slightly differs from that used for Cluster experimental data in Sect. 2. Presently, its lower bound (dashed line B) is defined where the B profile departs clearly (due to the foot region) from this linear fit. On the other hand, its upper bound (dashed line C) is defined where the B profile departs clearly (due to the proximity of the overshoot) from this linear fit. We do not use the overshoot as the upper bound since it is polluted by 535 reflected ions as these gyrate and penetrate downstream after being once reflected at the ramp (Leroy et al., 1982); then this would overestimate the ramp width. On the other hand, let us note that using the mean downstream B field as a reference (as done for experimental data analysis in Sect. 3) may be a source of inaccuracy when applied to the present numerical simulation results. Indeed, the size of the downstream region is not long enough to deduce a precise estimate of the averaged downstream B value. In the B profile of figure 10b, these lower and upper bounds are respectively 540 https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.  (c) and (d) respectively. Note that the width of "ramp-1" only is represented in panel (c) in order to avoid overwhelming the plot; the smalltime fluctuations in ramp-1 correspond to error bars of the measurements (which are based on y-averaged ̃t z profiles). located at = 3755 and 3745 which provides a ramp thickness ramp=10. Repeating the same measurement procedure at 555 different times within the self-reformation cycle allows to evidence that the ramp thickness stays always very narrow around a (time averaged) value of ramp= 3 / pe (Figure 10c); the width of the ramp is almost independent of time and of the strong fluctuations of the maximum B field amplitude at the overshoot (not shown herein), which is mainly driven by reflected ions (Leroy et al., 1982). Since one has to compare statistical results of experimental observations with simulations, it is important to clarify the following points : (i) considering a perpendicular shock represents a reference case in the sense the ramp is 560 expected to be the thinnest for 90°, as compared with oblique cases (lower than 90°), indeed, whistler wave emission (precursor) in oblique shocks tends to extend the thickness; (ii) herein, the concerned simulation is based on a relatively low i value (= 0.10) or equivalently a high ratio shock/ thi (= 31.66), which fits with the range of experimental conditions where data sets of shock crossing have been chosen ( Figure 7d); (iii) over a full self-reformation cyclic time period, the shock front covers a distance of Δ = 597 (or equivalently 9.8 / pi). In the case study shown on Figure 5c, the satellites inter-distances along 565 the normal are respectively 1.4 c/pi for C4-C2, 1.6 c/pi for C2-C1, and 2.6 c/pi for C1-C3 which makes the total distance for the shock crossing of 5.6 c/pi ; (iv) shock curvature effects and front rippling of the shock front have been neglected herein, in order to approach the experimental conditions where the inter-distance between satellites is small. The y-averaging of the simulated planar shock is in agreement with the hypothesis of considering locally a planar shock front in experimental data for estimating the relative speed of each satellite versus the shock velocity. Then one can extract at present the width of the fine 570 structures (foot and ramp) of the shock front mainly along the shock normal with a good accuracy. Moreover, in contrast with the methodology used for experimental data (Sect. 2), it is not necessary to use herein two different levels of identification ('visual' and 'automatic') for defining the upper and lower bounds of each microstructure (ramp and foot), since profiles have been already enough smoothed out by the y-averaging. Let us remind that this self-reformation process reveals to be a very robust mechanism since it has been evidenced in PIC 575 and hybrid simulations (Hellinger et al., 2002(Hellinger et al., , 2007Lembège and al., 2009;Yuan et al., 2009). It was observed originally in 1D PIC (Biskamp et Welter, 1972;Lembège et Dawson, 1987), later in 2D PIC (Lembège et Savoini, 1992) and in 3D PIC simulations (Shinohara et al., 2007). As will be discussed in Sect. 5.1, this self-reformation process persists well even in the presence of microinstabilities developing within the front for perpendicular shocks (Muschietti and Lembège, 2006). The https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. variation of the ramp width versus time reported in Figure 10c shows that three successive time ranges can be defined where 580 respectively the width decreases (range (i) defined by < 3456 ), very slightly decreases (range (ii) defined by 3456 < < 4400), and increases (range (iii) defined by > 4400) as represented by pink, blue and yellow areas respectively. Then, the width stays almost constant during a large time range (ii). Within this range, the width reaches a lower bound value around 3 c/ωpe, (dashed-dotted red line) which can be used as a reference scale in the comparison with experimental data instead of comparing exact scales measured at a given crossing time of the shock. In order to analyze the whole-time variation of the 585 ramp width in detail, it can be associated with that of the foot and with ion phase space as described in Sect. 4.4.2.

Results on the foot thickness
On the same, Figure 10d shows the time variation of the foot width issued from the same PIC simulation. The upstream edge of the foot (vertical dashed line A in Figure 10b) is defined as the location where the magnetic field has increased by 6.67% over its upstream value (Burgess et al., 1989;Yang et al., 2009Yang et al., , 2011. The value 6.67 % has been determined as corresponding 590 to the maximum amplitude in the magnetic field fluctuations measured upstream (which is sensitive to the particles number per cell). Then, the amplitude of the upstream magnetic field turbulence is below the value (1 + 6.67%) x Bo in our simulations.
The time variation of the foot width can be detailed with the help of the three time intervals defined for the ramp thickness in Sect. 4.4.1, and is associated to the ramp as follows. In range (i), the foot width increases in anticorrelation with the ramp variation; at that time, reflected ions are accelerated at the ramp and only initialize their gyromotion (not shown herein). The 595 range (ii) (where the ramp width slightly decreases) includes three successive phases for the foot dynamics: one (ii-a) where the foot width is constant (3456 < < 3984), another (ii-b) where the foot width strongly increases (3984 < < 4176), and a last one (ii-c) where the foot width is still increasing, but in addition a secondary foot can be identified upstream (blue curve for > 4176); a double foot appears as detailed below. Finally, in range (iii) ( > 4416), the foot width reaches a maximum value (Lfoot / ρci,us) = 0.6 which remains constant at later times while the ramp thickness starts increasing. 600 A detailed analysis of the ion phase space (not shown here) at these corresponding time intervals allows to interpret more precisely the time variation of both ramp and foot thicknesses. In range (i) (i.e. < 3456), the ramp (named ramp-1) is freshly formed from the steepening of the foot resulting from the previous self-reformation and the foot in panel (d) is only due to the old reflected ions (previous self-reformation) still located upstream of the new ramp-1 while accomplishing their full gyration (before being transmitted downstream); this steepening corresponds to the decrease of the ramp width (range (i) 605 in panel (c)). In range (ii-a), the ramp-1 starts to accelerate and reflect new upstream ions, which add to old reflected ions which did not reach the downstream region yet. Then, the local foot (named foot-1) is composed with two (old and new) reflected ion populations (reflected at two successive reformation cycles). During this acceleration phase, the energy of macroscopic fields at the ramp is transferred to newly reflected ions, which smooths the local steepening and the ramp-1 width slightly decreases only. Within the range (ii-b), all old reflected ions have reached the downstream region; only freshly reflected 610 ions are present upstream of the ramp, continue their gyromotion and contribute to the foot-1. Consequently, the foot-1 width https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. strongly increases as mainly carried by these newly reflected ions. Correspondingly, the width of the ramp-1 still slightly decreases. But, more and more new ions are reflected and accumulate in time, and a local new foot (named foot-2 ) builds up.
This ion accumulation is large enough so that the upstream edge of this new foot-2 steepens and becomes a new ramp (named ramp-2, not shown in panel (d)) which starts reflecting new ions; let us mention that the width of "ramp-2" cannot be 615 precisely identified before < 4368 by using the best linear fit (as in Figure 10b) since this part of the front microstructure is still mainly dominated by the growing foot (steepening is not strong enough in the upstream edge of "foot-2" within this time range). This new foot-2 (blue curve in panel d) co-exists with the old foot-1 (red curve) within the same time interval > 4176. Simultaneously, the amplitude of the new overshoot (not shown herein) associated to the new ramp-2 increases and the field gradient between the old and the new overshoot becomes weaker. Therefore, the ramp-1 width increases at later time 620 (time range (iii)). For > 4608, the old ramp-1 is totally absorbed in the shock profile (between the old and the new overshoots), and its width cannot be identified precisely.
Then, the shock front may acquire a multi-peak profile including traces of old and new successive ramps (not to be confused with emission of large amplitude whistler waves emitted from the main ramp). The width of the whole shock front is the summation of the ramp and the feet contribution and can be estimated from the upper edge of the new foot-2 to the 625 lower edge of the old ramp in the front. Present results provide maximum value measured at the final stage of self-reformation cycle ( shock,max= 176.3 = 2.93 / pi at = 4608). Let us stress that this value represents an upper bound of the shock front width and not a common measurement of the front width itself.

Comparison with some previous works
This section is dedicated to comparing the present results with those issued from previous works. The main emphasis is given 630 not only on the agreement/disagreement observed between the different results, but also -and even more importantly-on the strategy used in previous works for "extracting" the spatial width of the shock front microstructures. We will distinguish previous works dedicated to measurements of the ramp and of the foot thicknesses respectively.

Spatial measurements of the ramp thickness
Newburry and Russell (1996) 635 Most works based on the analysis of experimental data issued from ISEE 1 & 2 satellites indicate a scaling of the ramp width within the ion inertial length. In a very detailed analysis of a single shock crossing, Scudder et al. (1986) have reported the ramp thickness to be 0.2 c/pi , a measurement which was based on an exponential fit of the foot-ramp transition. By using the same fit technic, Russell and Greenstadt (1979) have measured a slightly thicker ramp of about 0.4 c/pi in an initial study of ISEE bow shock observations. Later on, Newburry and Russell (1996) has analysed seven nearly perpendicular shocks (Bn 640 >80°) and found a range of ramp widths (based on a linear fit to the ramp transition) from 0.5-0.8 c/pi. However, in one case https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. the ramp width was particularly thin 0.05 c/pi (= 2 c/pe). To the knowledge of the authors, it was the first time that a so thin ramp has been observed associated to the terrestrial shock. All the observations lead to the following statements: a) In contrast with the tentative explanations proposed at that time (summarized in Newburry et al., (1998)), this very thin ramp (2 c/pe) takes a full meaning herein and is in quite good agreement with the present statistical results. It represents 645 quite well a potential signature of the self-reformation process which allows the ramp to access to very small thickness (a few c/pe). b) one question raises up: why most previous works evidence mainly a ramp width scaling with a noticeable part of the ion inertial length (except the very narrow ramp mentioned above) in contrast with the present results? Several answers can be proposed: First, the method commonly used in previous works and based on exponential fitting to the foot-ramp region is not 650 appropriate to measure carefully the ramp and foot thickness respectively, since both fine structures were mixed together. This fitting methodology is a source of important errors. Second, the solar wind conditions corresponding to the set of shock crossings chosen at that time need to be verified carefully. Indeed, as i is relatively low (i = 0.17 as in Hada et al. (2003), or equivalently vsh/vthi equals several decades as in Scholer et al., (2003)), the reflected ions describe a focussed (narrow) beam during their gyration These maintain a well coherent gyromotion which forces them to accumulate at a certain distance from 655 the ramp and to initialize the foot. The growth of the foot separates clearly from the ramp and is characterized by a trapping loop in the ion phase space; self reformation takes place and the front is nonstationary.
In contrast, as i is relatively high (i = 0.35 as in Hada et al., 2003, or vsh/vthi equals a few units as in Scholer et al., 2003), the reflected ions rapidly diffuse as these start their gyromotion; these lose their large coherent motion and spread out within the whole front (no clear ion trapping loop). The foot grows up but does not separate from the ramp. Then, ramp and 660 foot are mixed, and a larger ramp thickness is expected as shown in PIC simulations of Lembège et al. (2013aLembège et al. ( , 2013b. Then, no self-reformation can take place and the shock front is stationary.

Newburry et al. (1998)
At the ISEE 1 & 2 time, the ramp was usually characterized by two different scales: (i) a large scale (or global ramp width) 665 within which the main transition from upstream to downstream magnetic fields takes place, and (ii) a thinner subramp scale which contains steep jumps in the magnetic field with amplitudes sometimes comparable to the overall change in the magnetic field at the shock. One conclusion (Newburry et al., 1998), was that both scales are characteristic of the quasi-stationary shock profile (which -at that time-was considered as stationary within the ion gyroperiod). More precisely, Newburry et al. (1998) have analysed a set of 20 shock crossings and found that the width of the overall ramp transition is typically within the range 670 0.5-1.5 c/pi which is quite large as compared with our present measurements. This difference can be easily explained by the fact that the authors have defined the ramp thickness based on a linear fit of the total magnetic field profile from the foot to the overshoot, which enlarges the estimate of the ramp thickness, as compared with the present method used in Sect. 3 and 4. https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.

28
Cluster data and present statistical results allow to correct this view within the frame of the intrinsic shock front selfreformation. The so called "large scale" of the ramp presents a large width variability (extending from a part of the ion inertia 675 length to a few electron inertia lengths as in Fig. 8a), which was not observed at that time for the reason mentioned above. But, simultaneously, the so called "small scales" can result from instability developing with the foot region as the ECDI which holds for "Electron Cyclotron Drift Instability" (Muschietti and Lembège, 2006) or MTSI Matsukyio and Schöler, 2003;Muschietti and Lembège, 2017) for instance, and can propagate (convection effects). In present work, the use of a linear fit applied within the ramp region only (as described in Sect. 3) allows to average these local small scale 680 fluctuations and to measure the ramp width with minimizing errors. Analysing small scale field fluctuations (as those also observed within the ramp in Figure 2) requires a further investigation which will be presented elsewhere.

Bale et al. 2003
From the spacecraft floating potential measured on board of the four Cluster satellites, Bale et al. (2003) have determined the 685 electron plasma density, and studied the macroscopic density transition scale at 98 crossings of the quasiperpendicular terrestrial bow shock. This first tentative for measuring the shock transition width from Cluster data was attractive but contains several limitations for the following reason: an hyperbolic tangent function has been fitted to each density transition in order to capture the main shock transition but no microstructures of the front -as foot, ramp, overshoot-have been identified and scaled. Such a fit partially includes parts of the shock front outside the ramp itself so that it cannot apply to determine the ramp 690 thickness with acceptable accuracy since both ramp and foot structures are not separated. This technique takes into account only the widest transition scale at the shock, which overestimates the ramp width and restricts drastically its application. Then, the conclusions stating that for high Mach number shocks (including the range MA=4-5 as that considered in the present work), only the convected gyroradius is the preferred scale for the shock density transition are strongly questionable since the real ramp is much thinner. This resulting scale corresponds to the scale of the foot rather than that of the ramp and no conclusive 695 statements can be obtained.
On the other hand, in the light of Cluster data, looking for an estimate of the shock transition presents nowadays a limited interest, since recent results show that terrestrial shock is rather nonstationary and a "unique" typical spatial scaling must be replaced by some interval or bounds of variation (lower bound/upper bound) within which the spatial scales of the fine structures can vary. 700 Figure 11 helps to better compare our results with those by Bale et al. (2003) which have been reproduced on Figure   11a, 11b and 11d. The same representation as on panels 11b and 11d has been used for panels 11c and 11e to plot the front thickness computed as the sum of the ramp and the foot thicknesses from our methodology versus the magnetosonic Mach number Mms. The pink area on panel 11b illustrates the range of Mms explored from the present study. Figure 11c reveals that in this range of Mms values no clear dependence considering the large error bars (as on above panel 11b) can be identified from 705 the present study. It would be obviously impossible here to conclude that the front thickness simply increases with Mms and conversely the best fit would conclude here that there is a trend to decrease. Our results are therefore in strong contradiction https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.

29
with conclusions of Bale et al. (2003). Moreover, if we consider another Mms range of similar extension as the one shown on panel 11c from our study and highlighted by the green area on panel 11b, a similar representation (not shown here) would conclude to exactly the opposite to what the dotted line on panel 11b shows, considering again the large error bars. This 710 illustrates again the impossibility to make a simple scaling of the front thickness by a single parameter as the magnetosonic Mach number. Panel 11e shows that the front thickness normalized by the downstream gyroradius is more or less consistent with the results of Bale et al. (2003) shown on Panel 11d. There is a trend there to decrease with Mms increase. As mentioned in the present study, the use what of the so-called 'density transition scale' by Bale et al. (2003) shown in Figure 11a is inappropriate since it mixes the time variability of the ramp and the ion foot. They also mention a stationary shock from their 715  Bale et al. (2003) for the shock crossing on 2000-12-15/01:26:15, while panels (c) and (e) are similar to panels (b) and (d) but correspond to results from the present analysis. Panel (c) has to be compared with an enlarged view of panel (b) within the magnetosonic Mach number range 2 < Mms < 6.5. In panel (a), the transition of the density from upstream (unshocked) to downstream (shocked) states for a magnetosonic Mach number Mms=3.5, Bn = 81° shock (see text), where Mms = Vsw,sc/cms with cms= (cs 2 + vA 2 ) 1/2 , cs is the sound speed, and Vsw,sc is the bulk velocity of solar wind measured in the frame of the spacecraft; 725 the green line is the hyperbolic tangent fit, and red vertical lines show the density transition scale. In panels (b) and (d), pink and green areas are superimposed to results of Bale et al. (2003) to indicate respectively the same Mms range as for our study and another Mms range (same width) for higher values for comparison. Panel (c) issued from present results shows the variation of the shock front thickness L(ramp+ foot) versus the Mms range corresponding to the pink area defined in panel (b) where L is normalized to the ion inertial length ⁄ defined upstream. Panel (e) issued from present results shows the variation of the shock front thickness L(ramp+ foot) versus Mms for the same range,

730
where L is normalized to (vsh,n /ci,2) in order to compare with panel (d), where vsh,n is the shock normal velocity defined in the frame of the solar wind and ci,2 is the proton gyrofrequency defined downstream.
https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. Krasnoselskikh et al. (2002) have proposed another mechanism as a possible source of nonstationarity of the shock front (so 735 called "gradient catastrophe"). In this mechanism, the transition to nonstationarity is characterized by the lack of the phase standing whistler wave train inside the shock front. This process takes place for oblique quasi-perpendicular shock (for a direction Bn not too far from 90°) as the Alfven Mach MA is above a certain critical Mach number given by Mnw = (1 / 2) ½ (Mi/ me) ½ cosBn; this condition can be equivalently expressed for a fixed MA as the direction Bn is below a critical angular value.

Krasnoselskikh et al. (2002)
In this case, the nonlinear steepening of the wave train cannot be balanced by the dispersion/dissipation effects. However, as 740 mentioned in Sect. 1, this theoretical mechanism has strong limitations since (i) it is based on 1D model only and (ii) it excludes any dissipative effect. Indeed, reflected ions (which play a key role in supercritical Mach regime) are neglected and their impact on the shock dynamics -and in particular on the temporal fluctuations of the shock front micro-structures which is our center of interest in the present paper-is totally absent. Because of the above strong limitations, this process is still a source of controversy. 745

Lobzin et al (2007)
Nevertheless, a further study made by Lobzin et al., (2007) has proposed the gradient catastrophe process as being at the origin of the shock front nonstationarity observed by Cluster on 24 January 2001. While the front nonstationarity is quite clear in the experimental data which show quite different shock front profiles of the magnetic field crossed by the 4 satellites of Cluster, 750 the proposed interpretation stays quite questionable for the following reasons: (i) the overall analysis is based on one shock crossing only (no statistical analysis has been performed), (ii) no procedure detailing the transfer from time to spatial scaling of the shock front (in particular focused separately on the ramp/foot structures) has been indicated so one ignores the percentages of errors underlying the spatial field profiles measurements which are not shown; (iii) no mapping and/or basic information have been provided on the 3D disposal and inter-distance of the satellites, so that one ignores whether the four 755 satellites cross the same shock front within the same or different successive self-reformation cycles; (iv) only one measurement of the ion density showing some fluctuations has been presented which is not enough for confirming that it could be due to the gradient catastrophe process. Moreover, note that such similar density fluctuations have been also measured quantitatively versus time in previous numerical simulations ( Lembège and Savoini, 1992;Yang et al., 2009) but have been attributed as being due to the self-reformation process fed by the accumulation of reflected ions and not to the gradient catastrophe process. 760 Similar self-reformations have been also obtained as supported by the MTSI for an oblique propagating shock (Comisel et al., 2011); (v) a deeper analysis shows that the upstream ion temperature used in this paper is not the correct one since the CIS instrument was in magnetospheric mode (see Sect. 2) so that the moment calculations provide a large overestimate of the solar wind proton core temperature. Indeed, the determination of the upstream ion thermal velocity requires a delicate procedure and is of central interest since the stationarity/nonstationarity of the shock front is strongly dependent on the ratio vsh/vthi When using the solar proton temperature from the ACE satellite (4.5 eV thermal energy), the ion beta found is 0.6 instead of the value of 2.0 reported by Lobzin et al. (2007).

Comisel et al (2011) 770
The same data of Cluster measured on 24 January 2001 has been re-examined in detail by Comisel et al. (2011) and reinterpreted quite differently by using PIC simulations in plasma conditions and Mach number regimes chosen to be appropriate with experimental conditions. Basically, the authors have shown that the self-reformation is quite well recovered but is not controlled by the gradient catastrophe process (dispersive effects) mentioned above but rather by the self-reformation mechanism fed up by the MTSI which is an indirect consequence of the ions reflection (dissipative effects which have been 775 neglected in Krasnoselskikh et al. (2002)). Indeed, for slightly oblique shocks (out of 90°), the MTSI develops between the incoming electrons and incoming ions, with a growth rate large enough to develop during the gyration of the reflected ions (Schöler and Matsukiyo, 2004). These authors have shown that for a nonrealistic mass ratio, the self-reformation is still controlled by the collection of reflected ions at the upstream edge of the foot as in Lembège and Savoini (1992), since the MTSI has a weak growth rate and has no impact on the nonstationary dynamics of the shock front. In contrast, for realistic 780 mass ratio, the growth rate is much larger and MTSI can be easily excited which leads to ion phase mixing and thermalization.
Then, the ion pressure increases considerably locally somewhere within the foot itself (and not at the edge of the foot) with associated building up and overtaking of the magnetic field. A new shock ramp builds up within the foot which initiates to reflect new incoming ions and a cyclic reformation is initialized. The cyclic period is lower than the one associated with the self-reformation due to the accumulation of reflected ions (as in Lembège et Savoini, (1992)), since the process takes place 785 before ions accomplish their full gyration. More precisely, by using the appropriate conditions of the 24th January 2001 shock crossing in the PIC simulations, Comisel et al. (2011) have shown that the waves triggered by this instability is of whistler type, have phase velocities directed downstream in the shock frame, and the associated Poynting flux (and wave group velocity) is also directed in the downstream direction. However, as these approach the density and magnetic field overshoot, the waves are refracted; as a consequence, the Poynting flux changes direction and is directed upstream in the shock frame. In summary, 790 there is no indication of a phase standing linear or nonlinear whistler precursor produced by the shock ramp itself in this process. This scenario does not support that proposed by Krasnoselskikh et al. (2002), where nonlinear whistler precursor is required to initialize the "gradient catastrophe" process as a source of the front nonstationarity. Hobara et al., 2010 795 A statistical study based on 77 crossings of the terrestrial bow shock (30 by THEMIS and 47 by Cluster) have been performed by Hobara et al. (2010) in order to analyze the ramp scale and the associated gradient width. Their study claims that the spatial scale of the ramp decreases as the shock Alfven Mach number increases. However, detailed comparison with our results is difficult (almost impossible) at present time because of important missing information: (i) No detailed information is provided on the procedure used for extracting the spatial widths of the ramp from the temporal experimental data as detailed in present Sect. 2; (ii) the duration of the ramp crossing has been defined as a time interval between the upstream edge of the ramp and the maximum of overshoot. However, let us note that the upstream edge is not defined precisely (see the caution to account for in our Section 3) and the use of the overshoot can lead to an overestimate of the ramp width as compared to the method used herein (Sect. 3) since it is partially polluted by reflected ions; (iii) statistical errors are indicated in their Figures 2 and 3 only for range of Mach numbers and not on individual measurements; (iv) no statistics are shown versus the direction Bn of the shock front normal; as a consequence, one ignores how/whether the measurements of the ramp /foot bounds are polluted by the whistler (either linear either nonlinear) whose emission strongly depends on the angle Bn; (v) no statistics/information are indicated on the ion i value (which has a strong impact on the shock front self-reformation) ; (vi) the authors use two different statistics: one based on the ramp spatial scale (their Figure 2), the other based on magnetic ramp spatial gradient scale (their Figure 3). Moreover, they used a simple theoretical argument (based on dispersive whistler) to 810 stress that the wavelength Lw of the linear/nonlinear whistler varies as Lw ~ 1 / MA. But, the fact that both statistics show that the width of the measured structure decreases as MA increases does not imply that dispersion determines the size of the magnetic ramp even for supercritical shocks (as discussed later on); (vii) at least, only a very limited number of measurements (whatever the MA value is within the concerned range MA=2-12 in Hobara et al., 2010) indicates that the ramp width succeeds to reach a few electron inertial lengths (equal or lower than 5 / pe, in their Figures 2 and 3) which is in strong contrast with our 815 statistical results of Figure 8a. The lack on statistics versus Bn and i (above points (iv) and (v)) and the absence of precise information on the methodology used in the statistics of Hobara et al. (2010) (as listed above) could explain this difference.
Moreover, Hobara et al. (2010) deduce that their measurements (in particular the decrease of the ramp thickness as MA increases) evidence that the major process responsible for the shock front nonstationarity is the "gradient catastrophe" process proposed by Krassnoselskikh et al. (2002), andLobzin et al. (2007). However, three main facts contradict this statement: (i) 820 arguments describing the limitations of the theoretical model (Krassnoselskikh et al., 2002), and the revised analysis of Comisel et al. (2011) as already described above in present section, (ii) Hobara et al. (2010) use the overshoot (observed in the set of their experimental measurement) as an upper bound to define the ramp width; but the overshoot of supercritical shocks is a clear signature of the reflected ions as these gyrate downstream after being once reflected at the ramp (Leroy et al., 1982).
Then, a contradiction occurs since dissipation effects (strongly carried by the reflected ions) have been fully neglected in the 825 model of "gradient catastrophe" which cannot support their statements; (iii) a similar variation of the ramp thickness as MA increases has been also observed in the self-reformation process in PIC simulations (not shown here), which means that this variation cannot be used as a unique signature of the "gradient catastrophe" process. At least, in order to support their statements, Hobara et al. (2010) remind that PIC simulations evidencing the self-reformation due to the accumulation of reflected ions at a foot distance from the ramp (as in Lembège and Savoini, (1992)) as well as due to the MTSI (as in Matsukyio 830 and Schöler, (2006) use low pe/ce ratio, which overestimates the electric field to magnetic field ratio (in practice, using a realistic pe/ce ratio is very difficult at present time for computational constraints). If were true, this argument would suggest that the electric field would be too large and that the self-reformation process would result from an artefact. However, no https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. quantitative result has been obtained until now to support such a statement. In contrast, let us remind that (i) the parametric analysis performed by Lembège et al. (2009) clearly evidences that both hybrid and PIC simulation retrieve the same self-835 reformation process of the shock front in similar plasma conditions and Mach regime, provided that the space grid is small enough (i.e. high resolution) in hybrid simulations as shown initially by Hellinger et al. (2002); this last result would suggest that the field gradients at the ramp are important and not the absolute fields amplitude; (ii) electron scales are neglected in common hybrid simulations (as those used in Hellinger et al., (2007); Lembège et al., (2009) and references therein) which are independent of the ratio pe/ce. This suggests that the ratio pe/ce has no direct impact on the self-reformation. Note, that the 840 ramp width may be estimated in terms of percentages of ion inertial length in hybrid simulations and both in terms of electron and ion inertial lengths in PIC simulations.
In summary, the measurements of Hobara et al. (2010) do not support the model describing the ramp as an evolving nonlinear whistler wave, neither the mechanism of "gradient catastrophe" process of Krassnoselskikh et al. (2002) as a source of the front nonstationarity. 845 Schwartz et al. (2011) Schwartz et al. (2011 have determined the scale of the electron temperature gradient via electron distributions measured in situ by the Cluster spacecraft and have identified two different scales; the Mach number is supercritical (MA=3.8) and the shock propagation is oblique (Bn =83°). First, the authors found that half of the electron heating coincides with a narrow layer 850 covering several electron inertia lengths (c/pe) and found an inflated electron distribution due to primarily the electron filling in/or the entrapment in regions of phase space that would be not accessible. Second, the remaining of the solar wind maximum is only present in the smoother increase of the field which precedes the ramp; this beam is totally destroyed by the time this electron scale ramp is crossed. The authors found that this scale is comparable to that deduced from the theoretical limiting case of a wave capable of phase standing in the incident flow and concluded that this suggests that supercritical shocks steepen 855 to this whistler limit since dissipation processes are insufficient to broaden the transition further.
However, several elements of information are missing: (i) the authors mentioned to have used a technique to convert the time series of data to distance along the shock normal similar to that of Schwartz (1998) i.e. the procedure to determine the normal and the velocity of the shock from multi-spacecraft measurements (as in Sect. 2) but did not provide detailed information on the procedure itself, associated errors estimate, or on the definition of the ramp itself (their Figure 1); (ii) the 860 results are limited to one shock crossing (one spacecraft only); some attempts to other shock crossings are mentioned but no statistical results have been shown or summarized; (iii) the observed results cannot lead to the conclusions that dissipation processes are insufficient to broaden the transition further, since once again the key dissipation processes driven by the (reflected) ions are not mentioned /analysed in this study. Of course, dispersive effects must be considered in the global balance of nonlinear steepening of the front by both dispersive/dissipative effects, but this does not mean that the shock front dynamics 865 is controlled by dispersive effects. Moreover, other sources of fields fluctuations may take place in the shock front (such as https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.
formation of multi-peaks (signatures of old/new self-reformation (Sec. 4.4.2), and/or whistler waves excited by the MTSI within the front) without invoking any dispersive effect. Walker et al. (2004) investigated short scale structures in the electric field that are observed during crossings of the quasi-870 perpendicular bow shock using also data from the Cluster satellites. An example is reproduced from this paper on Figure 12a.

Spatial measurements of the substructures
The structures observed at each Cluster satellite magnetic front exhibit large amplitudes, as high as 70mV/m and the authors argue they make a significant contribution to the overall change in potential at the shock front. They have shown that the scale size of these short-lived electric field structures is of the order of a few c/pe as shown on their statistics reproduced in Figure   12c. They also studied the relationships between the scale size and the upstream Mach number and Bn. They found that the 875 scale size of these structures decreases with increasing plasma i and as Bn ~90°. magnetic field magnitude; following panels display the amplitude of the electric field and Ex and Ey GSE components, respectively). Panel (b) reproduces the distribution of magnetic field ramp obtained from the present study already shown in Figure 8a to compare with the distribution of the electric field structures of Walker et al. (2004) shown on panel (c).
Interestingly, the comparison of distributions displayed on Figure 12b (present analysis) and Figure 12c (Walker et al., 2004) 885 shows that the electric field structures are consistent with substructures of the ramp since the typical scales found seem smaller.
Of course, to be fully conclusive, one would need to compare the size of the electric field structure and of the magnetic ramp https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.
for each shock crossing. However, both results seem very consistent. Recently, Dimmock et al. (2019) has evidenced electronscale field structures inside the shock ramp from a case study using the observations of only two Cluster spacecraft which were very close to each other and proposed that it is associated to the transition to nonstationarity. 890

Spatial measurements of the foot thickness
As compared with studies focused on the ramp width, previous works dedicated to the spatial measurements of the foot thickness are in relatively limited number. First observations of the foot have been made by Paul et al. (1965). Woods (1969; has proposed the reflected ions to account the foot structure and has estimated a foot width for a strictly perpendicular shock from the turning point distance (Eq. 4 of Sect. 4.3 which is recalled here): 895 assuming a specular reflection of ions, where ρci,sw is the upstream ion gyradius estimated from the bulk velocity of the solar wind (and not from the ion thermal velocity); the reflection is specular in the sense that the incident ion's velocity component along the shock normal is reversed at the shock, while the component perpendicular to the shock normal remains unchanged.
A good agreement has been found with ISEE 1 & 2 data by Paschman et al., (1982), and Sckopke et al., (1983). However, two 900 points need to be stressed out: (i) the measurement has been made from the ramp to the turning out point of reflected ions for five ISEE 1 & 2 crossings of the shock (the thickness of the foot should be comparable to this turnaround distance), (ii) the condition of specular reflection is not simple when applied to oblique shocks, as shown by Livesey et al. (1984)  where Vn is the acute angle defined between the incoming solar wind direction and the normal to the shock front. A good agreement has been obtained between the experimental measurements (distance from the ramp and to the turning point of reflected ions for several shock crossings made by ISEE 1 & 2) and values deduced from Eq. (1). Moreover, Moses et al. (1985) found that the shock velocity vsh in the frame of the spacecraft can be expressed as: where x = 0.68 ρci,us cos (Vn) sin 2 (Bn)/ Ωci Δt (5) "-" sign refers to a transition from upstream to downstream, the "+" sign refers to the inverted transition, Ωci is the upstream ion gyrofrequency and Δt is the foot traversal time. Later on, Gosling and Thomsen (1985) have derived a more general 915 expression for the turnaround distance dturning for arbitrary Bn which differs from (3) particularly for Bn < 60°, and which incorporates above Eq. (1) and (2) of Sect. 4.3. Initial steps of their derivation of dturning are closely similar to the formalism initially developed by Schwartz et al. (1983) and leads to (Gosling and Robson, 1985): https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.
Since present statistics concern quasi-perpendicular shock with shock normal direction range 74° < Bn < 90° (Figures 8 et 9 (Lembège and Savoini, 1992;Yang et al., 2009;Comisel et al., 2011); similarly, bursts of electron (with parallel kinetic energy) are formed for oblique shocks as long as the self-reformation persists for oblique shock until a certain critical angle ; (ii) the shock front variability has also a strong impact on the nature of the acceleration process SSA (shock surfing acceleration) versus SDA (shock drift acceleration) during their 935 reflection as shown by Yang et al. (2009).
In the same way, a super-upper bound value of the foot thickness can be obtained from the simplified Eq. 3 defined for a strictly perpendicular shock. The deduced value Lfoot/ρci,us = 0.68 (herein the convection bulk velocity of the solar wind is equal to the shock velocity since we are in the rest solar wind reference frame) has been reported in simulation results of Figure 10d (perpendicular shock). Eq. 3 provides a very good upper limit of the time variation of the foot thickness. 940 In summary, since the shock front is highly nonstationary, no fixed theoretical scale can be provided as a reference in contrast with statements of many previous works. Present results confirm that considering the whole shock front to determine the transition layer width can be a source of errors and misunderstanding. Instead, (i) a detailed and separate analysis of microstructures (ramp, foot -and later overshoot-) is necessary, and (ii) a comparison with a lower bound value (for the ramp width) and upper bound value (for the foot width) reveals to have a full meaning. 945 6 Conclusions 6.1 About the methodology used for extracting the spatial width of each fine structure: i) this works presents in detail the different steps necessary to determine carefully and without ambiguity the spatial thicknesses of each fine structure from the magnetic field measurements which are summarized as follows: https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License.

37
-First, an upstream time interval is defined in the time series for each of the four spacecraft and an average magnetic 950 field is determined. If all average fields agree in direction and magnitude within predefined limits, an unperturbed upstream field B0 is obtained for the global shock crossing by the four satellites; -Second, for each spacecraft, the entry times in the magnetic foot are automatically determined when the values of the components and magnitude exceed the upstream values by a certain value defined from the upstream standard deviations; then one validates by visual inspection that results are not polluted by any upstream disturbance or detached wave packet. 955 -Third, a downstream interval is defined to measure the downstream asymptotic value of the magnetic B field magnitude.
-Fourth, for each spacecraft, an automatic procedure is applied on the magnitude of the magnetic field to determine the entry in an exit times from the magnetic ramp. This procedure is initialized by preliminary estimates using the upstream and downstream asymptotic levels and based on a linear regression used only to better define the time interval where the 960 gradient is the strongest with the best correlation coefficient; -Fifth, the exit time from the first overshoot is estimated first by the magnetic field magnitude reaching the downstream asymptotic value for the closest time to the first undershoot and then adjusted by allowing a reasonable variation (for instance 5%) of the downstream asymptotic value (see Sect. 2, Fig. 2a); -Sixth, the interval variations of the entry/exit times for each microstructure are used to define error bars on them; 965 -Seventh, the values of both the normal to the shock and its velocity in each spacecraft frame is determined using the four-spacecraft timing method. The obtained normal vector n0 is then used to determine the general geometry of the shock characterized by the angle with the unperturbed upstream field vector B0 determined above. The error bars on comes from both the determination of B0 and n0 from the timing errors; -Eighth, the spatial profile along the normal direction n0 is derived for each spacecraft in units of upstream ion inertial 970 length and the microstructure thicknesses are derived in terms of appropriate physical lengths.
ii) the fitting technics based on exponential or hyperbolic function and used in previous works to match to match a shock front profile is questionable since these lead to mixing spatial scales of the foot and the ramp. Rather, our statistical results state that a linear best fit with appropriate identification of the upper and lower bounds of the temporal time ranges of the ramp and the foot is much more accurate and allows a precise error estimate on the spatial width of each microstructure. 975

About comparing experimental with numerical simulation/theoretical results:
Simulation of a strictly perpendicular case as presented in Sect. 4 can be considered herein as a reference case allowing to provide some bounds range within which statistics of experimental results may be inserted. More precisely, these provide the following information: the time variation of the ramp thickness is relatively limited (almost constant) and one can define a 'lower limit' (a few electron inertial lengths) for the ramp width. In contrast, the foot thickness is strongly varying in time and 980 one can define an 'upper limit' provided by Eq. 2 (which simplifies into Eq. 3 for a strictly perpendicular shock). Moreover, statistical analysis of measurements appear to be much more essential than the analysis of a single shock crossing only, in order to gather relevant information on the widths of the front microstructures. Since no shock crossing will be exactly identical one to each other even after a deep selection of data, we have to collect several shock crossings in slightly different conditions in terms of Mach number regime, of shock normal direction, of i plasma conditions and of orientation of 985 the shock crossing by the 4 satellites set. In addition, each shock crossing reveals some variability of fine structures thickness within the front. Then, no one to one correspondence in comparing numerical and experimental data for each crossing is possible; instead, the use of 'lower' and 'upper' limits mentioned above reveals to have a full meaning.

About the possible identification of the process responsible for the shock front nonstationary:
Several points need to be emphasized: 990 (i) All shocks profiles analysed herein exhibit a foot and overshoot which are clear signatures (in supercritical regime) of the gyration of reflected ions upstream of the ramp, and of these same reflected ions which succeed to penetrate downstream after being once reflected at the ramp. Then, dissipation effects carried by these reflected ions need to be fully considered.
Among the different mechanisms proposed to account for the shock front nonstationary, only two processes named 'selfreformation' include self consistently (without any approximation) the competition between the nonlinear, dissipative and 995 dispersive effects, and in particular the important role of reflected ions on the shock dynamics (and vice versa). Up to now, one still ignores whether the self-reformation is mainly controlled by the accumulation of reflected ions which is observed for perpendicular and quasi-perpendicular shocks (Lembège et Savoini, 1992) or by the MTSI-1/MTSI-2 (Schöler andMatsukyio, 2004: Muschietti andLembège , 2017) which is triggered only for oblique shocks (out of 90°); MTSI-1 and MTSI-2 refer to the instability excited in the foot region by the relative drift between the reflected ions/incoming electrons and the incoming 1000 ions/incoming electrons respectively. This open question will require a further analysis which is under investigation.
(ii) Among others, one typical signature of a front self-reformation is the anticorrelation between the time variation of the 'new' growing foot and the 'old' overshoot as illustrated in Figure 10a. Indeed, as the amplitude of the new foot increases, new incoming ions start to be reflected, which has an impact in the need for local dissipation at the overshoot located just behind. Then, a clear identification of a given self-reformation process would require an inter-satellites distance such that the 1005 four satellites cross the same shock front within its same reformation cycle. This requires a short inter-satellites distance for the quartet. This seems more easily addressed with the Magnetospheric MultiScale (MMS) mission data though the spacecraft separation is so small that low-scale fluctuations of the shock front (ripples) has a different impact on the normal determination.
These low-scale ripples have been identified in a single-case study by Johlander et al. (2016) and also indicated from the observed variations of the cross-shock potential in another single-case study by Hanson et al. (2019). 1010 (iii) We emphasize the importance of dimensionality effects in simulations studies for identifying what mechanisms can be dominant in the shock front nonstationarity. For a strictly perpendicular propagating shock, both 2D hybrid and PIC simulations have clearly evidenced the emission of nonlinear whistler waves which stay stationary with the ramp, and propagate obliquely with respect to both the shock normal and the upstream magnetostatic field Bo  https://doi.org/10.5194/angeo-2020-54 Preprint. Discussion started: 25 September 2020 c Author(s) 2020. CC BY 4.0 License. Hellinger et al., 2007). This emission holds only as the Bo field is laying within the simulation plane (and self-reformation is 1015 absent) and disappears as Bo is out of the simulation plane (and self-reformation is present). The source of this whistler emission has not been clearly confirmed yet; ion Weibel instability has been recently proposed as a possible mechanism (Burgess et al., 2016). In other words, the emission of these whistlers seems to be associated with the disappearance of any self-reformation process. But, this enigmatic observation reveals that 3D PIC simulation is necessary in which the Bo field is fully involved whatever its orientation is in the simulation box. Preliminary results of Shinohara et al. (2007) based on 3D PIC 1020 simulations show that both processes (emission of nonlinear whistler from the ramp and front self-reformation) can co-exist.
In other words, the time variation of the ramp/foot thickness associated to the self-reformation should persist in 3D simulation; a further quantitative study is also required in order to check the scales of the front microstructures.