Non-stationarity of the quasi-perpendicular bow shock : comparison between Cluster observations and simulations

We have performed full particle electromagnetic simulations of a quasi-perpendicular shock. The shock parameters have been chosen to be appropriate for the quasiperpendicular Earth’s bow shock observed by Cluster on 24 January 2001 ( Lobzin et al. , 2007). We have performed two simulations with different ion to electron mass ratio: run 1 with mi/me= 1840 and run 2 withmi/me= 100. In run 1 the growth rate of the modified two-stream instability (MTSI) is large enough to get excited during the reflection and upstream gyration of part of the incident solar wind ions. The waves due to the MTSI are on the whistler mode branch and have downstream directed phase velocities in the shock frame. The Poynting flux (and wave group velocity) far upstream in the foot is also directed in the downstream direction. However, in the density and magnetic field compression region of the overshoot the waves are refracted and the Poynting flux in the shock frame is directed upstream. The MTSI is suppressed in the low mass ratio run 2. The low mass ratio run shows more clearly the non-stationarity of the shock with a larger time scale of the order of an inverse ion gyrofrequency ( ci): the magnetic field profile flattens and steepens with a period of ∼ 1.5 ci . This non-stationarity is different from reformation seen in previous simulations of perpendicular or quasi-perpendicular shocks. Beginning with a sharp shock ramp the large electric field in the normal direction leads to high reflection rate of solar wind protons. As they propagate upstream, the ion bulk velocity decreases and the magnetic field increases in the foot, which results in a flattening of the magnetic field profile and in a decrease of the normal electric field. Subsequently the reflection rate decreases and the whole shock profile steepens again. Superimposed on this ’breathing’ behavior are in the realistic mass ratio case the waves due to the MTSI. The simulations lead Correspondence to: H. Comişel (comisel@spacescience.ro) us to a re-interpretation of the 24 January 2001 bow shock observations reported by Lobzin et al.(2007). It is suggested that the high frequency waves observed in the magnetic field data are due to the MTSI and are not related to a nonlinear phase standing whistler. Different profiles at the different spacecraft are due to the non-stationary behavior on the larger time scale.


