electrodynamic simulations

Abstract. We have performed 3 one-dimensional full particle electromagnetic simulations of a quasi-perpendicular shock with the same Alfven Mach number M A ~5, shock normal-magnetic field angle Θ Bn =87° and ion and electron beta (particle to magnetic field pressure) of 0.1. In the first run we used an ion to electron mass ratio close to the physical one ( m i /m e =1024). As expected from previous high mass ratio simulations the Modified Two-Stream instability develops in the foot of the shock, and the shock periodically reforms itself. We have then self-consistently included in the simulation 10% pickup protons distributed on a shell in velocity space as a third component. In a run with an unrealistically low mass ratios of 200 the shock still reforms itself; reformation is due to accumulation of specularly reflected particles at the upstream edge of the foot. In a third run including pickup protons we used a mass ratio of 1024. The shock reforms periodically as in the low mass ratio run with a somewhat smaller time constant. The specular reflection of pickup protons results in an increase of the shock potential some distance ahead of the shock foot and ramp. The minimum scale of the cross shock potential during reformation is about 7 electron inertial length λ e . We do not find any pickup proton acceleration in the ramp or downstream of the shock beyond the energy which specularly reflected ions gain by the motional electric field of the solar wind during their upstream gyration.


Introduction
It is generally accepted that interstellar neutral atoms, which are ionized in the heliosphere and picked up by the solar Correspondence to: M. Scholer (mbs@mpe.mpg.de)wind, are ultimately accelerated at the termination shock by first order Fermi acceleration, also known as diffusive shock acceleration, up to several 100 MeV nucleon −1 (Pesses et al., 1981) and constitute the so-called anomalous cosmic rays (ACRs).If pickup ions are sufficiently energetic to begin with, diffusive acceleration at a collisionless termination shock can accelerate them to the energies of the ACRs.The acceleration process has to be fast, i.e., in order for ACRs to be singly ionized the acceleration time for 10 MeV nucleon −1 is limited to ∼4.6 yr (Jokipii, 1992).Jokipii (1992) pointed out that such a short acceleration time can only be achieved by diffusive acceleration at the quasiperpendicular termination shock assuming weak (in the hardsphere sense) scattering perpendicular to the magnetic field.However, the requirement for the scattering to be weak requires that pickup ions are already energetic in order to be Fermi accelerated, since otherwise particles downstream of the shock are not capable of diffusing upstream (Jokipii, 1987), (Webb et al., 1995).Zank et al. (2001) derived as a condition for particles to be Fermi accelerated that the injection velocity v is larger than 3(U 1 /r)(1+η 2 ) 1/2 where U 1 is the upstream solar wind speed, r the shock compression ratio, and η=λ || /r g , the ratio of the particle's parallel mean free path λ || to the gyroradius r g .Since for weak scattering η 1, the energy at which particles are injected into a Fermi process is larger than the energy of the pickup ions, and this is the so-called injection problem. Lee et al. (1996) and Zank et al. (1996) have investigated the possibility that pickup ions are accelerated at a quasiperpendicular shock by so-called shock surfing, where part of the pickup ions incident on the shock are trapped between the electrostatic potential of the shock and the upstream Lorentz force.The efficiency of this mechanism strongly depends on the scale of the cross-shock potential: the surfing mechanism, also called multiply reflected ion (MRI) mechanism, is only efficient when the cross-shock potential scale is of the order of the electron inertial length.Zank et al. (2001) S. Matsukiyo et al.: Quasi-perpendicular shocks have proposed that the MRI mechanism works as a pickup ion injection mechanism into a first order Fermi acceleration mechanism at the termination shock.The MRI mechanism discriminates against heavier pickup ions: the reflection fraction of pickup ions decreases with increasing mass.This is at variance with the conclusion drawn from ACR observations in the outer heliosphere: the injection efficiency for ACRs seems to be higher for higher mass pickup ions (Cummings and Stone, 1996).To circumvent this problem Zank et al. (2001) have argued, that since according to quasi-linear resonant scattering theory the parallel mean free path is an increasing function of rigidity R=pc/M, η decreases with pickup ion mass M (p is particle momentum, Q is charge, c is speed of light).Since in the weak scattering limit the injection velocity for diffusive acceleration has to be larger than 3(U 1 /r)(1+η 2 ) 1/2 it follows that the injection velocity is for heaver ions at a lower value, where fluxes of MRI accelerated pickup ions are higher.This overcompensates for the lower acceleration efficiency of the MRI mechanism with increasing mass.
The crucial question is whether the scale of the crossshock potential at a quasi-perpendicular collisionless termination shock of the order of the electron inertial scale, a prerequisite for the MRI mechanism.Scholer et al. (2003) have performed full particle simulations of quasi-perpendicular shocks with the physical mass ratio and have found that low ion beta, β i , shocks periodically reform.The length scale of the potential during reformation cycle can, at best, get as small as about 4 electron inertial lengths λ e =c/ω pe (ω pe is the electron plasma frequency).They have therefore concluded that the MRI process is not a viable mechanism for pre-acceleration of pickup ions at the termination shock, although Lipatov and Zank (1999) found such an acceleration process in their finite electron mass hybrid shock simulations where pickup protons were included.We will comment on the Lipatov and Zank (1999) result at the end of the paper.An additional problem for the MRI mechanism is the existence of higher dimensional effects on the structure of the shock surface, like the shock ripples seen in the hybrid simulations by Lowe and Burgess (2003).
In the full particle simulations by Scholer et al. (2003) pickup protons were not included as a third component.Furthermore, they based their arguments on runs with rather low ion beta (magnetic field to particle pressure).However, the process of shock self-reformation and the associated shock transition scale depends critically on β i .The Voyager data show a strong deviation of the temperature in the inner as well in the outer heliosphere from an adiabatic temperature profile (Richardson and Smith, 2003).This is attributed in the inner heliosphere to stream interaction and shocks and in the outer heliosphere to turbulent heating by interstellar pickup ions.At ∼80 AU the temperature is about by a factor 60 above the adiabatic value.Taking a value of 1.5×10 4 K, and assuming a magnetic field strength of 0.1 nT and a density of 1.6×10 −3 cm −3 we obtain β i ∼0.083.The magnetic field strength is closer to 0.06 nT (Burlaga et al., 2003).According to the Wang et al. (2000) model for the evolution of the solar wind the density changes in the outer heliosphere synchronously with the magnetic field.Assuming thus a reduced density of 0.8×10 −3 results in β i ∼0.11.We assume that β i =0.1 is a representative value upstream of the heliospheric termination shock.We want to explore in the following the consequences of pickup protons for a quasiperpendicular shock.We assume in the present paper a magnetic field -shock normal angle of Bn =87 • .The exactly perpendicular shock is a singular case since in a 1-D simulation waves with a component of the wave vectors parallel to the magnetic field are excluded.For instance a 1-D shock simulation of an exactly perpendicular shock does not exhibit the Modified Two-Stream instability (MTSI) in the shock foot region (Scholer et al., 2004).
Recently Lee at al. (2005) have for the first time included pickup protons self-consistenlly in a full particle perpendicular shock simulation with parameters appropriate for the termination shock.Two important parameters enter a PIC (Particle-In-Cell) simulation: the mass ratio m i /m e and the ratio of electron plasma frequency to gyrofrequency ω pe / ce .Lee at al. (2005) used the large value of 20 for ω pe / ce , but had to compromise because of computational reasons by using the rather low mass ratio of m i /m e =20. Lee at al. (2005) found that a perpendicular shock with 10% pickup protons added also reforms.Due to the dynamics of the shock the reflected pickup protons are accelerated up to 20 times the solar wind proton energy.Since we found in previous PIC simulations that the use of the physical mass ratio is important for the dynamics of low β i quasi-perpendicular shocks we will present in the following results for high mass ratio simulations.However, in order to achieve reasonable run times and simulation domains a ratio of ω pe / ce =2 had to be used.

