1 Introduction
For perpendicular collisionless shocks above a certain critical Mach number there is no stable two-fluid solution based solely on resistivity. However, quasi-steady solutions are in fact possible, and such shocks observed, for supercritical shocks because some fraction of the incoming ions are specularly (or near-specularly) reflected at the shock. The gyration of the reflected ions as they return to the shock creates a foot (a gradual decrease in average velocity and increase in density and magnetic field) upstream of the ramp (where the flow quantities change most rapidly). At the ramp, and for a short distance downstream, the magnetic field increases above the downstream asymptotic value, forming what is termed the overshoot. The gyrating reflected ions provide the downstream temperature increase required by the shock conservation relations. The foot/ramp/overshoot structure of the supercritical perpendicular (and quasi-perpendicular) shock is well established from both observations (e.g. Sckopke et al. Reference Sckopke, Paschmann, Bame, Gosling and Russell1983; Scudder et al. Reference Scudder, Mangeney, Lacombe, Harvey, Aggson, Anderson, Gosling, Paschmann and Russell1986) and simulations (e.g. Leroy et al. Reference Leroy, Winske, Goodrich, Wu and Papadopoulos1982); see also Burgess & Scholer (Reference Burgess and Scholer2015).
Numerical simulations (Biskamp & Welter Reference Biskamp and Welter1972; Quest Reference Quest1985; Lembège & Dawson Reference Lembège and Dawson1987; Lembège & Savoini Reference Lembège and Savoini1992) revealed a strongly non-stationary behaviour of quasi-perpendicular shocks in the supercritical regime. Similar behaviour is also evidenced with parameters relevant to the heliospheric and astrophysical contexts (Shimada & Hoshino Reference Shimada and Hoshino2000; Schmitz, Chapman & Dendy Reference Schmitz, Chapman and Dendy2002). This non-stationarity is characterized by a periodic self-reformation of the shock front over ion time scales. In the case of strictly perpendicular shocks, the self-reformation process reveals itself to be very sensitive to Alfvén Mach number
$M_{A}$
and ion upstream
${\it\beta}_{i}$
(Hellinger, Trávníček & Matsumoto Reference Hellinger, Trávníček and Matsumoto2002; Hada et al.
Reference Hada, Onishi, Lembège and Savoini2003). Recent reviews (Hellinger Reference Hellinger2003; Lembège et al.
Reference Lembège, Giacalone, Scholer, Hada, Hoshino, Krasnoselskikh, Kucharek, Savoini and Terasawa2004) stress the importance of the self-reformation process for the shock properties. By using a 1-D full particle-in-cell (PIC) simulation code, Scholer, Shinohara & Matsukiyo (Reference Scholer, Shinohara and Matsukiyo2003) pointed out that the non-stationary self-reformation process switches off some acceleration mechanisms, for example the shock surfing mechanism, which has been proposed for time-stationary shock solutions (Sagdeev Reference Sagdeev and Leontovich1966; Lipatov & Zank Reference Lipatov and Zank1999). On the other hand, 1-D PIC simulations (Schmitz et al.
Reference Schmitz, Chapman and Dendy2002; Lee, Chapman & Dendy Reference Lee, Chapman and Dendy2004) show that the shock reformation leads to a strong energization of a portion of reflected ions during their subsequent interaction with the non-stationary shock.
The study of non-stationarity is dominated by results based on simulations and theory, which reflect the difficulty of such observation at heliospheric shocks. The difficulties are associated with interpreting differences between multipoint spacecraft measurements when only a few spatial scales can be observed and the fact that the conversion from time to space in the data relies on knowledge of the shock speed relative to the observation points. The separation between macrostructure and microstructure was already considered in an influential two-spacecraft study using International Sun-Earth Explorer (ISEE) data, where it was found that, by using inter-correlation with different averaging periods, the correlation between the two spacecrafts decreased strongly when the averaging period was approximately 0.15 times the upstream proton gyroperiod (Scudder et al. Reference Scudder, Mangeney, Lacombe, Harvey, Aggson, Anderson, Gosling, Paschmann and Russell1986). The inference is that fluctuations seen at short time scales are associated with non-stationarity, or the breakdown of the usual assumption of one-dimensionality. Using four-point Cluster data, Lobzin et al. (Reference Lobzin, Krasnoselskikh, Bosqued, Pinçon, Schwartz and Dunlop2007) presented evidence of shock reformation based on cross-correlation analysis, claiming that the macrostructure as seen by the different Cluster spacecraft could be explained by reformation. Also, the fraction of reflected ions was seen to vary strongly, and bursts of high-frequency waves were seen. The interpretation was that the reformation was driven by a fluctuating reflected fraction, and the generation of non-stationary whistler wave packets followed the reformation mechanism proposed by Krasnoselskikh et al. (Reference Krasnoselskikh, Lembège, Savoini and Lobzin2002).
Observations of reformation are usually in the context of some assumed model, which means that interpretations can vary. For example, Comişel et al. (Reference Comişel, Scholer, Soucek and Matsukiyo2011) performed 1-D PIC simulations and made comparisons with the same shock crossing studied by Lobzin et al. (Reference Lobzin, Krasnoselskikh, Bosqued, Pinçon, Schwartz and Dunlop2007). Based on the simulation results, they concluded that the high-frequency fluctuations were due to the modified two-stream instability (MTSI), rather than the whistler gradient catastrophe model of Krasnoselskikh et al. (Reference Krasnoselskikh, Lembège, Savoini and Lobzin2002). It could be noted that the use of a 1-D simulation introduces some strong assumptions about the microstructure of the shock, which could affect comparisons with observations.
Two-dimensional hybrid simulations were initially used to study the possible competition between Alfvén ion cyclotron (AIC) and mirror instabilities due to the perpendicular ion temperature anisotropy produced by the reflected-gyrating ions (Winske & Quest Reference Winske and Quest1988). It was found that the ‘ripples’ (surface fluctuations) at the ramp were most likely due to the AIC instability, although no exact agreement with linear homogeneous theory was found. Lowe & Burgess (Reference Lowe and Burgess2003), using similar simulations, found that the shock ramp fluctuations propagated across the shock at the Alfvén speed of the shock overshoot, confirming the importance of AIC-like waves. Both these studies used a 2-D simulation plane which was parallel to the shock coplanarity plane, i.e. the upstream magnetic field was in the simulation plane. When the 2-D simulation plane is arranged perpendicular to the shock coplanarity plane (i.e. the upstream magnetic field is perpendicular to the simulation plane) it is found that a new type of shock fluctuation dominates the shock (Burgess & Scholer Reference Burgess and Scholer2007). These fluctuations are associated with a modulation of the fraction of reflected ions, and propagate across the shock at a speed, and in a direction, determined by the reflected ions.
The question of non-stationarity at perpendicular and quasi-perpendicular shocks has also been addressed by studies using 2-D simulations. Hellinger et al. (Reference Hellinger, Trávníček, Lembège and Savoini2007) found that shock reformation could be suppressed due to the presence of oblique whistlers, which are not present in 1-D simulations. Lembège et al. (Reference Lembège, Savoini, Hellinger and Trávníček2009) reported full-particle PIC simulations with a reasonably high ion–electron mass ratio of 400, and showed that whistler waves in the foot can also suppress reformation. The 2-D hybrid simulations of Yuan et al. (Reference Yuan, Cairns, Trichtchenko, Rankin and Danskin2009), on the other hand, seemed to show that reformation could occur despite the presence of whistler waves. However, it has been shown that reformation can appear to be suppressed if quantities are averaged over a direction tangential to the shock surface, even though the time evolution at a fixed point on the shock surface does show oscillations attributed to reformation (Umeda et al. Reference Umeda, Kidani, Yamao, Matsukiyo and Yamazaki2010).
There are only a limited number of studies using three-dimensional shock simulations (e.g. Thomas Reference Thomas1989; Hellinger, Mangeney & Matthews Reference Hellinger, Mangeney and Matthews1996; Giacalone & Ellison Reference Giacalone and Ellison2000; Guo & Giacalone Reference Guo and Giacalone2013), so it is still an open question to what extent the behaviour and microstructure seen in 2-D simulations persists in three dimensions. There is also the possibility that new types of behaviour emerge in three dimensions. Thus in this paper we perform a series of 2-D and 3-D simulations to describe the microstructure see in hybrid simulations of the perpendicular shock. We restrict the study to the hybrid simulation method so as to concentrate on ion-scale behaviour, and to supercritical but moderate Mach number. The latter condition allows the different mechanisms to be studied without excessive nonlinear cross-coupling. In particular, the shocks we study are not dominated by reformation processes, due to the choice of Mach number, ion beta and the inclusion of a relatively strong resistivity. Although we also show some results to indicate the implications of relaxing these constraints.
The main aim of this paper is to document and explain the physics of microstructure, as seen in hybrid simulations, of moderate Mach number shocks. An important motivation is that the hybrid simulation technique is increasingly being used to model astrophysical shocks and particle acceleration (e.g. Burgess Reference Burgess2006; Caprioli & Spitkovsky Reference Caprioli and Spitkovsky2014), as well as the global interaction between the solar wind and magnetospheres (e.g. Trávníček et al.
Reference Trávníček, Hellinger, Schriver, Herčík, Slavin and Anderson2009), and therefore it is important to understand the physical mechanisms responsible for shock fluctuations and microstructure in such simulations. We are particularly interested in the structure seen in the shock ramp with high-resolution low-noise 3-D simulations, that might not be seen in coarser simulations. This aim is framed in the context of simulation studies. In the studies shown in this paper the smallest spatial resolution is
$0.1{\it\lambda}_{i}$
, where
${\it\lambda}_{i}$
is the ion inertial length. However, it is important to note that at real shocks the microstructure extends across all plasma scales, down to the Debye scale. In particular, there have been recent studies which indicate that the smallest magnetic field scales extend to the electron inertial length
${\it\lambda}_{e}$
(Hobara et al.
Reference Hobara, Balikhin, Krasnoselskikh, Gedalin and Yamagishi2010; Mazelle et al.
Reference Mazelle, Lembège, Morgenthaler, Meziane, Horbury, Génot, Lucek and Dandouras2010) and that electron heating can occur on scales of a few
${\it\lambda}_{e}$
. Such recent observations, as reviewed in Krasnoselskikh et al. (Reference Krasnoselskikh, Balikhin, Walker, Schwartz, Sundkvist, Lobzin, Gedalin, Bale, Mozer and Soucek2013), are used to infer the importance of electron dynamics for the shock structure (e.g. Krasnoselskikh et al.
Reference Krasnoselskikh, Vinogradova, Balikhin, Alleyne, Pardaens, Woolliscroft, Klimov, Petrukovich, Mier-Jedrzejowicz and Southwood1991; Sundkvist et al.
Reference Sundkvist, Krasnoselskikh, Bale, Schwartz, Soucek and Mozer2012). There is a coupling between the structure of the foot introduced by the presence of reflected-gyrating ions, whose scale is of the order of
$M_{A}{\it\lambda}_{i}$
, and the presence and behaviour of microstructure down to scales of order
${\it\lambda}_{e}$
. The full nature of this cross-scale coupling, and how it depends on shock parameters, are not completely known. A second aim of the paper is, by documenting the structure seen in hybrid simulations, to set the context for observational tests which can probe the validity of the hybrid method at ion scales, in the presence of such cross-scale processes.
A few comments on the limitations of simulations are appropriate. We use a hybrid simulation method in two and three (spatial) dimensions, where electrons are considered as a massless, charge neutralizing adiabatic fluid whereas ions are treated as macro-particles. The algorithm used is described in Matthews (Reference Matthews1994). The hybrid method neglects electron Landau and cyclotron resonance effects, although subgrid processes are included via a resistivity term in the field equations. The hybrid method has proved very successful for modelling ion dynamics in space plasmas, however it can be argued that it does not include the effects at the electron inertial scale. The importance of this neglect can be tested against observations. For example, high-resolution comparisons of hybrid simulations with spacecraft observations of the quasi-parallel shock have shown good agreement with the first stage of ion energization from the thermal population (Sundberg et al. Reference Sundberg, Haynes, Burgess and Mazelle2016). In principle, explicit full-particle PIC simulations should be able to model shock microstructure from ion to Debye scales. However, there are computational constraints which limit the size of the simulation and the time over which it can be run. Values are chosen for ratios of ion-to-electron masses and plasma-to-cyclotron frequencies which reduce the computational cost and make such simulations viable. However, it has been argued by Krasnoselskikh et al. (Reference Krasnoselskikh, Balikhin, Walker, Schwartz, Sundkvist, Lobzin, Gedalin, Bale, Mozer and Soucek2013) that the artificial values of the electron plasma-to-electron cyclotron frequency ratio, necessitated by the simulations (e.g. Comişel et al. Reference Comişel, Scholer, Soucek and Matsukiyo2011), produce electric fields in the simulations at least two orders of magnitude larger than those observed in heliospheric shocks. This underlines the importance of interpreting simulation results and the importance of validation against observations.
As the above discussion implies, an obvious question is to what extent does the microstructure seen in hybrid simulations correspond to a coarser version of what is seen in explicit PIC simulations. This question is still open, given the disparity in length and time scales possible with the two methods and the frequent use of non-realistic simulation parameters (although use of a realistic particle mass ratio is more common than in the past). PIC simulations have been used to investigate whistler structure generated by modified two-stream instabilities in the foot driven by drift between both incoming and outgoing ions and the electrons (e.g. Matsukiyo & Scholer Reference Matsukiyo and Scholer2003, Reference Matsukiyo and Scholer2006). The importance of whistler-like structure generated by these instabilities for electron acceleration has been investigated by Riquelme & Spitkovsky (Reference Riquelme and Spitkovsky2011). The Buneman instability is an alternative source of microstructure within the shock transition, although it is likely only to operate at relatively high Mach number shocks (e.g. Amano & Hoshino Reference Amano and Hoshino2007). The electron temperature anisotropy may also play a role in the foot and shock transition, but similarly, this is neglected in the current simulations since the hybrid code has no information on the electron particle distribution function.
In the majority of the simulations shown below, the 2-D simulation box is
$600\times 512$
with a spatial resolution
${\rm\Delta}x=0.1c/{\it\omega}_{pi}$
and
${\rm\Delta}y=0.2c/{\it\omega}_{pi}$
and there are 1024 particles per cell in the upstream region. The 3-D simulation box is
$512\times 128\times 128$
points with a spatial resolution
${\rm\Delta}x=0.1c/{\it\omega}_{pi}$
and
${\rm\Delta}y={\rm\Delta}z=0.2c/{\it\omega}_{pi}$
, and there are 256 particles per cell in the upstream region. The transverse size of the box (in the
$y$
and
$z$
directions) is large enough to allow many wavelengths of the fluctuations to be included, so that there is a good sampling of allowed propagation angles within the shock transition. The figures with spatial variations usually only show a small section of the simulation grid used for the Fourier space analysis.
The time step for the particle advance is
${\rm\Delta}t=0.005/{\it\Omega}_{i}$
whereas the magnetic field is advanced with
${\rm\Delta}t_{B}={\rm\Delta}t/10$
. In these definitions,
$c$
is the speed of light,
${\it\omega}_{pi}$
is the upstream ion plasma frequency,
${\it\lambda}_{i}$
is the ion inertial length and
${\it\Omega}_{i}$
is the upstream ion cyclotron frequency. Protons and electrons have initial (upstream) ratios between the particle and magnetic pressures
${\it\beta}_{i}=0.2$
and
${\it\beta}_{e}=0.01$
, respectively. In the generalized Ohm’s law we use a small resistivity
${\it\eta}=10^{-2}{\it\mu}_{0}v_{A}^{2}/{\it\Omega}_{i}$
; here
${\it\mu}_{0}$
is the magnetic permeability of vacuum and
$v_{A}$
is the upstream Alfvén velocity. For 2-D and 3-D simulations the plasma is streaming along the
$x$
axis (protons are continuously injected at
$x=0$
) and interacts with a piston at
$x_{max}$
which launches a shock travelling towards negative
$x$
.
In this paper we demonstrate the role of the ion Weibel instability as a driver of whistler-like fluctuations in the foot of the shock. The ion Weibel instability has received attention as crucial for the formation of shocks in unmagnetized (or low magnetization) plasmas (Spitkovsky Reference Spitkovsky2008), but here we show that it also plays a role in magnetized shocks as found at the Earth’s bow shock and in other heliospheric environments. Furthermore, we show that the 3-D microstructure is dominated, at ion scales, by a combination of the shock rippling and reflection fraction instability seen in 2-D simulations made in, and perpendicular to, the coplanarity plane. The result is a pattern of fluctuations propagating across the shock surface at an angle to both the magnetic field direction and the direction of motion of the reflected-gyrating ions. We also present results showing the consequences of relaxing the constraints of moderate Mach number and small ion beta.
2 Results: 2-D and 3-D hybrid simulations
2.1 2-D hybrid simulation in the coplanarity plane
We first consider a 2-D simulation of a perpendicular shock with
${\it\theta}_{Bn}=90^{\circ }$
and
$M_{A}=3.3$
, such that the upstream magnetic field is in the simulation plane, meaning that the simulation plane is also the coplanarity plane. Animations of the magnetic field fluctuations in the foot show a complex behaviour of the phase fronts. At the ramp and overshoot field-parallel propagating waves/structures are seen corresponding to the AIC modes previously studied (Winske & Quest Reference Winske and Quest1988; Hellinger & Mangeney Reference Hellinger and Mangeney1997; Lowe & Burgess Reference Lowe and Burgess2003). Also present are shorter-wavelength fluctuations with obliquely oriented phase fronts (neither field aligned nor aligned with the shock normal) which seem to have properties varying between propagating and phase standing. In figure 1 we show an overview of the shock fluctuations with a spectral analysis. The envelope of magnetic fluctuations is shown in figure 1(a,b), demonstrating, as earlier works have shown, that the amplitudes of fluctuations in all components of the field are comparable. This makes impossible any analysis based on simplistic expectations of mode polarization. (For example, the AIC is transverse and so would be expected to be seen only in the
$B_{x}$
and
$B_{z}$
components.)
2-D hybrid simulation in the coplanarity plane: (a) the average (see the text) profile of
$B_{y}$
(solid line) and the two extreme values of
$B_{y}$
(dotted lines); (b) average profiles of
$B_{x}$
and
$B_{z}$
(solid line), the dotted and dashed lines show the two extreme values of
$B_{z}$
and
$B_{x}$
, respectively. (c) The colour scale plot of
$B_{y}$
at
$t=40{\it\Omega}_{i}^{-1}$
as a function of
$x$
and
$y$
; (d) the colour scale plot of the fluctuating magnetic field
${\it\delta}B$
at
$t=40{\it\Omega}_{i}^{-1}$
as a function of
$x$
and
$k_{y}$
. (e) the colour scale plot of the fluctuating magnetic field
${\it\delta}B$
at
$t=40{\it\Omega}_{i}^{-1}$
as a function of
$x$
and
$k_{x}$
(using the wavelet transform). (f) The fluctuating magnetic field
${\it\delta}B$
between 20 and
$40{\it\Omega}_{i}^{-1}$
(within the shock front) as a function of
$k_{y}$
and
${\it\omega}$
in the average shock rest frame; here the mean shock profile has been removed.
A snapshot of the field magnitude fluctuations shows (figure 1
c) the complex pattern of phase fronts at one time (
$t_{m}=40{\it\Omega}_{i}^{-1}$
). Only a small section of the simulation is shown, whereas the
$k$
-space and
${\it\omega}{-}k$
analysis uses the full domain. The average profile in
$x$
has been removed. Spectral (wavevector) analysis of this pattern using a wavelet transform in
$x$
and Fourier transform in
$y$
is shown in figure 1(d,e). As expected, the field magnetic fluctuation power across wavevector spaces
$k_{x}$
and
$k_{y}$
peaks around the shock ramp at
$x\sim 8$
. The dependence of power with
$k_{x}$
, especially the spread of power at
$k_{x}\sim 5$
, is dominated by the wavelet response to the finite size of the foot. On the other hand, the dependence on
$k_{y}$
show peaks at small
$k_{y}$
(apparently close to zero) and at
$k_{y}\sim 2$
.
The latter structure is analysed in the
${\it\omega}{-}k_{y}$
power spectrum shown in figure 1(f). This is the result of a 2-D space–time transform of the variation of the field magnitude within the shock front as a function of
$y$
between
$t=20$
and
$40{\it\Omega}_{i}^{-1}$
, and in the shock rest frame with the mean shock profile removed. Recall that the power spectrum has 2-fold symmetry with
$P(-{\it\omega},-k_{y})=P({\it\omega},k_{y})$
. The spectrum consists of three components: the first is centred at
$k_{y}=0$
, with enhancements across a range of frequencies, and corresponds to fluctuations in the mean shock speed and mean shock profile. Such variations are a common feature of hybrid simulations, as noted in early 1-D studies (Leroy et al.
Reference Leroy, Winske, Goodrich, Wu and Papadopoulos1982). The second component has the signature of propagating waves with phase speed
${\it\omega}/k_{y}\approx \pm 2V_{A}$
. This component corresponds to the field-parallel fluctuations propagating at approximately the Alfvén speed of the shock ramp/overshoot as reported earlier (Lowe & Burgess Reference Lowe and Burgess2003). The third component in the power spectrum is an enhancement centred on
${\it\omega}=0$
and
$k_{y}\approx \pm 2$
, but also with the signature of propagating fluctuations with a slope
${\it\omega}/k_{y}\approx \pm 4V_{A}$
. Visually, the component can be described as a cross in
${\it\omega}{-}k_{y}$
shifted to non-zero
$k_{y}$
. This third component can be explained as the sidebands produced by the modulation of a propagating signal by a zero-frequency
$k_{y}\sim 2$
signal. Although the propagating signal might be tentatively identified with oblique whistler mode fluctuations (seen in figure 1
c), a source for a zero-frequency
$k_{y}\approx 2$
signal is not obvious.
We next identify the zero-frequency
$k_{y}\approx 2$
component in the field fluctuations and then explore the hypothesis that it is due to the ion Weibel instability. As discussed earlier, there are several coupled systems at work within the foot which can make it difficult to disentangle cause and effect. For example, an important factor is that at any point in the ramp the reflected fraction depends on the fields at that point. But the reflected fraction modulates the foot (upstream of the ramp) with a delay determined by the particle trajectories. Thus, considering spatial variation across the shock surface (i.e. in the
$y$
direction), fluctuations at the ramp, e.g. due to AIC-like fluctuations, can spatially modulate the density of reflected ions at the upstream edge of the foot. In order to eliminate this particular source of modulation of the foot, we carry out simulations which create, for a short time, the idealized case of uniform reflected-ion density across the shock front. This allows us to identify fluctuations driven by instabilities of the reflected-ion beam, without the complication of spatial variation of beam density. The method is to carry out a 2-D simulation, but to split the simulation domain into two parts. In the first part (nearest the wall which launches the shock) the shock is constrained to behave as if it were one-dimensional by, for each field component and at each
$x$
position, averaging over the
$y$
direction and then forcing the fields to take that average value for all
$y$
. Any variation with
$y$
(transverse to the shock normal) is suppressed and the shock structure depends only on
$x$
. However, in the second (left-hand) part of the simulation domain, the full 2-D system is allowed to evolve as normal. The shock is launched at the right-hand wall and travels to the left. Thus, the shock is initially one-dimensional with no variation across the
$y$
direction. However, the shock will gradually move into the left-hand part where 2-D structure can develop. It will be firstly the upstream edge of the foot, then the foot itself, followed by the ramp which progressively is allowed to develop 2-D structure. Of course, there is some loss of self-consistency at the boundary between 1-D and 2-D subdomains, but by following the evolution until the shock is fully two-dimensional it is apparent that any artefacts which have been introduced are limited to only a few cells.
2-D hybrid simulation in the coplanarity plane with a forced 1-D to 2-D transition: (a)
$B_{y}$
profile averaged over
$y$
direction, with 1-D and 2-D regions indicated; (b–d) colour scale plots of fluctuations in
$B_{x}$
,
$B_{y}$
and
$B_{z}$
(a–e) after removal of
$y$
-averaged profile; (e) power in component fluctuations as a function of
$x$
. Only a small region of the simulation domain, around the shock transition, is shown.
This technique of transitioning the shock foot from one- to two-dimensional is illustrated in figure 2 which shows the
$y$
-averaged field magnitude profile, the magnetic field component fluctuations as a function of
$x$
and
$y$
and the
$y$
-averaged profile of the power in the components. The
$y$
range in this simulation is twice that shown, and only a small region around the shock in the
$x$
direction is shown. The regions which are treated as 1-D and 2-D are marked. It can be seen that the shock has propagated far enough that the front edge of the foot is allowed to behave two-dimensionally, whereas the ramp at
$x=18.3$
is still constrained to be one-dimensional. The density of reflected ions at this time has been confirmed to be uniform in
$y$
across the foot. Fluctuations in all field components grow extremely rapidly almost immediately as the foot moves into the 2-D subdomain. The
$B_{x}$
component shows a wave pattern with wavevector in the
$y$
direction and a wavelength of approximately
$3{\it\lambda}_{i}$
(i.e.
$k{\it\lambda}_{i}\approx 2$
). The
$B_{y}$
and
$B_{z}$
components show slightly different behaviour, with a more complicated pattern which seems to develop further downstream of the front of the foot. The power in the
$B_{y}$
and
$B_{z}$
components is higher than in
$B_{x}$
, but all components develop in a similar way across the front edge of the foot.
Evolution of
$B_{x}$
as a function of
$y$
and time, for the shock experiment shown in figure 2, at a fixed distance upstream of the shock ramp. At early times the shock at this position transitions from 1-D to 2-D, and at later times the entire foot and ramp have fully developed 2-D structures.
The
$B_{x}$
component of the fluctuations at the front (upstream) edge of the foot seems to have the appropriate value of
$k_{y}\approx 2$
for the modulation seen in the
${\it\omega}{-}k_{y}$
power distribution in figure 1. But that modulation also requires
${\it\omega}\approx 0$
. To investigate the time behaviour, in figure 3 we plot the time evolution of the
$B_{x}$
fluctuations as a function of time and
$y$
, at a fixed distance upstream of the ramp in a frame in which the shock is at rest. The time range covers that corresponding from when the foot first enters the 2-D simulation subdomain until when all the foot and ramp are fully two dimensional. At the earliest time (after
${\it\Omega}_{i}t=16$
) a very clean monochromatic wave is seen with a wavelength of approximately
$3{\it\lambda}_{i}$
. There appears to be little propagation, so the frequency is zero, or close to zero. At slightly later times, after
${\it\Omega}_{i}t=19$
, the time evolution becomes more complex and more variable across
$y$
. It appears that the pattern is still dominated by fluctuations with
${\it\omega}\sim 0$
. The density of reflected ions at this position is uniform across
$y$
at early times, but strongly variable once the pattern in
$B_{x}$
becomes more variable. It can be noted that the time scale for a reflected ion to return to the shock is approximately
$1.9{\it\Omega}_{i}^{-1}$
and the flow transit time is approximately
${\it\Omega}_{i}^{-1}$
, so the time period taken for the reflected ion density to become non-uniform over
$y$
is consistent with modulation of the reflected fraction by the fluctuations seen at earlier times, once they convect into the ramp.
Evolution of
$B_{x}$
,
$B_{y}$
and
$B_{z}$
(a–c) as a function of
$y$
and time for the 2-D hybrid simulation shown in figure 2, at a fixed distance upstream of the shock ramp (as shown in g). Frequency – wavevector power distributions for
$B_{x}$
,
$B_{y}$
and
$B_{z}$
(d–f) for the fluctuations shown in (a–c).
We can further distinguish the behaviour of the fluctuation components in figure 4 which shows the space–time and also
${\it\omega}{-}k_{y}$
plots taken at a position a fixed distance upstream of the shock ramp in the shock rest frame, as shown in (g). The sampling position is very near the upstream edge of the shock, and the time period covers the emergence of the shock front part of the foot into the 2-D subdomain. The
$B_{x}$
fluctuations show a clear non-propagating signal, confirmed by the
${\it\omega}{-}k_{y}$
behaviour, with a strong signal at
$k{\it\lambda}_{i}\approx 2$
. The
$B_{y}$
and
$B_{z}$
components are, again, more complex, but the
${\it\omega}{-}k_{y}$
plot shows a mixture of zero frequency and propagating characteristics shifted in
$k_{y}$
.
Our interpretation of these results is as follows: at the very upstream edge of the foot, the reflected ions, near their turnaround point, can drive an instability with a zero frequency, relatively high
$k_{y}$
wave pattern and with a wavevector in the field-parallel
$y$
direction. The signal produced is most clear in
$B_{x}$
, but is coupled to
$B_{y}$
and
$B_{z}$
which produces propagating (non-zero frequency) waves. The propagating waves are essentially the fluctuations which would be identified as ‘whistler fluctuations’ since they are driven at high
$k_{y}$
and have corresponding relatively high phase speeds. Examination of their polarization shows that they are right handed, as expected for whistler-like fluctuations. The coupling between the field components is not straightforward, and it is complicated by the strong inhomogeneity in field and velocity as the plasma decelerates in the foot and ramp. In a shock with a fully 2-D developed structure the fraction of reflected ions is non-uniform over
$y$
due to fluctuations at the ramp, which in turn may be due to AIC fluctuations generated there, or fluctuations coupled to the foot. Despite the non-uniform density of the reflected-ion beam at the upstream edge of the shock, the zero-frequency high
$k_{y}$
instability can still operate. The wave pattern from this instability drives (propagating) whistler fluctuations, and it is this modulation by the zero-frequency component that is responsible for the power contribution shifted to
$k{\it\lambda}_{i}\approx \pm 2$
, as seen in the field magnitude at the ramp (figure 1
f).
An analogy is the pattern of surface fluctuations produced by a linear grid of fine wires perpendicular to a steady fluid flow, where the wires are fine enough to generate dispersive ripples. Very close to the grid (relative to the inter-wire spacing) the
${\it\omega}{-}k$
spectrum transverse to the flow will mainly show the spatial frequency of the grid spacing. Further downstream, the spectrum will show propagating modes associated with the ripples, but modulated by the spatial frequency of the grid. Far downstream, interaction between waves from many grid wires may produce a more complex, even turbulent, state.
In the 1-D/2-D simulations, the rapid appearance of the fluctuations at the upstream edge as it enters the 2-D subdomain indicates that they are driven by an instability of the reflected-ion beam. Consider specularly reflected ions at the perpendicular shock. At the point, where the reflected ions turn around, their
$x$
velocity is zero, and they form a beam in velocity space perpendicular to the magnetic field with
$v_{z}\approx \sqrt{3}M_{A}V_{A}$
. For this situation it is known that an instability due to a cross-field ion drift (Chang, Wong & Wu Reference Chang, Wong and Wu1990; Wu et al.
Reference Wu, Yoon, Ziebell, Chang and Wong1992) can drive an electromagnetic, zero-frequency wave with a high growth rate. The mechanism for this instability is essentially that of the Weibel instability (Weibel Reference Weibel1959), as described by Fried (Reference Fried1959). The instability mechanism is usually considered for electrons, but it is also viable for ions, and so in this context we follow Chang et al. (Reference Chang, Wong and Wu1990) who used the term ‘ion Weibel instability’.
For the homogeneous, cold ion beam case, with assumptions similar to the hybrid simulation (see appendix A), there is instability at zero frequency for all
$k$
, with a maximum growth rate for
$\boldsymbol{k}$
aligned field parallel. The growth rate
${\it\gamma}$
decreases at small
$k$
, and has a typical value of approximately
${\it\gamma}\sim 2.2{\it\Omega}_{i}$
, for beam density
$n_{b}=0.2n_{i}$
and cross-field velocity
$v_{b}=5V_{A}$
. An estimate of a typical value of
$k$
for which the growth rate is close to its maximum is approximately
$k_{m}\sim 2.2$
for the simulation shown here, which is close to the value seen in the results. The simplifications used in these estimates are discussed in the appendix A.
2.2 2-D hybrid simulation perpendicular to the coplanarity plane
We have presented evidence that, in 2-D hybrid simulations, some of the microstructure which is ‘whistler-like’ is due to driving from non-propagating waves associated with the ion Weibel instability at the upstream edge of the ion foot. It is an important question whether this scenario is valid in the case where the shock is allowed to have 3-D spatial structure. In this section, as an intermediate case, we consider a 2-D hybrid simulation, similar to that shown in figure 1, but with the simulation plane perpendicular to the shock coplanarity plane. This corresponds to the upstream magnetic field being arranged perpendicular to the simulation plane. This configuration was studied by Burgess & Scholer (Reference Burgess and Scholer2007) who found a nonlinear instability associated with the reflected-ion fraction. The pattern formed in this instability is carried by the specularly reflected ions, and thus travels with their speed in the cross-field direction.
2-D hybrid simulation with simulation plane perpendicular to the shock coplanarity plane: (a) the average (see the text) profile of
$B_{z}$
(solid line) and the two extreme values of
$B_{z}$
(dotted lines); (b) average profiles of
$B_{x}$
and
$B_{y}$
(solid line), the dotted and dashed lines show the two extreme values of
$B_{y}$
and
$B_{x}$
, respectively. (c) The colour scale plot of
$B_{z}$
at
$t=40{\it\Omega}_{i}^{-1}$
as a function of
$x$
and
$y$
; (d) the colour scale plot of the fluctuating magnetic field
${\it\delta}B$
at
$t=40{\it\Omega}_{i}^{-1}$
as a function of
$x$
and
$k_{y}$
. (e) The colour scale plot of the fluctuating magnetic field
${\it\delta}B$
at
$t=40{\it\Omega}_{i}^{-1}$
as a function of
$x$
and
$k_{x}$
(using the wavelet transform). (f) The fluctuating magnetic field
${\it\delta}B$
between 20 and
$40{\it\Omega}_{i}^{-1}$
(within the shock front) as a function of
$k_{y}$
and
${\it\omega}$
in the average shock rest frame; here the mean shock profile has been removed.
In figure 5 we show the results of a simulation with parameters similar to the shock shown in figure 1, where the upstream magnetic field is in the
$z$
direction and the simulation plane is the
$x{-}y$
plane. Comparing with figure 1, it is evident that the behaviour of fluctuations in the foot/ramp is very different compared to when fluctuations with a field-parallel orientation are allowed. The overall level of fluctuations is much lower (a,b), the
$k_{y}$
and
$k_{x}$
spectra of fluctuations are very different from the earlier case, with negligible power at high
$k$
. In the
$k_{y}$
spectrum there is only a single, strong, component at
$k_{y}{\it\lambda}_{i}\approx 0.75$
. From the snapshot of field magnitude (c) this is seen to take the same form as the nonlinear instability found by Burgess & Scholer (Reference Burgess and Scholer2007), i.e. a pattern of ‘tongues’ of field magnitude (and also reflected-ion density) in the foot. The
${\it\omega}{-}k_{y}$
power spectrum (f) shows a component spread over
${\it\omega}$
but at zero
$k_{y}$
. This probably corresponds to small variations in the shock speed, independent of
$y$
, which are not taken into account when translating into the shock frame assuming a constant shock speed. The only other component corresponds to peak power at
$k_{y}\approx 0.75$
and
${\it\omega}\approx 4.1$
, or a phase speed of 5.5
$V_{A}$
. The cross-field velocity of the reflected ions at their turnaround point is
$1.7M_{A}V_{A}\approx 5.6$
, so the phase velocity of the pattern seen in the field magnitude matches the cross-field velocity of the reflected ions at their turnaround point, in agreement with Burgess & Scholer (Reference Burgess and Scholer2007).
In the context of the current study, excluding fluctuations with a field-parallel orientation (by the choice of simulation plane) suppresses all of the mechanisms for microstructure discussed in § 2.1. Instead, a strong and monochromatic variation of the foot/ramp is found to be associated with the cross-field motion of the reflected ions.
2.3 3-D hybrid simulations
In order to determine which of the types of microstructure described above survive or dominate when full 3-D spatial dependence is allowed, we have performed a 3-D hybrid simulation of a shock with a similar Mach number (
$M_{A}=3.1$
) and otherwise identical parameters as for the earlier simulations. The upstream magnetic field is aligned with the
$y$
axis,
$\boldsymbol{B}=(0,B_{0},0)$
, so that the coplanarity plane is the
$x$
–
$y$
plane. Specularly reflected ions, from their sense of gyration, will have a
$-v_{z}$
velocity at their turnaround point.
The shock microstructure is shown in figure 6. In order to facilitate comparisons with the 2-D simulations, we refer to the 2-D simulations as ‘simulation plane parallel to shock coplanarity plane’ (
$\text{SP}\Vert \text{CP}$
, figure 1), and ‘simulation plane perpendicular to shock coplanarity plane (
$\text{SP}\bot \text{CP}$
, figure 5). Although we will not discuss the ion distributions, we note that Burgess & Scholer (Reference Burgess and Scholer2007) showed that the ion temperature anisotropy at the shock ramp is much reduced for the
$\text{SP}\Vert \text{CP}$
case compared to the
$\text{SP}\bot \text{CP}$
case. We find that for 3-D simulations the anisotropy is similar to the
$\text{SP}\Vert \text{CP}$
case, indicating the importance of field-parallel propagating fluctuations for isotropization.
3-D hybrid simulation of a perpendicular shock with
$M_{A}=3.1$
and
${\it\beta}_{i}=0.2$
: (a–c) show (solid) the average profiles and (dotted) the extreme values of (a)
$B_{y}$
, (b)
$B_{z}$
and (c)
$B_{x}$
as functions of
$x$
. (d–f) show colour scale 2-D cuts of
$B_{y}$
$t=30{\it\Omega}_{i}^{-1}$
(d) in the
$x$
–
$y$
plane, (e) in the
$x$
–
$z$
plane and (f) in the
$y$
–
$z$
plane. (g–i) show
${\it\delta}B$
at
$t=30{\it\Omega}_{i}^{-1}$
as a function of (g)
$x$
and
$k_{x}$
, (h)
$x$
and
$k_{y}$
and (i)
$x$
and
$k_{z}$
. (j–l) show the spectra of
${\it\delta}B$
between 10 and
$30{\it\Omega}_{i}^{-1}$
(within the shock front) calculated in the average shock rest frame with mean shock profile removed: (j)
${\it\delta}B$
as a function of
$k_{y}$
and
${\it\omega}$
, (k)
${\it\delta}B$
as a function of
$k_{z}$
and
${\it\omega}$
and (l)
${\it\delta}B$
as a function of
$k_{y}$
and
$k_{z}$
.
Firstly, the envelopes of the field components fluctuations (figure 6
a–c) are most similar to the
$\text{SP}\Vert \text{CP}$
case, but the
$B_{z}$
component (i.e. perpendicular to the coplanarity plane) has a much reduced level, with
${\rm\Delta}B_{z}/B_{0}\sim 1$
compared to
${\rm\Delta}B_{z}/B_{0}\sim 2$
. (d–f) Show colour scale 2-D cuts of
$B_{y}$
at
$t=30{\it\Omega}_{i}^{-1}$
in the
$x$
–
$y$
plane (d), in
$x$
–
$z$
plane (e) and in
$y$
–
$z$
plane (f). The
$x$
–
$y$
and
$x$
–
$z$
planes can be compared to figures 1(c) and 5(c), respectively, since these are parallel and perpendicular to the coplanarity plane. There are similarities in both cases. For the
$x$
–
$y$
plane (coplanarity plane) there is a pattern of oblique phase fronts, which seems to originate from the leading edge of the foot, where the wavelength is
${\sim}3{\it\lambda}_{i}$
. In the case of the
$x$
–
$z$
plane (perpendicular to coplanarity plane) there is (as for the 2-D
$\text{SP}\bot \text{CP}$
simulation) a pattern of ‘tongues’ of field increase, but the dominant wavelength is shorter than in the 2-D simulation case. The pattern in the
$y$
–
$z$
plane (essentially across the face of the shock) is remarkably regular, and the maxima form a centred square lattice with some dislocations.
The
$k$
-space representation for the magnetic field magnitude fluctuations are shown in (g–i). The
$k_{x}$
profile is similar to the
$\text{SP}\Vert \text{CP}$
simulation, but with reduced overall magnitude. The
$k_{y}$
profile shows enhancements in the foot between
$k_{y}\sim 1.5{-}2.7$
and a weaker component at higher values only in the foot. The
$k_{z}$
profile against
$x$
shows a somewhat similar pattern, but at lower amplitude.
The
${\it\omega}{-}k$
-space representations are shown in (e,f), taken at
$x=18.3$
near the ramp. The
${\it\omega}{-}k_{z}$
plot shows that the pattern is dominated by a
$z$
phase speed of
${\approx}-2.5V_{A}$
. This phase speed is in the same sense as the gyration velocity of the reflected ions at turnaround, but somewhat less than what would be expected from the velocity at the ion turnaround position due to specular reflection. The
${\it\omega}{-}k_{y}$
plot shows components again with a range of
$k_{y}$
but the phase speed is now
${\approx}\pm 2.5V_{A}$
. There is some evidence of a component with
$k_{y}{\it\lambda}_{i}\approx 1.75$
across a range of
${\it\omega}$
. This is consistent as a signature of the high-
$k$
(zero-frequency) modulation pattern seen in the
$\text{SP}\Vert \text{CP}$
simulation, although there is no indication of whistler-like fluctuations modulated by this signal. The
$k_{y}$
–
$k_{z}$
plot (l) shows a cruciform pattern as expected from the
${\it\omega}{-}k$
-space together with the
$k_{y}\sim 1.75$
(constant) component.
Combining the information from these plots, we infer that the real space pattern in
$y$
–
$z$
(f) comprises, across a range of
$k$
, two patterns propagating at right angles to each other with
$\boldsymbol{k}\approx 2.5(\pm \widehat{\boldsymbol{y}}-\widehat{\boldsymbol{z}})$
. Since the
$k_{y}$
(field-parallel) power indicates propagating modes in both directions, the obvious inference is that it is associated with the AIC-like ripples seen in the
$\text{SP}\Vert \text{CP}$
simulation, and the phase speed is similar, i.e. the Alfvén speed of the shock ramp or overshoot. Similarly, the
$k_{z}$
power shows a phase speed in a single direction, in the sense of gyration of reflected ions, and so it is natural to associate this with the pattern produced by the reflected ions in the
$\text{SP}\bot \text{CP}$
simulations. However, the phase speed is not exactly the same as in that case, being approximately
$-2.5V_{A}$
rather than the approximately
$-5.3V_{A}$
expected from specular reflection. Inspection of the ion velocity space distributions in the foot of the shock indicates that the lower
$v_{z}$
value of
$-2.5V_{A}$
corresponds to the gyrating ions on the part of their trajectory soon after reflection, rather than at the turnaround point. This is consistent with the pattern of fluctuations (i.e. in the
$y$
–
$z$
plane) being a combination of field-parallel propagation at the ramp/overshoot and motion of the reflected ions at the same location. The inference is that the currents responsible for what has been called the AIC-like fluctuations are carried by the reflected ions, and so the pattern of fluctuations propagates in the sense of the reflected ions.
3-D hybrid simulation for a quasi-perpendicular shock with
$M_{A}=3.0$
,
${\it\theta}_{Bn}=80^{\circ }$
and
${\it\beta}_{i}=0.2$
. Results are shown in the same format as figure 6.
We now show results from simulations which relax the constraints of perpendicular geometry, moderate Mach number and small ion beta. Firstly, we show that the type of structure described above continues to dominate for quasi-perpendicular geometry. Figure 7 shows the results from a simulation for a shock with
$M_{A}=3.0$
,
${\it\theta}_{Bn}=80^{\circ }$
and
${\it\beta}_{i}=0.2$
, where the format is the same as for figure 6. Comparing with figure 6, we see that the magnitude of fluctuations is similar (a–c), the spatial structure (d–f) is similar although in the
$x{-}y$
plane the phase planes have a preferred sense rather than being symmetric. Importantly, the pattern in the
$y{-}z$
plane is very similar, although slightly weaker. In Fourier space (g–l), the pattern of wave power is again similar. The
${\it\omega}{-}k_{z}$
plot shows the same strong asymmetry corresponding to propagation in the
$-z$
direction, in the direction of reflected-gyrating ions. The
${\it\omega}{-}k_{y}$
plot is dominated by field-aligned propagation, but due the quasi-perpendicular geometry there is an asymmetry between the positive and negative
$y$
directions. As before, there is a contribution with
${\it\omega}=0$
and
$k_{y}\approx 2.0$
, as well as propagating components with frequencies to at least
$6{\it\Omega}_{i}$
. In the
$k_{y}{-}k_{z}$
plane these components form a similar cruciform pattern as in figure 6.
3-D hybrid simulation for a perpendicular shock with
$M_{A}=5.5$
and
${\it\beta}_{i}=0.2$
. Results are shown in the same format as figure 6.
We next turn to the effect of increasing the Mach number: figure 8 shows the results from a simulation for a perpendicular shock with
$M_{A}=5.5$
and
${\it\beta}_{i}=0.2$
, again in the same format as figure 6. The fluctuation amplitudes (a–c) are now large, indeed, in the foot, much larger than the average, and at the ramp of the same order as the average. At this level of amplitude one expects strong nonlinear coupling between the various fluctuation components seen at more moderate Mach number, and this is seen in the spatial fluctuation plots (d–f). Both in the
$x{-}y$
and
$x{-}z$
planes, the patterns are irregular, while in the
$y{-}z$
plane it appears to be dominated by ripples with
$y$
aligned ‘ripples’ together with a more complex structure in the
$z$
direction. In the Fourier domain (j–l) the possibility of nonlinear coupling is supported by the broad spread of spectral power, which is very different from the previous results. However, some features remain visible. The signatures of non-stationarity at
$k_{y}=k_{z}=0$
and
${\it\omega}<{\it\Omega}_{i}$
are now prominent, and this seems to link to propagating waves in the
${\it\omega}{-}k_{y}$
plot with phase velocity
${\sim}2.2v_{A}$
. There are, however, fluctuations over a broad range of
$k$
and
${\it\omega}$
, including around
${\it\omega}=0$
and
$k_{y}\approx 2.0$
where peaks in the power were seen previously. Although the
$y{-}z$
plane spatial pattern is dominated by ripples propagating in the
$\pm y$
directions, the
${\it\omega}{-}k_{z}$
plot still shows the strong asymmetry indicating that the
$z$
structured fluctuations are carried in the direction of the reflected-gyrating ions, and with a higher speed than previously, which is expected from the high Mach number. The power in the
$k_{y}{-}k_{z}$
plane hardly shows separate components, and is dominated by power with a maximum around
$k_{z}=0$
. Evidence of coupling between the
$y$
-directed ripples and
$-z$
directed structuring can just be seen in the splitting of the power around
$k_{y}\approx 1.5$
.
We have also performed a simulation of a perpendicular shock with
$M_{A}=4.3$
and
${\it\beta}_{i}=0.2$
. The results are intermediate in form between those for
$M_{A}=3.1$
and
$M_{A}=5.5$
in the sense that there is a weak contribution across a range of
$k$
and
${\it\omega}$
(due to nonlinear coupling) and the separated components seen at lower Mach number are less clear. For example, the component seen with a small range of
$k_{y}$
between 1.5 and 2.0 in the
${\it\omega}{-}k_{y}$
plot is present, but spread over a wider range of
$k_{y}$
.
Finally, we turn to the effect of increasing ion beta. For moderate Mach number
$M_{A}=3.6$
and
${\it\beta}_{i}=1.0$
(not shown) it is remarkable that the fluctuation amplitudes are reduced to a level where there are only minimal fluctuations and the shock is approximately one-dimensional. On the scale of our previous plots, the only fluctuation corresponds to the weak non-stationarity of the shock ramp which is seen in quasi-steady 1-D hybrid simulations (Leroy et al.
Reference Leroy, Winske, Goodrich, Wu and Papadopoulos1982).
3-D hybrid simulation for a perpendicular shock with
$M_{A}=4.8$
and
${\it\beta}_{i}=1.0$
. Results are shown in the same format as figure 6.
At higher Mach number
$M_{A}=4.8$
, the effect of increasing ion beta to
${\it\beta}_{i}=1.0$
can be seen in figure 9. Relatively, the amplitudes of the fluctuations are reduced, and the patterns in both spatial data and Fourier space resume some, but not all, of the regularity seen at lower Mach number and lower ion beta. In the
$y{-}z$
plane there is some sign of the regular checker-board pattern seen in figure 6, but the dominant wave is fairly long-wavelength
$\pm y$
propagating ‘ripples’. This is seen in the
$k_{y}{-}k_{z}$
Fourier space. The importance of the long-wavelength ripples are also seen in the inner cruciform pattern in the
${\it\omega}{-}k_{y}$
plot. Although the spatial pattern in the
$z{-}x$
plane appears smooth, the Fourier space
${\it\omega}{-}k_{z}$
plot shows that, together with a strong signature of non-stationarity around
$k_{z}=0$
, there is also the asymmetric signature associated with the direction of motion of the reflected-gyrating ions.
3 Conclusions
We have carried out a series of high-resolution hybrid simulations to better understand the sources of microstructure and fluctuations in the foot and ramp of perpendicular shocks. We are here concerned with supercritical shocks where specularly reflected ions dominate the foot, ramp and overshoot structure, and provide the bulk of the ion heating. We have chosen parameters such that the shock Mach number is relatively small but still supercritical, so that nonlinear processes such as reformation are not important. Similarly the chosen parameters (perpendicular shock geometry and value of resistivity) do not favour the presence of dispersive whistlers generated at the shock ramp and then propagating upstream into the foot. The 3-D simulations we have carried out are at higher resolution in space and time than previous works, and also use a larger simulation domain. This has allowed us to use spectral analysis to characterize the microstructure fluctuations.
Using 3-D and 2-D simulations in different configurations has enabled us to separate the contributions from different mechanisms, and highlighted the importance of 3-D simulations to correctly model all types of fluctuations in the shock transition. The hybrid method simplifies the electron response, so that short length scale, high-frequency microstructure is not accurately modelled. On the other hand, full-particle simulations are currently not possible with realistic plasma parameters over the length and time scales achievable in hybrid simulations. This leads to hybrid simulations being more widely used to study other phenomena, such as particle acceleration or global systems, and therefore it is important to have a firm understanding of the mechanisms creating microstructure and fluctuations in hybrid shock simulations.
2-D simulations of a perpendicular shock with the simulation plane parallel to the shock coplanarity plane show field-parallel structure consisting of: (1) ripples at the shock ramp which propagate across the shock surface parallel to the magnetic field at approximately the Alfvén speed of the overshoot (Winske & Quest Reference Winske and Quest1988; Lowe & Burgess Reference Lowe and Burgess2003); and (2) whistler-like fluctuations with short wavelengths (Hellinger et al.
Reference Hellinger, Trávníček, Lembège and Savoini2007; Lembège et al.
Reference Lembège, Savoini, Hellinger and Trávníček2009). In addition, we have shown that the whistler-like fluctuations are modulated by a non-propagating mode at the upstream edge of the foot that we have associated with the ion Weibel instability. Using numerical experiments in which the ramp fluctuations are suppressed, we have shown that this instability is associated with the beam of specularly reflected ions at their turnaround point at the upstream edge of the ion foot. We note that this modulation can be seen, albeit at low resolution, in the spectra for
$B_{x}$
shown in Lowe & Burgess (Reference Lowe and Burgess2003) taken in the foot of a shock with
${\it\theta}_{Bn}=88^{\circ }$
and
$M_{A}=5.7$
.
For 2-D simulations where the simulation plane is taken perpendicular to the coplanarity plane, field-parallel structure is suppressed, and instead the microstructure is dominated by fluctuations associated with a spatio-temporal pattern generated by a modulation of the reflected ion fraction (Burgess & Scholer Reference Burgess and Scholer2007). These fluctuations are not linked to any homogeneous plasma wave mode, and are characterized most strongly by the fact that they propagate in the direction and at the speed of the reflected ions at their upstream turnaround position. The relationship between these fluctuations and oscillations of the maximum of the shock overshoot and potential in 1-D (Leroy et al. Reference Leroy, Winske, Goodrich, Wu and Papadopoulos1982), and studies of shock reformation in 2-D (Yuan et al. Reference Yuan, Cairns, Trichtchenko, Rankin and Danskin2009), have not yet been completely clarified.
In 3-D simulations there are signatures of several types of fluctuations seen in the 2-D simulations, indicating a complex, coupled system. There is a remnant of the zero-frequency modulation associated with the ion Weibel instability. However, the fluctuations are now dominated by a pattern of fluctuations generated by two counter-propagating waves inclined to the field direction. An analysis of the
${\it\omega}{-}k$
-space representation shows that the fluctuations have propagation velocity components field-parallel (and anti-parallel) and, additionally, in the sense of motion of the reflected-gyrating ions. The inference is that fluctuations are carried by the reflected ions and hence have a propagation velocity consistent with their cross-field motion. Thus the dominant microstructure seen in the 3-D simulations is the result of a coupling between the field-parallel fluctuations of the ‘coplanarity plane’ 2-D simulations and the reflected-ion fluctuations of the ‘perpendicular to coplanarity plane’ 2-D simulations. The identification of such fluctuations in spacecraft data will be challenging, but may be possible with closely spaced multi-spacecraft constellations, such as Cluster or Magnetospheric Multiscale (MMS).
We have chosen shock parameters suited to study microstructure at the supercritical shock, as modelled by the hybrid plasma simulation method, for a relatively simple situation with moderate Alfvén Mach number and relatively small ion beta. We have also presented results which indicate (1) the mechanisms discussed continue to operate when the shock geometry is quasi-perpendicular; (2) increasing Mach number leads to larger-amplitude fluctuations with strong nonlinear coupling leading to a more ‘turbulent’ signature with fluctuation energy spread in
${\it\omega}{-}k$
space; (3) increasing ion beta leads to a reduction in the fluctuation amplitudes so that at moderate Mach number the spatial structure becomes quasi-unidimensional, whereas as higher Mach number nonlinear coupling is reduced and some of the regularity in
${\it\omega}{-}k$
space seen at lower Mach number returns. These trends in parameter space deserve further investigation, just as the dependence of shock reformation on Mach number and ion beta were found to be crucial in 1-D hybrid simulations (Hellinger et al.
Reference Hellinger, Trávníček and Matsumoto2002). The results presented here will be important as the hybrid method is increasingly used to study, for example, particle acceleration and global shock interactions. As discussed in the Introduction, there are caveats associated with the use of the hybrid simulation to study microstructure. It will be important to test the results presented here against observations, but also against full-particle PIC simulations, in order to better understand the limitations of shock simulations.
Acknowledgements
This work was supported by the UK Science and Technology Facilities Council (STFC) grant ST/J001546/1. The research leading to the presented results has received funding from the European Commission’s Seventh Framework Programme FP7 under the grant agreement SHOCK (project no. 284515). This work used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National EInfrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National E-Infrastructure. This work was also supported by the projects RVO:67985815 and RVO:68378289 funded by the Czech Republic.
Appendix A. Ion Weibel instability
There is an extensive literature for instabilities driven by cross-field drifts. Chang et al. (Reference Chang, Wong and Wu1990) considered the situation of a cold ion streaming perpendicular to a magnetic field, and, by including electromagnetic ion effects, found a purely growing (zero-frequency) mode for wavevectors parallel, or at small oblique angles, to the magnetic field. Their analysis assumed unmagnetized ions, so that an ion trajectory was treated as a straight line. One of the first applications of this instability was in the study of current disruption during substorm expansions in the geomagnetic tail (Lui et al.
Reference Lui, Chang, Mankofsky, Wong and Winske1991). Wu et al. (Reference Wu, Yoon, Ziebell, Chang and Wong1992) subsequently studied the instability more generally, and found that, even when relaxing the condition of unmagnetized ions and small
${\it\beta}_{i}$
, the zero-frequency instability persists, and in contrast to other cross-field drift instabilities (such as the modified two-stream instability), the growth rate remains high even for high
${\it\beta}_{i}$
. The mechanism of the instability is essentially that of the Weibel instability (Weibel Reference Weibel1959) as described by Fried (Reference Fried1959) for counter-streaming electrons. For this reason Chang et al. (Reference Chang, Wong and Wu1990) dubbed it the ion Weibel instability (IWI). However, as Wu et al. (Reference Wu, Yoon, Ziebell, Chang and Wong1992) note, for the case discussed here there is a net equilibrium current and the instability is due to a cross-field drift between electrons and ions, whereas the Weibel instability is usually driven by an electron temperature anisotropy in the low field limit.
A simple derivation of the growth rate can be obtained as follows. Consider the situation of a plasma consisting of a core population of ions and electrons, together with a cross-field ion beam. We take the magnetic field uniform in the
$y$
direction
$\boldsymbol{B}=B_{0}\widehat{\boldsymbol{y}}$
. At the front (upstream) edge of the shock foot, where the specularly reflected ions reach the turning point in their gyration, they form a single, cross-field beam which we take to have velocity
$v_{b}\widehat{\boldsymbol{z}}$
and density
$n_{b}$
. The core ion density is
$n_{i}$
which is equal to that of the upstream flow; the electron density adjusts for charge neutrality. The density of reflected ions is typically of the order of 0.2–
$0.1n_{i}$
. The foot is strongly inhomogeneous due to the flow deceleration, the presence of outgoing and returning reflected ions and the increases of the magnetic field and shock potential. Nevertheless, at its upstream edge, the flow is not yet strongly disturbed, so we can use the upstream flow frame (where the core populations are stationary) for the instability analysis. We assume perturbations
${\it\delta}B_{x}$
and
${\it\delta}E_{z}$
(i.e. linearly polarized) which only have
$y$
and
$t$
dependence, so that they are coupled via the Faraday law
$\partial _{y}{\it\delta}E_{z}=-\partial _{t}{\it\delta}B_{x}$
. By linearizing the beam ion equation of motion, the conservation of number applied to the beam ions and Ampere’s law it is possible to derive a consistent set of equations. With the usual harmonic ansatz, the condition for non-trivial solutions is found to be
${\it\omega}^{2}=-F(k)$
where the function
$F(k)$
is positive for all
$k$
. This corresponds to a zero real frequency (i.e. non-propagating) mode with growth rate given by
$$\begin{eqnarray}{\it\gamma}=\left(\frac{n_{b}}{n_{i}}\right)^{1/2}\left(\frac{v_{b}}{V_{A}}\right)\left(\frac{(k{\it\lambda}_{i})^{2}}{(k{\it\lambda}_{i})^{2}+n_{b}/n_{i}}\right)^{1/2}{\it\Omega}_{i}.\end{eqnarray}$$
Here
$V_{A}$
and
${\it\lambda}_{i}$
are those of the upstream flow. A simplified derivation of the Wiebel instability, without explicit inclusion of the beam density, is presented in Matsukiyo & Scholer (Reference Matsukiyo and Scholer2006).
This simple derivation ignores the response from the core ion and electrons, so that the mode is carried entirely by the cross-field beam. Nevertheless, the dependence of the growth rate on
$k$
is similar to that of the more complete analysis (figure 1, Chang et al.
Reference Chang, Wong and Wu1990). The growth rate
${\it\gamma}\rightarrow 0$
as
$k\rightarrow 0$
, and asymptotes as
$k$
increases to
$$\begin{eqnarray}{\it\gamma}\sim \left(\frac{n_{b}}{n_{i}}\right)^{1/2}\frac{v_{b}}{V_{A}}{\it\Omega}_{i},\end{eqnarray}$$
which is similar to the expression from Chang et al., except for the factor depending on the beam density. (For the comparison, recall that
${\it\omega}_{pi}/c={\it\Omega}_{i}/V_{A}$
.) For the straight line orbit approximation to be valid,
${\it\gamma}>{\it\Omega}_{i}$
, which for
$n_{b}/n_{i}\sim 0.2$
and
$v_{b}/V_{A}\sim 5$
is barely satisfied. However, as shown by Wu et al. (Reference Wu, Yoon, Ziebell, Chang and Wong1992), the instability persists even for magnetized ions. A typical scale for the instability can be estimated from the value
$k_{m}$
at which
${\it\gamma}$
achieves
$1/\sqrt{2}$
of its asymptotic value, giving
$k_{m}{\it\lambda}_{i}\sim 3.2{-}2.2$
for
$n_{b}/n_{i}\sim 0.1{-}0.2$
. This is approximately in the range found in the simulations.
This estimate of the characteristics of the instability depends on some drastic simplifications. The frame used for the analysis is the upstream flow frame, rather than the plasma frame which decelerates through the shock foot. Effects of inhomogeneity in the shock foot are neglected, and is also the fact that the reflected ions will have already had their orbits perturbed before they reach the upstream edge of the foot (at the turnaround point where their
$x$
component of velocity is zero). Also, in the simulations there is evidence of oblique oriented phase fronts and coupling to the other magnetic field components, indicating that probably inhomogeneity, growth of the IWI at oblique angles and coupling to propagating oblique whistlers may play a role.