Introduction
Already the very early kinetic particle-in-cell simulations of collisionless shocks have revealed that quasi-perpendicular shocks can be intrinsically instationary and can exhibit cyclic reformation.This process is due to the fact that above a critical Mach number part of the incoming ions are specularly reflected at the shock front by the cross-shock potential and the magnetic field increase.The specularly reflected ions gyrate in the upstream magnetic field and lead to an extended foot in the magnetic field profile upstream of the shock ramp.During their upstream gyration the ions are accelerated by the motional electric field in the direction perpendicular to the magnetic field and parallel to the shock plane and return eventually to the shock with a velocity nearly twice the original encounter velocity.This is considered as the first step in the ion heating process by the shock.At shocks below a critical Mach number given by M w = 1 2 m i and m e are the ion and electron mass, respectively.Above this critical Mach number a nonlinear whistler may exist.Biskamp and Welter (1972) performed PIC simulations of a shock with Alfvén Mach number M A = 5 and Bn = 45 • and found that the shock was nonstationary.They explained their result in terms of a nonlinear instability between incoming and reflected ion beam if a (spatially) periodic electric field due to the almost phase-standing whistler train is present.The free energy for this instability is the velocity between incoming and reflected ions.The Biskamp and Welter (1972) simulations were done for relatively small m i /m e mass ratios of 64, 128, and 256.For the m i /m e = 64 and 128 runs the shock was supercritical relative to the whistler critical Mach number (the whistler critical Mach number in the m i /m e = 64 and m i /m e = 128 runs was M w = 2.8 and M w = 4, respectively).Scholer and Burgess (2007) have recently performed full particle simulations of oblique shocks with the physical ion to electron mass ratio and have confirmed the Biskamp and Welter (1972) whistler induced reformation (WIR) scenario for the lower Bn regime ( Bn below ≈ 80 • ): in less oblique shocks above the whistler critical Mach number the whistler amplitude in the foot grows, leading to vortices of the incoming ions and of the reflected ions in velocity phase space and eventually to phase mixing and reformation.
The situation is different for more oblique shocks or for an exactly perpendicular shock.From a simulation of an exactly perpendicular shock Lembège and Dawson (1987) have concluded that the shock periodically reforms itself with a time scale of the order of the inverse ion gyrofrequency ci due to accumulation of reflected ions at the upstream edge of the foot.When the density of the reflected ions is large the compression of the magnetic field at the upstream edge of the foot where the reflected ions are turned around can ignite the emergence of a new shock ramp.Hada et al. (2003) have developed a semi-analytical model for this shock reformation process, which is based on the coupling between incoming and reflected ions.They determined a critical fraction of reflected to incoming ions beyond which no stationary solution for the coupling process exists.In other words, for a high reflection rate the coupling becomes so strong that the shock is nonstationary and reforms.As pointed out by Hada et al. (2003) the proposed model applies for a large Mach number range, but is restricted to the case where the reflected ions can be described by a mono-energetic ion population, which is approximately true for low ion beta.As shown by Schmitz et al. (2002), Scholer et al. (2003), and Hada et al. (2003) the reformation process disappears in PIC simulations of higher beta shocks (β i of the order of 1).Scholer et al. (2003) suggested that the important quantity determining whether reformation occurs or not is actually not the ion beta, but the difference between the upstream plasma velocity and the ion thermal velocity: when this difference is large the incoming and reflected ion beam interact at the upstream edge of the foot, and this interaction starts a reformation process.When this difference is small, as in β i > 0.4 simulations of medium Mach number shocks, the incoming and reflected ion beams are overlapping in velocity space in the region close to the shock ramp.The interaction then occurs smoothly over the whole foot and a stationary profile results.
The simulations by Lembège and Dawson (1987) and Hada et al. (2003) were performed for exactly perpendicular shocks.However, as seen from an evaluation of the linear theory for parameters appropriate to the foot of quasiparallel shocks by Matsukiyo and Scholer (2003) and from simulations by Scholer et al. (2003) the modified two-stream instability (MTSI) can get excited and will modify the reformation and heating process for shocks.As the specularly reflected ions propagate upstream away from the shock and constitute a foot, the incoming electrons are decelerated in order to achieve zero electrical current in the shock normal direction.In a not exactly perpendicular configuration, i.e., when a component of wave vectors parallel to the magnetic field is allowed for, the resulting difference between incoming ions and incoming electrons leads to the excitation of the MTSI.The waves generated by this instability are on the whistler mode branch and in the shock frame the k vectors are directed toward the shock.Investigation of this process requires simulations with realistic ion to electron mass ratio, since the growth rate of the MTSI in units of the inverse ion gyrofrequency depends strongly on m i /m e (Matsukiyo and Scholer, 2003): for a mass ratio below ∼ 400 the growth rate is smaller than the ion gyroperiod and the instability can not get excited within one gyromotion of the reflected ions.In the low m i /m e ratio simulations by Biskamp and Welter (1972) the MTSI could thus not get excited.Simulations with a realistic mass ratio have been performed by Scholer and Matsukiyo (2004) and it has been shown by these authors that under certain conditions the MTSI in the foot leads to ion phase mixing and eventually in shock reformation.
A fourth process for shock reformation has been proposed by Krasnoselskikh et al. (2002).Krasnoselskikh et al. (2002) have suggested that above a nonlinear whistler critical Mach number given by the nonlinear upstream whistler train becomes unstable to a gradient catastrophy.This is supposedly due to the fact that the nonlinear steepening of the wave train cannot be balanced anymore by the effects of dispersion and dissipation.The wave train then becomes instable with respect to overturning and reformation results.The analysis by Krasnoselskikh et al. (2002) neglects the modification due to the specularly reflected ions.In this respect the whistler induced reformation by Biskamp and Welter (1972) is different from the Krasnoselskikh et al. (2002) catastrophy model: in the WIR mechanism the free energy source for the instability is the velocity difference between reflected and incoming ions.It should, however, be noted that although the Biskamp and Welter (1972) shock with M A = 5 was supercritical in the m i /m e = 128 run with respect to the whistler critical Mach number (M w ∼ 4), it was subcritical with respect to the nonlinear whistler critical Mach number (M nw ∼ 5.7).
Only for the low mass ratio m i /m e = 64 run does the shock Mach number exceed both the linear (2.8) and the nonlinear whistler critical Mach number (4.0).This with respect to the nonlinear whistler critical Mach number supercritical case could be expected to exhibit the gradient catastrophy; the simulation result of only one whistler train becoming nonlinear and turning into a new ramp may indeed be interpreted in terms of the catastrophy model.The nonlinear whistler catastrophy mechanism was supported by full particle simulations with an ion to electron mass ratio of 200 (Krasnoselskikh et al., 2002).However, considering for example a Bn = 80 • shock such a low mass ratio leads to the unrealistic result that the first critical Mach number (above which the shock reflects particles in order to achieve dissipation) is larger than both the whistler critical Mach number and the nonlinear whistler critical Mach number.In other words a realistic mass ratio has to be used in order to investigate at large values of Bn the change of the shock behavior when increasing the Mach number from the whistler critical Mach number through the nonlinear whistler critical Mach number.
We should like to address briefly the recent discussion on reformation in one-dimensional (1-D) and 2-D simulations.Hellinger et al. (2007) claimed on the basis of 2-D hybrid and full particle simulations of exactly perpendicular shocks that shock reformation is suppressed by an oblique whistler in the foot which is not seen in 1-D simulations.Lembège et al. (2009) also reported full particle simulations of an exactly perpendicular shock with a mass ratio of 400 and showed that whistler waves are excited in the foot which inhibit shock reformation.This is seemingly at variance with 2-D hybrid simulations by Yuan et al. (2009) who showed that for not exactly perpendicular shocks reformation still prevails, although the reformation periodicity is slightly changed due to the whistler.This discrepancy has recently been addressed by Umeda et al. (2010).These authors have shown that the shock magnetic field averaged over the tangential direction does indeed not exhibit reformation cycles.However, the time evolution of the local shock magnetic field does also in 2-D simulations exhibit periodic oscillations.
From the four-spacecraft Cluster mission there is some evidence for shock nonstationarity at Earth's quasiperpendicular bow shock.The magnetic field observations at the four spacecraft which were several hundred kilometers apart revealed structures in the foot of the quasiperpendicular bow shock (Horbury et al., 2001).The highly localized magnetic field activity in the foot was not phasestanding but was convected into the shock.During other quasi-perpendicular bow shock crossings Horbury et al. (2002) found vastly different magnetic field profiles at the four spacecraft so that they actually had to discard these ex-amples for a determination of the bow shock orientation and motion.Mazelle et al. (2009) have recently presented evidence for reformation of the quasi-perpendicular bow shock by a statistical analysis of Cluster bow shock crossings.Lobzin et al. (2007) investigated in detail a Cluster bow shock crossing on 24 January 2001.This bow shock was quasi-perpendicular with Bn ≈ 81 • and an Alfven Mach number of M A ≈ 10.The magnetic field profiles at the four spacecraft differed considerably from each other.Oscillations in the frequency range from 3 to 8 Hz were superimposed on the large scale structure and are interpreted as a whistler wave train nested in the shock.After removing fluctuations with frequencies higher than 2 Hz by Fourierfiltering one or two large amplitude structures with a characteristic time of about 2 s survived, and it was concluded from cross-correlation studies that the structures were temporal and not spatial.This has been interpreted by Lobzin et al. (2007) as an indication for whistler induced reformation: the observed peaks in the overshoot region are a part of nonstationary whistler wave packets in the Krasnoselskikh et al. (2002) reformation mechanism.Additional support for this interpretation comes from the temporal behavior of the reflected ions.The count number of reflected ions exhibits periodic variations with a period of about 8 s.This has been interpreted by Lobzin et al. (2007) in terms of the Krasnoselskikh et al. (2002) mechanism: when the leading wave train attains a large enough amplitude a new population of reflected ions appears upstream of the precursor.
We report in the following on one-dimensional full particle electromagnetic simulations of a quasi-perpendicular shock with parameters appropriate for the 24 January 2001 Cluster bow shock crossing.In short, we can reproduce at artificial spacecraft time series like those observed by Cluster, but we can not confirm the interpretation in terms of a nonlinear whistler induced reformation mechanism.Quite to the contrary, the simulated shock is not reforming, but exhibits a "breathing-like" temporal behavior with a periodically steeper and flatter ramp.The high frequency waves in the foot are due to the modified two-stream instability and the differences in the large amplitude structures at different spacecraft are due to the differences in spatial sampling of the periodically changing shock profile.

