First observations of irregular surface of interplanetary shocks at ion scales by Cluster

We present the first observational evidence of the irregular surface of interplanetary (IP) shocks by using multi-spacecraft observations of the Cluster mission. In total we discuss observations of four IP shocks that exhibit moderate Alfv\'enic Mach numbers (M$_A\leq$6.5). Three of them are high-$\beta$ shocks with upstream $\beta$ = 2.2--3.7. During the times when these shocks were observed, the Cluster spacecraft formed constellations with inter-spacecraft separations ranging from less than one upstream ion inertial length (d$_i$) up to 100~d$_i$. Expressed in kilometers, the distances ranged between 38~km and $\sim$10$^4$~km. We show that magnetic field profiles and the local shock normals of observed shocks are very similar when the spacecraft are of the order of one d$_i$ apart, but are strikingly different when the distances increase to ten or more d$_i$. We interpret these differences to be due to the irregular surface of IP shocks and discuss possible causes for such irregularity. We strengthen our interpretation by comparing observed shock profiles with profiles of simulated shocks. The latter had similar characteristics (M$_A$, $\theta_{BN}$, upstream ion $\beta$) as observed shocks and the profiles were obtained at separations across the simulation domain equivalent to the Cluster inter-spacecraft distances.


Introduction
Interplanetary (IP) shocks are ubiquitous in the heliosphere and are mostly associated with interplanetary coronal mass ejections [ICME, e.g. Sheeley et al., 1985] and stream interaction regions [SIR, e.g. Gosling & Pizzo, 1999]. Since the solar wind (SW) is a collisionless plasma, the heliospheric shocks are also collisionless. Marshall [1955] realized that purely resistive dissipation mechanisms cannot sustain collisionless shocks when their Mach number exceeds some critical value M c . Sagdeev [1966] suggested that particle acceleration and reflection at the shock surface could be an efficient mechanism to dissipate the kinetic energy of the incoming SW.
Numerical simulations show that one of the characteristics of supercritical collisionless shocks is the nonstationarity or reformation of their surface which leads to its irregularity. Depending on the shock properties (Alfvénic Mach number M A , the angle θ BN between the upstream B-field and the shock normal and the ratio β of the upstream thermal to magnetic pressure), the nonstationarity may be due to the self-reformation of the shock surface [e.g., Biskamp & Welter, 1972, Leroy et al., 1982, the whistler mode waves emitted in the foot/ramp [Hellinger et al., 2007, Lembge et al., 2009, Scholer & Burgess, 2007 or the upstream ultra-low frequency (ULF) waves and steepened foreshock structures [e.g., Burgess, 1989, Krauss-Varban et al., 2008, Krauss-Varban & Omidi, 1991, Schwartz & Burgess, 1991.
Self-reformation is favoured by large M A and small upstream β and produces shock rippling at small spatial scales along the shock surface ( 1 ion inertial length, d i ). Large amplitude whistler waves emitted in the foot cause rippling at spatial scales of a few d i . Finally, irregularities due to upstream ULF waves occur at spatial scales of several tens of d i .
Shock surface rippling was first studied by comparing with 2D hybrid (kinetic ions, masseless fluid electrons) simulations by Winske & Quest [1988] and later by, for example, Lowe & Burgess [2003] and Ofman & Gedalin [2013]. The latter showed that the difference between the global shock normal (determined by the far upstream and far downstream states) and the local normal (determined by the local direction of the fastest variation of the magnetic field) due to rippling can be as large as 40 • . Burgess [2006a,b] studied simulated high Mach number (M A =9.7), almost perpendicular (θ BN 90 • ) shocks and showed that the B-field profiles observed by virtual probes varied even when the probes were closely separated (≤2.5 d i ) and that the shock surface ripples cause large variations of the local θ BN .
Larger-scale irregularity of shock surface was reproduced by Krauss-Varban et al. [2008] who performed local hybrid simulations of a shock with M A =4.7 and θ BN =50 • . They observed that in the case of large-scale shocks, such as IP shocks, even a very small amount of reflected ions generate upstream compressional waves that bend the B-field lines and change the local θ BN . The portions of the shock that become more parallel eject more protons back upstream and further enhance the compressional waves there. These regions were found to travel along the shock surface leading to irregularities with a wavelength of 100 d i or more.
Although the particle acceleration mechanisms at Earth's bow-shock and at IP shocks should be the same, the IP shocks at 1 au have much larger curvature radii (∼0.5 a.u.) and the majority exhibit smaller magnetosonic Mach numbers [up to four, e.g. Blanco-Cano et al., 2016]. This is reflected in the way the ions are accelerated at IP shocks. For example, Kajdič et al. [2017] reported first observations of field-aligned ions upstream of an IP shock with much higher energies that those observed at the Earth's bowshock [e.g., Gosling et al., 1978].
Here we use the multi-spacecraft capabilities of the Cluster mission to study the structure of IP shock fronts. At the times of the shocks the four probes were separated between several tens to several thousands of kilometers, which corresponds to less than one to several tens of upstream ion inertial lengths (d i ). All the shocks exhibit moderate Alfvénic Mach numbers (M A ≤6.5). Three of the shocks also exhibit upstream ion β values (ratio between ion thermal and magnetic pressures) larger than unity. This is different from the studies of the rippling of the Earth's bow-shock surface where β was either larger than 1 but M A were also very large [e.g., Moullard et al., 2006] or β was much smaller than 1 and M A was moderate [e.g., Johlander et al., 2016, Yang et al., 2018.
We compare shock profiles, local shock normals and geometries at each spacecraft and show that these vary even at small (<10 d i ) spacecraft separations. We attribute these differences to irregular IP shock surface and further strengthen our case with 2D hybrid simulations.

Observations: Shock profiles
In this section we compare B-field profiles of four shocks observed on 17 January 2001, 5 April 2010, 26 February 2012 and 8 March 2012. We use data obtained by the Flux Gate Magnetometers [FGM, Balogh et al., 1997, 2001 onboard the four identical Cluster [Escoubet et al., 1997] probes. Figure 1 shows Cluster FGM data with time resolution of 22 vectors per second at four different times during which the shocks were observed. Different colors correspond to different spacecraft with black, blue, green and red lines representing the Cluster 1 (C1), 2 (C2), 3 (C3) and 4 (C4) data, respectively. The IP shocks were selected from the Catalog of IP shocks observed in the Earth's neighborhood by multiple spacecraft available at http://usuarios.geofisica.unam.mx/primoz/IPShocks.html. Basic shock parameters are exhibited on each panel: the range of angles θ BN , the shock's Alfvénic (M A ) and magnetosonic (M ms ) Mach numbers (if all the data are available, otherwise only M A is provided) and the upstream ion β. If the data were missing in the catalog, values were obtained from the Comprehensive database of interplanetary shocks at www.ipshocks.fi.
During these times the Cluster probes formed different constellations.
Separations of pairs of spacecraft are provided in Table 1.

17 January 2001 shock
On 17 January 2001 a fast forward shock (panel a in Figure 1) was detected at 16:27:48 UT by the C1 and C3 spacecraft. C2 detected the shock about a second earlier, while C4 observed the initial B-field increase at the same time as C1 and C3, but the increase was much more gradual in its data. We see from Table 1 that in this case the spacecraft separations were between 530 km (7.4 d i ) and 1300 km (18.1 d i ). C2 was the most sunward spacecraft, therefore detecting the shock first. It was followed by C3 and C4 while C1 was closest to the Earth. Although C3 and C4 had the most similar X GSE coordinates it was C1 and C3 that detected the shock ramp simultaneously. The spacecraft positions and the local shock orientation in the Y GSE -Z GSE plane probably account for this.
We find that shock profiles observed by the four spacecraft are not the same ( Figure 1a). C2 (blue) observed a ramp, an overshoot and possibly a foot. A small peak just upstream of the shock ramp could be a small whistler wave precursor. The B-field magnitude of the overshoot was 7.2 nT, which is more than in the other three cases. C1 and C3 (black and green) observed the shock ramp simultaneously but the B-field magnitudes were slightly different (7.0 nT and 6.5 nT, respectively). Just after the ramp, both spacecraft observed an overshoot, a short lasting undershoot and another peak with similar B-field magnitude value as the first peak. Finally, the shock ramp observed by C4 (red) is not as steep as in the other three cases. After the overshoot there is a strong dip (undershoot) after which the B-field magnitude rises again. Just before the shock ramp, C1 (black) observes a small peak similar to that observed by C2. C4 does not observe any foot, however it detects a distinct compressive structure which lasts ∼2 s peaking at 16:27:46 UT. Similar structures are observed by C1 at 16:27:35.5 UT and 16:27:42.5 UT.

5 April 2010 shock
This shock was detected at 08:24:59 UT by C1 (Figure 1b). C2 detected it about two seconds earlier while C3 and C4 observed it simultaneously some 4 seconds later. The positions of the spacecraft were such that C3 and C4 were separated by only 200 km (1.5 d i ), so they observed almost identical shock profiles. These two probes observe a well defined ramp and overshoot followed by two dips and peaks. The latter could be old overshoots from previous reformation cycles or compressive waves. C3 detected a short lasting (∼1 s) whistler precursor with frequency of ∼3 Hz which can barely be distinguished in C4 data. C2 (blue) observed a steep ramp and an overshoot followed by a very short lasting dip and another increase of B. Upstream of the ramp C2 observed short lasting whistler precursor. Further downstream there is a deep dip in the C2 data, similar to that observed by C3 and C4. C1 (black) observed a much more gradual shock transition, typical of quasi-parallel (θ BN ≤45 • ) shocks.

26 February 2012 shock
C1 observed this shock at 21:38:18 UT ( Figure 1c) followed by C3 and C4 spacecraft that observed it simultaneously ∼3 s later while C2 observed it ∼5 s later. C3 and C4 were separated by only 38 km (0.37 d i ) while the other separations were much larger. It is interesting that while all θ BN values (13 • -39 • ) indicate that this is a quasi-parallel shock, all profiles resemble those of quasi-perpendicular (θ BN >45 • ) shocks with steep ramps and overshoots. C1 and C2 observed a high frequency whistler precursor just upstream of the shock. These waves exhibit small amplitudes (0.1 nT -0.2 nT) and frequencies of ∼2 Hz. In the case of C1 the precursor lasts for ∼1.5 s while in C2 data it lasts for ∼7 s. C2, C3 and C4 profiles show three well defined peaks separated by two dips. Although these peaks could be some compressive waves they could also be overshoots from previous reformation cycles (see Section 3). In the case of C2 the first peak is actually separated from the shock ramp and exhibits higher B-field magnitude than the actual overshoot. C1 profile also shows these features but with much smaller amplitudes.

8 March 2012 shock
Our last shock was observed by C1 at ∼11:02:43 UT (Figure 1d) then by C2 about 3 seconds later and lastly by C3 and C4 spacecraft simultaneously. C3 and C4 were separated by 55 km (0.92 d i ) while the other separations were much larger (Table 1). C3 and C4 profiles are thus almost identical. All spacecraft observe whistler mode precursors upstream of the shock ramp but their appearance varies from spacecraft to spacecraft. In C1 data the whistlers lasted for ∼7 s featuring at least three wave-trains, while the other three spacecraft observed one wave-train during ∼5 s time interval. Frequencies of these waves were ∼2 Hz and their amplitudes increased with proximity to the shock ramp. The shock ramp was steepest in C1 data, while the other profiles show whistler waves inside the shock ramps.

Observations: Local shock normals
The θ BN angles shown in Figure 1 were obtained using the magnetic coplanarity theorem. This requires averaging of upstream and downstream fields during chosen time intervals (but exclude the shock transition) which are then used to calculate the shock normal and θ BN . Thus one obtains some time-averaged values. When multiple inter-spacecraft separations are small ( 100 d i ), one would expect the shock normals calculated this way to coincide within the margin of error. This is because nonstationarity is quasicyclic so local shock normals and θ BN vary in time around some average value (see Section 4). In order to study shock surface irregularity, we need local shock normals at the times when the shocks were observed by each spacecraft and see how they vary as a function of inter-spacecraft separation.
For this we use a novel one-spacecraft method based on the shock normal coordinate system (SNCS). The latter contains three perpendicular axes, n, l and m. The n-axis is parallel to the shock normal, the l-axis contains a projection of the upstream B-field on the shock plane, while the m-axis completes the right-hand coordinate system. In this coordinate system only B l component changes from upstream to downstream. This is of course strictly true only for MHD shocks. In the case of collisionless shocks there exist an out-of-plane component of the magnetic field produced in the shocks's foot and overshoot. However the largest variation of the B-field is still produced due to the shock ramp itself and it occurs in the l direction.
In order to find the SNCS using given interval, we first smooth the B-field data by using a 4-second sliding window in order to remove the upstream whistlers. We then perform minimum variance analysis [MVA, Sonnerup & Scheible, 1998] of the B-field across the shock and postulate that the direction of maximum variance gives us the l-direction. We also obtain two more vectors, perpendicular to l. We then rotate one of them around the l-axis and calculate the absolute value of the mean of the B-field projection along it. Once this value reaches its minimum close to 0, we take the corresponding vector to point along the m-axis and the remaining vector has to point along n.
This method is not without errors. Their main sources are the intrinsic error of the MVA and the choice of time intervals used for the MVA. Details on how the errors were estimated and examples of B-field profiles of shocks in the SNCS are provided in the Appendix A and in the supplement repository available at https://zenodo.org/record/2587992. Figure 2 shows angles between pairs of normals (θ N N ) as a function of spacecraft separation in units of upstream d i . The results are summarized in Table 1. It can be seen that θ N N angles match very well at small spacecraft separations. As long as the spacecraft are 5 d i apart, their normals are within ≤4 • . As the spacecraft separation increases to ∼100 d i , the θ N N angles grow to ∼30 • . Table 1 also shows local θ BN angles. We estimate the errors of θ BN to be ∼ ±5 • . These stem mostly from errors of normal directions and the errors due to selection of time intervals used for calculating the upstream B-field. It can be seen that the average θ BN angles vary substantially from probe to probe.

Simulations results
In order to further strengthen our case in favor of IP shock irregularity at ion scales, we use the 2D HYPSI numerical code [Burgess & Scholer, 2015, Gingell et al., 2017 to carry out two simulations of collisionless shocks.
We use a grid of N x ×N y = 1000×800 cell with cell size of 0.5 d i in both directions. The SW is injected at the left boundary with velocity of 3.0 V A and 4.5 V A resulting in M A ∼4.4 and 6, respectively, for the shock. The upstream β values were 2.4 and 0.2, respectively, which is close to the observed values. Density and B-field are normalized to their upstream values. Initially there were 100 simulation particles per cell. At the right boundary the particles are reflected while periodic boundary conditions are applied in the y direction. The upstream B-field lies in the simulation plane at an angle of 50 • with respect to the x-axis. This is also the value of the average θ BN which we denominate as θ BN,0 .
In Figure 3a we show simulation results for high-β, low M A shock at time t=222.5 Ω −1 i (Ω i is the upstream proton gyrofrequency) with x = [-25,25] d i and y = [40, 120] d i . Results for high-M A shock can be seen in the Appendix B and in the supplement repository. x=0 is the average shock position obtained from the averaged (in y) B-field profile. The colors represent the B-field magnitude. The white curve marks the shock front determined by B≥2.2 and then smoothed with eleven-cell (5.5 d i ) wide runing window. Finally, we calculate the local shock normals (blue arrows).
In the upstream region, there are compressive structures shaded in darkerorange shades. These features are convected towards the shock surface and are the primary cause of its reformation. To show this, four animations are available in the supplement repository. Two of those animations correspond to the Figure 3a). The long animation begins at t=102.5 Ω −1 , ends at t=327.5 Ω −1 and is 30 s long. The shorter animation is 5 s long and shows the shock evolution durint t=220 Ω −1 and t=232.5 Ω −1 . Another two animations correspond to the Figure 7a), have durations of 5 and 36 s and show shock evolution during time intervals of 210-232 Ω −1 and 52.5-317.5 Ω −1 , respectively. In both cases, we can observe compressive upstream magnetic structures being convected towards the shock front. When they get really close and start touching the shock, their amplitudes rise sharply. They merge with the shock front and their upstream edges become new shock ramps. Figure 3b shows four B-field profiles for y=50 d i (black), 60 d i (red), 65 d i (cyan) and 100 d i (magenta). The black profile exhibits a gradual rise and several downstream peaks and dips. The red profile shows upstream whistlers and a steep ramp. The cyan profile does not show any whistlers, only a steep ramp. The magenta profile exhibits a more gradual rise and a strong downstream dip.
We can see in Figure 3a that the upstream variations seen in black, red and cyan profiles are due to upstream structures that have arrived close to the shock and/or the whistler precursors. The multiple dowstream peaks may be due to old overshoots from previous reformation cycles (red profile). The large dips seen in black and magenta profiles coincide with downstream regions with B values similar to those in the upstream region. These regions, although downstream of the shock, have not yet been fully compressed. Figure 3c shows the distribution of θ BN angles at time t=222.5 Ω i . The average and median (µ) θ BN angles are close to θ BN,0 , but the local θ BN s have values anywhere between ∼10 • and ∼90 • . Figure 3d shows the time evolution of the θ BN for the point on the shock surface at y=100 d i . We see that θ BN oscillates around θ BN,0 . The Fourier spectrum of the θ BN variation (Figure 3e) reveals the presence of several periods. Figure 3f shows θ N N angles for all pairs of normals on panel a) as a function of distance. These angles increase and decrease due to the shock irregularities as seen on panel a. There is no real tendency between θ N N and the distance, although θ N N tend to be smaller at small distances.
The high-M A , low-β run (see Appendix B and the supplement repository) differs from the one discussed here mainly in that upstream and downstream compressive structures and the shock exhibit larger amplitudes, the standard deviation of the θ BN distribution is larger, the periods in the θ BN spectra are shorter, and the shock is more structured especially at smaller separations ( 10 d i ), so there is even less correlation between θ N N and the distance.

Discussion and conclusions
In this work we present the first direct observational evidence for an irregular surface of IP shocks at ion scales. We show four case studies (Figure 1) that were observed by the four Cluster probes with inter-spacecraft separations between 38 km and ∼10 4 km (0.37 d i -92 d i ). We show that B-field shock profiles vary from probe to probe. When the spacecraft were 1 d i apart (Table 1), the profiles are very similar. When spacecraft are separated by more than 10 d i the shock profiles differ significantly. We attribute these differences to be associated with an irregular shock front.
We further strengthen our argument by calculating local shock normals at each spacecraft for which we design a new one-spacecraft analysis method (see Section 3). We plot the angle between pairs of single spacecraft normals, θ N N , as a function of the distance between the probes. On average these angles tend to be smaller when the spacecraft are 5 d i apart and then increase as the distance increases up to ∼100 d i .
We also calculate θ BN angles at each spacecraft and show that these can be very different at different points on its surface (Table 1). This should be taken into account in any interpretation of data from IP shocks.
Our findings fit well with the 2D hybrid simulation results which show that: • Shock profiles and local θ BN may vary significantly at separations of the order of 5 d i and more.
• At any given time, different locations on the shock surface exhibit values of θ BN anywhere between 10 • and 90 • .
• The geometry at the particular point on the shock surface varies with time and the location on the shock surface.
• The cause of irregularities of observed shocks may be upstream magnetic structures that are convected towards the shocks. These can arise due to small amount of backstreaming ions reflected by the shocks.
The fact that irregularities of simulated shock surfaces may occur at quite small spatial scales is not reflected in our observations (Figure 2), possibly due to low number of our case studies.
In the supplement repository and Appendix B we show that fluctuations in the ULF frequency range (0.01-0.1 Hz) exist upstream of all four shocks (Figure 1), although in the case of February 2012 shock their compressive component is weak. These fluctuations could be responsible for irregular structure of observed shocks.
There is a possibility that the upstream B-field fluctuations have already been present in the upstream SW and the shocks just caught up with them. In the supplement repository and Appendix B we show ion spectra for the January 2001 and April 2010 shocks in Figure 2. The ion particle energy flux peaks at shock transitions, suggesting the ions are accelerated by the IP shocks (for the other two shocks the data were not available). In the case of the April 2010 shock part of these ions and ULF fluctuations may have actually come from the Earth's bow shock, as the data suggest that the Cluster probes have entered and exited it on several occasions before the IP shock arrival. These excursions are marked by increased suprathermal ion fluxes (green trace in the middle panel of the Figure 2b). The last excursion into the foreshock occured ∼20 minutes prior to the IP shock arrival and may have lasted some minutes after the shock was observed.
Acknowledgments Authors acknowledge ClWeb and Cluster Science Archive teams for easy access and visualization of the data. PK's work was supported by PAPIIT grant IA101118. XBC's work was supported by PAPIIT IN-105218 and Conacyt 255203 grants. DT acknowledges support of a studentship funded

Appendix A
In this section we explain how the shock normals were calculated and how the errors of θ N N and θ BN were estimated. There are two main sources of errors. The first is the error of the MVA method itself which depends on the number of measurement points and the calculated eigenvalues [Sonnerup & Scheible, 1998]: Here λ 2 , λ 3 and M are the intermediate and minimum eigenvalues and the number of measurement points, respectively. This error is stated in the Figure 4.
The second source of errors comes from determining time intervals which are used for the MVA. These intervals need to include the shock transition but also some upstream and downstream regions. One needs to select the intervals carefully so not to include large B-field rotations that are not associated with shocks and could affect the the determination of the direction of maximum variance. We select the time intervals by hand. We repeat the process for each shock and spacecraft ten times. We then proceed to calculate angles between pairs of normals from different spacecraft (θ N N ) and calculate the the average angles and the error of the mean. Next we sum this error with θ Err in order to estimate the total error of our method. The latter is shown in Table 1 and in Figure 2 in form of error bars.
The θ BN errors stem from the MVA method and the selection of the upstream time intervals over which we calculate the average B-field direction.
They are typically ∼15 seconds long. After repeating this selection for 10 times we estimate the errors to be ∼5 • .

Appendix B
Here we present a) wavelet spectra of magnetic field fluctuations observed by Cluster 1 spacecraft upstream of the four interplanetary (IP) shocks (Section 7.1); b) ion spectrograms and energy fluxes around two of the shocks for which the data were available (Section 7.2); and c) Simulation results of our high-M A , low-β run (Section 7.3). Panels ii) and iii) exhibit B and B x,GSE wavelet spectra, respectively. We can see that compressive and/or transverse B-field fluctuation in the frequency range 10 −2 -10 −1 Hz are present upstream of all four shocks. In general, there is more power in the transverse component of these fluctuations than in the compressive component. Figure 6 shows magnetic field, particle spectrogram and particle energy fluxes at time when a) 17 January 2001 and b) 5 April 2010 shocks were observed. In both cases the suprathermal ion energy fluxes in units (in units of keV/(s·cm −2 ·sr·keV)) start increasing before the shock arrival and peak at shock transition, suggesting they are accelerated by the shocks. The suprathermal ion energy flux (and magnetic ULF fluctuations) in the case of the 5 April 2010 could partially arrive from the Earth's bow-shock, since the ion spectrogram suggests that prior to and possibly during the shock encounter, the Earh's foreshock has been observed intermittently.  The colors represent the B-field magnitude. The white curve marks the shock front and blue arrows show the directions of local shock normals. Figure 7b shows four B-field profiles for y=46 d i (black), 50 d i (red), 97 d i (cyan) and 115 d i (magenta). Animations for this run can be found at http://usuarios.geofisica.unam.mx/primoz/IPShockRipplingSupplement/ and are titled BfieldLowBeta.avi, BfieldHighLowShort.avi. Figure 7c shows the distribution of θ BN angles at time t=112.5 Ω i . Figure 7d shows the time evolution of the θ BN for the point on the shock surface at y=97 d i . Figure 7e shows θ N N angles for pairs of normals shown on panel a).   The dotted horizontal lines indicate zero value. All 160 profiles may be seen in the supplement located on-line at http://usuarios.geofisica.unam.mx/primoz/IPShockRipplingSupplement/SNCS.pdf