The escape model for Galactic cosmic rays

The escape model explains the cosmic ray (CR) knee by energy-dependent CR leakage from the Milky Way, with an excellent fit to all existing data. We test this model calculating the trajectories of individual CRs in the Galactic magnetic field. We find that the CR escape time τesc(E) exhibits a knee-like structure around E/Z = few × 1015 eV for small coherence lengths and strengths of the turbulent magnetic field. The resulting intensities for different groups of nuclei are consistent with the ones determined by KASCADE and KASCADE-Grande, using simple power-laws as injection spectra. The transition from Galactic to extragalactic CRs happens in this model at low energies and is terminated below ≈ 3 × 1018 eV. The intermediate energy region up to the ankle is populated by CRs accelerated in starburst galaxies. This model provides a good fit to ln(A) data, while the estimated CR dipole anisotropy is close to, or below, upper limits in the energy range 1017 – 1018 eV. The phase of the dipole is expected to change between 1 × 1017 and 3 × 1018 eV.


Introduction
The all-particle cosmic ray (CR) energy spectrum is a nearly featureless power-law between ∼ 10 10 eV and ∼ 10 20 eV, with only a few breaks in its spectral index. The two most prominent ones are the knee at E k ≈ 4 PeV, and the ankle at E a ≈ 4 EeV. Another related open question is where the transition from Galactic to extragalactic CRs occurs. For the knee, two main explanations currently remain possible. First, it may be the signature of the maximum energy to which Galactic CR sources can accelerate protons [1]. Second, the knee could be caused by a change in the energy dependence of the CR diffusion coefficient and thence confinement time in the Galaxy [2,3,4], if the CR Larmor radius is the order of the coherence length l c of the turbulent Galactic magnetic field (GMF) at E k . In Ref. [4], we have studied this possibilitywhich we denote as the "escape model"-by propagating individual CRs in recent GMF models. We showed that the escape model is viable and can explain the individual fluxes of CR groups as measured by KASCADE and KASCADE-Grande. Moreover, our estimate for the dipole anisotropy in this model was consistent within uncertainties with observations.
Here we extend our previous study [4] and present preliminary results for the entire energy region between 300 GeV/Z and the ankle [5]. In particular, a more detailed study of the Galactic CR primary composition between E k and E a is presented here and compared to observations. We show that any remaining heavy nuclei flux in the sub-ankle region would be dominated by only one or few local sources. We use limits on the iron fraction at > ∼ 7 × 10 17 eV determined by the Auger collaboration together with ln(A) measurements to constrain the transition energy between Galactic and extragalactic CRs, deducing R max = E max,Fe /(26e) ∼ 10 17 V as the maximal rigidity R max to which Galactic CR sources are able to accelerate CRs. The recovery of the proton and helium spectra above E/Z ∼ 10 16 eV in the KASCADE-Grande data is maninly explained by the specific shape of the escape rate τ esc (E) discovered in [4]. We show also that observational constraints from anisotropy limits are compatible with the escape model. A natural extension of the escape model to other normal galaxies suggests that the intermediate energy region up to ankle is composed of CRs accelerated in starburst galaxies. In both cases, observational constraints from anisotropy and composition measurements are compatible with the escape model.