Simulation
The shock is produced by the so-called injection method: a high-speed plasma consisting of electrons and ions is injected from the left hand boundary of a one-dimensional simulation system and travels toward positive x.The plasma carries a uniform magnetic field which has a B z and a B x component.At the right hand boundary the particles are specularly reflected.A shock then propagates in the −x direction, i.e., the simulation system is the downstream rest frame, and the shock normal is the x-axis.Furthermore, the simulations are  done in the so-called normal incidence frame where the upstream bulk velocity is parallel to the shock normal.Initially there are 100 particles for each species in a computational cell; the simulation box consists of 40 000 grid cells, each of the size of one Debye length λ D .CGS Gaussian units are used in the code.In the following, time will be given in units of the inverse of the ion cyclotron frequency ci , distances in units of the electron inertial length λ e = c/ω pe (c = velocity of light, ω pe = plasma frequency), the velocity in units of the upstream Alfvén speed v A , magnetic field and the density in units of their upstream values B 0 and n 0 , respectively.The electric field unit is B 0 ; thus the potential is given in units of B 0 λ e (there is a typo concerning this unit in Scholer et al., 2003, andScholer andMatsukiyo, 2004).We use Bn = 81 • and inject the plasma with Alfvén Mach number M A0 = 6.0 which leads to a M A ∼ 9 shock.The ion beta during the 24 January 2001 bow shock crossing was given by Lobzin et al. (2007) as β i ≈ 2.0.This quantity has been recently reanalyzed and is closer to ∼ 0.6 (C.Mazelle, private communication, 2010).We have used in the simulations β i = 0.6 and β e = 1.7.Two important parameters enter a PIC simulation: (1) the mass ratio m i /m e and (2) the ratio of the electron plasma frequency to gyrofrequency ω pe / ce .The latter parameter is in the solar wind near the Earth's orbit 100-200.It is not possible at present to have both parameters realistic values.We present mainly results for the physical mass ratio and for ω pe / ce = 8.Note that this value is considerably larger than what has been used by, for example, Scholer and Matsukiyo (2004) (ω pe / ce = 2), but smaller than the value used by, for example, Shimada and Hoshino (2005) (ω pe / ce = 20).Note, however, that the latter authors had to compromise by using a mass ratio of 20.In addition, we will present a low mass ratio (m i /m e = 100,ω pe / ce = 8) simulation where the modified-two-stream instability is sup- pressed.For both m i /m e ratios used the shock is supercritical with respect to the linear and to the nonlinear whistler critical Mach number.

