ULF Wave Transmission Across Collisionless Shocks: 2.5D Local Hybrid Simulations
Abstract
We study the interaction of upstream ultralow frequency (ULF) waves with collisionless shocks by analyzing the outputs of 11 2D local hybrid simulation runs. Our simulated shocks have Alfvénic Mach numbers between 4.29 and 7.42 and their 
angles are 
, 
, 
, and 
. The ULF wave foreshocks develop upstream of all of them. The wavelength and the amplitude of the upstream waves exhibit a complex dependence on the shock's 
and 
. The wavelength positively correlates with both parameters, with the dependence on 
being much stronger. The amplitude of the ULF waves is proportional to the product of the reflected beam velocity and density, which also depend on 
and 
. The interaction of the ULF waves with the shock causes large-scale (several tens of upstream ion inertial lengths) shock rippling. The properties of the shock ripples are related to the ULF wave properties, namely their wavelength and amplitude. In turn, the ripples have a large impact on the ULF wave transmission across the shock because they change local shock properties (
, strength), so that different sections of the same ULF wavefront encounter shock with different characteristics. Downstream fluctuations do not resemble the upstream waves in terms the wavefront extension, orientation or their wavelength. However, some features are conserved in the Fourier spectra of downstream compressive waves that present a bump or flattening at wavelengths approximately corresponding to those of the upstream ULF waves. In the transverse downstream spectra, these features are weaker.
Key Points
-
Upstream ultralow frequency (ULF) waves are not simply transmitted across collisionless shocks
-
Some spectral features of compressive upstream ULF waves are conserved in the spectra of downstream fluctuations
-
There is some correlation between the properties of the shock ripples and the foreshock ULF waves
Plain Language Summary
We address the problem of what happens to upstream ultralow frequency (ULF) waves as they reach shocks and are carried into the downstream region. We do this by analyzing the results of 11 2.5D local hybrid simulations of collisionless shocks. We find that the waves are not simply transmitted into the downstream region but that their identity is largely destroyed, as the downstream fluctuations do not resemble the upstream waves neither in wavelengths nor in appearance. However, some features observed in the Fourier spectra of upstream ULF waves are conserved in the Fourier spectra of compressive downstream fluctuations.
1 Introduction
The Earth's magnetosheath (MSH, Lucek et al., 2005) is a region sandwiched between the bow-shock (BS) and the magnetopause of Earth. Its existence was first predicted by Kellogg (1962). MSH is highly turbulent, perturbed by different low-frequency (
proton gyrofrequency, 
0.01–0.1 Hz, Schwartz et al., 1996; Song & Russell, 1997) fluctuations. Already at early times various authors (Barnes, 1970; Fairfield, 1976; Fairfield & Ness, 1970; Greenstadt et al., 1970; McKenzie & Westphal, 1969) suggested that the magnetic field fluctuations in the MSH may be generated in various ways: in the magnetosheath itself, at the BS or in the foreshock region (Eastwood et al., 2005).
There is a general consensus (e.g., Schwartz et al., 1996, and the references therein) that in the region of the MSH that lies downstream of the quasi-perpendicular BS (
, where 
is the angle between the upstream magnetic field and the local shock normal), also called quasiperpendicular magnetosheath, the dominant wave modes are the Alfvén Ion-Cyclotron (AIC) and mirror mode (MM) waves which form either behind the BS or deeper in the MSH. The AIC instability is prevalent in the regions, where the plasma 
(ratio between plasma thermal and magnetic pressures) is low (
1), while MM instability occurs in plasmas with 
1 (e.g., Gary & Winske, 1993; Gary et al., 1976, 1993; Samsonov et al., 2007; Schwartz et al., 1996; Sckopke et al., 1990; Song et al., 1992; Tsurutani et al., 1982). In the quasiparallel magnetosheath (
), the wave activity is much higher. It is thought that fluctuations there are either transmitted foreshock waves or that they can be formed at the BS (i.e., Czaykowska et al., 2001; Fairfield & Ness, 1970; Luhmann et al., 1986). In the following sections, we briefly summarize the state of the art regarding wave activity in the MSH.
1.1 Waves in Quasiparallel Magnetosheath
The indirect evidence for the association of the MSH fluctuations with the bow-shock geometry came from observing that these fluctuations are not uniformly distributed in the MSH. Dawn-dusk asymmetry in the long-term average level of fluctuations, consistent with the idea that strong waves are associated with the quasiparallel shock, was reported already by Fairfield and Ness (1970). The reason behind the dawn-dusk asymmetry lies in the average orientation of the interplanetary magnetic field (IMF) in the ecliptic, which follows a nominal Parker spiral. It is more probable for the BS on the dawn side to exhibit quasiparallel geometry, while on the dusk side the quasiperpendicular configuration is more common. Luhmann et al. (1986) were the first to search for a relationship between the spatial distribution of the MSH fluctuations and the IMF orientation and acknowledge the quasiparallel BS to be an important contributor to the fluctuations in the dayside MSH.
It is believed that upstream ultralow frequency (ULF, periods between 10 and 100 s) waves from the foreshock may strongly influence the formation of downstream fluctuations. They may be partially transmitted into the downstream region and mode converted into downstream Alfvénic turbulence. Another source of downstream waves at quasiparallel shocks is an interface instability that arises due to the interaction of incident ions and partially thermalized plasma at the shocks. These waves are only present near the shock transition region (Krauss-Varban & Omidi, 1991; Krauss-Varban, 1995; Scholer et al., 1997).
Czaykowska et al. (2001) performed a study of fluctuations during 132 bow-shock crossings. They found that in the vicinity of quasiparallel shocks the magnetic power was strongly enhanced, which the authors attributed either to wave generation at the shock or an amplification of convecting upstream waves at the shock interface. Du et al. (2008) analyzed magnetosheath fluctuations in the 4–240 s range. They observed large-amplitude compressional and transverse fluctuations downstream of quasi-parallel shocks. Dimmock et al. (2014) studied magnetic field fluctuations in the dayside MSH in the frequency range from 0.1 to 2 Hz. The authors showed a tendency of such fluctuations to exhibit higher amplitudes in the dawn flank magnetosheath and close to the magnetopause during southward IMF. A dawn-dusk asymmetry of Pc3 velocity fluctuations in the dayside MSH was shown to exist by Dimmock et al. (2016). In general, larger amplitudes were observed downstream of the quasi-parallel BS and during the times of fast solar wind (SW).
1.2 Waves in the Quasiperpendicular Magnetosheath
It was suggested by, for example, Omidi et al. (1994) and McKean et al. (1995), that AIC and MM waves grow at or near quasi-perpendicular shocks and are then convected away by the sheath plasma. The fact that AIC and MM fluctuations are often observed downstream of quasi-perpendicular shocks is not surprising. These shocks heat the plasma, preferentially in the perpendicular direction with respect to the local magnetic field (Winske & Quest, 1988), thereby enhancing the temperature anisotropy which is needed for the growth of these two types of waves.
McKean et al. (1995) found that as the short wavelength waves propagate downstream in the MSH, they are heavily damped so that in these regions only wave modes with longer wavelengths survive. The authors suggested that this is probably due to wave-particle scattering which results in gyrotropic and approximately bi-Maxwellian ion distributions.
It was shown by several authors that the temperature anisotropy 
/
is inversely correlated with the parallel proton beta, 
(Anderson et al., 1994; Fuselier et al., 1994; Gary & Winske, 1993; Peter Gary et al., 1995; Phan et al., 1994). Fuselier et al. (1994) argued that this inverse correlation is a consequence of pitch angle scattering of ions by electromagnetic ion cyclotron waves, which regulate the anisotropy and restore a marginally stable plasma. Chaston et al. (2013) showed that broadband kinetic Alfvén waves heat magnetosheath ions. The authors observed that the energy density of ions correlates well with the energy density of these waves. The heating occurs predominantly in the direction perpendicular to the local magnetic field and is limited by the threshold condition for anisotropy instability.
Hubert et al. (1998) used data from ISEE 1 and 2 spacecraft to study the nature of low-frequency waves during a crossing of the Earth's MSH downstream of a quasiperpendicular (
= 
) bow-shock. The authors observed a region of purely AIC waves in a 0.3 
thick layer adjacent to the shock, followed by a region 2
thick where AIC and MM waves coexisted and finally a pure MM region. The authors thus argued that the dominant wave mode is controlled by the depth in the MSH. Czaykowska et al. (2001) found that downstream of quasi-perpendicular shocks the magnetic wave activity is significantly enhanced compared to upstream. During times of low 
the authors found left-hand polarized AIC waves in the MSH. Strong correlation between the temperature anisotropy and the intensity of these waves was also found. Downstream of some highly oblique shocks, MM waves were also observed.
Dimmock et al. (2015) showed that temperature asymmetry 
1 is favored on the dusk (quasiperpendicular) side and this is reflected in a similar asymmetry of mirror mode activity. 
/
decreases with increasing SW Alfvénic Mach number (
), whereas mirror mode occurrence exhibits the opposite trend. The dawn-dusk asymmetry diminishes with increasing 
. Also, during the transition from low to moderate 
the authors observed a shift in the data from MM dips to MM peaks.
Soucek et al. (2015) presented a statistical study of the spatial distribution of MM and AIC waves in the MSH as a function of relevant plasma parameters, such as ion temperature anisotropy and ion 
. The authors showed a strong dependence of the two plasma parameters and the occurence of the waves on the shock's 
and its 
. It was found that the AIC waves occur almost exclusively in the plasma that is stable to MM instability. MM were found to occur in the quasi-perpendicular MSH in correlation with high 
, while lower 
favor the AIC waves. Both are rarely observed in the quasiparallel MSH. AIC are observed behind the shock in the MSH flanks and less so near the subsolar MSH.
1.3 Wave Transmission Across the Bow-Shock
Although plenty of works have been published on magnetosheath fluctuations, questions such as what happens to the foreshock waves as they cross the BS and how the properties of downstream fluctuations relate to those of the foreshock waves, are still not completely answered.
It has been shown that the upstream ULF waves are a mixture of Alfvén and magnetosonic waves (e.g., Eastwood et al., 2003; Hoppe & Russell, 1983; Sentman et al., 1981). Some of the early works (Asséo & Berthomieu, 1970; McKenzie, 1970; McKenzie & Westphal, 1969, 1970) concluded that the magnetosonic and Alfvén waves are strongly amplified on passage through the shock. On the other hand, it was found by McKenzie and Bornatici (1974), that the waves also impact the shock: Alfvén waves diminish the shock compression ratio while sound waves either enhance or diminish the compression ratio depending on whether the incident wave vector is parallel or antiparallel to the upstream flow direction. Asséo and Berthomieu (1970) showed that when an incident wave strikes a shock front, it is refracted and gives rise to five other hydromagnetic waves. Hassam (1978) concluded that small amplitude Alfvén waves may perturb the shock surface and give rise to transmitted waves consisting of a fast magnetosonic wave, forward and backward slow magnetosonic waves and an entropy wave. Whang et al. (1987) showed that the number of downstream waves excited due to transmission of the upstream waves depends on the angle of incidence: when this angle is lower than some critical value, six diverging downstream waves are excited. When the angle of incidence is larger than a critical value, the number of downstream excited waves may be smaller than six. Lu et al. (2009) used 2D hybrid simulations in order to study the interaction of Alfvén waves with a quasi-perpendicular shock. The authors found that the Alfvén waves are transmitted through the shock and that their amplitude is enhanced 10–30 times. The authors also found that the shock ripples form due to the upstream Alfvén waves.
There have been some works showing indirect evidence in favor of the upstream ULF wave transmission into the magnetosheath and even into the magnetosphere. It has been known since roughly 50 years now that a subset of the magnetospheric ULF waves, dayside Pc3-4 pulsations, is enhanced during radial IMF configurations (Troitskaya & Plyasova-Bakunina, 1971). This has been explained in terms of the foreshock ULF waves being the source of the magnetospheric waves in the same frequency range. Several studies have observed simultaneously waves of similar frequencies in the foreshock region and in the magnetosphere (e.g., Clausen et al., 2009; Engebretson et al., 1991; Francia et al., 2012; Greenstadt et al., 1983; Lin et al., 1991; Russell et al., 1983; Takahashi et al., 2016, 2021; Villante et al., 2011). Some local (Krauss-Varban & Omidi, 1991) and global (Shi et al., 2013, 2017) hybrid simulations also favor upstream ULF wave transmission into the magnetosheath and the magnetosphere. On the other hand, Narita and Glassmeier (2005), Narita and Glassmeier (2006), and Narita et al. (2006) do not favor the wave transmission across the bow-shock. Their arguments are that from the viewpoint of the dispersion relation and other wave properties, such as propagation angle and polarization, the foreshock waves are not transmitted into the magnetosheath.
Spacecraft observations alone cannot provide a complete answer to the question of what happens to the waves as they cross the shock, since such data are spatially limited. However, local and global hybrid simulations of upstream wave-shock interaction may provide us with some clues on how this interaction occurs. Local simulations (with self-consistent foreshock) are complementary to global simulations with complex foreshock which links regions of shock with different 
. In this work, we try to answer this question by analyzing the outputs of local 2.5D hybrid (kinetic ions, fluid electrons) simulations of collisionless shocks and the corresponding upstream and downstream regions. In total we show results from 11 local runs with Alfvénic Mach numbers 
ranging between 4.29 and 7.42 and the 
between 
and 
. Thus, we simulate quasiparallel and marginally quasi-perpendicular shocks. These shock properties are comparable to those of the Earth's bow-shock. In the case of the latter, the typical 
at its nose ranges between 6 and 7 (Winterhalter & Kivelson, 1988) while it is lower toward the flanks.
In the remainder of this paper, we first describe the simulation setup in Section 2. Next, we discuss the simulation results in Section 3. In this section, we begin by analyzing the degree of the shock rippling at different shocks. We emphasize this process before studying the downstream ULF wave properties since, as we show later on, shock rippling is crucial for understanding, how the ULF waves are affected by the shocks. This process results from the fact that the shock reformation is not in phase all over the shock surface. Thus, the shock properties (
, compression ratio, etc.) change with the location on the shock surface, so different sections of a single upstream ULF wave encounter the shock with different properties, which then affects the wave transition. In Section 4, we analyze an example of a ULF wave interacting with a shock. Finally, in Section 5 we discuss the results of this work and present the conclusions.
2 Simulation Setup
Local hybrid simulations were performed with the 2.5D HYPSI numerical code (Burgess & Scholer, 2015; Gingell et al., 2017; Trotta & Burgess, 2019). In these simulations, the electrons are treated as a massless, charge-neutralizing fluid, while ions, in our case pure protons, are treated kinetically. We perform a series of simulations with a grid of 
= 1000 
800 cells whose dimensions are of 0.5 
(
is the upstream ion inertial length) in both directions, respectively. The SW is injected from the left along the 
-axis with inflow velocities of 3.3 
, 4.5 
, and 5.5 
(
is the initial Alfvén speed). The upstream magnetic field lies in the XY plane. Its initial magnitude is one and its orientation is such that four different shock geometries are produced with 
= 
, 
, 
, and 
. The upstream ion and electron 
(ratio between thermal and magnetic pressures) are 0.5. Initially there were 100 particles per cell. The simulation is periodic in 
direction, while the right boundary acts as a reflective wall. Time is measured in units of the inverse of proton gyrofrequency (
). The simulation timestep is 0.005 
and the outputs are produced every 2.5 
. In total, we performed 11 simulations with initial conditions and the resulting shock Mach numbers that are summarized in Table 1. As the incident protons are reflected at the right boundary, a shock is created which then starts to propagate toward the left. Throughout the paper, the physical quantities, such as velocity and dynamic pressure, are calculated in the shock-rest frame. The running times of the simulations range between and 250–325 
.
  |
 ( ) |
 |
 |
