{'fraction_non_alphanumeric': 0.03482342774089641, 'fraction_numerical': 0.02083906194759791, 'mean_word_length': 4.196468667074913, 'pattern_counts': {'":': 1, '<': 4, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cosmic Ray Anisotropy as Signature for the Transition from Galactic to Extragalactic Cosmic Rays

We constrain the energy at which the transition from Galactic to extragalactic cosmic rays occurs by computing the anisotropy at Earth of cosmic rays emitted by Galactic sources. Since the diffusion approximation starts to loose its validity for E/Z>~ 10^(16-17) eV, we propagate individual cosmic rays (CRs) using Galactic magnetic field models and taking into account both their regular and turbulent components. The turbulent field is generated on a nested grid which allows spatial resolution down to fractions of a parsec. Assuming sufficiently frequent Galactic CR sources, the dipole amplitude computed for a mostly light or intermediate primary composition exceeds the dipole bounds measured by the Auger collaboration around E ~ 10^18 eV. Therefore, a transition at the ankle or above would require a heavy composition or a rather extreme Galactic magnetic field with strength>~ 10 muG. Moreover, the fast rising proton contribution suggested by KASCADE-Grande data between 10^17 eV and 10^18 eV should be of extragalactic origin. In case heavy nuclei dominate the flux at E>~ 10^18 eV, the transition energy can be close to the ankle, if Galactic CRs are produced by sufficiently frequent transients as e.g. magnetars.


Introduction
The question at which energy the transition from Galactic to extragalactic cosmic rays (CRs) takes place is one of the major unresolved issues of cosmic ray physics. Two promising possibilities are to associate the transition with one of the two evident features of the cosmic ray spectrum: The second knee around E ≃ 5 × 10 17 eV or the ankle at E ≃ 3 × 10 18 eV. Since the chemical composition of galactic and extragalactic CRs should differ in general, both because of propagation effects and of the different nature of their sources, the transition may be detected experimentally studying the chemical composition of CRs as function of energy.
In the case of a transition around the second knee, Galactic CR sources such as e.g. supernova remnants would accelerate CRs up to the rigidity-dependent knee, which is close to 10 17 eV for iron. If the extragalactic CR flux dominating at higher energies would consist mainly of protons, the ankle could be explained as a dip in the extragalactic CR spectrum due to the pair-production losses of protons on cosmic microwave background (CMB) photons p+γ CMB → p+e + +e − [1]. Below ∼ 10 17−18 eV, the extragalactic CR flux may be suppressed because of CR propagation in extragalactic magnetic fields [2,3]. On the other hand, the scenario of Ref. [4] would favour a transition at the ankle. The composition of the CR flux at high energies is the subject of current debate due to the facts that hadronic physics must be extrapolated from lower energies and that the complex experimental analyses for different experiments are not yet completely reconciled. The scenario of Ref. [1] is supported by the composition measurements of HiRes [5] and the first results of the Telescope Array [6], which are consistent with a light composition around the ankle and above. On the other hand, recent results from the Pierre Auger Observatory [7,8] indicate a composition becoming heavier with increasing energy above the ankle, and the Yakutsk EAS array muon data suggests a non negligible fraction of heavy nuclei above ≃ 10 19 eV [9]. Moreover, the measurements of the KASCADE-Grande [10] collaboration are consistent with a dominantly heavy composition up to 10 18 eV. However, the KASCADE-Grande data indicate a fast rising proton contribution above 10 17 eV.
Thus at present the experimental data on the CR composition do not allow us yet to determine the transition energy between Galactic and extragalactic CRs. In this paper we suggest to use instead experimental limits on the anisotropy of the arrival directions of UHECRs to constrain the maximal contribution of Galactic CRs at E 10 18 eV. At energies below 10 17 eV, the diffusive propagation of Galactic cosmic rays and their resulting anisotropy at Earth was studied in details in Refs. [11,12].
Since the propagation of CRs in the Galactic magnetic field (GMF) is not longer diffusive at E 10 17 eV, we directly propagate UHECRs in the GMF using the numerical code developed in Refs. [13,14]. We present also a way to generate the turbulent field on a nested grid without limitation on its spatial resolution. This method allows us to include magnetic field fluctuations spanning the required large dynamical range of scales, from negligible compared to the CR Larmor radii up to 300 pc. As main result of this work we show that the existing limits on CR anisotropies strongly restrict the contribution of the CNO element group to the Galactic CR component above E 1 EeV, while the contribution of iron is restricted above E 3 EeV.
Details of the method to generate turbulent magnetic fields are discussed in the Section 2. In Section 3, we review the GMF models used and discuss how the CR anisotropy is calculated. Results of numerical simulations are presented in the Sections 4 and 5 for anisotropies and the spectrum of UHECR.

Modeling Turbulent Magnetic Fields
We adopt in this section a convenient way to generate turbulent magnetic fields on nested grids which allows to include a large dynamic range of spatial scales contributing to the turbulence.
A turbulent magnetic field B satisfies B(r) = 0 and B(r) 2 ≡ B 2 rms > 0. Let us denote k the modulus of wave vectors and α the spectral index of the field: α = 5/3, 3/2 and 1 respectively for Kolmogorov, Kraichnan and Bohm spectra. The power spectrum of the field satisfies P(k) ∝ k −α , and the amplitudes of its Fourier modes are |B(k)| 2 ∝ k −α−2 . The spectral index of the turbulent Galactic magnetic field is poorly constrained. While α ≃ 1 appears hardly plausible, both Kolmogorov and Kraichnan spectra could be allowed by the data. To study the dependence of our results on the spectral index, we present below computations for α = 5/3 and 3/2, as examples. Wave vector moduli satisfy 2π/L max ≤ k = |k| ≤ 2π/L min , where L min and L max are respectively the minimal and the maximal variation scales in the turbulent field. In practice, L min corresponds to the damping scale of the field, which could be as low as an astronomical unit. We choose here L min = 1 AU. For α = 5/3 and 3/2, the value of L min does not noticeably affect the results, because the larger α is, the more the energy is concentrated in the modes with large spatial variations. We take L max = 100 − 300 pc. The correlation length L c of the field, defined as in [15], is equal to As discussed in Refs. [14,16], there are two main numerical methods to generate turbulent magnetic fields. First, they can be generated as a superposition of plane waves as in Ref. [17] and computed in any point of the space. Second, values of the field can be pre-computed with the Fast Fourier Transform on a three dimensional cubic grid, which is periodically repeated in space. The value of the field can be extrapolated to any position from these values. Computing the individual trajectories of millions of cosmic rays with rigidities as low as E/Z ∼ 3 × 10 16 eV is achievable within reasonable computing times only with the second method. The number of vertices on such cubic grids is N 3 . N ∼ 256 − 512 is typically the limit above which the grid cannot be loaded in a 2 gigabyte RAM memory. The ratio L max /L min is limited by N /2, when L max equals the size of the cubic box. Moreover, we take L max /L min to be smaller than N /2, by at least a factor of a few. This ensures that the modes with the largest spatial variations ≃ L max have a few oscillations within the box size. Otherwise, the generated turbulent field can be highly anisotropic. Cosmic rays which diffuse in turbulent magnetic fields are mostly sensitive to modes with wave numbers k close to ∼ 2π/r L , where r L is their Larmor radius. For E/Z = 10 18 /26 eV and a field of strength 6 µG, it is equal to r L ≃ 7 pc. In the numerical simulations, one can disregard modes with 2π/k ≪ r L because they have a negligible influence on the particle trajectories. On the contrary, modes with 2π/k ∈ [r L , L max ] which isotropize cosmic rays in a non trivial way have to be taken into account. Therefore, instead of using 2π/k ∈ [L min , L max ] = [1 AU, 100−300 pc], we truncate the minimal scale of spatial variations for the generated field and restrict ourselves to [L ′ min , L max ] with L ′ min sufficiently small compared to r L . However, L max /L ′ min is still too large to fit in one magnetic field grid of reasonable size. To solve this issue, we use the method of nested grids explained in the following.
Let us assume that B(r) is the sum of N + 1 components: In practice, N = 2 is sufficient for this work. B 0 (r) contains all Fourier modes with 2π/k ∈ [L min , L ′ min ], and the fields B i (r) (1 ≤ i ≤ N ) respectively contain the modes with 2π/k ∈ [L i , L i+1 ], where L 1 = L ′ min and L N +1 = L max . The ratios L i /L i+1 are all chosen to be smaller than N /2 by a factor of a few.
The root mean square (rms) strength of the total field, B rms , satisfies [18] which yields the rms amplitude of B i (r), B rms,i . The turbulent field B turb generated for the computations is equal to the sum of the N components . B turb is equal to the total turbulent field B(r) after subtracting the modes with spatial variation scales smaller than L ′ min ≪ r L . Each B i (r) is generated on a cubic grid of lateral size N L i /2. Each grid is periodically repeated in physical space. The B i (r) with large i contain the modes with large spatial variation scales and the B i (r) with small i, the modes with small variation scales. In any space point, the magnetic field from the large (respectively small) resolution grid is evaluated as the 8-point linear interpolation of the values on vertices of the large (respectively small) scale resolution grid.
We have verified that we recover with this code the results found by the earlier studies of Refs. [16,20]. As an example, we present in the appendix our computations of the CR diffusion coefficient for pure magnetic turbulence, as well as the parallel and perpendicular diffusion coefficients for turbulence superimposed to a regular field. Our results are found to be in very good agreement with those of these previous studies.
To summarize, we take in the following L min = 1 AU for the normalisation of the turbulent field strength, so that the rms for the total field with modes satisfying 2π/k ∈ [L min = 1 AU, L max ] would be B rms . In all computations for CRs with rigidities E/Z < 1 EeV, we take N = 2 and set L ′ min , the actual minimal scale of fluctuations in the generated field, to 1 pc. The intermediate scale L 2 between the N = 2 grids is 20 pc. In practice, for large rigidities E/Z ≥ 1 EeV, r L ≥ 180 pc and we can drop the smaller scale grid and only use the larger one : N = 1 and L ′ min = 20 pc. We use the standard Runge Kutta method with adaptative step size of Ref. [19]. Removing the smaller scale grid increases the step size of the integrator and allows us to reduce computing time for particles with E/Z ≥ 1 EeV.

Galactic Magnetic Field models and Method to compute the Anisotropy
The Galactic magnetic field (GMF) can be regarded as the sum of a regular component (large scale variations) and a turbulent component (small scale variations).
We described in the previous Section a method to generate numerically the turbulent component. The spatial profile of the rms strength of the turbulent field, B rms (r, z), is poorly constrained. Therefore, we use two different types of profiles as examples. First, we take a model with an exponentially decaying field strength in the Galactic halo [21]. We will refer to it as the "Profile 1": where r is the Galactocentric radius and z the distance to the Galactic plane. The parameter z 0 denotes the scale height of the random field into the z−direction. We will take z 0 = where B 0 is defined as the value of B rms close to the Sun. Second, we also consider a constant rms strength within a box of size r ≤ 20 kpc and |z| ≤ z 0 ("Profile 2"): Although this profile is very likely less realistic than the previous one, we test it because it corresponds to the profile used in the usual "leaky-box approximation". The global geometry of the regular GMF is still poorly known. The Faraday rotation measures (RM) for extragalactic sources suggest that it is made of at least two different components, in the disk and in the halo, with different geometries [22,23]. The field in the disk is believed to be symmetric with respect to the Galactic plane, while the field in the halo is believed to be antisymmetric [22,23]. The RM at high latitudes show the existence of a toroidal field in the halo, on each side of the Galactic plane. This field is counter clockwise in the Northern halo and clockwise in the Southern halo, as seen from the Galactic North pole. Several analytical models have been proposed to describe the regular GMF. As shown in Refs. [22,24], presently no theoretical GMF model can fit all experimental data. However, Ref. [23] presents the two first models that fit reasonably well all extragalactic RM data in most regions of the sky, which represents a significant improvement of our knowledge of the GMF. Since it is impossible to distinguish between bisymmetric (BSS) and antisymmetric (ASS) geometries for the disk field, the authors of Ref. [23] propose two different benchmark models for the regular GMF. Below, we will refer to them as the "PTKN-BSS" and "PTKN-ASS" models. They contain disk and toroidal contributions. Let us use Galactocentric cylindrical coordinates (r, θ, z), where r = x 2 + y 2 1/2 , and Cartesian coordinates x, y and z. The Earth is assumed to be at (x = 0, y = r ⊙ = 8.5 kpc, z = 0), where θ is set to zero at the position of the Earth, and increases clockwise, as seen from the Galactic North pole. The components in cylindrical coordinates of the disk field strength, B r and B θ , are defined as with p = −5 • and p = −6 • respectively for the ASS and BSS models. For the ASS disk field, while, for the BSS disk field, where with B T0 , r T0 , h T chosen as in Table 1. For |z| ≤ h T , w T = 0.25 kpc, and for |z| > h T , w T = 0.4 kpc. The strength of the halo field decays towards the Galactic center, for r < r T0 . For most of the following computations, we use the PTKN-BSS model as an example. We test the dependence of our results on the regular GMF by also using the PTKN-ASS model, the "ASS+RING" model of Ref. [25] (which we will refer to as "Sun08" in the following), and the Prouza and Smida (PS) model [26,27] with the parameters given in Ref. [13].
The Galactic center [28][29][30], some types of supernovae [31], magnetars [32][33][34][35] or GRBs [36][37][38][39][40][41][42] have been discussed as potential Galactic sources able to accelerate CRs to ultra-high energies. The spatial extension of the region containing Galactic CR sources is better constrained than the GMF parameters. Sources are expected to be distributed in the Galactic disk, within ≃ ±(200-500) pc from the Galactic plane z = 0 [43,44]. The Galactocentric 6 kpc 6 kpc 6 kpc 5 kpc h T 1.3 kpc 1.3 kpc 1.5 kpc 1.5 kpc Table 1. Values for the Northern (z > 0) and Southern (z < 0) halo parameters B T0 , r T0 and h T , in the ASS and BSS versions of the PTKN model. radius r up to which the source region extends is less constrained. Since there should not be a significant number of sources with r > 20 kpc [43], we will take in the following, for most cases, r = 20 kpc as the limit of the source region. At sufficiently low rigidities, the Larmor radius of cosmic rays is smaller than the coherence length of the turbulent GMF. Previous studies that predicted the amplitude of the cosmic ray anisotropy at Earth assumed CRs are diffusive. While the diffusion approximation is justified for rigidities smaller than E/Z ∼ 10 17 eV, it starts to fail in the rigidity range investigated in this work: E/Z ≥ (10 18 /26) eV. Between these rigidities, one typically expects a transition from the diffusive regime to the ballistic regime for CR propagation. The transition does not happen abruptly at a given rigidity, which leads to non-trivial modes of CR propagation. This can have a non-trivial impact on the anisotropy of Galactic CRs at Earth. Therefore, we propagate in this work individual cosmic rays in models of the GMF. Figure 1 shows trajectories of anti-iron nuclei with energies 10 18 eV ≤ E ≤ 10 20 eV backtraced in one GMF model from the Earth, located at (x = 0, y = 8.5 kpc, z = 0). It shows the variety of CR propagation types in the transition from "purely" diffusive (here at 10 18 eV) to "purely"ballistic (here above > ∼ (2 − 4) × 10 19 eV). This energy range is shifted when the magnetic field parameters are changed. The regular GMF used for these plots is the PTKN-BSS model. The turbulent component has a Kolmogorov spectrum with L min = 1 AU, L max = 200 pc, the profile 1 with z 0 = 2 kpc, and a strength set to B rms = 4 µG.
The left panel of Figure 1 display the trajectories of 1, 4 and 8 ×10 18 eV iron antinuclei. Values of spatial coordinates on the axes are given in kilo-parsecs. For these GMF parameters, the Larmor radius of the 10 18 eV nuclei is smaller than the correlation length L c ≃ 40 pc of the turbulent component. The trajectory of this cosmic ray resembles a random walk, see the red line. The green (4×10 18 eV) and blue (8×10 18 eV) trajectories, respectively, correspond to diffusion in the regimes when r L ≃ L c and r L > L c . The trajectory of the 8 × 10 18 eV iron anti-nucleus is still confined in the Galactic plane for an extended time. On the right panel of Figure 1, one can see that this anti-nucleus goes back and forth in the disk. It propagates especially along the regular field lines which are locally approximately oriented along the x axis.
The right panel of Figure 1 shows the trajectories of 1, 4 and 10×10 19 eV iron anti-nuclei. At 10 19 eV (red line), the CR is still strongly deflected before escaping the Galaxy. If one sums up all deflections along its trajectory, it exceeds 360 • . This iron anti-nucleus is weakly deflected over distances up to ∼ 1 kpc. It is strongly deflected only locally, when it reaches regions with stronger turbulent magnetic field fluctuations. At 10 20 eV, the trajectories are fully ballistic, see the blue line. At such energies, one expects that iron nuclei suffer deflections of the order of ∼ 20 • − 40 • before escaping the Galaxy. The 4 × 10 19 eV anti-nuclei are not diffusive any more. However, they still experience large deflections. For instance, there is a big wiggle on the 40 EeV particle trajectory, see the green line in the right panel.
In principle, the best way to compute the anisotropy of Galactic CR at Earth is to use forward tracking. One should inject cosmic rays in the Galaxy at the source locations and only record the momenta of CR which cross a sphere around the Earth. The radius of this sphere should be small compared to the CR Larmor radius. This is, however, not feasible within reasonable computing times for the lowest rigidities we study. Therefore, we use a method first proposed in Ref. [45], and reused in more recent works such as Ref. [46]. It consists in backtracking anti-particles with random initial momenta from the Earth to outside the Galaxy, and to record for each one the total path length in the source region. This corresponds to assuming a continuous and homogeneous source distribution inside the source region. Except for the Galactic center, potential Galactic CR sources should be transient. The method we use here assumes the existence of sufficiently frequent transient sources in the Galactic plane, so that the continuous source distribution hypothesis is fulfilled. For rare transient sources whose periodicity start to be comparable to the CR confinement time in the Galaxy, the current anisotropy may strongly differ from the average anisotropy [47]. Ref. [47] shows that if GRBs were to be the sources of Galactic CRs up to the ankle, strong variations of the Galactic CR flux and anisotropy should be expected on time scales of a hundred Myr. This would mean that CR anisotropy limits may be compatible with Pierre Auger measurements if we live in atypical times, when it drops below ≃ 2% at EeV energies [48]. Since one can hardly go beyond this statement for rare transients such as GRBs, we restrict our work to the case of sufficiently frequent transients. We will discuss below in more detail the domain of validity of this approximation, see Section 5.
We count the length of particle trajectories contained in the source region without any weighting depending on the position inside the region. This corresponds to assuming that the sources are homogeneously distributed inside −200 pc ≤ z ≤ 200 pc and r ≤ 20 kpc. In practice, taking a more realistic source distribution with, for example, a source density decreasing with |z|, would only increase the Galactic CR anisotropy at Earth. The values presented below can be regarded as lower limits.
Finally, we compute the amplitude of the dipole of CR anisotropies and compare it to the upper limits presented by the Pierre Auger Collaboration, for energies E ≥ 1 EeV [48]. To do so, we associate to each of the N backtracked cosmic rays a vector v(θ, φ) whose direction on the sky (θ, φ) corresponds to the initial CR direction at Earth. Its length |v| corresponds to the trajectory path length in the source region. The dipole direction is given by the sum of all vectors i v i . In spherical coordinates (r, θ, φ), v = L (1 + D cos θ) u r for a dipole of amplitude D and direction u z , with L = 1 N i |v i | being the average length of vectors and θ the angle between u r and u z . Since the vectors v are isotropically distributed, If higher order multipoles are present, Eq. (3.9) is unchanged because only the dipole has a non-zero contribution to the sum i v i .

Anisotropy of Galactic Cosmic Rays predicted at Earth
In this Section, we compute the anisotropy at Earth of cosmic rays emitted by sources distributed in the Galactic plane, and compare it to the upper limits on the anisotropy as measured by the Pierre Auger Collaboration [48]. Within a given Galactic magnetic field model, if the predicted anisotropy of Galactic cosmic rays exceeds these limits at a given energy, a sufficiently large contribution of extragalactic CRs is required above that energy. Extragalactic CRs of energies below the GZK threshold around ≃ 4 × 10 19 eV [49,50] tend to have anisotropies at the percent level or below because either a large number of sources at cosmological distances can contribute to the flux or, in case of relatively strong deflections in extragalactic magnetic fields, such deflections tend to wash out anisotropies over the long path lengths propagated over the typical energy loss distance. Resulting anisotropies of extragalactic cosmic rays can, therefore, easily be below the current upper limits. If sources at cosmological distances dominate the extragalactic CR flux, the motion of the Sun with respect to the CMB frame would induce an anisotropy of the extragalactic flux of ≃ 0.6% [51], due to the Compton-Getting effect [52]. In subsection 4.1, we present our predictions for the Galactic CR dipolar anisotropy and discuss its implications on the energy at which the transition from Galactic to extragalactic CR should occur. Figs. 2-5 show computations for this anisotropy. Assuming in a first approximation that the extragalactic flux is isotropic, one can deduce from these figures the maximum contribution of Galactic CR sources to the total flux. For instance, at the transition from Galactic to extragalatic CRs, half of the CR flux is of Galactic origin and half of extragalactic origin. Therefore, the exact transition energy must be in the energy range where half of the dipole amplitude of Figs. 2-5 does not overshoot the experimental upper limits on it. In 4.2, we study in more detail the anisotropy predicted by Galactic protons below the ankle and its dependence on the GMF parameters.

Dipole Amplitude and predicted Energy of the Transition
We assume that the sources of Galactic CRs are located in a cylinder along the Galactic disc with height z max and radius r max . Varying the extension of the source region from r max = 15 kpc to r max = 20 kpc, we verified that for rigidities E/Z 3 × 10 18 eV, the dipole amplitude does not change by more than 3%. For higher rigidities, the difference rarely exceed (5 − 10)%. In practice, the conclusions below will not change for r ≤ 15 kpc or r ≤ 20 kpc. In the following we assume r ≤ 20 kpc. Larger extensions of the source region along the z−direction would reduce the predicted anisotropy. We verified that, as long as the source region is less extended than |z| 500 pc, no strong modification of the results and conclusions would arise. We assume |z| = 200 pc for most of the figures below.
For each of the following plots we backtrack 10 4 particles. We find that this induces an error ≤ ±3% on the predictions of the dipole amplitude. Reducing further this error to ≤ ±1% could in principle be achieved by backtracking ≃ 10 times more CRs, but it is in practice impossible due to computing time reasons. Figure 2 presents the simulated predictions for the dipolar amplitude at Earth for several different Galactic CR primaries: protons, helium, carbon (representative for the CNO group), silicon and iron. CR sources are assumed to be located in the thin disk with −200 pc ≤ z ≤ 200 pc and r ≤ 20 kpc. The PTKN-BSS model is assumed for the regular GMF, and the turbulent field parameters are taken to be B 0 = 4 µG, z 0 = 2 kpc, L min = 1 AU and L max = 200 pc. In the left panel, we use profile 1 given in Eq. (3.1) for the turbulent field, whereas in the right panel we take profile 2 from Eq. (3.3). We include in these figures measurements of, and 99% C.L. upper limits on, the dipole amplitude in right ascension from the Pierre Auger Collaboration, as indicated in the captions. We assume that the true dipole vector does not lie in the equatorial plane, so that these upper limits in right ascension do not overconstrain significantly the total amplitude. In Fig. 2, the predictions for the dipole amplitude for silicon and iron primaries may appear to exceed the Auger upper limits in the whole energy range we consider. However, since error bars on our computations are ≃ ±3%, some parts of these silicon and iron lines are compatible with the Auger limits at low energies. For instance, one cannot exclude the dipole amplitude for 1-few EeV Galactic iron to be below the Auger limits. In the next two figures, 3 and 4, the red dashed lines represent the lower error bars on the red shaded bands. These bands represent the allowed ranges for the dipole amplitude with Galactic iron primaries, for different turbulent GMF parameters. In the energy ranges where this dashed line falls below Auger upper limits, a pure flux of Galactic iron CRs cannot be excluded on anisotropy grounds. In practice, we expect that with ten times larger statistics, predicted dipole amplitudes at the lowest rigidities (E/Z < ∼ few EeV/26), which are already smaller than a few percent, would fall further below the Auger limits. We have 10 sets of 10 4 CRs for different turbulent GMF parameters. Their individual dipole directions look random at the lowest rigidities and if one adds up these ten sets, the resulting anisotropy for iron at E = 10 18 eV falls to ≈ 0.6%, which is below the Auger limit. This gives an idea of what should be expected for ten times larger statistics. Therefore, below the ankle, both silicon and iron of Galactic origin are compatible with Pierre Auger Observatory limits. For the set of GMF parameters assumed here, in case of a predominantly heavy composition below the ankle and sufficiently frequent transient sources, CRs may still be of Galactic origin up to the ankle. Iron nuclei of Galactic origin up to ≃ 10 EeV cannot currently be ruled out from the point of view of the CR anisotropy.
Depending on the composition at E ≃ 10 18 eV, this has an important implication for the transition energy between Galactic and extragalactic CRs: For the set of GMF parameters assumed here, if the CR primary composition is predominantly light (p, He) or intermediate (C, N, O) at these energies, the predicted anisotropy at Earth would be larger than the 99% C.L. upper limits from the Pierre Auger experiment if these nuclei were of Galactic origin, as seen in Fig. 2. This implies that if the composition at E ≃ 10 18 eV is measured to be light or intermediate, scenarios in which the transition from Galactic to extragalactic CRs occurs at the ankle are strongly disfavoured, at least for a wide range of GMF parameters. We investigate below the ranges of parameters for which this conclusion would be valid. Figure 2 also shows that the conclusions above do not strongly depend on the turbulent field profile. For profile 2 (right panel), the predicted dipolar anisotropy grows slightly more slowly with energy than for profile 1 (left panel). This is expected because for profile 1 the gradient of the turbulent field tends to drive CRs towards larger z in the Galactic halo slightly faster than for the constant field of profile 2. Since predicted anisotropies are not very different for the two profiles, we will mostly focus on profile 1 in the following.
We have also tested the dependence of these results on the regular GMF model. For rigidities E/Z 3 EeV, the dipole amplitude and direction depends on the regular GMF model. However, the change of the dipole amplitude is too small to affect significantly our findings. Moreover, the PTKN-BSS model which we use in all the following figures is one of the models with the lowest dipole amplitudes among those tested. Figures 3 and 4 present the dependence of the previous results on the other turbulent GMF parameters, which are poorly constrained: the index α of the fluctuation spectrum ( Fig. 3 -left panel), the maximal length of field fluctuations L max (Fig. 3, right panel), the field strength normalization B 0 (Fig. 4, left panel) and the scale height z 0 (Fig. 4 -right  panel). The shaded areas of the filled cures represent the relative change of the results for p, C and Fe primaries when varying separately the four above parameters and keeping all other parameters at the values in Fig. 2. The red dashed line is computed as the lower boundary of shaded areas for iron primaries minus 3% from our statistical uncertainties.
At low rigidities E/Z 4 EeV, the dipole amplitude grows with E/Z as expected. At larger rigidities, CRs start to enter the ballistic regime and higher order multipoles start to make a significant contribution to the total Galactic CR anisotropy at Earth. As seen in Figs. 3 and 4, the dipole amplitude may then become smaller and/or vary with E/Z. Dipole amplitude (percent) Energy (x 10 18 eV) Fe C p Figure 3. Dependence of the predicted dipole amplitude on the turbulent field spectral index and on its maximum spatial variation scale. For comparison, Pierre Auger data [48] are shown in blue as in   Thus a decrease of the dipole amplitude at high energies does not necessarily imply that the distribution of CR arrival directions at Earth becomes more isotropic. The widths of the filled curves indicate that results mostly vary with the turbulent GMF strength and its maximum spatial variation scales. Results are less sensitive to the spectral index of the field and no strong difference in the dipole amplitude is found between the fields with Kolmogorov and Kraichnan spectra. The amplitudes at E/Z 3 EeV/26 are slightly larger for the Kraichnan spectrum because less power is present in the large length scale modes, which are relevant at such rigidities, than for the Kolmogorov spectrum. Results are only marginally sensitive to the extension z 0 of the turbulent field into the halo, see right panel of Fig. 4. This is due to the fact that CRs which escape the source region and propagate to large z in the halo rarely come back to the source region.
Results are mostly sensitive to L max and B 0 . The maximal length scale of the turbulence L max determines up to which energy CRs still scatter on the turbulent field inhomogeneities. For larger L max , CRs can be diffusive up to larger energies, which therefore reduces their expected anisotropy at Earth. One can see in the right panel of Figure 3 that for L max = 300 pc and B 0 = 4 µG, the dipole amplitude below E ≃ 15 EeV for a pure iron composition may be compatible with the current 99% C.L. upper limits from the Pierre Auger Observatory.
The dependence of our results on the turbulent field strength B 0 is very strong, see the left panel of Fig. 4. The upper parts of the shaded areas correspond to B 0 = 2 µG and the lower to B 0 = 8 µG. For B 0 = 8 µG, the anisotropy at Earth of iron primaries is a priori compatible with the Pierre Auger upper limits up to E ≃ 20 EeV. For B 0 = 2 µG, the dipole amplitude starts to overshoot the Pierre Auger limits around E ≃ 3 × 10 18 eV, while for B 0 = 4 µG, the amplitude starts to exceed the Pierre Auger limits around 10 EeV. For all cases, a light or intermediate composition at E ≃ 10 18 eV would exceed the Pierre Auger upper limits and rule out the ankle as the transition from Galactic to extragalactic cosmic rays.

Dipole Amplitude at E/Z = 10 18 eV
In this section we demonstrate that for any reasonable combination of turbulent GMF parameters, at E ∼ 10 18 eV the dipolar anisotropy predicted by light primaries of Galactic origin is always larger than the observational limit from the Pierre Auger experiment. Therefore, having a reliable composition measurements at such energies is crucial for knowing if the ankle can or cannot be the signature of the transition from Galactic to extragalactic CRs. Low energy extensions such as HEAT [53] and AMIGA [54] can solve this important question. when the scale height z 0 is allowed to vary in the range 1 to 8 kpc (shaded areas). Both turbulent field profiles are tested, and a ±500 pc width source region is also tested. Both for Kolmogorov (left panel) and Kraichnan spectra (right panel), predicted dipole amplitudes are above 10%, considerably higher than the ∼ 2% upper limit from the Pierre Auger experiment at such energies. Therefore, the ankle cannot be the signature of the transition from Galactic to extragalactic CRs, if the contribution of Galactic protons to the CR flux at 1 EeV is larger than ≃ 20%. Such a scenario would be consistent with the dip model [1]. A CR flux at 1 EeV consisting mostly of light nuclei cannot be of Galactic origin, except in the very unlikely case of B 0 ≫ 10 µG. Figure 6 presents histograms of the relative fraction of CRs spending a certain time in the source region, for rigidities ranging from 1 EeV/26 to 17.8 EeV/26. We use the PTKN-BSS model for the regular GMF component. The turbulent GMF strength is set to B 0 = 4 µG, its extension in the halo to z 0 = 2 kpc. We take a Kolmogorov spectrum with maximum spatial variation scale L max = 200 pc. With such parameters, the average time spent in the source region for cosmic rays with E/Z = 1 EeV/26 is ≃ 0.5 Myr. The CR escape times from the magnetized region of the Galaxy defined as −10 kpc z 10 kpc and r 20 kpc, are found to be ≃ 5 times larger than the times spent in the source region for the turbulent GMF profile 1. CRs which escape the Galactic thin disk containing CR sources still stay a non-negligible amount of time in the halo compared to the time spent in the source region.

Energy Spectrum of Galactic Cosmic Rays and Sources contributing at Earth
For 1 EeV iron nuclei (red curve in Fig. 6 -left panel), only a few percent of CRs stay more than 1 Myr in the source region. This implies that rare transient sources such as gamma ray bursts (GRBs) are very unlikely to be sources of Galactic CRs in the sub-ankle region, even if heavy nuclei were able to escape such sources. For 10 EeV iron nuclei, the average time spent in the source region is ≃ 0.06 Myr, which is nearly ten times smaller than at 1 EeV. We found that the average time spent in the source region is approximately proportional to 1/E -or slightly softer-in the rigidity range E/Z ∈ [1 EeV/26, 20 EeV/26]. This implies that for a source injection spectrum proportional to E α , the spectrum at Earth approximately goes as E α−1 . Interestingly, this result is compatible with predictions from diffusion in the regime where the Larmor radius is larger than the coherence length of the turbulent field: The distance traveled by CRs within a given amount of time is proportional to the square root of the diffusion coefficient, and the diffusion coefficient is proportional to E 2 when Larmor radii are larger than ≈ L c [20]. We found no significant change to this conclusion when varying turbulent GMF parameters in the ranges tested in the previous section. However, the time spent in the source region for B 0 = 8 µG is ≃ 30% smaller than for B 0 = 4 µG at E/Z ∼ 10 18 eV. For stronger turbulent fields, the larger turbulent field gradient towards z tends to make CRs leave the source region faster.
We show in Fig. 7, which regions of the Galactic plane are passed most by CRs back-traced from the Earth for E/Z = 10 18 eV/26 (upper left panel), 3 × 10 18 eV/26 (upper right), 10 19 eV/26 (lower left) and 3 × 10 19 eV/26 (lower right). This equivalently shows which parts of the source region contribute most to the Galactic CR flux detected at Earth. The color code presents the fraction of all particles backtraced from the Earth which pass cubes of 200 pc lateral size located in the Galactic disk. Here, CRs are not counted more than once. Location of Earth is marked by the bright spot at (x, y) = (0, 8.5 kpc). One can clearly see the shapes of the spiral arms present in the model of the regular GMF model. CRs indeed diffuse or propagate faster along the regular field direction. The Galactic center region appears strongly demagnified, which makes it unlikely to significantly contribute to the observed fluxes. This is due to the stronger field in the disk and especially towards the Galactic bulge. For the forward tracking point of view, CRs potentially emitted by the Galactic center would escape from the Galactic thin disk which contains the Earth, and propagate towards larger z, before reaching Earth. At E/Z = 3 × 10 19 eV/26 (lower right panel), the CRs start to be in the ballistic regime and only the region within a few hundreds of parsecs from Earth could contain sources. At these rigidities, CRs cannot have Galactic origin any more because anisotropies would exceed the Pierre Auger upper limits on the dipole anisotropy, as discussed in the previous section. When the confinement time of CRs in the source region (≃ 0.5 Myr for 1 EeV iron nuclei) starts to be comparable to the period between two potential Galactic UHECR sources, the continuous source distribution approximation used in this paper breaks down. For GRBs, this happens at E/Z (0.1−1) EeV/26. In this case, the expected Galactic CR flux at Earth and its anisotropy would strongly vary on time scales of several Myr, see for instance Ref. [47]. The CR anisotropy at Earth may then be substantially smaller or larger than those computed for a continuous source distribution. Potential sources such as magnetars are expected to have a larger rate of ∼ 1 per 1000 years. For such rates, ∼ (0.5 Myr)/(1000 yr) ∼ several hundreds of sources would contribute to the flux observed at E/Z = 10 18 eV/26. This larger number of sources substantially reduces the fluctuations in time of the Galactic CR flux and anisotropy at Earth. Therefore, in this case, the anisotropy does not significantly differ from the values presented in this work, and the continuous source distribution approximation is valid.
Let us now estimate more quantitatively when the continuous source distribution approximation breaks down. The average number of sources that would contribute to the Galactic CR flux at one given rigidity can be estimated as the average time spent by CRs in the source region multiplied by the source rate R. We plot in Fig. 8 (left panel) this estimate of the number of contributing sources, for three different rates R = 10, 1, 0.01 kyr −1 and for CR rigidities in the range E/Z = (1 − 32) EeV/26. We assume here for the turbulent GMF a Kolmogorov spectrum with L max = 200 pc, strength B 0 = 4 µG and scale height z 0 = 2 kpc. For sources with a rate R = 1 kyr −1 (green line) similar to that expected for magnetars, the average number of contributing sources stays above 100, up to the ankle for iron nuclei. It decreases from ∼ 500 for 1 EeV iron to only ∼ 10 at 32 EeV. As shown below, ∼ 10 sources is too small for the continuous source distribution approximation to be valid. With a ten times larger source rate R = 10 kyr −1 (red line), there would be 100 sources contributing for all the explored rigidity range, but 10 kyr −1 is of the order of the Galactic supernovae rate which looks very unlikely for extreme CR accelerators. The blue line corresponds to rarer transients R = 0.01 kyr −1 . At the ankle, only one source would contribute in average. Sources such as Galactic GRBs with R ∼ 1 Myr −1 cannot be described by the continuous source distribution approximation. However, such sources are unlikely to be responsible for the sub-ankle CR flux if it were to be of Galactic origin. Indeed, one can for example hardly match the bumpy CR spectrum resulting from rare Galactic transients to the observed smooth power law spectrum [47]. We provide in Fig. 8 (right panel) an estimate of the minimum rate R of Galactic sources that would be required to maintain relative fluctuations of the flux σ(F )/ F below five given thresholds (≈ 5, 10, 25, 50, 100%) on ≫ R −1 time scales. Since we follow individual CR trajectories, we cannot directly provide σ(F )/ F for computing time reasons. We can however estimate its value: The panels of Fig. 7 also give an estimate of the total flux that would be received in any point of the Galactic disk from one source located at (x, y) = (0, 8.5 kpc). Looking for the total flux received at the Earth position from N sources -with same power-located in N random positions in the Galactic disk is then roughly equivalent to putting N observers in N random locations in the disk and summing up the total flux they receive from the single source located at (x, y) = (0, 8.5 kpc). We compute the flux F in the latter way for 10 4 different configurations, using the computations of Fig. 8 (left panel) for N . This yields the estimate for σ(F )/ F that is used in Fig. 8 (right panel). To maintain σ(F )/ F below ≈ 5% (red solid line) at the ankle for Galactic iron primaries, sources with rates comparable with that of Galactic supernovae R ≈ 30 kyr −1 would be needed. For rates R ∼ 1 kyr −1 , σ(F )/ F ≈ 10% for 1 EeV iron nuclei and remains 25% below the ankle. In this case, the continuous source distribution approximation is still valid. For energies E (10 − 20) EeV, it quickly starts to break down: For 20 EeV (resp. 30 EeV) iron nuclei, σ(F )/ F ≈ 50% (resp. ≈ 100%). For Galactic source rates R 0.01 kyr −1 , σ(F )/ F exceeds 100% above ∼ 1 EeV and the anisotropy measured at Earth is expected to significantly differ from the averaged values presented in the previous section.

Conclusions and Perspectives
In this work we studied the consistency of a transition from Galactic to extragalactic CRs with existing anisotropy limits as a function of energy above E = 10 18 eV. The diffusion approx-imation predicts a dipole anisotropy δ = −3D ij ∇ j ln(n) increasing with energy, since both the diffusion tensor D ij and the relative CR gradient ∇ j ln(n) increase with energy. However, this approximation becomes unreliable at O(E/Z) ∼ 10 16 eV, and therefore we studied the propagation of CRs in the Galactic magnetic field directly by backtracking trajectories. We simulated the turbulent magnetic field on nested grids which allows one to include turbulent field modes B(k) with arbitrary small wave-lengths. For the regular Galactic magnetic field we used up-to-date models from Ref. [23]. Because the global structure of the GMF is still rather uncertain, we studied the dependence of the resulting anisotropy on the magnetic field parameters such as its strength B 0 , scale height z 0 , correlation length L c and exponent α of its power-spectrum. We also examined the dependence of our results on the width and height of the disk in which sources are located.
The main results of this study are presented in the Figs. 3-5. They show that the anisotropy mostly depends on the amplitude B 0 of the magnetic field in the disk. As our main conclusion from this study, we found that existing anisotropy limits are not compatible with light (proton) and intermediate (CNO) nuclei of Galactic origin as dominant contribution to the CR flux above 1 EeV. By contrast, Galactic iron nuclei as CR primaries are consistent with the existing limits even up to 10-20 EeV, if the strength of the turbulent field is as large as B rms ∼ 8 µG. This finding implies that determining the chemical composition of the CR flux around 10 18 eV settles also the question of the transition energy between Galactic and extragalactic component: As light nuclei at this energy are not sufficiently isotropized, they have to be extragalactic. Therefore the fast increasing proton contribution indicated by the KASCADE-Grande collaboration between 10 17 eV and 10 18 eV suggests the beginning of an extragalactic component.
We also studied qualitatively the dependence of the anisotropy on the effective density of sources, see Figs. 6-8. The average escape time of iron nuclei with 10 EeV energy from the Galaxy is ∼ 10 5 yr. Assuming for magnetars a rate of 10 −3 /yr, the effective density of magnetars as sources of CR at 10 EeV is ∼ 100/Galaxy. Thus magnetars satisfy the anisotropy constraint and can be natural candidates for the sources of the high-energy end of the Galactic CR flux in the scenario where the transition from Galactic to extragalactic cosmic rays occurs at the ankle, provided they are able to accelerate iron up to few ×10 18 eV.
In summary, we conclude that models with a transition from Galactic to extragalactic cosmic rays around the ankle are consistent with the existing anisotropy limits if the composition of Galactic cosmic rays at E 10 18 eV is dominated by heavy nuclei. In contrast, if the chemical composition at these energies turns out to be light or intermediate, a transition at the ankle would be very strongly disfavoured.

A Validitation of the nested grid method
In this appendix, we verify the validity of the new method we propose to generate turbulent magnetic fields on nested grids, see Section 2. We have reproduced the earlier results of Refs. [16,20], and present below our computations for the diffusion coefficient versus the Casse et al. ones [20].
In Fig. 9, we present the numerical results for the parallel (right panel) and perpendicular (left panel) diffusion coefficients for 1 PeV to 600 PeV CR protons diffusing in a field containing both a regular and a turbulent component. Let us denote B reg the strength of the regular component, and B rms the root mean square strength of the turbulent one. Two levels of turbulence are used η = 0.1 and η = 0.46, with η = B 2 rms /(B 2 reg + B 2 rms ). The magnetic field strength is set to 4 µG. For the turbulent component, we take a Kolmogorov spectrum (α = 5/3) with L max = 150 pc and L ′ min = 0.1 pc. Therefore, L ′ min is smaller than the Larmor radius r L for all energies. We average over a few turbulent field configurations and propagate 2000 protons.
Red symbols in Fig. 9 represent our results with the nested grid code introduced in Section 2, and green symbols represent the results of Figs. 4 and 5 of Ref. [20] adapted to the values we use here for L max and the magnetic field strength. The uncertainties of our values and those of Casse et al. [20] can be estimated from the fluctuations from one point to another compared to averaged behaviour of the diffusion coefficients. Both for η = 0.1 and for η = 0.46, our results reproduce very well those of Casse et al.. We did not report in Fig. 9 the results of Casse et al. below E ≃ 4 PeV because they correspond to computations for CRs with Larmor radius smaller than the minimum size of turbulent magnetic field fluctuations in their grid. This stresses one of the advantages of our new method: There is no lower limitation on L ′ min /L max because one can always add another grid with smaller spacing and smaller scales of the magnetic field fluctuations. Therefore one can safely explore low rigidities.
In Fig. 10, we test our code for the case of 100 TeV to 1 EeV CR protons diffusing in a purely turbulent field, i.e. without any regular field, B reg = 0. The parameters are the same as for Fig. 9, except for E = 100 − 300 TeV where we take L ′ min = 0.01 pc to ensure that L ′ min is smaller than r L . For computing time reasons, we keep L ′ min = 0.1 pc for E ≥ 1 PeV. We use B rms = 4 µG and α = 5/3.
Red crosses in Fig. 10 correspond to our computations and green ones to those of Fig. 4 of Casse et al. [20] adapted to our values for L max and B rms . One can see that also in this case the results are in very good agreement. As predicted theoretically, the diffusion coefficient is proportional to E 1/3 (respectively to E 2 ) at low (respectively high) rigidities. The computations of Casse et al. for pure turbulence were done with a field generated superposing Fourier modes. This enabled them to check such a wide rigidity range. However, such computations are significantly slower than ours, as discussed in Section 2.
For the computations at E ≥ 1 PeV in Fig. 10, we take two grids and the intermediate scale L 2 between the two grids is 5 pc, which corresponds here to the Larmor radius of E ≈ 10 PeV protons. We have computed twice more values in the range 3 − 30 PeV than at other energies so as to check that this scale does not induce any artificial imprint in the results. As can be seen in the figure, the diffusion coefficient behaviour in this energy range is smooth and not affected by that intermediate scale.   Figure 10. Computations of the diffusion coefficient for CR protons with energies E = 100 TeV to 1 EeV, in pure isotropic magnetic turbulence (no regular field); α = 5/3, B rms = 4 µG and L max = 150 pc for the turbulent field parameters. Red crosses for our results and green ones for the Casse et al. results [20].