Realistic mass ratio, low ω pe / ce shock simulation
We report now results from the m i /m e = 1840, ω pe / ce = 8 shock simulation.Figure 1 shows the color-coded magnetic field B z component in the t − x plane.Time runs from 5 to 13 −1 ce and the scale on the x-axis is in electron inertial lengths λ e .The left hand side of the x-axis is arbitrarily set to zero.The magnetic field has been transformed from the simulation frame (downstream rest frame) into the shock frame by determining an average shock velocity over the whole time period.The magnetic field exhibits temporal variations on two different time scales: on a small time scale there are small wave-length waves originating upstream with a downstream directed phase velocity, and on a larger time scale of ∼ 1.8 −1 ci the whole shock structure changes periodically.This periodic change manifests itself (1) in a small forward and retreating motion of the shock ramp, (2) in a periodic excitation of the waves upstream, and (3) in a periodic change of the wave amplitudes behind the shock (maximum magnetic field intensity indicated by red color).This can be more clearly seen from Fig. 2 which shows the time evolution of the magnetic field in higher resolution over a smaller time period.The upstream waves are excited far upstream when the ramp magnetic profile is rather flat, for instance at t ci = 9.3 ; conversely, during times of a large downstream magnetic field (at t ci = 10.3)upstream waves exist only very close to the ramp.Furthermore, we note that the phase velocity of the waves increases considerably in the ramp and in the downstream region as compared to the upstream phase velocity.
Figures 3 and 4 show, from top to bottom, the magnetic field B z component, the ion density n i , and ion v ix −x phase space at two specific times (t ci = 9.5 and t ci = 10.3, respectively).At t ci = 9.5 (Fig. 3) there is no clearly defined sharp ramp, rather an extended foot exists with specularly reflected ions and strong magnetic field wave activity with a wavelength of ∼ 12λ e .Further upstream the waves lead to local decelerations of the ions and a simultaneous density increase; the amplitude of these structures grows toward the shock and results in vortices where the incoming ions and the reflected ions are phase-mixed.These vortices can be considered as small-scale reformation cycles.At t ci = 10.3 (Fig. 4) the shock has a considerably shorter foot and a sharper ramp.Ions just begin to be reflected from the ramp and wave activity in the foot is low; most of the (now) large amplitude waves are in the overshoot; these waves originate from the foot region at earlier times.The shock changes periodically between two states: a phase 1 with a rather well defined ramp and overshoot with large amplitude waves, which are remnants from previous small amplitude waves in the foot, and a phase 2 with an extended foot with specularly reflected ions leading to the excitation of smallamplitude waves in the foot.During phase 1 a part of the incoming ions are specularly reflected; this reflection ceases during phase 2 since the shock is less sharp and, as will be shown below, the maximum of the electric field in the shock normal direction, E x , is considerably smaller than in phase 1.

Linear instability analysis
The upstream waves are most probably not related to a shock produced upstream whistler, but due to the modified twostream instability (MTSI).During phase 1, part of the incoming ions are specularly reflected, propagate upstream away from the shock, and constitute a foot.In this region the incoming electrons are decelerated in order to achieve zero electrical current in the shock normal direction.In a not exactly perpendicular configuration, i.e., when a component of wave vectors parallel to the magnetic field is allowed for, the resulting difference between incoming ions and incoming electrons leads to the excitation of the MTSI.The waves generated by this instability are on the whistler mode branch and slightly modify the shock density and magnetic field profile.We have solved the linear dispersion relation for a situation appropriate to the conditions in the foot region during period 2. For the instability analysis it is assumed that there are three components: solar wind ions, reflected ions and solar wind electrons, which are all represented by drifting Maxwellians.Since we consider waves with much higher frequencies than the ion gyrofrequency the ions are assumed to be unmagnetized.Details of the linear instability analysis can be found in Matsukiyo and Scholer (2003).As will be shown below, the density of specularly reflected ions is highly variable between E ∼ 15% and ∼ 30%.For the linear analysis it has been assumed that 25% of the ions are specularly reflected, i.e., the analysis is performed for times of high reflected ion density, and the ion and the electron beta are assumed to be the same as far upstream.As can be seen from Fig. 3 the solar wind velocity decreases in the foot and the magnetic field B z increases.The latter means that the magnetic field angle relative to the reflected ion and incoming ion beam increases from the nominal value of Bn = 81 • .We have therefore calculated the growth rate for two different solar wind velocities and for a number of Bn values.Figure 5 shows the growth rate γ / i of the MTSI as function of k for different angles Bn .The upper panel shows the result for a solar wind velocity of M A = 9 and the lower panel shows the growth rate for M A = 5.Although the ion beta in this case is rather large, the growth rate is positive.The proton bulk speed in the foot is close to M A = 5.Assuming in the foot a value of Bn = 86 the linear analysis results in a k value of maximum growth of kc/ω pe ∼ 0.5.This corresponds to a wavelength of ∼ 12λ e which is in accordance with the wavelength found in the foot in the shock simulation.Conditions in the foot vary as a function of distance and, in particular, as will be shown later, as a function of time.The instability analysis should thus be done for a large variety of parameters appropriate to an individual location and time.However, our aim is to demonstrate that the MTSI growth rate is for parameters within the parameter range positive, and that the instability can in principle grow within one ion gyroperiod of the specularly reflected ions.Additional support for the MTSI origin of the waves in the foot comes from the shock directed phase velocity (see Fig. 2) and from a Poynting flux analysis described in the following subsection.
In addition to the MTSI based on the interaction between incoming electrons and incoming ions discussed above a second MTSI based on the interaction between incoming electrons and reflected ions is expected, which is also on the whistler mode branch.In a 2-D simulation of a set-up appropriate for the situation in the foot of a shock both instabilities have actually been found (Matsukiyo and Scholer, 2006).These authors have termed the MTSI based on the interaction between incoming electrons and incoming ions MTSI-1 and the MTSI between incoming electrons and reflected ions MTSI-2.The waves with maximum growth of the MTSI-1 have wave vectors almost perpendicular to the magnetic field and can thus be observed in 1-D shock simulations.However, the waves with maximum growth of the MTSI-2 propagate oblique to the magnetic field (e.g., Gary et al., 1987) and will thus not be seen in 1-D simulations, which allow only for k vectors in the shock normal direction, i.e. almost perpendicular to the magnetic field.The occurrence of the MTSI-2 can change the picture obtained from 1-D simulations; however, the linear stability analysis by Gary et al. (1987) has shown that the MTSI-2 is only dominant in the low electron beta case (β e < 0.5), while the present β e ∼ 1.7.We thus expect the MTSI-2 to be of minor importance.