|---|---|---|---|
| 3.3 | 15, 45, 50 | 0.5 | 4.29, 4.48, 4.57 |
| 4.5 | 15, 30, 45, 50 | 0.5 | 5.31, 5.82, 5.94, 6.05 |
| 5.5 | 15, 30, 45, 50 | 0.5 | 6.97, 7.07, 7.17, 7.42 |
3 Results
3.1 Overview
Figures 1 and 2 show simulation outputs for all our runs at times when the bow-shock was located at 
in the simulation domain. The upstream regions lie to the left of the shocks. Figure 1 exhibits maps of 
, while Figure 2 shows 
B
. The units of magnetic field are normalized to the initial 
B
value.
Simulation results of 11 local hybrid runs at times when the average shock location was x 
200 
. Colors represent the 
component and the two black lines mark the locations along which the Fourier spectra in Figures 4g–4k are calculated. 
increases from left to right and 
from top to bottom.
200 
. Colors represent the B magnitude and the two black lines mark the locations along which the Fourier spectra in Figures The different panels show shocks with different properties. The 
increases to the right, while the inflow velocity and thereby the 
of the shocks increases from top to bottom. We also observe that 
increases slightly with the increasing 
. This is due to the fact that for the same upstream SW properties (density, velocity, temperature), the quasiperpendicular shocks tend to propagate faster than the quasiparallel ones.
The colorscales in both figures are centered on the initial upstream values. Black lines mark the cuts along which the wave properties are studied in Section 3.3 and were selected so as to include as many wavefronts as possible. We first look at the 
plots shown in Figure 1. We can see that the wavelength and the orientation of the wave fronts of the upstream fluctuations varies from panel to panel. 
fluctuations are by definition transverse, which means that they propagate at small angles with respect to the magnetic field direction and their wavefronts are approximately perpendicular to B. The wavefronts in the first column (panels a, d, h), where the magnetic field makes an angle of 
with respect to the 
-axis, are oriented almost perpendicular to 
. On the panels in the last column, where nominal 
= 
the upstream wavefronts are strongly inclined with respect to the 
-axis.
The properties (wavelength and the amplitude) of the upstream waves change with 
(from top to bottom) and 
(from left to right), which is addressed in more detail in the Discussion section.
In the case of the 
B
fluctuations (Figure 2), the upstream wave fronts are much more aligned with the magnetic field. The analysis of sequential outputs reveals that these waves propagate obliquely to B. It was shown by Omidi (2007) and Blanco-Cano et al. (2011) that these fluctuations are compressive fast magnetosonic waves. We can see in Figure 2 that in the case of waves on panels (b), (c), (g), and (k), the wavefronts differ from other runs in that they do not seem to be always aligned with the magnetic field but many of them appear like more localized magnetic field enhancements.
Downstream of the shocks the fluctuations exhibit a very different appearance. In all cases a filamentary structure can be observed. The orientation of these filaments is different from that of the upstream waves and changes with the increasing nominal 
: as this angle increases, the orientation of the filaments forms an increasingly larger angle with respect to the direction of the nominal shock normal.
3.2 Statistical Analysis of the Shock Ripples
to 
100s 
. On the other hand, the simulated shocks exhibit surfaces, which are rippled and variable in time. Both facts influence importantly the wave transmission across the shocks. We show a close-up of one of the shocks in Figure 3. Here, colors represent the magnetic field magnitude, the vertical magenta line marks the average shock location (here, the coordinate system is shifted along the 
-axis so that this location is at 
= 0. Note that the Figure only shows a portion of the simulation domain, which extends further in the 
direction. There are other ripples, which result in the average shock position being centered on zero in this plot.) and the thick black line represents the shock surface. The latter was calculated the following way:
-
First, we average the magnetic field profile in
direction at the time t = 225.5 
. Then we determine the 
-position of the shock to be where the averaged B-magnitude first reaches the value of 2.5. This value was selected since it coincides with the shock transition in the averaged shock profiles. -
Afterward, we determine the position of the shock surface in the original (not averaged) magnetic field data. For each
-coordinate we search for the first 
-coordinate at which the B-magnitude reached the value of 2.35. The search was done in the vicinity of the shock position calculated from the averaged shock profile. The value of 2.35 was selected by eye after many attempts to approximate each individual shock surface with a single curve. The reason why the search for this value was made near the position of the shock determined from the average B-field profile, was to avoid the upstream compressive structures, such as shocklets, to be identified as parts of the shock. -
Next, we smooth the calculated
-coordinates using the Savitzky-Golay filter with width of 21 points and of the third order. This calculated surface can be seen in Figure 3 as a thick black line. We can see that the curve closely matches the magnetic field jump. -
Finally, we average all
-values of this black line to obtain an averaged shock location. This location is different from that calculated from the averaged magnetic field profile and this is the position marked by the purple vertical line in Figure 3.
An excerpt from our run e with the thick black curve representing the shock surface and the vertical purple line marking the average shock location. Here we use a coordinate system that is shifted along 
, so that the average shock position is at 
= 0.
In order to establish how rippled the shocks are, we show three diagnostics that describe three properties of the ripples: the amplitudes, the length scales, and how steep their sides are.
In Figures 4a and 4b, we show two histograms that exhibit the distributions of distances (
) for every point on the shock surfaces from the corresponding shock's average location. The panel on the left exhibits the histogram for the shock from the run a (from now on, we will refer to different runs with the letters that correspond to the panels in Figures 1 and 2), while the one on the right is for the shock from the run i (histograms for all the runs can be seen in Figure A1). The most extreme values in each histogram, whether negative or positive, provide information about the maximum amplitudes of the ripples. The width of the histograms (represented by 
) also provides information on the ripples: small 
means that the majority of the ripples have smaller amplitudes, while large 
point to ripples with larger amplitude. We can see on panel (4a), that the bulk of the 
x values is below 10 
. The distribution is strongly peaked at around 0 and the standard deviation is 
= 3.94 
, the smallest among all runs. The value of 
for the run i (4b) is 8.76 
, the largest of all. Figure 5a exhibits 
values for all the runs. The horizontal dashed line marks the value of 
= 6.0. We divide the simulated shocks into those with 
6.0 (from runs a, b, f, h, j, and k) and those with 
6.1 (from runs c, d, e, g, i). This division is arbitrary but it will help us with further analysis.
(a, b) Histograms of distances 
of the points on the shock surface from the shock's mean location. (c, d) 
angles for each pair of normals calculated for all points on the shock surfaces. (e, f) Power spectra of the curves representing the shock surfaces. Left panels represent data for the weakly rippled shock from the run a, while right panels are for a strongly rippled shock from run i.
Values of (a) 
, (b) inte, (c) 
, and (d) maximum and median 
for all the runs.
Figures 4c and 4d show power spectra of the curves describing the surfaces of the shocks from runs a and i. An example of such a curve is shown in Figure 3. The idea is to check whether the ripples exhibit specific wavelengths. Only two spectra are exhibited in this figure, however spectra for all shocks are presented in Figure A2. All of them are continuous and their power at wavelengths above 40 
(vertical blue lines) is enhanced compared to the trend at smaller wavelengths (red dashed lines). On the panels (4c) and (d) we state a quantity called “integrated power” which was obtained by first integrating the spectra and then subtracting the power that can be attributed to small-scale (below 40 
) background (red dashed lines on panels 4c and d). We can see that this value is 16.5 for the run a (smallest among all) and 50.4 for run i (largest of all). Again, based on Figure 5b, which shows integrated power values for all shocks, we can divide the shocks into two groups - runs a, b, c, f, g, h, j and k exhibit integrated power
26 (horizontal dashed line), while the runs d, e, and i exhibit the integrated power
26.
Figures 4e and 4f show another statistic from the simulations. For each point on the shock surface, we calculate the local normal at that position. Then we calculate the angle between every possible pair of normals, 
, and plot this as a function of distance between the points (plots for all the runs are exhibited in Figure A3). The maximum angles observed for each case are a proxy of how strongly the shocks are rippled.
We can make a thought experiment in order to have a clearer picture about what the plots in Figures 4e and 4f tell us. If the shock surface could be described with a single sinusoidal curve and if we took one point on a crest of one of the ripples, and plotted 
values as a function of a distance from that point, we would clearly observe 
peaks and dips at regular intervals corresponding to the half of the wavelength of that sinusoidal function. 
would exhibit values of 0 at the locations of the crests (minima or maxima since the shock normals by definition always point toward upstream). If we used more than one sinusoidal function to describe the shock surface, the picture would become more complicated, we would observe many irregular peaks at intervals that would exhibit some quasi periodicity. We would also get a range of 
values at each distance. If the typical wavelength of the ripples was small, we would get many thin peaks, if it was large, fewer and thicker peaks would be produced. 
values close to 
would only be obtained between points lying on the sides of steep ripples, pointing in opposite 
directions.
We can see from Figures 2a, 2i, 4e, and 4f that there seems to be a correlation between the upstream compressive waves and the peaks: in the case of the run i, the wavefronts are many and exhibit small widths, and this is likely why Figure 4e exhibits many peaks which tend to be narrower. On the contrary, wider and less numerous upstream waves, such as those in Figure 2i, are associated with wider and fewer peaks (Figure 4f).
We show three numbers on the panels (4e) and (f): the maximum and median 
and 
, which is the fraction of the area on the panels delimited by the black points. Once again, based on Figures 5c and 5d, we can divide the runs in two groups: the runs a, b, f, h, j, and k exhibit median 
values 
and 
0.42 (horizontal dashed lines on both panels), while the the runs c, d, e, g, i exhibit higher values of 
and 
.
We can see that the shocks from the runs a, b, f, h, j, and k all form a group that consistently exhibits smaller 
, integrated power, 
and 
values. Just for illustrative purposes, we call these shocks weakly rippled. The shocks from the runs d, e, i always exhibit larger values, so we refer to them as strongly rippled. The shocks from the runs c and g exhibit higher values of 
, 
, and 
, but smaller values of the integrated power. Hence, we denominate them as intermediately rippled.
3.3 Analysis of Upstream and Downstream Wave Properties
3.3.1 Fourier Spectra
Figures 4g, 4h, 4i, and 4j show Fourier spectra for upstream (black curve) and downstream (red) waves. Spectra from all runs can be seen in Figures A4 and A5. These are instantaneous spectra obtained at simulation times indicated in Figures 1 and 2 along the cuts marked with black lines. These cuts were chosen so that they are approximately perpendicular to the wave fronts and they include the maximum number of the waves. This way the units on the 
-axis in Fourier spectra represent true wavelengths, not some projections along the cuts and the power of the foreshock in the Fourier spectra waves is maximized. The spectra were obtained by first windowing the original data with the Tukey (tapered cosine) function (Harris, 1978) and then applying a 13-point running average filter on the FFT power spectra.
We can see that the downstream spectra exhibit higher power than upstream spectra. The upstream spectra of transverse (
) and compressive (
B
) fluctuations are rather flat at wavelengths 
10–60 
. The exact value depends on the run. We shall refer to these intervals as “ULF range” (shaded in blue) since the simulated upstream waves are analogous to the foreshock ULF waves. Parts of the upstream spectra at wavelengths above the ULF range tend to increase smoothly without strong features.
Normally, we expect the foreshock ULF waves to produce a bump in the Fourier power spectra (e.g., Greenstadt et al., 1995; Hoppe & Russell, 1983; Hoppe et al., 1982; Le & Russell, 1992, 1994). However in these studies the data contain many wavefronts. In our case, as can be seen in Figures 1 and 2, the upstream cuts include 
10 wavefronts, in some cases as few as three (see for example panels 1h and 2k). Due to the fact that the number of wavefronts is low, the small peaks in Fourier spectra do not merge into a single bump.
We can observe certain features in the upstream spectra that are also present downstream. In the case of the weakly rippled shock from the run a (4g and i) the flattening can be clearly observed in the 
B
and 
spectra. In the case of the strongly rippled shock from the run i, this flattening is much less obvious in the 
spectrum than in the 
B
spectrum. In general (see Figures A4 and A5), flattening in the ULF range can be observed in the case of the spectra shown on panels (a), (b), (f), (h), (j), and (k), which are classified as weakly rippled. In the case of the panel (g), the downstream spectrum does not exhibit flattening but several peaks. In the case of the runs (d) and (i) the downstream spectra are rather featureless. It thus seems that flattening in the ULF range is better observed downstream of weakly rippled shocks.
In general, the downstream spectra are more turbulence-like than the upstream spectra meaning that the non-turbulent features (bumps and flattening) are less prominent. The latter holds even more for the transverse, 
spectra in Figures 4i and 4j.
3.3.2 Wavelengths and Amplitudes of Upstream Waves
In previous sections, we imply that the upstream wave properties determine how they are going to be affected during shock crossing. Their wavelengths and amplitudes determine the properties of the shock ripples, which in turn affect the transmitted wave properties. Here, we take a look on how the wavelength and the amplitude of the waves depend on 
and 
.
and 
, however the correlation with these two parameters is not obvious at the first sight. Watanabe and Terasawa (1984) have shown that the frequency of the upstream ULF waves in the spacecraft frame, 
, depends on the upstream solar wind speed 
, the reflected ion beam speed, 
, the wave phase speed 
, and angles 
and 
, between the wave vector 
and the upstream SW velocity and B-field directions, respectively. These authors derived the expression for the ratio between 
and the proton gyrofrequency, 
:
(1)
, one needs to know the origin of the upstream ULF waves. It was first proposed by Fairfield (1969) that these waves are formed due to the electromagnetic ion beam-cyclotron instability excited by the reflected, upstream propagating field-aligned ion beams. This instability generates low frequency, right-hand mode waves. In order for them to be excited, a cyclotron resonance condition must be fulfilled:
(2)
and 
are the wave frequency and the beam velocity in the rest frame of the upstream plasma. The beam velocity in the spacecraft frame can be calculated by following Schwartz et al. (1983) calculations for magnetic moment-conserving reflection at the shock:
(3)
. Their propagation direction and thus their wave vectors 
are approximately parallel to the upstream magnetic field, thus yielding 
and 
. In our simulations, the “spacecraft rest frame” corresponds to the shock normal incidence frame in which the upstream SW velocity and the shock normal anti-parallel. Thus the acute angle 
(between the SW velocity and shock normal) is 0 and 
. This, plus the expression (Equation 3) may be introduced into the expression (Equation 1) to obtain
(4)
, we obtain:
(5)
, we arrive to the expression for the wavelength of the upstream waves:
(6)This wavelength does not depend on the selected frame of reference and thus corresponds to the wavelength measured in the rest frame of the simulation domain. The numeric factor 
is the same for all our simulations, so we can define the normalized wavelength as 
divided by this factor. This quantity is shown in Figure 6 for all our models. We can see that the quantity that most strongly correlates with the normalized 
is the 
. The normalized 
for models (h–k) (
= 5.5 
) varies between 
0.8 and 
1.7, so by a factor of 
2.1. Maintaining the 
= 
(models c, g, and k) and varying the 
by varying the inflow velocity 
, produces normalized 
values between 1.43 and 1.70, so only by a factor of 
1.2. By visual inspection of Figure 1, we can observe that the wavelength of the upstream waves indeed increases strongly with 
(from left to right), while its increase due to 
(from top to bottom) is somewhat smaller, although still significant.
Predicted normalized wavelength 
for different models. 
angles are stated for each model.
(7)We know from Equation 3 that for the acute angle 
= 0, the 
increases with increasing 
. Similar conclusion was reached by Burgess (1987). These authors, by employing 1D hybrid simulations, showed that (a) the field-aligned beam velocity increases and (b) its density decreases strongly with increasing 
for 
. This interplay between 
and 
is probably the reason why the amplitudes of the upstream waves exhibit a complicated dependence on the shock's 
, 
in Figures 1 and 2.
4 Case Study
The spectra presented in the previous section suggest that as the upstream waves cross the shock their properties are at least partially modified. Here we look at a sequence of outputs produced by the run e with 
= 7.08 and 
= 
in order to analyze what happens to the waves as they encounter the shock. This run was chosen since it produces a strongly rippled shock. Figures 7 and 8 exhibit a sequence of excerpts of outputs produced at times indicated on each panel. The former Figure shows the magnetic field magnitude while the latter exhibits the 
component.
Excerpts from the outputs of our run i at five different moments. The colors represent the B magnitude.
Excerpts from the outputs of our run i at five different moments. The colors represent the 
component.
We first look at Figure 7. It is quite obvious that the shock front is rippled and its appearance changes between panels. On panel (a) we observe a compressive structure in the upstream region, shaded with the light blue color, located approximately between y = (−20, 30) 
. This structure, the equivalent of steepened ULF waves or compressive structures such as shocklets, is already interacting with the shock at y 
(−0, 20) 
. Short wavelength fluctuations, equivalent of the whistler waves at real shocks, are present close to this interaction region.
At this location, the shock is locally concave (yellow arrows), i.e.that is, this portion of the surface is located slightly to the right compared to the rest of it. By t = 225 
(panel b) the region of the most intense interaction between the compressive upstream structure and the shock expanded to y 
(−19, 35) 
. The elevated magnetic field magnitude in this section, marked with the white trace between y 
(−3, 30), means that the upstream structure had steepened close to the shock. This interaction region is also permeated by small wavelength fluctuations, equivalent to the whistler mode waves commonly observed upstream of collisionless shocks. By t = 230 
(panel c) almost all of the upstream structure had passed into downstream. The shock is now much more locally concave than during previous times. At y 
−20 
the shock is convex and there is another upstream structure approaching it. We can see in panels (d) and (e) that the story now repeats on somewhat smaller scale as the shock surface between y 
(−20, −40) 
becomes concave.
We can see on panels (a), (b), and (c) that as different sections of the arriving wave interact with the shock, the concave section of the shock surface propagates along the shock front. In this case, it moves downward.
The transverse upstream fluctuations approaching the shock are shown in Figure 8. The colors represent the 
component. It can be seen that these fluctuations do not cause the shock rippling but are affected by it as different portions of a single wavefront reach the shock at different times and locations along the shock surface. Since the shock properties, such as local normal orientation and the shock strength, vary with location, the interaction of different portions of upstream waves is different and this leads to fragmentation of waves as they cross the shock.
The shock transition region (dark red trace at and immediately behind the shock transition in Figures 7 and 8) seems fairly complex. Regions of different magnetic field strength are mixed and tend to be oriented vertically, except behind the ripples, where they look much more disturbed. The wavefronts do not seem to conserve their shapes in these regions.
Farther downstream the magnetic field magnitude is lower than in the shock transition layer. Locations of magnetic field magnitude of different strength (blue, white, orange and red traces) are mixed together and exhibit filamentary structure with their orientations becoming increasingly vertical further downstream. They do not look similar to the shapes of the upstream waves. They exhibit much smaller 
, their sizes are smaller and their orientation is different. This is true for 
B
and 
fluctuations.
5 Discussion and Conclusions
We performed 11 2.5D local hybrid simulations of collisionless shocks in order to study the interaction of upstream ULF waves with the shocks. The shock's Alfvénic Mach numbers range between 4.29 and 7.42, while their nominal angles 
= 
, 
, 
, and 
. We divide our study in two parts: we characterize the degree of rippling of the simulated shocks in the first part and the properties of the magnetic fluctuations in the second.
We acknowledge that our shocks are much simpler than real ones, such as the bow-shock of Earth. For one, in the case of the bow-shocks, the upstream ULF waves may be a mix of waves that are produced upstream of the sections of the bow-shock that have different nominal geometries and Mach numbers and produce waves with different properties. Also, due to the draping of the IMF lines around the magnetopause the regions deep in the magnetosheath may be magnetically connected to the regions of the bow-shock that are much farther away than is the case in our simulations. Another difference is the propagation direction of the downstream waves. For example, Narita and Glassmeier (2006) and Narita et al. (2006) showed that magnetosheath low frequency waves propagate toward the flank magnetosheath for small zenith angles (close to the Sun-Earth line), while they propagate toward the magnetopause for large zenith angles. This ordering will not be observed in local simulations. Finally, in the case of planetary bow-shocks, the waves have much more time to grow in the corresponding foreshock regions. They also develop over larger distances from the bow-shock before being convected back to the shock.
Past works have shown that besides the ULF waves, there are also transient structures in the foreshock, such as foreshock cavitons (e.g., Blanco-Cano et al., 2009, 2011; Kajdič et al., 2011, 2013) and spontaneous hot flow anomalies (SHFAs, Omidi et al., 2013, Zhang et al., 2013). The latter form when the foreshock cavitons approach the bow shock which results in further local ion heating and plasma and magnetic field depletion. Omidi et al. (2014) showed that SHFAs form for 
5. These structures contribute to the shock rippling. To our best knowledge, the relative importance of SHFAs versus ULF waves has not yet been quantified. Kajdič et al. (2013) showed that spacecraft on average observe 
2 cavitons per each 24 h spent in the foreshock. Hybrid simulations (e.g., Omidi, 2007) indicate that for each caviton detected by the spacecraft several more are formed that go undetected but that also contribute to shock rippling. On the other hand, the observations reported by Kajdič et al. (2013) show that for each foreshock caviton, tens or even hundreds of compressive ULF waveforms are detected. This may indicate that the ULF waves influence shock rippling more than the SHFAs. Here we do not specifically search for the cavitons and SHFAs, however eight of our 11 shocks have 
5. The visual inspection of Figure 2 does not reveal any obvious candidates for SHFAs. The apparent absence of these structures may be due to the fact that in our simulated shocks the upstream region is quite limited in size, so that the compressive ULF waves have less time to grow than is the case for the real bow-shock. This results in lesser steepening of these waves. Since according to Omidi (2007) the foreshock cavitons form due to nonlinear interaction of compressive and transverse ULF waves, less wave steepening may result in a smaller number of SHFAs in our simulations.
In Section 3.3.2, we showed that the wavelength of the upstream ULF waves depends on 
and 
and that the dependence on the latter is much stronger. We also showed that the amplitude of these waves depends on the product of the reflected beam density and velocity, both of which also depend on the same two shock parameters. Hence, the ULF wave properties are a complex function of 
and 
, as can be seen in Figures 1 and 2. If ripple properties depend on the properties of the upstream waves, we may expect similarly complex dependence of the degree of rippling on 
and 
.
We determine how rippled the shocks are by studying how certain characteristics of the ripples (amplitudes, length scales, and the steepness of the ripple sides) vary between runs. We characterize the ripple amplitudes with the dispersion 
of the distributions of distances of points 
x on shock surfaces from the mean shock locations (Figures 4a, 4b, and 5a). The spatial scales of the ripples are studied with the Fourier spectra of the curves that we use to describe the surfaces (Figures 4c, 4d, and 5b). Finally, the 
plots describe the steepness of the sides of the ripples (Figures 4e, 4f, 5c, and 5d). Based on these statistics we divide the shocks into three groups: the weakly rippled shocks with small 
, mean 
, 
and Fourier power and the strongly rippled shocks for which these quantities are larger. The runs c and g are considered intermediately rippled since they exhibit low values of integrated power (Figure A2).
There is a correlation between the upstream waves and the shock ripples. Many short-wavelength waves produce many smaller-scale ripples. These exhibit histograms with small standard deviations (
) and many narrow peaks in 
plots (see panels 1a, 2a, 4a, and 5a). On the other hand few long wavelength upstream waves tend to produce less large-scale ripples and result in fewer thicker peaks in 
plots (compare panels 1i, 2i, 4b, and 5a). Large (small) amplitude waves tend produce shock ripples with larger (smaller) amplitudes and histograms with larger 
.
We also compare properties of upstream and downstream fluctuations. Most upstream spectra exhibit flattening or enhancements in the ULF range. These features are conserved in the spectra of the downstream compressive fluctuations. Flattening is better conserved in the case of weakly rippled shocks. Downstream transverse fluctuations however exhibit spectra that are more turbulence-like with flattening and/or enhancements conserved to a much lesser. This is probably due to the fact that downstream fluctuations are mostly compressive.
By examining a case study (Figures 7 and 8), we see that it is the compressive waves that locally bend the shock surface, creating concave ripples. However, since compressive fluctuations are not aligned with the shock surface, different portions of these waves interact with different sections of the shock at different times. As a consequence, the ripples follow this interaction region and thus travel along the shock surface. This also means that different portions of the waves encounter the shock with different properties (magnetic field jump, local normal orientation, etc.). This results in the loss of the wave's identity. Downstream fluctuations do not resemble the upstream waves in wavelengths, wavefront extension, and orientation nor amplitude. This is the reason why the transverse downstream Fourier spectra Figures 4i and 4j are much more similar to turbulent spectra. Even in the case of the compressive downstream spectra Figures 4g and 4h the features they exhibit are less prominent than in the case of upstream spectra.
The fact that flattening and bumps in the ULF range are conserved may be important in the context of ground magnetic pulsations, such as those in Pc2, Pc3 and Pc4 (100–200 mHz, 22–100 mHz, and 7–22 mHz, respectively) bands. Their appearance is known to be associated with the IMF orientation and their origin is commonly attributed to the foreshock ULF waves since both exhibit similar frequencies (e.g., Clausen et al., 2009; Engebretson et al., 1987; Verõ et al., 1998). The present work shows that some spectral features of compressive fluctuations in the ULF range may indeed survive in the magnetosheath and perhaps affect the magnetopause where they could act as pressure pulses exciting magnetospheric Pc fluctuations. Although local runs provide valuable information on the nature of the interaction between the shocks and the upstream ULF waves, further investigation with global hybrid runs is required in order to study the transmission of such fluctuations from the bow-shock all the way to the magnetopause.
Acknowledgments
The authors acknowledge the International Space Science Institute for supporting the ISSI team 448 on “Global study of the transmission of foreshock ULF waves into the magnetosheath and the magnetosphere” led by L. Turc and M. Palmroth, in which this study was initiated. P. Kempf's work was supported by DGAPA/PAPIIT grant IN105620. The work of N.Takahashi was supported by Grants-in-Aid for Scientific Research of Japan Society for the Promotion of Science (grant number: 16H06286, 18KK0099). K. Takahashi’s work was made supported by the NASA Grant NNX17AD34G. The work of L. Turc is supported by the Academy of Finland (grant number 322544). E. Kilpua acknowledges the ERC under the European Union's Horizon 2020 Research and Innovation Programme Project SolMAG 724391, and Academy of Finland Project 310445.” X. B. Cano’s work was supported by UNAM DGAPA PAPIIT IN-105218-3 grant. A. P. Dimmock received financial support from the Swedish National Space Agency (Grant #2020-00111) and the EU Horizon 2020 project SHARP, SHocks: structure, AcceleRation, dissiPation # 101004131. Y. P. Kempf’s work was supported by the ERC Consolidator grant no. 682068-PRESTISSIMO. L. Preisser's work is supported by the Austrian Science Fund (FWF): P 33285-N.
Appendix A
Here we show the outputs produced by all our models at simulation times when the shocks were located at 
. The outputs are stacked so that the 
increases from left to right. The models in the top row are produced with the inflow velocity 
, those in the middle row with 
while in the bottom row the 
. The 
increases mostly from top to bottom, although is also increases slightly with 
(from left to right). Figures A1 to A5 show: Histograms for distributions of x-coordinates along the shock surface with respect to its average position, Fourier spectra of curves representing the shock surfaces, 
angles between the normals of all pairs of points along the shock surface as a function of distance between the points, upstream (black) and downstream (red) Fourier spectra of 
fluctuations at times shown in Figure 2 and upstream (black) and downstream (red) Fourier spectra of 
fluctuations at times shown in Figure 1, respectively.
Histograms for distributions of x-coordinates along the shock surface with respect to its average position. Shaded panels represent what we classify as weakly rippled shocks.
Fourier spectra of curves representing the shock surfaces. Vertical blue line marks the location of 
= 40 
, while the red line marks the trends at 
40 
. The integrated power is calculated by calculating the total power of the spectra at 
40 
and subtracting from it the power delimited by the red line in the same wavelength range. Shaded panels represent what we classify as weakly rippled shocks.
angles between the normals of all pairs of points along the shock surface as a function of distance between the points. The maximum and median 
and the value S, which is the fraction of the surface delimited by black dots, are also shown. Shaded panels represent what we classify as weakly rippled shocks.
Upstream (black) and downstream (red) Fourier spectra of 
B
fluctuations at times shown in Figure 2. Blue shading indicates the ultralow frequency (ULF) range. Panels shaded in pink represent what we classify as weakly rippled shocks.
Upstream (black) and downstream (red) Fourier spectra of 
fluctuations at times shown in Figure 1. Blue shading indicates the ultralow frequency (ULF) range. Panels shaded in pink represent what we classify as weakly rippled shocks.
Open Research
Data Availability Statement
The simulation outputs from which Figures 2 and 1 were created are available at https://zenodo.org/record/5156128.
