Spectral intensities of Anomalous Cosmic Rays derived from the injection rate at the solar wind termination shock

We investigate the well accepted idea that anomalous cosmic ray (ACR) ions originate from suprathermal halo ions by means of diffusive shock acceleration. Usually the seed ions are taken to be the interplanetary pick-up ions, but here we want to enlarge this idea by taking in addition as seed candidates also into account normal solar wind ions which according to our recent calculations are reflected from the termination shock into the upstream solar wind flow and then are also picked-up as suprathermal ions. We start out from an ideally planar shock approximation and fix the ACR spectrum by absolute spectral intensities and maximum ACR energies, taking a precalculated fraction of the suprathermal ion flow as ACR injection rate. Comparison of our calculated spectral intensities with ACR measurements near 94 AU, shows that satisfying data fits only can be achieved, if about three percent of the suprathermal ions swept into the shock structure enter into the Fermi-1 acceleration process. We also show that the spectral slope of the ACR spectra is decreasing and the spectral intensities are increasing with increasing shock compression ratios. As maximum energies available from an ideally operating diffusive Fermi-1 acceleration process we find, depending on the shock compression ratio, ion energies ranging from a few MeV up to 10 3 MeV. Compared to observations this seems to be a little on the high side and may point to the fact that injection into Fermi-1 acceleration near the termination shock is occurring only sporadically due to variable upstream magnetic field orientations with respect to the shock normal vector, i.e. due to variations of the obliquity of the local shock surface with respect to the local upstream magnetic field. Correspondence to: H.-J. Fahr (hfahr@astro.uni-bonn.de)


Introduction: The ACR seed population
Since 1974 anomalous cosmic ray particles (ACRs) are considered as originating from accelerated pick-up ions (PUIs) which undergo Fermi-1 acceleration processes in the region close to the solar wind termination shock (Fisk et al., 1974;Pesses et al., 1981;Potgieter and Moraal, 1988;Cummings and Stone, 1996;Fichtner, 2001;Decker et al., 2006;Mc-Donald et al., 2006;Cummings and Stone, 2007, etc.).As shown by Christian et al. (1988) for instance ACR protons are especially seen at some phases of low solar activity on top of the modulated galactic cosmic ray spectrum in the spectral region between 10 and 100 MeV, amounting in this region (∼ 20 AU) to peak differential flux intensities of the order of 1 to 3 ACRs/(m 2 s MeV sr/nucl).Pick-up ions are thought to be the ideal seed population for ACRs, because in the solar wind frame they are initially injected into an unstable suprathermal ring distribution, and due to pitch angle scattering and preacceleration by Fermi-2 energy diffusion (Isenberg, 1987;Chalov and Fahr, 2000) finally arrive at the termination shock as ion species which, with good efficiency, can get magnetically or electrically reflected from the shock structure (see Chalov andFahr, 1996, 2000) and then undergo further Fermi-1 scattering processes.While the injection efficiencies are not yet precisely predictable (Kucharek and Scholer, 1995;Chalov and Fahr, 1996;Le Roux and Fichtner, 1997;Schwadron and McComas, 2003;Kucharek et al., 2006;Fahr et al., 2008), there exists a kind of tacit understanding that about 5 to 10 percent of the PUIs passing over the shock may enter diffusive acceleration up to ACR energies of the order of 10 to 100 MeV.
Serious hints towards PUIs as really being the seed of ACRs always came from considerations of the elemental abundances.The elemental abundances that appear in the local interstellar medium (LISM) should also reasonably well characterise the elemental abundances of PUIs -besides of heliosheath filtering effects (see Rucinski et al., 1993) -and originate from LISM neutral atoms penetrating into the inner heliosphere.The elemental composition indeed appears to be fairly consistent with that of the neutral component of the LISM (see Cummings andStone, 1990, 1995).One should, however, bear in mind that the ionisation degree of the LISM and filtration functions for the neutral LISM species in fact even nowadays are not well known.This is because the ionisation state of the LISM neither is characterised as a thermodynamical equilibrium nor perhaps even a quasistationary state (see Frisch, 1990;Reynolds, 1990;Frisch and Slavin, 2006;Breitschwerdt and de Avillez, 2006).This uncertainty in the neutral composition provides us with some leeway to explore other candidates for the ACR seed population in addition to PUIs.
In this respect it recently became interesting to us to perhaps also consider normal solar wind protons as a possible seed population at least for ACR protons.In recent papers by Verscharen and Fahr (2008a) and Fahr and Verscharen (2008) we found that, at least for the case of quasi-parallel termination shock conditions, a specific percentage of the normal solar wind ions passing over the shock will become reflected back to the upstream side of the shock due to nonlinear interactions of ions with shock-generated electrostatic turbulence.These reflected ions are also, analogously to PUIs, picked-up by the upstream solar wind and are again swept back into the shock structure after pitchangle redistribution, the fraction of reflected over all solar wind ions being about 15 percent and the most probable injection energy being about 4 times the one for PUIs.In the following we shall study the ACR spectrum as originating from these reflected solar wind ions and compare it with the one expected from PUI-generated ACRs.