Poynting flux determination
As the waves approach the shock and travel into the region of compressed magnetic field and density their amplitude and the phase velocity increases.In the overshoot, the phase velocity is still directed downstream.We have determined the Poynting flux direction in these waves for two different time periods.Figure 6 shows B z magnetic field profiles stacked in time.We have followed four magnetic field wave crests from the foot region into the shock overshoot.These wave crests are indicated at two different times by four heavy dots.The Poynting flux parallel to the shock normal was calculated using continuous wavelet transform (CWT) to resolve the local Poynting vector corresponding to the short wavepackets in the foot and ramp region of the shock.Electric and magnetic field components E i (x) and B i (x) were transformed in the spatial coordinate x using a complex Morlet wavelet for a set of 78 linearly spaced wavenumbers kc/ω pe from 0.2 to 1.8 of inverted grid step.The CWT was applied to a snapshot of 8192 data points centered at the shock ramp, but to avoid edge effects, only the 3000 points in the center of the interval were used in the analysis.
From the wavelet transformed field vectors Ẽ(x,k) and B(x,k) the Poynting flux for each x and k is calculated as where the asterisk denotes a complex conjugate.Since the electromagnetic fields were transformed from the simulation frame into an average shock rest frame, the sign of the component of the Poynting vector parallel to the shock normal indicates the direction of propagation of the wavepacket at the corresponding position and wave number relative to the shock front.In our coordinate system, positive values correspond to downstream propagation and vice versa.The result for the two time periods under investigation is shown in Figs.7 and 8. Shown is color coded the wavenumber versus x.Also shown is a trace of the B z profile during the respective time.The color coding indicates the direction of the Poynting flux (see color bar to the right).The four wave crests indicated in Fig. 6 by solid dots have been indicated on each of the magnetic field traces.Figure 7 for t ci = 9.2 shows that the wavenumber of the wavetrain indicated by the solid dots is kλ e ≈ 0.5, corresponding to a wavelength of about 12λ e .According to the color coding (red/yellow color) the Poynting flux is directed downstream.When the wave packet moves into the high density, high magnetic field region the group velocity changes and is upstream directed, as can be seen from the blue region in Fig. 8 at t ci = 9.5.This can be interpreted in terms of refraction of the waves when they propagate into the high density,  high magnetic field compression region of the shock.The change of the Poynting flux direction/group velocity during wave propagation from upstream in the foot into the ramp and overshoot shows that an interpretation of the wave characteristics from measurements at one particular time can be misleading: the wave packet with the largest amplitude, i.e., the waves in the compression region have an upstream directed group velocity and can be misinterpreted as shock produced whistlers, whereas the waves actually originate in the upstream in the foot with downstream directed group velocity.

