© European Geosciences Union 2005

Abstract. Acceleration processes associated with the heliospheric termination shock may provide a source of anomalous cosmic rays (ACRs). Recent kinetic simulations of supercritical, quasi-perpendicular shocks have yielded time varying shock solutions that cyclically reform on the spatio-temporal scales of the incoming protons. Whether a shock solution is stationary or reforming depends upon the plasma parameters which, for the termination shock, are ill defined but believed to be within the time-dependent regime. Here we present results from high phase space resolution particle-in-cell simulations for a three-component plasma (solar wind protons, electrons and pickup protons) appropriate for the termination shock. We find reforming shock solutions which generate suprathermal populations for both proton components, with the pickup ions reaching energies of order twenty times the solar wind inflow energy. This suprathermal "injection" population is required as a seed population for subsequent acceleration at the shock which can in turn generate ACRs.


Introduction
With the approach of the Voyager spacecraft towards the heliopause and in particular the termination shock (Burlaga et al., 2003;Gurnett et al., 2003), together with possible crossings thereof (Krimigis et al., 2003;McDonald et al., 2003), the problem of the acceleration of anomalous cosmic rays (ACRs) (Fisk et al., 1974;Pesses et al., 1981) in the outer heliosphere is at the forefront of current research (see, for example, Kucharek and Scholer, 1995;Zank et al., 1996;Rice et al., 2000;Zank et al., 2001;Lever et al., 2001).This provides an opportunity for comparing theoretical results with observations.
Correspondence to: R. E. Lee (leer@astro.warwick.ac.uk)Kinetic simulations are a useful tool to investigate, selfconsistently, local acceleration processes at the shock.Particle-in-cell (PIC) simulations of quasi-perpendicular shocks (see, for example, Biskamp and Welter, 1972;Lembège and Dawson, 1987;Lembège and Savoini, 1992;Shimada and Hoshino, 2000;Schmitz et al., 2002a;Hada et al., 2003;Lee et al., 2004a), have found that, for a certain range of parameters, the shock solutions are not static, rather reforming on the gyro scales of the incoming bulk protons.Hybrid simulations, on the other hand, typically generate static solutions (see, for example, Burgess et al., 1989;Kucharek and Scholer, 1995;Lipatov and Zank, 1999), although reforming solutions are recovered for a range of values of resistivity (Quest, 1985;Hellinger et al., 2002).The time dependent electromagnetic fields at reforming shocks have been found to accelerate a fraction of the bulk proton population to suprathermal energies (Lee et al., 2004a), and this is a possible cosmic ray injection mechanism at supernova remnants.Parameters relevant to the heliospheric termination shock, although ill defined, are consistent with those expected to yield reforming solutions in a quasiperpendicular geometry.The use of a quasi-perpendicular geometry is appropriate if we make the assumption that Parker's solar wind model is still valid in the outer heliosphere.Following Pesses et al. (1981) one concludes that the shock is close to perpendicular (the angle between magnetic field and shock normal >80 • ) at heliographic latitudes less than 75 • .Pioneer 11 measurements show the existence of two broad distributions of shock propagation angle, centred on 90 • and 270 • , with half widths of 25 • (Smith, 1993).
Importantly, the heliospheric termination shock is generated in the presence of pickup ions originating from the transheliopause interstellar medium.Two questions then arise, which we address here: first, whether the reforming shock solutions found previously persist with the self consistent inclusion of a pickup ion population; and second, to what extent are both solar wind and pickup protons accelerated in the resulting time dependent electromagnetic fields.
We use a 1.5-D (1x, 3v) relativistic electromagnetic PIC code to simulate the structure and evolution of a supercritical, collisionless, perpendicular magnetosonic shock.Simulations of reforming shocks in quasi-perpendicular and pureperpendicular reforming shocks have been performed, for example, by Scholer et al. (2003), where it is shown that the additional wavemodes supported by oblique geometry do not strongly affect the overall structure and dynamics of the reforming shock.In a PIC simulation, the distribution functions of all particle species are represented by computational super-particles, the phase space locations of which are evolved via the Lorentz force law, whilst the electromagnetic fields are defined on a spatial grid and evolved via the full Maxwell equations (see, for example, Birdsall and Langdon, 1991).The simulation setup is that described in Lee et al. (2004a), with the addition of pickup protons, the results of which are presented for the first time here.We find reforming shock solutions in this 1.5-D geometry, where scalar and all three components of vector quantities are functions of one spatial coordinate and time.Reforming shock solutions have also been found in higher dimensions (Lembège and Savoini, 1992).
To accurately reproduce the effects of pickup protons, we simulate three species within our code: electrons, solar wind protons, and pickup protons.Particles are injected on the left-hand side of the simulation box with a drift speed v inj .The simulation is conducted in the downstream rest frame, so if v 1 is the upstream and v 2 the downstream flow speed derived from the Rankine-Hugoniot relations, v inj can be defined from the Alfvénic Mach number, M A by: v inj =M A v A −v 2 .For solar wind protons and electrons, each individual particle has its initial speed modified by an additional Maxwellian random velocity, chosen from a thermal distribution characterised by a v thm for each species.Pickup protons are born with similar values of perpendicular velocity and arbitrary gyrophase, centred on the drift speed.The pickup protons are modelled by a population with velocities in the form of a spherical shell, centred on v inj , radius v 1 .On the particle injection boundary, the magnetic field (B z,1 ) is constant and the electric field (E y,1 ) is calculated self consistently.We use the piston method to generate the shock, and shock-following algorithms, to obtain a sufficiently long run time, as described in Schmitz et al. (2002a).
Here results are presented from simulations of a perpendicular shock propagating into a collisionless plasma, at an Alfvénic Mach number (M A ) of 8.The upstream ratio of plasma thermal pressure to magnetic field pressure for the solar wind thermal protons is taken as β i =0.2, and for the electrons β e =0.5.The pickup proton effective temperature is characterised by v 1 due to their initial shell distribution (after Kucharek and Scholer, 1995;Zank et al., 1996;Ellison et al., 1999).To enable ion and electron time scales to be captured within the same simulation, with reasonable computational overheads, it is necessary to reduce the simulation mass ratio for ions and electrons to M R =m i /m e =20, and the ratio of electron plasma frequency to electron cyclotron frequency ω pe /ω ce =20 (in common with Shimada and Hoshino, 2000;Schmitz et al., 2002a,b;Lee et al., 2004a).Simulations at higher mass ratios have been conducted by, for example, Lembège and Savoini (2002) at M R =400, and Scholer et al. (2003) at M R =1840; however, to achieve reasonable runtimes and simulation domains, a ratio of ω pe /ω ce =[2, √ 8] was used, for example, in the latter.In our simulations a lower mass ratio is chosen to allow for sufficient resolution of the ion phase space to study the ion dynamics.Scholer et al. (2003) demonstrate that reformation is a low β process, and is not an artifact of unrealistic mass ratios, and that the differences that occur between simulations at M R =20 and M R =1840 are in the nature of the instabilities in the foot region.Previous work suggests (see, Lee et al. 2004a, Lee et al. 2004b) that the reflected ion dynamics are insensitive to these instabilities, and importantly we resolve the same length scales for the foot region (∼λ ci ) and shock ramp (of the order of c/ω pe ), thus we expect our conclusions to hold at M R =1840.