Galactic magnetic field models and CR confinement in the Galaxy
An important constraint on CR propagation models comes from ratios of stable primaries and secondaries produced by CR interactions on gas in the Galactic disk. In particular, the B/C ratio has been recently measured by the AMS-02 experiment up to 670 GeV/nucleon [6]. In order to take advantage of the high-quality and the large energy range of the B/C data from AMS-02, we use now these data to derive the grammage traversed by CRs as function of their energies in two simple models. In Fig. 1, the two sets of grammages deduced are shown with magenta and blue error-bars. Note that the error-bars take into account only the statistical and systematic errors of the AMS-02 measurement, while uncertainties in the cross sections or deficiencies of our approximations are not accounted for. The latter can be estimated by the differences between the results from the two approximations.
In order to compare these measured values to those predicted in the escape model, we inject N cosmic rays at z = 0 in the Galaxy and follow their trajectories x i (t) until they reach the edge of the Galaxy. Then we calculate the average grammage X = N −1 c N i=1 dt ρ(x i (t)) summing up the density along the trajectories of individual CRs. Since the grammage X(E) ∝ E −δ scales as the confinement time τ (E) ∝ E −δ , this quantity serves also as an indicator for changes in the CR intensity induced by a variation of the CR leakage rate.
In Fig. 1, we compare the grammage calculated from simulated CR trajectories with the grammage deduced from B/C measurements. We use the JF model for the Galactic magnetic field [7], choosing as the maximal length of the fluctuations L max = 5l c = 10 pc, and consider two values of its root mean square (rms) strength, the original one suggested in [7] (β = 1) and a second one rescaling it to one tenth of its original value (β = 1/10). (Similiar results are obtained in the GMF model of Ref. [8].) Because of the large energy reach of the AMS-02 data, the extrapolation required from the lowest energy of our numerical calculations, E = 10 14 eV, to the measurements has shrunken to two orders of magnitude. Using the JF model with β = 1 as proposed in [7] would require a constant power spectrum of magnetic field fluctuations, P(k) ∝ k −α with α = 0, in the intermediate energy range. Such a power-spectrum is difficult to reconcile with the theoretical understanding of turbulence. Moreover, the CR spectrum is very close to a power-law above ≃ 200 GV. This implies that if D(E) would become significantly flatter beyond TeV energies (e.g. changing from D(E) ∝ E 1/3 to ∝ E 0 ), then the injection spectrum of sources has to have the exact opposite change of slope (e.g. respectively from ∝ E 2.4 to ∝ E 2.7 ). Alternatively, a change in the source density should compensate the change in D(E) such that the observed CR intensity remains a nearly featureless power-law [9]. Although such a conspiracy cannot be excluded, it appears to us as a not very attractive option.
We choose a Kolmogorov power-spectrum P(k) ∝ k −5/3 as the theoretical model with the smallest slope α, and have therefore to reduce B rms by a scaling factor β < 1. Next we examine how the shape of the grammage X as function of energy E/Z depends on the two parameters l c and β. In [4], we discovered a specific shape of X(E) that leads not only to a knee-like feature but reproduced also the recovery of the proton and helium spectra above E/Z ∼ 10 16 eV, visible in the KASCADE-Grande data. It is clear that a too strong turbulent field, β ∼ 1, results in knee-like feature at too high energy. Compensating a relatively strong turbulent field by decreasing the coherence length tapers off both the knee-like feature and the recovery. As a consequence, the allowed range of turbulent field strengths and coherence lengths is correlated and very restricted, l c ≃ (1 − 10) pc and β ≃ 1/10 − 1/8.  shown with the data from CREAM [11], KASCADE [16], KASCADE-Grande [17], TIBET [18], and Auger [19] experiments.