Low mass ratio (100), low ω pe / ce shock simulation
In order to investigate the periodically steepening and flattening of the shock ramp we have performed a simulation with a mass ratio m i /m e = 100.According to linear theory at this low mass ratio the MTSI should not get excited.Figure 9 shows the color-coded magnetic field B z component in the t −x plane over a time period of about about 8 −1 ci .As can be seen, upstream and downstream small wavelength waves are indeed absent (there is some wave activity far downstream), which shows that the MTSI is not excited.There is also no indication of a nonlinear whistler precursor.However, the periodically flattening and steepening of the shock ramp still prevails and is thus independent from the occurrence of the MTSI.However, the periodic flattening and steepening is different from reformation as described by, e.g., Hada et al. (2003) and Scholer et al. (2003): a new shock ramp does not build up at the upstream edge of the foot, but rather emerges again out of the flatter profile.This non-stationarity is better described by a "breathing" behavior of the shock with a periodically changing steepness of the ramp.This can be seen in Fig. 10  downstream rest frame, whereas Fig. 9 shows the color coded magnetic field in the shock frame.There is no new shock ramp emerging near the upstream edge of the foot, rather the existing ramp steepens and flattens periodically.Comparing Figs. 6 and 10 for the same shock parameters but for different ion to electron mass ratio also demonstrates drastically the importance of realistic mass ratio simulations for a comparison of simulations with in situ measurements.The lower mass ratio run allows determination of the ramp scale.We have defined the ramp scale L r as the distance between the position of the maximum field strength in the overshoot and the beginning of the foot.The latter is defined by the intersection of the magnetic field gradient in ramp and foot.At times of a steep ramp the length scale is L r ≈ 6λ e and at times of a flat ramp the lengths scale is L r ≥ 28λ e , i.e., the scale length of the ramp exhibits temporal variations by a factor 4-5.The foot exhibits similar temporal variations, so that the scale of the total shock transition varies by about a factor 4. Figure 11 shows in the top panel the average reflected ion density n ref upstream of the ramp.The density is determined by counting reflected ions from the maximum of the magnitude of the electric field x-component, |E x |, within a distance of 60 electron inertial lengths upstream.Ions in this region are defined as being reflected when they have an upstream directed x-component of the velocity larger than 5M A .By this criterium we eliminate ions which have been reflected at earlier times and are gyrating back into the shock.The reflected ion density in the region upstream of the shock (as defined by the maximum of the upstream directed normal electric field component) changes periodically by a factor three with a period of ∼ 1.6 −1 ci .The bottom panel of Fig. 11 shows the electric field component E x at the position where |E x | has the maximum value.The electric field changes periodically with the same periodicity as the reflected ion intensity.Comparing the bottom and top panel of Fig. 11 shows that when the upstream directed (negative) normal electric field component increases (larger negative value) the reflected ion density begins to increase, and vice versa, after a minimum of the upstream directed x-component of the electric field the reflected ion density begins to decrease.The electric field in the code given in units of B 0 in statV/(cm G) can be converted into the SI unit mV m −1 by multiplying with the conversion factor ς ≈ 300B up , where B up is the upstream magnetic field in nT.With B up ∼ 4 nT this leads to electric fields of |E x | ∼ 240 mV m −1 , which is rather large.This is due to the unrealistically small value of ω pe / ce used in the simulations.Note that the electric field is proportional to the upstream magnetic field.However, since ω pe / ce is constant, an increase in the upstream magnetic field results in a corresponding decrease of the electron inertial length, so that the potential stays the same.direction leads to high reflection rate of solar wind protons.
As they propagate upstream the ion bulk velocity decreases and the magnetic field increases in the foot, which results in a flattening of the magnetic field profile and at the same time in a decrease of the normal electric field.Subsequently the reflection rate decreases and the shock steepens again.Since this process is connected with the upstream propagation and gyration of the reflected ions the natural time scale is the inverse ion gyrofrequency.We note that a shock scenario with a periodic change between a steep profile and a flat profile has been proposed by Krasnoselskikh (1985).In this model the steep profile causes wave breaking, and a subsequent flat profile is due to dispersive broadening.

Comparison with Cluster observations
For a comparison of in situ spacecraft data with the simulation we have constructed artificial time series of magnetic field data: two probes separated by some distance along the shock normal move with constant speed from upstream to downstream through the high mass ratio 1-D simulation.Due to the small total time for which we can follow the shock development the distance between the two probes is rather small (∼ 50c/ω pe ).However, due to the periodicity of the shock profile the result for this small separation distance can be taken as a good representation for any probe separation.
Note that the separation used here is about the thickness of the shock ramp.The velocity of the probe relative to the shock is assumed to be ∼ 0.8M A . Figure 12 shows in the top panel the magnetic field B z component recorded by probe 1 (S/C 1).The profile exhibits a rather complicated foot, a ramp with overshoot, and a second downstream maximum.
The temporal foot profile is due to the fast changing spatial profile during the crossing of the probe.Superimposed on the large scale profile are the small wavelength waves due to the MTSI.The temporal profile shown in the top panel has then been high-pass filtered; the result is shown in the middle panel.The bottom panel shows the high-pass filtered profile obtained by the second probe S/C 2. This probe has passed the shock during a time with a flatter ramp and smaller overshoot.The two probes different amplitudes of the two peaks which are separated by about 3 − 4 −1 ci .This result compares favorably well with the temporal magnetic field profiles shown by Lobzin et al. (2007) for a Cluster bow shock crossing on 24 January 2001 (see their Fig. 1).Lobzin et al. (2007) furthermore demonstrated that during this event the intensity of reflected ions changed periodically with a periodicity of 8 s.We have already seen that in the low mass ratio simulation the intensity of reflected ions changes periodically with ∼ 1.5 −1 ci .Figure 13 shows as in Fig. 11 in the top panel the average reflected ion density in the foot and in the bottom panel the maximum of the upstream directed (negative) normal electric field.The sampling region for the reflected ions is now 150λ e as compared to 60λ e in the mass ratio 100 case.Note that the ion and electron inertial lengths are related by λ i = λ e √ (m e /m i ). Figure 13 shows that the density of reflected ions and the normal electric field exhibit periodic variations with a periodicity of ∼ 2 −1 ci .During the 24 January 2001 Cluster bow shock crossing the ion gyroperiod was ∼15 s (Lobzin et al., 2007).Thus the unit time in the simulation of one inverse ion gyrofrequency corresponds to 2.5 s and the periodicity of the reflected ions is according to the simulation ∼5 s in real time which compares favorably well with the 8 s periodicity reported by Lobzin et al. (2007).