Predicted ACR spectral intensities
According to the shock acceleration theory by Gleeson and Urch (1973); Krymskii (1977); Axford et al. (1977); Blandford and Ostriker (1978);Drury (1983), the spectral profile of the shock accelerated ions can be analytically calculated, at least for ideally planar shocks.The absolute value of the upstream distribution function f 1 (hereinafter, the index 1 indicates the upstream and the index 2 the downstream value) has, however, not been specified in this theory and, thus, needs to be fixed for the prevailing shock conditions.For that purpose the absolute spectral intensity calibration factor C in this function has to be fixed such that flux continuity relations at the shock are fulfilled expressing the fact that the total diffusive outflow of the ACR particle fluxes to the upstream and downstream sides of the termination shock, i.e. the sum of the upstream and downstream streamings (see Gleeson and Axford, 1968), respectively, has to be identical with the flux of incoming particles serving as injection seed to start the Fermi-1 process (see Fahr, 1990).
The two-dimensional surface divergence of the differential streaming has to balance the incoming injected particles since the injection at the shock represents a particle source: where U indicates the proton bulk velocity, n 1 the upstream solar wind ion density, and ε 1 is the fraction of the flux of inflowing solar wind ions which are reflected by the shock structure and injected into the Fermi-1 acceleration process, thus serving as seed of the ACR population.The situation is sketched in Fig. 1.
To calculate the particle streaming, the so-called Compton-Getting function C CG (p), developed by Gleeson and Axford (1968), in view of the low-momentum weight, here is approximated with its subrelativistic low-momentum value for the ACR particles.According to Eq. (3.5) from Gleeson (1969), the streaming can be written in the form with the spatial diffusion coefficient κ, kinetic energy T , and for the subrelativistic case.The Compton-Getting factor in this case is given by Hence, the above mentioned requirement of particle conservation can be formulated in the form (see Fahr, 1990;Scherer et al., 2006): Evaluating the above equation with the expression for f 1,2 given in Eq. ( 22) of Scherer et al. (2006) by where the power index q = 3s/(s − 1) is used with the shock compression ratio s = U 1 /U 2 .H (x) is the step function with H (x) = 1 for x ≥ 0 and H (x) = 0 for x ≤ 0. This then leads to the following input-output balance relation: and furthermore to the following result when reminding that at x = 0 (upstream shock position) the upstream and downstream distribution functions, as seen from Eq. ( 5) above, can be considered as identical, i.e. f 1 (x 0 ) = f 2 (x 0 ) = f 0 : This, when evaluating the above expression with the distribution function given by Eq. ( 5), finally leads to: where x max = p max /p 0 and x 0 = 1 are the upper and lower integration borders of the above integral.With the evaluation of the integral and a further re-ordering this then yields the following expression for the spectral intensity calibration factor: The dependence on U 1 in the above formula is hidden in the value of p 0 = p 0 (U 1 ) for the critical momentum of the particle injection into the shock acceleration.In order to inject particles into the diffusive acceleration process, it is necessary that these particles are not simply convected over the electric potential wall of the termination shock, but become reflected at this wall, at least for the first time, and then enter a continuous bouncing process between upstream and downstream flow regimes (see Chalov and Fahr, 1996).
The reflectivity value ε 1 for quasiparallel MHD shocks has recently been determined in the papers by Fahr and Verscharen (2008) and Verscharen and Fahr (2008b) to amount to ε 1 0.03.The associated injection momentum p 0 p ref v , connected with these shock-reflected ions is found to be p 0 2mU 1 .
For most probable compression values of about s 3, for example, one then obtains This points, depending on the value of x max , to a perhaps weak dependence of the absolute spectral calibration factor C on the upper momentum border p max = p 0 x max .This is why in Sect. 3 we shall aim at an estimation of this border.