Fluxes of Galactic CR groups
For the calculation of the CR flux at Earth, we adapt a procedure based on the use of templates as decribed in detail in Ref. [5]. Results for nuclei with charge Z are deduced from the calculations for protons by shifting the energy by a factor Z. We then interpolate the resulting CR nuclei fluxes to the same energies as for protons. At energies below Z × 100 TeV, we assume that diffusion in the Kolmogorov turbulence shifts the injection power-law ∝ E −α by 1/3 to the spectrum ∝ E −α−1/3 observed. This is indeed what we observe in our simulations in the energy range ≃ Z × (100 − 300) TeV. Therefore we assume that the CR spectrum of protons released by sources follows a power-law spectra ∝ E −2.4 ; the maximal energy E p of protons will be fixed later by considering constraints form the resulting dipole anisotropy and the observed nuclear composition of CRs. For all other nuclei, we use power-law spectra with either ∝ E −2.17 or ∝ E −2.22 and maximal energy ZE p . These power-law indices are chosen so as to fit the direct observations from CREAM at low energy. We fix the density of sources by normalizing the flux found in our simulations to the observed one at 100 TeV. On average, we require 130 sources per 100 kyr for a total energy per source of E tot = 10 50 erg, so as to fit the observed CR spectra. In Fig. 3, we plot the CR nuclei fluxes, multiplied by E 2.5 , as a function of energy. In the upper left and upper right panels of Fig. 3, we show the proton and helium fluxes, both for turbulent fields with L max = 10 pc and with L max = 25 pc. We plot orange lines ∝ E −2.4−1/3 (upper left panel) and ∝ E −2.22−1/3 (upper right panel), which represent the slopes expected theoretically at Earth, for our injection spectra with α = −2.4 and −2.22 as power indices and Kolmogorov turbulence. Note that the slopes of the injection spectra required for nuclei, α ≃ 2.2, coincide with the naive expectations from diffusive shock acceleration. Only the proton injection spectra requires a somewhat harder slope, α = 2.4, than expected.
(10 − 30) PeV region, where errorbars of both experiments are relatively small. In contrast, the helium flux from KASCADE-Grande is below the one measured by KASCADE . This behavior is explained by the insufficient discrimination power between protons and helium in the KASCADE-Grande experiment [13]. For this study, we then choose to reduce the proton flux of KASCADE-Grande by 40%, and add this difference to the helium flux, in same energy bins. By doing so, the CR fluxes of KASCADE-Grande and KASCADE experiments become consistent with each other. In the lower left panel of Fig. 3, we plot the CNO flux, which predominantly consists of carbon and oxygen. We calculate the carbon and oxygen fluxes by normalizing them to the CREAM fluxes interpolated to higher energies with power-laws, and then sum them up. The CREAM flux in this figure is the sum of its carbon and oxygen fluxes, where we use carbon energy bins for the binning, and interpolate the oxygen flux to these bins before summing up. KASCADE and KASCADE-Grande measurements of the CNO flux are directly compared to our fluxes. In the lower right panel of Fig. 3, we show the flux of heavy nuclei, which is dominated by Mg, Si and Fe nuclei. The combined 'heavy nuclei flux' measured by the KASCADE experiment is smooth and agrees well with a simple power-law extrapolation of the CREAM flux to higher energies. It also agrees with the KASCADE-Grande flux. As for the other components, the model presented in this work fits well the heavy nuclei flux too. As can be seen in Fig. 3 (upper right), the CR proton flux follows a power-law from 300 GeV up to about 1 PeV. It then changes to a steeper slope at the knee, and recovers at ≃ 10 PeV to a  flatter power-law with index α ≃ 2.5. Similar 'knee-like' cutoffs, shifted by factors Z in energy, are visible in the fluxes of all groups of CR nuclei-see the other panels of Fig. 3. These plots demonstrate that the escape model fits very well all these observations. As discussed previously in [4], the knee is due, in this model, to a change in behaviour with energy of the CR diffusion coefficient. The energy of the knee corresponds to the energy at which the Larmor radius of CR protons is of the order of the coherence length of the turbulent magnetic field (l c = L max /5 for a Kolmogorov spectrum). For the field strengths we consider in this paper, B rms ≃ 0.3 µG (or β = 0.1/8) close to the solar position, we find in our calculations a change in the slope of the CR flux at about 1 PeV, as observed in the proton data. Let us note that, in this model, the flatter part of the CR proton flux above ≃ 10 PeV is dominated by recent nearby sources. This is due to the fact that the confinement time of CR protons in the Galaxy quickly drops with energy beyond the energy of the knee. While the fluxes from the few nearby sources are obviously subdominent below the knee, they tend to explain most of the measured flux above 10 PeV. As should be expected in such a scenario, fluctuations with time (on ≃ (10 − 100) kyr time scales) of the CR proton flux at Earth are large above 10 PeV-what explains the large error bars of the theoretical prediction at high energies.