Simulation
The shock is produced by the so-called injection method: a high-speed plasma consisting of solar wind electrons, solar wind protons, and pickup protons 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 proton species as well as for the electrons in a computational cell.As in earlier work (Liewer et al., 1993), (Kucharek and Scholer, 1995)  the solar wind, which neglects adiabatic deceleration and velocity space diffusion.Since the shock moves in the simulation system with speed M s (in units of the Alfvén speed) from the right toward the left the shock Mach number is M A =M o +M s where M o is the injection velocity of the solar wind protons.In order to initialize and subsequently inject from the left hand boundary pickup protons for a termination shock simulation M s has to be known: the pickup protons have to be injected on a sphere in velocity space centered at M o with radius M A .In the simulations the injection velocity M o is assumed to be 3.5.This produces a shock which moves with M s ∼2.3 in the simulation frame to the left hand side, so that the shock Mach number is M A ∼5.8 shock.
The size of a cell is one Debye length λ D .In the following, time will be given in units of the inverse of the proton cyclotron frequency ci , distances in units of the electron inertial length λ e , 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 potential e is given in units of cB 0 /λ e .As outlined in the Introduction, we will investigate a low beta case and assume β i =β e =0.1.Because of computational reasons the parameter τ =(ω pe / ce ) 2 is set to 4. We will return to the problem of small τ in the Discussion section.
2.1 High mass ratio (m i /m e =1024), no pickup protons We will first discuss a high mass ratio m i /m e =1024 simulation run without the addition of pickup protons.Figure 1 shows the magnetic field component B z stacked in time; time runs from bottom to top, and beginning and end is indicated at the y-axis.As can be seen the magnetic field profile is not steady, but the foot and ramp is highly structured.The ramp reforms itself on a time scale of about 2 −1 ci , although it has to be noted that in the high mass ratio case we can follow the shock development only over three reformation cycles.How- ever, in contrast to earlier lower ion β i simulations presented in Scholer et al. (2003), there exists already an extended foot immediately after the new ramp has build up.Figures 2 and  3 show from top to bottom the magnetic field B z component, the ion density n i and the ion v ix −x phase space plot at two particular times of the simulation (t ci =5.73 and t ci =6.0, respectively).In Fig. 2 the shock ramp is at ∼390λ e , there is an extended foot region with specularly reflected ions in which the incoming solar wind protons are decelerated and the magnetic field is increased.The ion phase space plot in the bottom panel of Fig. 2 shows structure in the phase space distribution of the incoming ions.These ions interact with the incoming electrons: due to the high density of reflected ions the incoming electrons are decelerated relative to the incoming ions in order to achieve zero current in the shock normal direction.Due to the relative velocity between solar wind ions and solar wind electrons the solar wind ion beam mode can interact with the whistler mode, as shown in more detail by Matsukiyo and Scholer (2003), which results in the MTSI.Phase mixing of the reflected ions and the incoming ions due to the MTSI turbulence leads to solar wind ion thermalization in the foot region.Eventually a new shock ramp emerges at the upstream edge of the foot, as can be seen from Fig. 3 at a somewhat later time.Similar simulations have been discussed in detail in Scholer et al. (2003).
2.2 Low mass ratio (m i /m e =200), pickup protons included We have included self-consistently 10% pickup protons on a shell in velocity space.As described above the particle number of solar wind protons and pickup protons per cell is identical; the relative contribution of pickup protons is scaled down by assuming the appropriate mass and charge.First we present results from a low mass ratio (m i /m e =200) run. Lee  (2005) found in a mass ratio m i /m e =20 run that selfreformation indeed occurs.They have argued that the shock dynamics during reformation is important for the acceleration of the specularly reflected pickup protons in the foot region.Figure 4 shows the magnetic field stacked in time over an extended time period (longer run times can be achieved in this low mass ratio case).One can see repeated reformation cycles with a period of ∼1.5 −1 ci .The development of the magnetic field during one of these cycles can be seen in more detail from the lower panel of Fig. 4. Reformation is due to accumulation of specularly reflected solar wind protons at the upstream edge of the solar wind proton foot as discussed by Lembège and Savoini (1992) and Hada et al. (2003).In this lower mass ratio run the MTSI is absent.Figure 4 also shows that an extended second foot exists in front of the solar wind proton foot due to the reflected pickup protons.Self-reformation is independent of the mass ratio and the presence of (up to 10%) pickup ions, although the details of processes in the shock foot region during reformation strongly depend on these.Since Lee at al. (2005) also found reformation in simulations with ω pe / ce =20 it seems that reformation is also independent of the ratio of speed of light to Alfvén speed (because ω pe / ce =(c/v A )(m e /m i ) 1/2 ).