Transformation into spectral ACR energy fluxes
To compare our calculations with observational data it is appropriate to express the distribution function f 0 (p) given in Eq. ( 5) by a spectral energy flux function (E).For that purpose we use the relation E = p 2 /2m valid for subrelativistic ions and first find which delivers the differential energy distribution in the form where the approximation 1 ≈ 1 is used and the s-dependence of the prefactor is merged into the function The differential energy flux 0 per steradian (calculated with the particle velocity v(E) according to Fig. 2. ACR spectra for different compression ratios.A higher compression ratio leads to a higher and flatter power-law spectrum.
This flux then should be expressed in dimensions The termination shock is located at a distance of about 100 AU, which leads to a particle density of 5 × 10 −4 cm −3 .The upstream bulk velocity is assumed to amount to 400 km/s.

The Fermi-1 acceleration time
The Fermi-1 acceleration period τ acc (p max , p 0 ) needed to process ions from the injection threshold p 0 up to p max is given by Drury (1983) in the following form: where κ 1,2 and U 1,2 denote the spatial diffusion coefficients and the bulk flow velocities upstream and downstream of the shock, respectively.We assume that turbulence levels and diffusion coefficients are identical upstream and downstream of the shock and that these diffusion coefficients can be represented in the form (see e.g.McDonald et al., 1992;Scherer et al., 1998).For a more detailed approach to the spatial diffusion coefficient, it could be of interest beyond the scope of this work to use quasilinear theories (e.g.Forman et al., 1974;Matthaeus et al., 1990;Bieber et al., 1996) or even nonlinear theories (e.g.Forman, 1977;Bieber and Matthaeus, 1997).With our above assumption, one then obtains here x n dx x (18) evaluating with n = 1 (subcritical diffusion for ions with and evaluating with n = 2 (supercritical diffusion for ions with E ≥ 1.0 GeV) to In literature, the threshold between subcritical and supercritical diffusion at E ≈ 1 GeV is sometimes given in terms of the rigidity R = pc/e (Scherer et al., 1998).An ion energy of 1 GeV corresponds to a rigidity of about 0.4 GV.

Fermi-1 differential momentum gain and loss per time
We consider the differential energy gain due to ongoing Fermi-1 ion processing at some maximum momentum border p max and its balancing for the stationary case by leakage of such ions from the shock structure due to particle streamings into upstream and downstream directions.The temporal change of the ion distribution function (∂f/∂t) in the region of the accelerating TS structure is caused by particle processing in p-space (∂f/∂t) F and losses (∂f/∂t) L of these particles by the 1-dimensional configuration space divergence of the particle streaming.This net change, i.e. the leakage of these particles from the acceleration region, is therefore given by

Leakage rate
The leakage rate per momentum space differential volume p 2 dp is determined by the divergence of the particle streamings S corresponding to Eq. ( 6): with L sh as the extent of the leakage region.At the shock position the distribution functions are furtheron identical (f 1 = f 2 = f 0 at x = 0): Taking U 2 in terms of the compression ratio leads to For subrelativistic ions (i.e.kinetic energy T = p 2 /2m 0 ) the coefficient C CG is given by Forman and Gleeson (1975) in the form where the function α(T ) is defined by using the Lorentz factor γ = 1/ 1 − (v/c) 2 .Expressing the kinetic energy in terms of the momentum then leads to We neglect terms higher than third order in (v/c) and find The distribution function f 0 at the shock position (x = 0) is given by (cf.Eq. ( 5)).This spectrum can be inserted into the Eq. ( 27) for the Compton-Getting factor.Taking the partial derivative then leads to Reordering the terms yields Thus, the Compton-Getting factor is given as a function of the momentum and the compression ratio of the shock and can now be used to evaluate the leakage term.

Fermi-1 production rate
The corresponding differential gains per momentum space volume p 2 dp on the other hand are due to the momentumspace divergence of the Fermi-1 accelerative drift flow and thus is given by the following expression Here according to Lagage and Cesarsky (1983) the Fermi-1 accelerative momentum drift can be represented by the expression Insertion of the above relation into Eq.( 32) at p = p max then leads to Assuming the upstream and downstream diffusion coefficients to be essentially identical and reminding that the stationary distribution function at the shock is obtained as f (p, x = 0) = f 0 (p) = C • (p/p 0 ) −q (see Eq. ( 5)), then leads to We take κ = κ 0 v 0 c (p/p 0 ) 2 from Eq. ( 17) because we assume a shock which accelerates the ions to the supercritical diffusion regime.Hence, the evaluation of the derivative delivers Combining and arranging the terms yields Now, balancing losses and gains at p = p max then leads to the balance relation which can be written in the form To estimate the characteristic scale of the leakage region we assume that the following relation is valid, in order to have detailed equlilibrium established at low momenta (C CG = C CG,0 ): This suggests that the needed dimension for the Fermi-1 accelerative region may be suggested to be which then for higher momenta leads us to After insertion of the p-dependent Compton-Getting function derived in Eq. ( 31) the achievable maximum momentum for the Fermi-1 process is found as Introducing as an example s = 3, which means that q = 3s/(s − 1) = 9/2, then delivers meaning that p max 375p 0 or • E 0 (remember: v 0 2U 1 8 × 10 7 cm).This means that the maximum energy amounts to about 470 MeV, which shows that the magnitude of the supercritical energy range can be reached.The maximum values for the achievable particle energy in dependence of the shock strength (i.e. the compression ratio) is shown in Fig. 3. Obviously, there is only a narrow range for possible svalues.For s < 1 + √ 3 2.73 the argument of the square root becomes negative (with an unphysical exception at a certain range of compression ratios).Therefore, there is a cut-off at this position.On the other hand, the singularity at q = 4 leads to a pole of p max at s = 4.At higher compression ratios the argument of the square root is furtheron negative.
We have to distinguish two different kinds of diffusion, i.e. the subcritical and the supercritical.Up to an ion energy of about 1.0 GeV the diffusion coefficient depends linearly on the ion momentum.This branch is called subcritical.Above this threshold, this relation becomes nonlinear and therefore the diffusion coefficient in Eq. ( 17) has a changed slope in its momentum dependence.
Additionally, one has to clarify that due to our choice of n = 2 in Eq. ( 17) this result is only valid for a shock which reaches the conditions of supercritical diffusion (E > 1.0 GeV).For smaller compression ratios and, therefore, at smaller energies this assumption is no longer valid.The calculation with n = 1 for the subcritical diffusion leads straightforward to which has a real solution for (1/2 + √ 13/4) < s < 4.This result is additionally plotted to Fig. 3. Altogether, this means that from s 4 down to compression ratios of 1/2 + √ 13/4 2.3 our approach describes an effective acceleration mechanism.Evidently, there is an intermediate region where neither n = 1 nor n = 2 are valid because there both cases would lead to incorrect maximum energy.Within this gap another slope for the dependence of κ on p must be taken.
With these values the total spectrum can be calculated by use of Eq. ( 10): According to Eq. ( 15) the spectra for different compression ratios are shown in Fig. 2.

Conclusions
In Fig. 2 we have shown theoretically calculated absolute ACR intensity spectra for protons at a solar wind termination shock position of 94 AU, thereby starting from precalculated ion injection efficiencies.This position has meanwhile turned out also to be just the position where Voyager-1, Astrophys.Space Sci.Trans., 5, 21-30, 2009 www.astrophys-space-sci-trans.net/5/21/2009/ at 16 December 2004, has crossed the upwind solar wind shock (see Fisk, 2005;Stone et al., 2005;Decker et al., 2005) and hence a comparison of measured and calculated ACR data seems to be in place here.While the spectral slope of the measured and calculated ACR proton spectra seems to point to most probable compression ratios at the shock of 2.6 ≤ s ≤ 3.0, the measured spectral intensities appear to be lower by a factor of about 10 −1 in comparison to our theoretical intensities calculated for the above range of compression ratios.The question here may thus arise why this deviation by one order of magnitude occurs.
One evident reason emanating from our calculations is connected with the fact that the efficiency value ε 1 0.03 which we have applied in our calculations is theoretically justified for quasiparallel MHD shocks, but is not simply transferable or applicable to quasi-perpendicular shock configurations.The MHD shock configuration that was met by Voyager-1 when it crossed the upwind shock surface under an ecliptic latitude of about ϑ = 30 • must most probably have been much closer to the quasiperpendicular, rather than the quasiparallel one.This means that for a reasonable comparison with Voyager-1 data, injection rates ε ⊥ ≤ ε ε 1 have to be applied to satisfactorily well describe the injection situation for this actual configuration.
The injection efficiency ε 1 used in our present calculations is taken from the relative number of solar wind ions running into the quasi-parallel shock structure and becoming reflected therein due to nonlinear interactions with shockgenerated electrostatic turbulences (see Verscharen and Fahr, 2008a).This process of ion reflection from the shock structure studied by us is only operating with the calculated efficiency, if the upstream magnetic fields are essentially parallel to the shock normal.If the upstream magnetic field would have a tilt α with respect to the shock normal, then at best a reduced amount of this efficiency ε 1 would prevail and also the injection energy would be reduced to values of the order of E 0 (α) 2mU 2 1 cos 2 α.It is, however, very interesting to keep in mind that the case of the quasi-parallel shock will be prevailing for some specific fraction of time even in the upwind shock region, at least in regions which are touched by the neutral sheath of the interplanetary magnetic field (HCS = heliospheric current sheath).Above and below the HCS the magnetic sector polarity changes, i.e. the magnetic field orientation changes by 180 degrees from below to above, also passing through just that magnetic field configuration when the tilt angle α vanishes.According to Voyager data (Burlaga et al., 2003) a sector structure near the ecliptic lasts approximately between 8 and 10 days, while the change to the opposite polarity occurs in about one day.During this day the quasi-parallel configuration at the termination shock prevails, meaning that for about (1/10) of the time even in the upwind part of the termination shock, where usually a quasi-perpendicular shock is realised, the actually prevailing shock configuration will be of the quasi-parallel type, and thus for that fraction of the time our calculated injection efficiency ε 1 should apply.If for the other period of time (nearly quasi-perpendicular type) we neglect the associated injection efficiency, then with a temporal weight of (1/10) for quasi-parallelism to occur, we would just come to the right injection efficiency bringing us to a complete fit of spectral ACR intensities measured with Voyager (Decker et al., 2005).Additionally, the temporal nonsteadiness of the magnetic field geometry at the shock can explain lower observed maximum energies than in our calculation.If the time that is needed to accelerate particles up to the maximum energy E max is longer than the time scale on that the shock can be assumed to be quasi-parallel due to the changing sector structure, the gained maximum energy of the accelerated particles will be lower then the ideally calculated E max .However, this effect cannot be quantified at this moment since the time-dependence of the Fermi-1 acceleration process is not adequately understood at this time.The quasi-parallelity of the shock in our model is needed for the efficient injection and not for the Fermi-1 acceleration itself.The later acceleration can continue also under less favorable magnetic field orientations.
We should mention here that injection to the ACR population up to now was considered to be solely due to PUIs becoming reflected at the shock (see Kucharek and Scholer, 1995;Chalov andFahr, 1996, 2000;McComas and Schwadron, 2006;Fahr et al., 2008).The latter injection efficiency in that case also is highly sensitive to the magnetic tilt angle α realised at the shock, and at angles α ≥ 65 • practically drops to ε 0. In tilt ranges 0 ≤ α ≤ 60 • on the other hand the efficiency for PUI injection turns out to be roughly α-independent with a value of about ε 0.3.Taking again into account the HCS crossing statistics as already mentioned above would then reduce the effective PUI injection efficiency in the upwind part of the termination shock to a value of ε pui 0.07 (see Fahr et al., 2008).To now compare the PUI-induced injection with the one discussed in this paper, i.e. due to reflected solar wind ions, one should compare injection rates β i rather than bare injection efficiencies.For the PUI induced injection one thus should take into account the PUI abundance ξ pui at the shock and take as a rate: which should be compared with the solar wind ion injection rate for quasi-parallel shock configurations which evaluates as shown above to an average number of revealing the fact that the ratio between the two different injection rates is given by:  Fahr and Ruciński, 1999) leads us to the result that the ratio β i,pui /β i, 0.7, meaning that the contributions from injected PUIs and from reflected solar wind ions to the ACR population appear to be about equal.
The theoretical values of maximum ion energies E max originating from the Fermi-1 diffusive acceleration process as it has been studied in this paper are displayed in Fig. 3.This figure shows a strong dependence of E max on the shock compression ratio s and also shows that for compression ratio s ≥ 3 energies of E max (s ≥ 3) of the order of 1000 MeV can be expected under an ideal operation of the shock, i.e. for an ideally planar 1D shock without precursor.Voyager data (see Stone et al., 2005) seem, however, to give indications that maximum energies only of the order of 3.5 MeV can be ascribed to Fermi-1-accelerated ions near the shock.The reason for relatively low Fermi-1 maximum energies seen by Voyager-1 can in first place be ascribed to the fact that relatively low compression ratios of s ≤ 3 have been found from the observed spectral indices.For low compression values, however, also in view of our results shown in Fig. 3, the resulting maximum energies are dropping strongly.
In second place it must also be recognised that the operation mode of the solar wind termination shock in accelerating ions is not the one of an ideally planar shock as studied by Drury (1983) or Malkov and Drury (2001).In a realistic geometry of the solar wind termination shock problem, in addition to a balance between gains due to Fermi-1 acceleration and losses due to diffusive leakage from the shock structure, adiabatic energy losses have to be taken into account which then definitely will reduce the maximum energies presented in our calculation here.The fact that for compression ratios smaller than s = 2.3 our theoretical approach does not deliver reasonable values for E max is connected with the fact that in Eq. ( 28) we have made use of the moderately relativistic expression for the Compton-Getting function.Without keeping the energy-dependence in this function C CG = C CG (E) we cannot find a solution for E max within the frame of our balance consideration.The estimation of E max that we aim at with our approach presented here thus is only reasonable and logically justifyable, if energies E max are found for which an energy-variable Compton-Getting function C CG = C CG (E max ) like the one used in Eq. ( 28) can be adopted.
A comparison with the very recent observations by VOYAGER-2 during the termination shock transit is not conveniently possible due to a different magnetic field orientation.The absolute spectral intensity of the low energy ACR spectra (0.5 MeV -1.5 MeV), though reasonably well represented by our theoretical values, does essentially depend on the local termination shock conditions, i.e. the inclination angle between the local upstream magnetic field and the shock normal.At the VOYAGER-2 crossing this latter angle appeared to be 83 • and, therefore, the shock to be quasi-perpendicular at the transit (Richardson et al., 2008;Burlaga et al., 2008), while the injection rates discussed here in this paper essentially apply to the case of quasi-parallel configurations.
In the VOYAGER-1/2 shock crossings it became evident that the dominant portion of the shock-dissipated upstream kinetic energy reappears downstream in the form of thermal energy of suprathermal ions (i.e.pick-up ions) (Richardson et al., 2008).This does not, however, imply that much higher injection efficiencies than taken into account in our present approach (i.e.ε 1 0.03) are indicated by these most recent data.It only means that hot upstream ions are heated more efficiently than cool ones when passing over the shock, which can nicely be understood as an indication that the magnetic moment of the ions is conserved at the shock passage (see Fahr and Chalov, 2008).
In our considerations, we only take solar wind hydrogen ions into account.Heavier ACR particles cannot be discussed in line with our treatment of these protons because their reflectivities at the shock are different and have not yet been calculated by us.Protons, however, can be considered clearly as the dominant ion species of the upstream solar wind.

Fig. 1 .
Fig. 1.Schematic sketch of the acceleration at the termination shock structure.From the left hand side a certain number of particles are injected into the Fermi-1 acceleration in the shock layer.They are the source for the streamings S 1 and S 2 .
The achievable maximum ion energy in dependence of the shock compression ratio.Only the solutions for supercritical, physically relevant s-values is shown.The energy diverges at s = 4.If the maximum ion energy is smaller than 1 GeV, the subcritical diffusion is used and, therefore, the calculations for n = 1 are additionally shown.In the intermediate range the slope of the diffusion coefficient is not known.