Transition from Galactic to extragalactic CRs
Determining at which energy E * the CR flux starts to be dominated by extragalactic sources is one of the most important unsolved problems in CR physics. While in the 'dip model' of Refs. [14] the transition energy is as low as E * ≈ a few × 10 17 eV, in other models (e.g. [15]) it is expected to occur at E * ≈ a few × 10 18 eV. In our model, the maximum rigidity R max to which Galactic sources are able to accelerate CRs is a free parameter. The energy of the transition E * in our model can be determined by using additional observational constraints. One possibility is to constrain the maximum contribution of Galactic sources to the total CR flux by using the observational limits on the anisotropy of the CR flux. Another way to determine E * , is to study the elemental composition of primary CRs and use the fact that the composition of Galactic and extragalactic CRs should in principle differ from one another.
We start with the latter method and restrict our discussion to the case that the Auger iron limit is obeyed. As a first step, we derive the all-particle CR flux summing up all CR groups and compare it to the experimental data of KASCADE [16], KASCADE-Grande [17], TIBET [18], and Auger [19]. Then we deduce the extragalactic flux for a given R max by subtracting the predicted total Galactic flux from the measured total CR flux. The resulting extragalactic flux is shown in the right panel of Fig. 1 with a red dashed line for R max = 1.0 × 10 17 V. Next, we have to fix the nuclear composition of the extragalactic CR flux. As a first approximation, we can assume that its composition is constant in a sufficient small energy interval around E * . In contrast, the Galactic CR composition is strongly energy dependent between the knee and the cutoff of the Galactic flux. Thus, we expect that an observable like the average of the logarithmic mass number, ln(A), will be quickly changing for E < ∼ E * , while being approximately constant for energies slightly above E * .
In the left panel of Fig. 4, we plot measurements of ln(A) from several experiments, together with the values of ln(A) calculated within the 'escape model' studied here. The points for KASCADE have been computed by converting the flux measurements given in [12] into ln(A) values. The most striking feature, namely the peak in ln(A) around 5 × 10 16 eV, is clearly visible in all data sets, although its exact position and strength depend on the experiment. Our model reproduces the trend in the data very well. At higher energies, the composition becomes lighter because of the 'flattening' of the escape time at such energies, see Section 2. For the value of R max we consider here, extragalactic CRs start to contribute to the observed flux at ≈ 10 17 eV. Consequently, above this energy, the exact value and shape of ln(A) depends on the assumed composition of the extragalactic flux. In blue, we show ln(A) for an extragalactic flux made of protons only, in magenta for a mix of 50% p and 50% Fe, and in red for a mix of 60% p, 25% He, and 15% N. Independently of the composition chosen for the extragalactic component, we can identify the energy where ln(A) stops to decrease with the maximum energy E max,Fe to which Galactic sources can accelerate iron. It corresponds to the rigidity R max = E max,Fe /26e. This transition is clearly visible in the PAO data, around the ankle, and allows us to determine the maximum rigidity as R max ≃ 1 × 10 17 V. For the case of a mixture of 60% p, 25% He, and 15% N (red curve), we obtain a good agreement with the ln(A) data from PAO up to 2 × 10 18 eV. While this choice of composition is not unique, it is consistent with the results from the recent composition study published in [22]. In particular, Ref. [22] found the fraction of iron to be below 20% above 6 × 10 17 eV. This agrees well with the results of our model, where the Galactic flux at 6 × 10 17 eV consists purely of iron but contributes to only 15% of the total CR flux.
In addition to fitting the above observables, we still have to verify that the model presented here is also consistent with the existing upper limits on the CR anisotropy. In the diffusion approximation, the CR dipole anisotropy d is given by d = 3D ∇ ln(n)/c. Following the same procedure as in [4], we compute the average anisotropy and derive the energy dependence of D(E) from the escape rate as calculated previously, setting D(E/Z) ∝ 1/τ esc (E/Z). We fix the proportionality constant by requiring that the dipole amplitude d = k f k d k equals the dipole componentd observed by the EAS-TOP collaboration at E = 1.1 × 10 14 eV [23]. Here, k labels the groups of nuclei we consider in the Galactic flux plus an extragalactic component. The latter has a dipole amplitude which is independent of its composition and which we set equal to 0.6%, as expected for the extragalactic Compton-Getting effect [27]. The factor f k corresponds to the fraction the component k contributes to the total CR flux, and d k ∝ 1/τ esc (E/Z) to their individual dipole.
In the right panel of Fig. 4, we show the resulting dipole amplitude d as a function of energy E. As expected, the amplitude raises below the knee as E 1/3 , while it increases approximately as E 0.7 until 1 × 10 17 eV. At higher energies, the dipole amplitude decreases, which is due to the facts that the Galactic composition becomes heavier and that the extragalactic contribution grows. We also plot the values ofd observed by IceCube [24], as well as the 99% C.L. upper limits on d ⊥ from the Pierre Auger Observatory [25]. Comparing our estimate for the dipole amplitude with the upper limits in the energy range 10 17 − 10 18 eV, we should take into account that the approximation d ∝ 1/τ esc (E/Z) starts to break down above E/Z 10 17 eV, which leads to a sizeable error. We conclude therefore that our prediction is marginally consistent with these limits. The Pierre Auger Observatory should however be able to reach a detection of the dipole anisotropy. Let us also note that the escape model predicts that the phase of the dipole amplitude varies strongly in the energy range between 1 × 10 17 and 3 × 10 18 eV: This corresponds to the range where the transition from Galactic to extragalactic CRs lies. Such a picture is supported by current observations of the phase of the dipole, see References [23,24,25]. In summary, there are two reasons for having an early transition, from predominantly Galactic to predominantly extragalactic CRs, at E ≈ a few × 10 17 eV. First, the limits on the observed dipole anisotropy requires either a very heavy Galactic composition or a predominantly extragalactic contribution at E > ∼ 10 18 eV [28,25]. The former possibility is however strongly disfavored by the recent composition measurements from the Auger collaboration [21,22]. Second, identifying the energy where ln(A) stops to decrease with the maximum energy to which Galactic sources can accelerate iron, E max,Fe ≈ 3 × 10 18 eV, suggests that the maximal rigidity reached in Galactic sources satisfies R max = E max,Fe /(26e) ∼ 10 17 V. Note that the cutoff in the Galactic CR spectrum would be at lower energies, if no nearby active CR source exists [5].