High mass ratio (m i /m e =1024), pickup protons included
We now increase the ion to electron mass ratio to 1024 and include 10% pickup protons.Initially there is a total number of 2.4×10 6 particles in the system, but particles are continuously added so that by the end of the run the total number of particles has doubled.In the upper panel of Fig. 5 is shown the stacked magnetic field profile; the lower panel exhibits time-stacked profiles of the cross shock potential = E x dx.Reformation clearly occurs also in this case with a somewhat smaller cyclic time period of ∼1.7 −1 ci and is even more pronounced than in the case without pickup ions.The lower panel shows that the potential increases a considerable distance in front of the reflected solar wind proton generated foot.This is due to reflected pickup ions and pickup ions which are first transmitted and then re-enter the upstream region due to their large gyroradius.During their gyration in the upstream field these pickup ions contribute to an ion bulk velocity v y in the direction perpendicular to the field and to the shock normal.As pointed out by Lee at al. (2005) the electric field in the shock normal, x-direction can be approximated by The pickup ion bulk velocity in the y-direction results in an increase of the potential far ahead of the specularly reflected solar wind proton foot.In the foot the magnetic field rapidly increases (see top panel of Fig. 5) leading to a further increase of the cross shock potential .Figure 6 shows the magnetic field profile and the cross shock potential during two instances in time.At t −1 ci =4.9 there is a large foot in the magnetic field followed by a steep ramp.At this time more than 80% of the total potential change occurs across the foot.The potential decreases slightly before the ramp and has only a small increase in the ramp.At t −1 ci =5.5 a new ramp has formed at the upstream edge of the former foot and a very small foot in the magnetic field exists.At this time ∼1/3 of the potential occurs in the extended region upstream due to reflected and gyrating pickup protons, ∼1/3 of the potential occurs in the small new foot region and ∼1/3 in the ramp.This is contrary to the synchronous behavior of and B when no pickup ions are present.The scale of the ramp potential during this time period is ∼6−7λ e .
Figure 7 shows reduced distribution functions of pickup protons f pi (v x ), f pi (v y ), and f pi (v z ) in the region downstream of the shock at time t=7.1 −1 ci close to the end of the run.The distribution is averaged over a distance of 200λ e ; the total number of particles used for this averaged distribution is about 2×10 6 .Pickup protons gain maximum velocities of about 15 v A .The downstream distribution function does not exhibit a high energy tail, in contrast to what has been seen in the low mass ratio, high ω pe / ce ratio simulations by Lee at al. (2005).Also no further high energy tail in the distribution function of reflected pickup ions is either seen upstream or far downstream (not shown here).
In Fig. 8 we present the energy of the pickup ions normalized to the solar wind energy in in the downstream rest frame as a function of distance x at ci t=7.1.Also shown in the top panel for reference is the profile of the magnetic field B z component.In the lower panel we have plotted the normalized energy of each individual pickup ion in the system.As can be seen, the maximum pickup ion energy is below pi / in ≈17, i.e., at this time there are no ions with higher energy in the system.In the Appendix we give a simple analytic derivation of the maximum energy obtained by a pickup ion during the gyration in the foot.The analytically obtained value of pi / in =14.1 is close to the value obtained in the simulation.