Summary and discussion
We have performed full particle simulations of a high Mach number quasi-perpendicular collisionless shock.The shock parameters, i.e., the ion beta (β i = 0.6) , the electron beta (β e = 1.7), the Mach number (M A = 9), and the magnetic field-shock normal angle ( Bn = 81 • ), have been chosen to be similar to the ones observed by Cluster during a bow shock crossing on 24 January 2001.Two simulations have been performed: run 1 with the physical ion to electron mass ratio and run 2 with an ion to electron mass ratio of 100.The results can be summarized as follows.
High frequency waves in the foot and ramp with wavelengths of ∼ 12λ e are excited and have a downstream directed phase velocity.Based on linear theory these waves are identified as being due to the modified two-stream instability (MTSI).In the shock frame the group velocity of small amplitude waves in the foot is directed downstream; as the wave amplitudes increase when they propagate into the high density compression region of the ramp and overshoot the waves are refracted and the group velocity in the shock frame is directed upstream.There is no indication of a phase standing linear or nonlinear whistler precursor produced by the shock ramp itself.
In addition to the high frequency waves in the foot and ramp the shock exhibits a nonstationarity on a time scale of 1 − 2 ci .This nonstationarity is different from shock reformation where a new ramp emerges at the upstream edge of the foot due to accumulation of reflected ions (Hada et al., 2003).The nonstationarity found in the present high Mach number, higher ion beta quasi-perpendicular shock is better described by "breathing" behavior of the shock with a periodically flatter and steeper ramp.In the m i /m e = 100 run the MTSI is suppressed.This allows particularly well to analyze the nonstationary behavior of the large scale shock profile.Specular reflected ions result in a deceleration of the solar wind in the foot, an increase of the magnetic field in the foot and to a flattening of the shock ramp.Combined with the resulting decrease of the normal electric field this leads to a decrease of the ion reflection rate.After the reflected ions gyrate and move downstream the shock ramp steepens again and the cycle repeats itself.The periodic flattening and steepening of the magnetic field profile is different from the reformation process where a new shock emerges at the upstream edge of the foot.
The question arises why the present result of a breathing behavior of the shock front differs from either the classical picture of shock reformation by reflected ion accumulation at the upstream edge of the foot (Hada et al., 2003) or from reformation by phase mixing of reflected and incoming ions due to the MTSI as proposed by Scholer and Matsukiyo (2004).Responsible for the differences in shock behavior is clearly the variety of shock parameters, as Bn , β i,e , and M A .First, we note that in 1-D simulations of exactly perpendicular shocks the MTSI can, per design, not occur, and it has to be seen whether in 2-D simulations of exactly perpendicular shock the MTSI occurs and would lead to reformation by subsequent phase mixing.Lembège and Savoini (1992) have actually performed 2-D simulations of exactly perpendicular shocks; however, the low ion to electron mass ratio in these simulations does not allow the excitation of the MTSI, so that no conclusion can be made.Why then does the MTSI in the present case not result in phase mixing as in the Scholer and Matsukiyo (2004) simulations?The MTSI growth rate depends on the reflected ion density, on the ion beta β i , and on Mach number.The Scholer and Matsukiyo (2004) simulations were done for a rather small β i of 0.1 and a shock with Mach number M A = 4.5.In this case the reflection rate is high (≈ 50%) and the growth rate for the MTSI is large and of the order of γ ∼ 20 ci (Scholer et al., 2003).
The large growth rate leads to large amplitude waves already well upstream in the foot and to subsequent phase mixing.At M A ≈ 4−5 and β i ≥ 0.4 Scholer and Matsukiyo (2004) found that reformation ceases and the shock is smooth and steady.In the present case β i = 0.6, reflection rate is lower, and at a low Mach number the MTSI would be suppressed.However due to the large Mach number of M A ∼ 9 the MTSI growth rate is positive and of the order of γ ∼ 5 ci (Fig. 5), i.e. by a factor 4 smaller.As can be seen from Fig. 3 a large region in the foot well upstream of the ramp is actually void of MTSI produced waves.The MTSI is only excited closer to the ramp and only here wave amplitudes grow to significant values.This may prevent phase mixing and subsequent reformation.To summarize, the present case is a high ion beta case, which at lower Mach number would result in a steady shock without the excitation of the MTSI.However, due to the large Mach number the MTSI is excited, although with a small growth rate.
The simulation can explain many features of the Lobzin et al. (2007) observations during the 24 January 2001 Cluster bow shock crossing.It is suggested that the high frequency waves are due to the MTSI.At small wave amplitudes further upstream we predict that the Poynting flux is directed downstream.This is not necessarily so when the waves move into high density, high magnetic field compression region in the ramp and overshoot.Here wave refraction can lead to an upstream directed group velocity.In situ wave observations with upstream directed group velocity can not necessarily been interpreted as being related to an almost phase standing whistler radiated from the shock ramp.The nonstationary "breathing" behavior on ion time scales can lead to differences in the temporal profiles of spacecraft separated by more than ∼ 50λ e or more than about one ion inertial length.This nonstationarity also lead to a bursty behavior of the ion reflection rate on time scales of several inverse ion gyrofrequencies.We thus conclude that the observations by Lobzin et al. (2007) do not support the reformation mechanism of a nonlinear whistler catastrophe proposed by Krasnoselskikh et al. (2002), but can easily be interpreted on the basis of the excitation of the MTSI in a nonstationary shock.
A few caveats are in order.In the present simulations the plasma is still strongly magnetized since a value of ω pe / ce = 8 has been used compared to a value of 230 during the 24 January 2001 event (V.V. Krasnoselskikh, private communication, 2009).The low value of ω pe / ce used here may result in strong overestimation of the role of the electric field for particle dynamics.The present simulations are onedimensional in space.The one-dimensionality may overemphasize the amplitude of the MTSI waves since a rather limited region of k space is available.Furthermore if a higher spatial dimensionality is allowed, other instabilities with different k vectors relative to the magnetic field may get excited (Matsukiyo and Scholer, 2006), (Lembège et al., 2009).As with all PIC simulations there is the problem with the low particle number per cell used in the simulations.This may cause an artificially high noise from which the fluctuations can grow, and the conditions for the amplitudes to grow to a significant amplitude are easer to be satisfied.On the other hand the real solar wind may include physical noise other than a thermal noise which is not included in the simulations and which may act as seed fluctuations for instabilities.