Results
The key phenomenology found in our PIC simulations arises from the non-time-stationary nature of the reforming shock.This is shown in Figs. 1 and 2; these plot magnetic field strength and E x dx as functions of x and time for two simulations, on the left, for solar wind protons only and on the right, with the addition of a pickup proton population at 10% of the density of the background.As with all figures in this paper, they are presented in a frame where the downstream plasma is at rest, using data collected from a time segment of the simulation when the shock has formed and is propagating independently of the boundary conditions.Units are normalised to upstream parameters for the solar wind protons, thus λ ci is the upstream solar wind proton cyclotron radius, ω ci is the upstream solar wind proton cyclotron frequency, and E inj,SW =m i v 2 inj /2 is the solar wind proton injection energy.The shock can be seen to move upstream, from the right to the left of the simulation box, over time, with peaks in magnetic field (the shock ramp) recurring, and the potential well expanding and collapsing, on the time scale of the local proton cyclotron period.Over the course of each cycle, a stationary shock ramp forms at the turnaround point of protons in the foot region.Then a new foot region extends into the upstream region as protons reflect from this new shock front.A new ramp (magnetic field peak) then starts to form at the upstream edge of this foot region, and another cycle begins.The last of the original foot region protons gyrate into the downstream plasma as the new shock starts to reflect a new foot region population.The shock thus propagates upstream (left) in a stepwise fashion.Comparison of a shock simulation containing no pickup ions (panel A in Figs. 1 and  2, and see also, Lembège and Savoini, 1992;Shimada and Hoshino, 2000;Schmitz et al., 2002a;Lee et al., 2004a) to a shock simulation containing 10 % pickup protons (panel B  in Figs. 1 and 2), shows that the reformation spatio-temporal scales are similar.
Figure 3 shows a snapshot of the spatial distributions of xvelocity and kinetic energy for both proton species, together with E x dx, and perpendicular magnetic field, at the beginning of a reformation cycle (t=34.5×2πω−1 ci ).The temporal evolution of both proton species in phase space over a reformation cycle is shown in Fig. 4, together with the corresponding scalar E x dx and perpendicular magnetic field in Fig. 5. From this figure we see that the addition of 10% pickup protons does, however, lead to an extended foot region, observable by the slow decline into the shock's potential structure starting ∼6λ ci upstream, and in the downstream acts to suppress magnetic field fluctuations over a distance of ∼5λ ci downstream of the shock front.
In Fig. 3 protons from both species flow in from the lefthand, injection boundary.Solar wind protons continue until they reach the shock front (currently situated at ∼75λ ci ).At this stage of the reformation cycle, the shock front is almost stationary (see Fig. 1), so that a fraction of the solar wind protons specularly reflect, in the downstream rest frame, from the ramp in the magnetic field and its accompanying potential, with v x →−v x (panel 5).The reflected solar wind pro- tons move upstream to form a relatively energetic population, as shown in panel 4, gaining energies dispersed from 1−6×E inj,SW (as seen in Lee et al., 2004a) as they cross the shock front and move into the downstream region.
Pickup protons also flow in until they reach the shock front, at which stage they, too, are reflected (panel 3).The presence of pickup protons leads to an extended foot region with two components: the reforming foot region associated with the solar wind protons close to the shock, plus a weak foot region extending up to ∼6λ ci upstream.The relative spatial scales of these foot region populations may be understood in terms of specular reflection at the shock ramp.If the protons simply reflect specularly, that is, reverse their component of velocity normal to the shock front (v n ) as they encounter the ramp, then the length-scale for each population is L foot =|v n |/ω ci .This is governed by v n =v inj + thm sin(ω ci t+φ): for thermal protons this implies v n v inj because v thm,SW v inj ; whereas for the pickup protons v thm,PI v inj , so that 0<v n <2v inj .Thus, after specular reflection, the pickup protons create an extended foot region by moving further back into the upstream region, with a corresponding larger excursion in y.The evolution of the phase space of both proton species over the course of a reformation cycle is shown in Fig. 4; this is accompanied by Fig. 5 which shows the magnetic field and E x dx at the same times.Visible in E x dx in Fig. 5 is the weak potential feature, associated with the pickup protons, that extends ∼6λ ci upstream of the potential feature associated with the solar wind protons.This extended potential changes little over a reformation cycle, particularly when compared to the changes visible in the potential associated with reflecting solar wind protons.The maximum upstream positions for reflected pickup protons and solar wind protons are indicated by arrows on Fig. 4 and on the E x dx panels in Fig. 5, to highlight this point.The presence of an extended foot region is also visible to a lesser degree in the magnetic field profiles; these also clearly show the nontimestationary nature of the shock ramp.When it is possible to do so (panels A, C and D) we observe a shock ramp thickness of the order of a few electron inertial lengths (in common with Lembège and Savoini, 1992;Scholer et al., 2003).This value of the shock thickness can also be seen in Fig. 6 which shows magnetic field strength and E x dx profiles for a simulation in the absence of pickup protons, at four times over the course of a reformation cycle comparable to those shown in Fig. 5. Comparing Figs. 5 and 6 shows E x dx and magnetic field structure associated with the bulk solar wind protons, that are qualitatively similar to that arising in the presence of pickup protons.From the magnetic field profiles it can be seen that the field is B z /B z,1 ∼3, on average, in the downstream region.Taking the average time between minima in the shock's foot region value of E x dx, the reformation time is ∼1.5×2π ω −1 ci in the simulations with both 0% and 10% pickup protons, as can be seen in Fig. 2. In terms of the average field just downstream of the overshoot, the reformation time is then ∼0.5 local ion gyroperiods.
The differences in temporal evolution for the two species can be seen when the counts of reflected protons are compared.We count the number of protons upstream of the minima in E x dx with, in the case of solar wind protons, velocities directed away from the shock front and in the case of the pickup protons, those protons with velocities outside of the injected shell in velocity space.The time variation of these reflected populations is shown in Fig. 7.In contrast to the solar wind protons, the pickup proton reflected ion count (and consequently foot region structure) is stable over a reformation cycle, with slight changes coinciding with times of changing position for the minima in E x dx as the reforming shock advances.The solar wind proton reflection count, on the other hand, varies by a factor of ∼2.5 over the course of a reformation cycle.
To understand the energy gain processes, it is useful to consider the evolution of F•vdt, and its components, calculated along proton trajectories.Specifically, in Fig. 8 we examine a solar wind proton that reaches ∼6E inj,SW , and a pickup proton that journeys ∼8λ ci upstream after reflection, before returning and reaching high energy on crossing the shock front.To relate these changes to the shock dynamics, the proton trajectories are superimposed upon the spatiotemporally evolving E x dx in the bottom panel of Fig. 8.The corresponding local E x dx experienced by these two protons is plotted in the middle panel of Fig. 8. Detailed analysis of the energetics of the solar wind protons (Lee et al., 2004a) suggests that the electromagnetic fields relevant to acceleration to suprathermal energies are on ion scales.As discussed in Lee et al. (2004a,b), for example, the fluid equations treated appropriately lead to defining the bulk potential in terms of integral E xidx , then removes features, such as electron phase space holes, but maintains ion scale features, such as the foot region potential structure.The rate of change of energy eE•v integrated over the trajectories of energetic ions with the electric field from Eq. ( 1) is then found to track that calculated using the full electric field from the PIC simulation.Previous PIC simulations for a range of values of m i /m e (see, for example, Lembège and Savoini, 1992;Scholer et al., 2003) show a variety of kinetic instabilities in the foot region, and in our results, features, possibly related to the Buneman instability, are visible in the foot region of the solar wind ion phase space in Figs. 3 and 4.These features occur on elec-  tron spatio-temporal scales, and here we simply note that, as shown in Lee et al. (2004b), the ion dynamics will, from the above argument, be insensitive to structures on electron kinetic scales.
The proton trajectories in Fig. 8 display both similarities and differences during their shock interactions.Both particles encounter the shock when the magnetic field is at its strongest (Fig. 1), and the potential at its narrowest (panel 5).They are reflected off the shock front, moving upstream into the foot region gaining energy (panels 1 and 2) in the local time dependent electric field.Their gyration and the propagation of the shock then brings them back towards the shock front, where their energies are too high for a second reflection, and they enter the downstream region (panel 5).
The x-components of F•vdt for both protons are similar (red, panels 1 and 2 of Fig. 8).Initially, as the protons move towards the shock, there is an x-energy loss of roughly  equal magnitude (∼E inj,SW ).After reflection this is followed by an x-energy gain, also of roughly equal magnitude, as the protons move out into the foot region.The magnitudes are similar because the protons experience similar potential structures.However, there is a slight disparity because the pickup proton interacts with the shock front and transfers some energy from its x-component to its y-component (at t∼30.3×2πω−1 ci ), whereas this process is absent from the solar wind proton trajectory which remains within the foot region (t=30.3 to 32.8×2πω −1 ci in panel 5).Examining the y-components of energy (blue, panels 1 and 2 of Fig. 8), the origin of the disparity in final energies becomes clear.The pickup proton moves into the shock at a higher velocity, and is therefore reflected at that higher velocity back into the upstream region (compare panels 3 and 5 of Fig. 3).Consequently, the distance travelled by pickup protons in the foot region, upstream of the shock, is greater than that for the solar wind protons.It follows that the drift in the y-position is greater (panel 4), so that the gain in the y-component of energy, through E y motional drift, is greater.Examining E x dx at the particle positions (panel 3), we see that the overall potential structures around the times of yenergy gain are similar for both particles.Energisation occurs during a phase when the local potential is approximately constant: again, as the pickup proton travels further in this state, and as the y-energy gain is roughly linear, it gains more energy.As the particle passes through the foot region for a second time, the energisation finishes and the particle crosses the shock front.Consequently, solar wind protons are seen at energies up to and beyond ∼6E inj,SW and pickup protons at up to ∼20E inj,SW .The self-consistent inclusion of pickup ions leads to additional downstream magnetic field structures, that can lead to further acceleration of the pickup protons, as indicated by the scattering around x 83λ ci in the pickup proton energies in Fig. 3.
To illustrate changes in energisation mechanisms that occur due to the presence of pickup ions, the normalized distribution functions of the proton energies, f (E), for each species are plotted on semi-log axes verses energy, normalized to the solar wind proton injection energy, in Figs. 9 and 10.The distribution functions shown here were collected over ten upstream proton gyroperiods from a region starting at x 12λ ci downstream of the shock front, and extending for x 5λ ci .These distances are sufficiently far downstream to be in a region of strongly fluctuating fields and over spatio-temporal scales sufficiently large as to average over the small-scale variations present within the downstream populations.For comparison the distribution functions in the upstream region are also shown.These space and time averaged distributions capture the extended tail of f (E) which cannot be seen in the snapshots of phase space (Figs. 3  and 4).
Figure 9 shows the solar wind proton distribution functions.Comparing the distribution functions of the simulations containing no pickup protons to one containing 10% pickup protons shows only minor differences.The peak in the downstream distribution functions is at zero, as the plots refer to the downstream rest frame.The distribution functions decay exponentially, with changes in exponent at ∼2E inj,SW , the energy gained by protons that specularly reflect form the stationary shock front (see, Lee et al., 2004a), and at ∼5E inj,SW , close to the maximum energy observed for protons as they leave the shock front and enter the downstream region.Energisation can also be seen to be weaker in the presence of pickup protons.
Examining the pickup proton downstream distribution function (Fig. 10, blue), we see a peak at around 2E inj,SW , the maximum energy of the inflowing pickup protons.The distribution function varies slowly on this log axis up to ∼17E inj,SW , this value is the maximum energy observed for pickup protons leaving the shock front and entering the downstream (Fig. 3, panel 3).At higher energies the distribution decays exponentially, which is suggestive of acceleration in the downstream region.