Conclusions
We have shown in [4,5] that the knee can be entirely explained by energy-dependent CR leakage from the Milky Way, with an excellent fit to all existing data from E/Z ∼ 300 GeV to 100 PeV. In particular, all deviations from a single power-law behavior that are observed in the CR intensity of individual CR groups in this energy range are consistently explained by rigidity-dependent CR escape. This model requires small coherence lengths of the turbulent field and relatively small turbulent magnetic fields. If these two conditions are fullfilled, then the CR escape time τ esc (E) exhibits a knee-like structure around E/Z = few × 10 15 eV together with a recovery around E/Z ≃ 10 16 eV.
In this presentation, we have considered only the case that a nearby CR source is present (for the opposite case see [5]). Then the maximal rigidity R max to which Galactic CR sources are able to accelerate CRs can be determined by identifying ZeR max with the energy where ln(A) derived from PAO measurements stops to decrease. This allowed us to deduce R max = E max,Fe /(26e) ∼ 10 17 V. As a result, the flux ratio of Galactic and extragalactic sources is in this case 1:1 at E * ≈ 2 × 10 17 eV, dropping to 0:1 at 2 × 10 18 eV. Since the transition from Galactic to extragalactic CRs happens in this model at rather low energies, the estimated CR dipole anisotropy is consistent within uncertainties with upper limits in the energy range 10 17 − 10 18 eV, while it reproduces the measurements at lower energies from EAS-TOP and IceCube. The dipole phase is expected to change between 1×10 17 and 3×10 18 eV, i.e. the energy range of the transition from Galactic to extragalactoc CRs. Such a behavior corresponds to the one observed, providing thus additional evidence for a transition from Galactic to extragalactic CRs in this energy region.
The ankle is interpreted in the escape model as the transition to a new, second extragalactic source class, as e.g. active galactic nuclei or gamma-ray bursts. Moreover, this model suggests that the CR knee in starburst galaxies is shifted by up to two orders of magnitude to higher energies [29]. Therefore, the CR flux in the intermediate energy region up to ankle should be composed mainly of CRs accelerated in starbust galaxies.
Let us note also that the weakness of the turbulent GMF in the "escape model" would have important consequence for the search of UHECR sources: UHECRs from a single source would be mainly deflected by the regular component of the GMF, while the spread of their arrival directions due to the turbulent GMF should be small. As a result, the search for UHECR sources, at least in the case of protons or light nuclei, should be more easily than thought before. Even for heavier nuclei, the deflections in the regular field of the Galaxy can be traced back in those patches of the sky with small turbulent fields [30]. Weaker magnetic fields will also simplify the search of nuclei sources using the methods discussed in Refs. [31,32]. Thus the results of this work are an additional motivation for future searches of UHECR sources, performed by future all-sky missions as e.g. JEM-EUSO [33] and KLYPVE [34].