Fig. 1 .
Fig. 1.Color coded magnetic field Bz component in the t − x plane.The magnetic field has been shifted into an average frame standing with the shock.

Fig. 1 .
Fig. 1.Color coded magnetic field B z component in the t −x plane.The magnetic field has been shifted into an average frame standing with the shock.

Fig. 3 .Fig. 3 .
Fig. 3. From top to bottom: magnetic field component Bz, ion density, and ion vix − x phase space plots in the foot and the ramp for the mi/me = 1840 run at tΩci = 9.5.

Fig. 4 .Fig. 4 .
Fig. 4. From top to bottom: magnetic field component Bz, ion density, and ion vix − x phase space plots in the foot and the ramp for the mi/me = 1840 run at tΩci = 10.3.

AFig. 5 .
Fig.5.Growth rate γ / ci versus k in units of ω pe /c for various values of the angle between the ion beams and the magnetic field.In the upper panel a velocity of M A = 9 has been assumed and in the lower panel a value of M A = 5 was assumed.

Fig. 6 .Fig. 6 .
Fig. 6.Stacked profiles of the magnetic field Bz component.Four wave crests are indicated by filled circles at two different times.

Fig. 7 .Fig. 7 .
Fig. 7. Color coded Poynting flux S as a function of wavenumber and x shock in the shock ramp at tΩci = 9.2.The trace of the magnetic field Bz component is shown for reference.The wave train indicated by filled dots is followed in time.

Fig. 9 .
Fig. 9. Color coded magnetic field Bz component in the t − x plane for the ion to electron mass ratio 100 run.The magnetic field has been shifted into an average frame standing with the shock.

Fig. 9 .
Fig. 9. Color coded magnetic field B z component in the t −x plane for the ion to electron mass ratio 100 run.The magnetic field has been shifted into an average frame standing with the shock.

Fig. 10 .Fig. 10 .
Fig. 10.Stacked profiles of the magnetic field Bz component for the ion to electron mass ratio 100 run.

Fig. 11 . 27 Fig. 11 .
Fig. 11.Number of reflected ions within a distance of 50λe upstream of the maximum of the cross shock potential.

Fig. 12 .Fig. 12 .
Fig. 12. From top to bottom: magnetic field Bz profile recorded by an artificial spacectraft (S/C 1) moving in the simulation from upstream to downstream, same profile after high-pass filtering, high pass-filtered magnetic field profile recorded by second spacecraft (S/C 2) moving with the same velocity but displaced by ∼ 50λe relative to S/C 1.

Fig. 13 .Fig. 13 .
Fig. 13.Same as Figure 11 for the ion to electron mass ratio 1840 case.For details see text.