Conclusions
We have performed PIC simulations for a quasiperpendicular shock in a parameter range where the shock solutions are reforming and which is appropriate for the heliospheric termination shock.We include a shell distribution of 10% pickup protons, centred in velocity space on v inj , with a radius v 1 25×v thm,SW .
1.The addition of 10% pickup protons does not modify the solar wind proton phase space and dynamics of the reforming shock found previously (Lembège and Savoini, 1992;Shimada and Hoshino, 2000;Schmitz et al., 2002a;Lee et al., 2004a).
2. The dynamics of the energised solar wind proton population proceeds as in simulations where pickup protons are absent.Typical energies reach ∼6E inj,SW and beyond, however, energisation is weaker in the presence of pickup protons.
3. A subset of the pickup protons reflect off the shock to form a weak foot region which stands upstream of the shock and is quasi-stationary.
4. The reflected pickup protons gain E∼20×E inj,SW due to this dynamics.This acts as a suprathermal population which could then act as a seed population for subsequent acceleration at the shock, which can in turn generate ACRs.
5. There is some evidence of further acceleration of pickup protons downstream of the shock ramp.
The energisation seen in our simulations takes place in a broad volume spanning the shock, its extended upstream foot, and downstream.Both the degree of energisation, and the extent of the region in which it occurs, are substantial.This implies that some of these effects may be observable, depending upon instrumental availability and resolution, and may assist the interpretation of heliopause data.
Finally, our simulations, whilst generating a suprathermal population of protons, do not directly generate ACRs.The reforming shock solutions, under the restrictions of these simulations, provide an "injection" or seed population for other mechanisms, such as diffusive Fermi acceleration (Bell, 1978;Malkov and Drury, 2001).Importantly, the reforming shock solutions have an effective reflection probability for the protons which differs from that of a static solution (Lee et al., 2004b).Whether or not a solar wind proton is reflected depends critically upon when it encounters the shock ramp, rather than where it is located in the phase space (Burgess et al., 1989, compare).This may have implications for processes such as shock surfing (as conjectured by Scholer et al., 2003) and multiple reflection (as found for static shock solutions by Zank et al., 1996;Lipatov and Zank, 1999).In this context Zank et al. (1996) have noted that the shock ramp thickness is a critical parameter.In common with previous simulations (see, for example, Shimada and Hoshino, 2000;Schmitz et al., 2002a) and with the assumption of (Zank et al., 1996), the shock ramp width in the simulation presented here is ∼2c/ω pe .However, the extension of these PIC simulations to higher dimensions, whilst retaining reforming shock solutions (Lembège and Savoini, 1992), will allow wavemodes with k vectors oblique to the magnetic field rather than strictly perpendicular, as in our 1.5-D simulations.The occurrence of multiple reflection thus remains an open question.