Conclusions
We have presented in this paper 3 full particle simulations of a quasi-perpendicular shock with the same Alfvén Mach number M A and shock normal -magnetic field angle Bn but different ion to electron mass ratios, and with or without contribution of pickup protons.The parameters are such that they might be characteristic for the heliospheric termination shock.The results of the present study can be summarized as follows.
1. Without the addition of pickup protons in a high ion to electron mass ratio simulation the Modified Two-Stream instability (MTSI) occurs in the foot of the shock and leads to important changes of the shock structure as compared to low mass ratio simulations.
2. In a low mass ratio (m i /m e =200) simulation reformation still occurs on about the same time scale (∼2 −1 ci , but the magnetic field in the foot is smooth and the MTSI is absent.Reformation is due to accumulation of specularly reflected ions at the upstream edge of the foot.A feedback effect leads to a magnetic field increase, further deceleration of incoming ions and density increase, and eventually to a new shock ramp.This mechanism has been described by a semi-analytical model for exactly perpendicular shocks by Hada et al. (2003) .
3. In the high mass ratio simulation without pickup ions and with 10% pickup protons included the MTSI occurs.Reformation is still due to accumulation of specularly reflected ions at the upstream edge of the foot as described above.
4. A large part of the cross shock potential drop occurs already in the upstream region, where reflected and/or escaping pickup protons gyrate in the upstream magnetic field.
5. The cross shock potential in the ramp has scales varying between ∼6 and ∼40λ e during a reformation cycle.
Small scale lengths occur when there is a clear separation between foot and a steepened up ramp.However, during such times only one third of the total potential change occurs across the ramp.Approximately one third of the potential drop occurs over the extended upstream region with gyrating reflected pickup protons and one third occurs in the foot associated with the reflected solar wind protons.
6. Reduced distribution functions of pickup ions do not indicate that there exists any acceleration process beyond the energy gain of reflected pickup protons during upstream gyration by the solar wind convection electric field.
Some of the results of this paper are at variance with previous work.Lipatov and Zank (1999) found in finite electron mass hybrid simulations by including pickup protons selfconsistently that these ions were accelerated to higher energies by the shock surfing process.However, as shown in the review by Lembège et al. (2004), large differences can occur in the pickup proton spectra obtained from such simulations depending on whether the electrons are treated adiabatically or by implicitly solving the electron energy equation.
In the latter case considerable artificial electron heating occurs through the shock ramp due to resistivity.This leads to a sharp increase of the electron pressure and, in turn, of the potential, through the shock ramp.In their PIC simulations of perpendicular shocks with pickup protons included Lee at al. (2005) did not find any acceleration of the pickup protons due to shock surfing.However, these authors found acceleration of pickup protons up to 20 times the solar wind proton energy, which they attributed to the shock dynamics during self-reformation.This is not seen in the present simulation either.
There are a number of reasons for the lack of pickup ion acceleration found here as compared to the Lee at al. (2005) result.The lack of acceleration may either be due to (1) the high, realistic mass ratio used here, (2) the small, unrealistic value of ω pe / ce used here, or (3) the combination of both.At a small value of ω pe / ce electrostatic effects are suppressed and strong electric fields possibly occurring on electron scales are reduced or absent.This is definitively a drawback of the present high mass ratio simulations.On the other hand, ion dynamics is expected to be independent of electric fields occurring on electron scales.Thus we conclude that small scale electric fields present in high ω pe / ce simulations have no effect on ion dynamics.However, in low mass ratio, high ω pe / ce simulations it might well be that the low mass protons can still be accelerated in small scale electric field structures.Such an acceleration process would disappear in large ω pe / ce simulations if the physical mass ratio were used.
Other limitations of the present simulations are the limited system size and run time.Lee at al. (2005) found a high energy tail in the downstream pickup proton distribution, which they attributed to acceleration downstream of the shock.Since we have rather short run times we would probably not observe such a stochastic acceleration mechanism.Furthermore, the downstream region is still rather small, so that our particle statistics is poor.But any process dealing with the dynamics downstream of a quasi-perpendicular shock is anyway beyond what can be inferred from 1-D shock simulations.The downstream large scale dynamics in a 2-or 3-dimensional system is dominated by waves propagating parallel to the magnetic field, i.e., almost perpendicular to the x-axis.Such waves are not allowed for in 1-D simulations.The same holds in the ramp region: instabilities with k vectors less oblique to the magnetic field may have larger growth rates and can modify the reformation process.In particular the modified two-stream instability between solar wind electrons and reflected ions with k vector components parallel to the magnetic field can occur in a multi-dimensional spatial system (Gary et al., 1987) and plays an important role in the foot region (Matsukiyo and Scholer, 2006).
How can then the injection problem at quasi-perpendicular shocks be attacked?As recently pointed out by Giacalone (2005) there might not really exist an injection problem.In 2-D hybrid simulations of quasi-perpendicular shocks Giacalone (2005) found that part of the thermal protons are reflected and move upstream along magnetic field lines that are multiply connected to other locations on the shock surface.The latter is due to long-wave lengths fluctuations superimposed on the upstream magnetic field.Such a reflection and acceleration mechanism is also possible, or even more likely, for the pickup protons.The proposal by Scholer (2004) that fluctuations with a length scale of the order of 0.1 AU lead to a locally quasi-parallel shock, which preferentially reflects and accelerates pickup ions, is on the same line.It is likely that in the case of ion acceleration large-scale higher dimensional effects are more significant than the small scale dynamics of the shock.It would nevertheless be very important to verify whether a micro-scale ion acceleration process like the one reported by Lee at al. (2004) exists.Computer resources to do this with mass ratios sufficiently high (≈200) to clearly separate ion and electron spatial and temporal scales and with ω pe / ce ratios sufficiently high (≈10) so that the speed of light and the Alfvén speed are sufficiently far apart are now in reach.

Fig. 1 .
Fig. 1.Magnetic field B z component stacked in time for the ion to electron mass ratio 1024 run.Time runs from bottom to top.

Fig. 2 .Fig. 3 .
Fig. 2.From top to bottom: magnetic field B z component, ion density n i , and ion v ix phase space versus shock normal direction x for the ion to electron mass ratio 1024 run.

SFig. 4 .
Fig. 4. Upper panel: Magnetic field B z component stacked in time for the ion to electron mass ratio 200 run.Lowe panel: Same for a single reflection cycle.

Fig. 5 .
Fig. 5. Upper panel: Magnetic field B z component stacked in time for the ion to electron mass ratio 1024 run with pickup protons included.Lower panel: Same for the cross shock potential .

Fig. 6 .
Fig. 6.Magnetic field B z component and cross shock potential versus shock normal direction x for two different times during the reformation cycle.

SFig. 7 .
Fig. 7. Reduced pickup proton distribution functions downstream of the shock averaged over 200λ e .From top to bottom is shown f pi (v x ), f pi (v y ), and f pi (v z ).

Fig. 8 .
Fig. 8. Magnetic field B z component (top panel) and normalized pickup ion energy pi / in (bottom panel) versus shock normal coordinate x.