Fig. 1 .
Fig. 1.Evolution of perpendicular magnetic field strength, over time (vertical axis), and space (horizontal axis), for a simulation in the absence of pickup ions (panel A) and a simulation containing 10% pickup ions (panel B).

Fig. 3 .
Fig. 3. Spatial cross section of the simulation in the shock region at t=34.5×2πω −1 ci .The top panel (panel 1) shows perpendicular magnetic field strength, B z , normalised to the upstream value, B z,1 .Panel 2 shows φ= E x dx normalised to the solar wind proton injection energy E inj,SW .Panels 3 and 5 show the energies of the pickup protons and the solar wind protons, respectively, again normalised to E inj,SW .Panels 4 and 6 show the pickup and solar wind proton x-velocities, normalised to the injection velocity v inj .

Fig. 4 .
Fig. 4. Phase space (v x /v inj vs. x/λ ci ) for solar wind protons (lower panel of each pair) and pickup protons (upper panels) in the shock region.Four pairs of panels at equal time intervals of 0.4×2πω −1ci , starting at time t=33.3×2πω−1 ci (lowest two panels).The maximum upstream position for the reflected ions in each species are indicated for each snapshot.

Fig. 5 .
Fig.5.Four pairs of panels, corresponding to the spatio-temporal coordinates of the panels in Fig.4, showing perpendicular magnetic field strength (upper panel of each pair), and φ= E x dx (lower panels) in the shock region.As in Fig.4, the maximum upstream positions of reflected protons for each species are indicated, with the furthest upstream arrow relating to the pickup ion species.

Fig. 6 .
Fig. 6.Four pairs of panels, showing perpendicular magnetic field strength (upper panel of each pair), and φ= E x dx (lower panels) in the shock region for a simulation in the absence of pickup ions.The panels are separated by equal time intervals of 0.4×2πω −1 ci , starting at time t=33.3×2πω −1 ci (lowest two panels).

Fig. 7 .
Fig. 7. Reflection rates over a series of reformation cycles.Solar wind protons are in black, and pickup protons in blue.Counts are of particles upstream of the minima in E x dx associated with the shock front, with velocities, in the case of solar wind ions directed in the direction of shock propagation, or for pickup ions, outside of the initial injected spherical distribution.Pickup ion counts are normalised to solar wind counts.

Fig. 8 .
Fig. 8. Top panel (panel 1) shows the energy of a pickup proton, calculated from F•vdt and its components (x−red, y−blue, z−green, total−black), panel 2 the energy of a solar wind proton.Panel 3 shows E x dx at the position of both protons (solar wind−black, pickup−purple), panel 4 the y-positions, calculated from y= v y dt, and panel 5 the x-positions overlayed onto the shock potential on ion scales (plotted as E x,i dx).

Fig. 9 .
Fig.9.Solar wind proton energy distribution functions, kinetic energy normalised to E inj,SW vs. normalised frequency on a log 10 scale.Results are shown both in the absence of pickup protons (black) and with 10% pickup protons (blue).For comparison, the upstream distribution functions are also shown, in the absence of pickup protons in red and with 10% in magenta.

Fig. 10 .
Fig. 10.Pickup proton energy distribution functions, kinetic energy normalised to E inj,SW vs. normalised frequency on a log 10 scale.Colours are as in Fig. 9.