© ESO 2018

1. Introduction

The sources of Galactic cosmic rays (GCRs) remain elusive in spite of decades of intense observational and theoretical efforts. Supernova remnants (SNR; Blandford & Ostriker 1978; Krymsky et al. 1979; Meyer et al. 1997) and superbubbles (Higdon et al. 1998; Binns et al. 2005; Bykov & Fleishman 1992; Parizot et al. 2004) have long been acknowledged as promising candidates, based on energy considerations, isotopic composition arguments and a detailed understanding of the characteristics of particle acceleration. Several issues remain outstanding, however, including the ²²Ne signature of GCRs and the maximum energy levels that can be accounted for (Lagage & Cesarsky 1983). Furthermore, while there is no doubt that these astrophysical environments do accelerate particles, as shown by the high-energy radiation that they generate (Koyama et al. 1995), many questions remain about the magnitude of their actual contribution to the locally observed GCRs. In addition, new observations of unexpected structures in the low-energy GCR spectrum and composition (Adriani et al. 2011; Aguilar et al. 2015) raise questions about the respective contributions of different sources in different energy ranges.

In this context, growing evidence for episodes of intense dynamical activity in the inner regions of the Galaxy (Acero et al. 2016; Abramowski et al. 2016) justifies an evaluation of their potential contributions to GCRs and implications for the characteristics of high-energy particle acceleration (Cheng et al. 2012; Tibolla & Blandford 2018). Indeed, a total energy release of up to 10⁵⁷ ergs has been proposed (Guo & Mathews 2012), which is enough to compete with the average SNR power in the entire Galaxy if the repetition time is of the order of 10⁷ years. From a study of how so-called Fermi bubbles interact with the Milky Way hot gas halo, Miller & Bregman (2016) have estimated that the average energy injection rate is in a 1–7  ×  10⁴² erg range, which exceeds the kinetic power due to SN explosions in the interstellar medium. These results have motivated us to investigate whether the particles accelerated during these episodes may account for a significant fraction of the GCRs, at Earth and/or elsewhere in the Galaxy, or, conversely, what constraints can be derived about Galactic transport of these particles if their contribution is negligible. In the following, we will refer to these particles as Galactic centre cosmic Rays (GCCRs).

In this paper, we have made a first attempt to address these important questions by studying the contribution of a continuous source of energetic particles at the centre of the Galaxy to the local GCRs. Our calculations rely on a simplified propagation model similar to that which is used in generic studies of GCR phenomenology (Ginzburg 2013; Strong & Moskalenko 1998; Taillet & Maurin 2003; Bringmann & Salati 2007; Boudaud et al. 2015b; Giesen et al. 2015; Genolini et al. 2015). This model includes energy-dependent diffusion and advection in a Galactic wind, energy losses and particle re-acceleration, and is described in Sect. 2. The formalism and resolution scheme are presented in Sect. 3 and results are shown in Sect. 4. A summary and discussion of the results are proposed in Sect. 5.

Fig. 1. Geometry of the diffusive volume.

2. Model description

For the present exploratory study, we used a classical, simplified, two-zone model of the Galaxy, consisting of a cylindrical homogeneous disk with radius R = 20 kpc and half-thickness h = 100 pc, surrounded by a cylindrical magnetic halo with the same axis and radius, R, and with a half-thickness L ≫ h (see e.g. Taillet & Maurin 2003). Cosmic-ray transport is treated through a standard diffusion equation (see e.g. Strong & Moskalenko 1998; Taillet & Maurin 2003), including a convective term corresponding to a wind in the z direction orthogonal to the Galactic disk, dragging cosmic rays away from the disk with a velocity V _wind. Diffusion is assumed to be isotropic and homogeneous over the entire “confinement volume”. Diffusion coefficient D depends on the magnetic rigidity of the particles (see e.g. Bringmann & Salati 2007):

(1)

where β = v/c, R_GV is the rigidity in units of GV, D₀ is a constant expressed in kpc² Myr⁻¹, and δ depends on the type of turbulence underlying the diffusion process. Here, we allow for both Kolmogorov and Kraichnan turbulence spectra, corresponding to δ = 1/3 and δ = 1/2, respectively.

It is convenient to use the kinetic energy of the particles, K, as the primary variable. The particle spectral density at energy K, Ψ(K, r , t) is written as a function of time t and position with respect to the Galactic centre r. Neglecting the spallation of heavier species, it must satisfy the following diffusion equation:

(2)

where Q(K, r, t) is a source term to be specified, Γ_ISM is the rate of “catastrophic losses” of the particles (decay or destruction due to interactions with ambient matter in the interstellar medium), b_loss(K) is the rate of energy loss, and χ(K) is an effective diffusion coefficient in energy space, associated with a re-acceleration process accompanying the diffusion of particles in space. In its simplest form, the latter is traditionally modeled as a Fermi second-order acceleration process related to the ambient magnetic turbulence, which is assumed to be isotropic, and can be expressed in terms of a single parameter, namely the Alfvén speed, v_A. In this model, the diffusion in energy space is thus inseparable from that in geometric space, and the coefficient χ(K) can in practice be related to D(K) by (see Maurin et al. 2001):

(3)

where E_tot = K + mc².

Energy loss occurs due to ionization and adiabatic processes, such that b^loss = b^ion + b^adiab. We further consider that the interstellar matter only fills the “infinitely thin” disk (with respect to other dimensions) and use cylindrical co-ordinates centred on the galactic axis. Assuming cylindrical symmetry, the energy loss term in Eq. (2) may be written as

(4)

where r is the galactocentric radius: r = r u_r + z u_z. If the “discofugal” wind reaches its nominal velocity V_wind at the top and bottom of the disk, the overall effect of adiabatic losses can be written as:

(5)

It has been also argued by Ptuskin et al. (1997) that re-acceleration may be mostly efficient within the disk, because it would not work in the adopted wind model. We do not believe that this should be necessarily the case, since some degree of re-acceleration should also accompany diffusion, even in the case of a non-isotropic turbulence. However, this assumption was also adopted by Maurin et al. (2001), Taillet & Maurin (2003), Donato et al. (2004) and Giesen et al. (2015) in their study of GCR propagation, and we shall use it here for the sake of simplicity and to allow a direct comparison of our results. We thus write:

(6)

Likewise, the source term can be written Q(K, r, t)=Q(K, r, z = 0, t)×2h δ(z) if the sources are distributed in the Galactic disk, and

(7)

for a central source, where f(t) is a function of time allowing for time changes of the particle injection rate. For dN/dK, we take a power law in momentum (in some relevant energy range), with logarithmic index α: dN/dp = N₀(p/p₀)^−α. The normalization of this injection term is directly related to the total injection power from that source:

(8)

Boundary conditions follow from the standard assumption that the magnetic halo has finite dimensions, such that diffusion only confines cosmic rays within a limited volume. Outside this “volume”, particles escape freely with a velocity close to c, resulting in a practically vanishing density. We thus impose:

(9)

3. Resolution scheme

3.1. Discrete Fourier-Bessel expansion

With these boundary conditions, Eq. (2) can be solved in cylindrical co-ordinates by expanding the cosmic-ray spectral density as a series of Bessel functions (see e.g. Morse & Feshbach 1953):

(10)

where {u_i}_i ≥ 1 are the zeros, ranked by increasing order, of the zeroth-order Bessel function of the first kind, J₀(x). J₀(u_ir/R) is an eigenfunction of the Laplacian operator appearing in the diffusion equation, with eigenvalue . Thus:

(11)

Applying the Fourier-Bessel transform, , to Eq. (2), one derives equations for each Bessel “coefficient” of order i, P_i(K, z, t):

(12)

where Q_i(K, z, t) is the ith order Bessel coefficient of the source term, which writes, for a source term as in Eq. (7):

(13)

By symmetry, Ψ(K, r, z, t) is expected to be an even function of z, and so will be P_i(K, z, t), ∀i.

3.2. Steady-state solution

In this paper, we concentrate on steady-state solutions, such that ∂P_i/∂t = 0 in Eq. (12). The resolution follows standard procedures, as in Taillet & Maurin (2003) for example. We first solve Eq. (12) in the halo, i.e. outside the disk: z ≠ 0, where both the right-hand side of Eq. (12) and the energy loss/re-acceleration term vanish. Restricting ourselves to z >  0, where V_wind is constant, we rewrite the equation as follows (after dividing by −D(K)):

(14)

Integrating with the boundary conditions P_i(z = L)=0 yields:

(15)

where S_i = [(V_wind/2D(K))² + (u_i/R)²]^1/2.

To obtain the solution in the disk plane, P_i(K, 0), one integrates Eq. (12) through the disk between [ − ϵ, ϵ] and take the limit as ϵ → 0. This gives:

(16)

where f(t)=1 s⁻¹ in steady-state conditions. For non-radioactive primary particles, Γ_ISM = n_ISMσ_D(K)v, where v is the particle velocity and σ_D(K) the total destruction cross section for that particle due to interactions in the ISM of homogeneous density n_ISM.

Equation (16) is then solved by discretizing in energy space, as detailed in the Appendix, providing P_i(K, 0) and thus P_i(K, z) for all i, from which the particle spectral density is finally obtained at all positions (r, z) by summing over a sufficiently large number of terms (see Eq. (10)). In practice, the truncated series , is slowly oscillating when N increases, and we found that accelerated convergence is ensured by computing the average of a large enough number of terms, beyond a certain order. Here, we used for most calculations .

3.3. Secondary particles

The above approach can easily be extended to compute the distribution of secondary particles, produced in flight by the interactions of the primary cosmic rays with the ISM, through the standard spallation process. The formalism remains the same with a source term appropriate for spallation. Noting σ_{I + T → II}(K′,K) the differential cross section for the production of secondary nuclei II at energy K, by the interaction of a primary cosmic ray of type I of energy K′ with a target nucleus of type T with density , and summing over all spallation channels, the source term for nuclei S is given by:

(17)

For our present purposes, we approximate this source term by assuming that spallation products in the cosmic rays keep the same energy per nucleon as its energetic progenitor, and consider only protons in the ISM, with an average density n_ISM ≃ 1.3 cm⁻³, with cross sections taken from Silberberg & Tsao (1973).

This leads to the Fourier-Bessel coefficient of order i for the source term of secondary nuclei of type II:

(18)

where A_I and A_II are the atomic mass numbers of the primary and secondary particles, respectively.

3.4. Solar modulation

To fully describe CR transport to Earth, one has to include the influence of the Sun for the very last part of their flight. Close to Earth, CRs penetrate the Sun’s sphere of influence and are subjected to a phenomenon called Solar Modulation (“Smod”). The solar wind and associated magnetic field significantly reduce the kinetic energy of low energy CR (T ≲ 10 GeV/n) and prevent these from reaching our planet. This effect can be effectively described by a Fisk potential Φ_F in the “force field approximation”. The flux in the local interstellar environment (LIS) is modulated to obtain the flux on Earth Φ_⊕(K) as follows (Gleeson & Axford 1968; Boudaud et al. 2015b)

(19)

For Pamela data, conservative estimates of Φ_F are 0.1 GV <  Φ_F <  1.1 GV and for AMS 0 GV <  Φ_F <  2 GV.

4. Results

4.1. Typical expectations

To evaluate contributions from episodes of intense activity at the Galactic centre, we first evaluate the power that needs to be injected in GCCRs in these episodes to obtain local CR fluxes comparable to those observed at Earth. We first use benchmark values for the propagation parameters, taken from models that reproduce the observed secondary-to-primary ratios with a distribution of primary sources corresponding to supernovae remnants (Case & Bhattacharya 1998; Maurin et al. 2001). These models are called MIN, MED and MAX (for minimum, medium and maximum, referring primarily to the thickness of the halo) and are referenced in Table 1. To obtain the required power injected as CRs at the Galactic centre (see Eq. (8)), we specify the injection rate dN/dK by specifying values for the logarithmic index α and the normalization coefficient N₀ in the momentum power law. We choose values between 2.2 and 2.4 (Achterberg et al. 2001) for the former, and adjust the value of the latter numerically so that the calculated flux matches the data at high energy (here K = 600 GeV). Thus, coefficient N₀ is just a scaling factor which does not affect trends in the results. The data are taken from the Cosmic Ray Database (Maurin et al. 2014 and references therein).

We find that for the MIN model the required injected power is much larger than that available, due to the small halo size, L (see Sect. 4.2). However, recent studies based on synchrotron radio emission (Di Bernardo et al. 2013; Fornengo et al. 2014) but also positrons (Boudaud et al. 2015a; Lavalle et al. 2014) and anti-protons (Giesen et al. 2015), do not support the thin halo of MIN models. For the other two benchmark models (MED and MAX), a fraction of the global power budget from high energy events at the Galactic centre (Guo & Mathews 2012; Miller & Bregman 2016) is sufficient to match the observed flux.

Table 1 gives a list of the parameters of the different models and the corresponding injection power, as well as the parameters and injection power for the models of Sect. 4.2. Figure 2 shows the cosmic-ray spectral densities at Earth obtained from simulations for the three benchmark models, assuming the indicated injection power, compared to the observed flux.

Fig. 2. H protons flux Φp at Earth, rescaled by K2, for the three standard models (at the obtained injection powers) against kinetic energy K. Data from Pamela and AMS02 are also displayed.

From Fig. 2, it appears that with only 10% of the matching injection power (downward shift of the curves by one decade), Galactic centre bursts can still be expected to generate features in the observed flux. With the parameters of the MED or MAX models, a CR injection power of ∼1–6  ×  10⁴⁰ erg s⁻¹, i.e. merely a few percent at most of the total energy injection rate of 1–7  ×  10⁴² erg s⁻¹ of events leading to Fermi Bubbles (Guo & Mathews 2012; Miller & Bregman 2016), is sufficient to make a significant contribution to the locally observed GCR fluxes, at least in some energy range. We note that these total injection powers in GCRs are close to those that are derived from models involving a classical source distribution (Strong et al. 2010).

From this simple estimate, one may conclude that a global model of the local GCRs should not be limited to standard sources distributed throughout the Galaxy and should also include the contribution of GCCRs. For there are only two alternatives: the contribution of these cosmic rays is negligible or significant. In the former case, one must understand why the high energy sources that are present in the Galaxy centre do not contribute much to the local GCRs. This could be due to a very inefficient conversion of the available power into energetic particles or to an unexpectedly thin halo, but this would likely conflict with previous arguments (Di Bernardo et al. 2013; Fornengo et al. 2014; Boudaud et al. 2015a; Lavalle et al. 2014). If the contribution of GCCRs is significant, the “standard parameters” derived for the propagation of GCRs should be modified accordingly. As shown below, GCCRs might even account for the majority of the observed GCRs for some particular parameter choices, which emphasizes further that common results on GCR propagation rely on specific assumptions regarding the acceleration and transport of energetic particles from high-power events in the central part of the Galaxy.

4.2. Influence of parameters

4.2.1. Energy losses and re-acceleration

Energy losses and re-acceleration play a role in the low energy range of the GCRs spectrum. The typical timescales for energy losses, τ_loss = −K/b_loss(K), and diffusive re-acceleration, τ_DR = K²/χ(K), can be compared with the typical confinement timescale either derived from a one dimensional slab model or from the Leaky Box model (e.g. Boudaud et al. 2015b). In these models involving a homogeneous distribution of sources, both energy losses and re-acceleration become negligible at kinetic energies above of a few GeV.

In the case of a central source, however, the relevant timescale is the typical time needed to reach the Earth from the central source by diffusion τ_Diff, which is of the order of r_⊙²/2D(K), where r_⊙ ≈ 8.5 kpc the distance of the Solar-system to the Galactic centre. These timescales are represented in Fig. 3 (right). We find that energy losses are significant up to ∼10² GeV. This can also be verified with simulations by turning on and off the energy losses and diffusive re-acceleration. Figure 3 (left) shows the resulting fluxes for the Kr2.4 model (see Table 1), with arbitrary normalization. Also shown is the relative difference ΔΦ/Φ = |Φ_loss − Φ_no−loss|/Φ_loss, which can be seen to drop below 10% above ∼150 GeV in this configuration. At higher energies, parameters V_wind and v_A have essentially no incidence on the cosmic ray flux, whose level is then controlled by the other parameters.

Fig. 3. Left panel: H protons flux Φp at Earth, for the Kr2.4 model with and without energy losses and re-acceleration against kinetic energy (K). Solar modulation was turned off in both cases, such that ΦF = 0. Right panel: timescales ratios against kinetic energy for the Kr2.4 model.

4.2.2. Halo size

As mentioned in Sect. 4.1, the halo size has a large influence on the power that needs to be injected to obtain CR fluxes comparable with the ones observed. More particles indeed escape the halo before arriving in the Earth’s vicinity as the confinement box gets smaller. This can be seen within a purely diffusive model, which is a good approximation of our complete model at high energy (see Sect. 4.2.1). Figure 4 shows the corresponding density profiles in the Galactic disk as a function of the galactocentric radius r, for various sizes of the halo, L, and for a fixed diffusion coefficient D (and thus energy K). As can be seen, at a distance r_⊙ ≈ 8.5 kpc, the size of the halo becomes critical around L ≤ 2 kpc: for lower values of L, the local spectral density Ψ drops significantly below from that obtained for L = 5 kpc, L = 10 kpc or for an infinite box case. This explains why a very large injection power is needed to compete with locally observed GCR flux in benchmark model MIN with a thin halo.

Fig. 4. Density profiles for a various set of halo size L for a purely diffusive equation. The profile for an infinite box (in both direction) is also displayed and is labeled 1/r as Ψ = N0/4πDr in this case.

4.2.3. Degeneracy in the choice of parameters

By solving Eq. (16) without energy losses and re-acceleration (which only influence the low-energy part of the spectrum), one may write Ψ(K, r, z = 0)=N₀K^−α/D(K)×g(r, R, h, L), where g is a function of spatial variables only. Consequently, the power injected as GCCRs, 𝒫_inj, required to produce CR fluxes comparable to the observed ones at Earth is (in our model) a function of the halo size L, the diffusive coefficient D and the spectral index α. However, there is some degeneracy in the choice of parameters, similar to that expressed by the three benchmark models. More precisely, a negligible influence of the activity at the Galactic centre on the observed GCRs can be obtained by invoking either an unusually thin halo or a particularly large diffusion coefficient for example (since Ψ ∝ 1/D(K) at high energy). Conversely, if a source at the Galactic centre turns out to have a significant impact on the observed fluxes, some degeneracy must be accepted at this stage between L and D.

With this caveat in mind, we investigated whether it could even be possible to account for the entire GCR spectrum with a single central source. Figures 5 and 6 show flux values that are obtained for the Kr2.4 model, i.e. a model with a source spectrum in p^−2.4 and a Kraichnan-like turbulence scaling for the energy diffusion coefficient (see above). With the parameters given in Table 1, we find that an injection power 𝒫_inj = 2.8 × 10⁴¹ erg s⁻¹ allows an excellent fit to the observed flux. Such a power appears accessible in principle, given the average power range of 1–7  ×  10⁴² erg s⁻¹ available from high energy events at the Galactic centre.

Fig. 5. H protons flux Φp at Earth, rescaled by K2, for the Kr2.4 model (see Table 1) against kinetic energy K.

Fig. 6. Carbon flux ΦC at Earth, rescaled by (K/n)2, for the Kr2.4 model (see Table 1) against kinetic energy per nucleon K/n.

Other sets of parameters can also reproduce the data for reasonable values of the injection power. An example is shown on Fig. 7, with a model using a Kolmogorov-like diffusion coefficient (δ = 1/3) and a halo size L = 10 kpc. All the other parameters are also displayed.

Fig. 7. H protons flux Φp at Earth, rescaled by K2, for a Kolmogorov model, against kinetic energy K. Propagation parameters are displayed and we obtain an injection power 𝒫inj = 6 × 1041 erg s−1.

In addition to reproducing the flux of primary GCR nuclei, it is instructive to investigate whether the observed ratios of secondary nuclei to primary nuclei (so called secondary-to-primary ratios) and secondary radioactive nuclei to secondary stable nuclei ratios could also be matched by a central source of cosmic rays. This is addressed in the next sections.

4.3. Secondary-to-primary ratios

From the secondary-to-primary ratio (II/I), one classically constrains the diffusion coefficients parameters (Genolini et al. 2015) using the high energy part of the ratio. Indeed, in that energy range where energy losses can be neglected, the ratio (II/I) scales as D_I/D_II² where D_I and D_II are respectively the diffusion coefficient of species I and II at the same given kinetic energy per nucleon K/n (see Sect. 3.3).

In our case, the energy domain in which energy losses are inefficient and the diffusion coefficient can be estimated directly corresponds to K/n ≳ 10² GeV, where the observational data become scarse and have larger uncertainties. We can nevertheless estimate the normalisation coefficient, D₀, to be of the order of 0.07–0.1 kpc² Myr⁻¹, should the GCCR Boron-to-Carbon (B/C) ratio at Earth be matching the observed one at high energy (see Fig. 8).

Fig. 8. Boron to Carbon (B/C) ratio, at high energy (K ≥ 200 GeV), at Earth, vs. the kinetic energy per nucleon K/n, for L = 10 kpc, δ = 0.5 (Kraichnan turbulence) and two values of the parameter D0.

Once D₀ is fixed in this way, the B/C ratio of the GCCRs at Earth can be obtained over the entire energy range, under the assumption that the other parameters are chosen, say, to reproduce the observed primary spectra (e.g. protons and carbon nuclei). The results are displayed on Fig. 9, where we show the B/C ratio as a fuction of kinetic energy per nucleon, for the models Kr2.4 and Kol2.55, corresponding to a Kraichnan-like and Kolmogorov-like energy scaling of D(K), respectively (see Table 1). As can be seen, if one adopts the source spectral indices, resp. 2.4 and 2.55, which lead to primary spectra matching the observed ones, the resulting evolution of the B/C ratio with energy is very different from the observed one in the case of a Kolmogorov-like scaling (Kol2.55 model), but remarkably close to it in the case of a Kraichnan-like scaling (Kr2.4 model).

Fig. 9. Boron to Carbon (B/C) ratio at Earth against kinetic energy per nucleon K/n for the Kr2.4 (left panel) and Kol2.55 (right panel) models (see Table 1).

4.4. “Cosmic-ray clocks”

Another important ingredient of GCR phenomenology is related to the so-called cosmic-ray clocks, i.e. radioactive secondary nuclei with half-lives comparable with the GCR dynamical timescales, which can therefore provide some information about the timescale(s) of cosmic-ray transport in the Galaxy (typical age, confinement or homogenization time, depending on the model). The key observable, in the respect, is the radioactive-to-stable isotopic ratio for appropriate secondary nuclei. A classical example is the ¹⁰Be/⁹Be ratio: both isotopes are produced exclusively by (inverse) spallation during CR transport, and ¹⁰Be has a half-life of τ_1/2 ≃ 1.4 Myr, while ⁹Be is stable.

Qualitatively, it is easy to see how such cosmic-ray clocks can be used in standard GCR phenomenological models to constrain, for instance, the size of the magnetic halo. At low energy, the GCR confinement time is large enough (and the relativistic time dilation effect weak enough) that most ¹⁰Be nuclei produced in flight decay before escaping the Galaxy. The resulting ¹⁰Be/⁹Be ratio is thus expected to be “low”, at a level that depends on the size of the halo, L: for thicker halos, more radioactive nuclei will be able to decay before escaping, thereby reducing the radioactive-to-stable ratio (see e.g. Strong & Moskalenko 1998). At higher energy, the confinement time is reduced and the time dilation effect increases, so that for sufficiently large energies, the observed ratio is essentially expected to match the production ratio, taken to be the same as in Strong & Moskalenko (1998) for example. The lack of precise measurements over a sufficient large energy range, however, prevents an accurate determination of the halo size. Thus, the L/D degeneracy cannot totally broken, as illustrated by benchmark models that are still in use. For a more detailed description, see the review of Strong et al. (2007) and references therein.

The above behaviour of radioactive-to-stable ratios is of course expected to be similarly observed in the case of GCCRs. It would be misleading to expect that the resulting ¹⁰Be/⁹Be ratio be significantly lower in the case of GCCRs, because the typical propagation time from the central source to the Earth vicinity is longer than the typical age of cosmic-rays reaching us from more nearby sources in the case of distributed GCR source scenarios (spanning the Galactic disk). However, both Be isotopes are secondary, and even in the case of GCCRs, most of those who are actually observed on Earth are still produced by spallation in our local part of the Galactic disk, at a rate and level that essentially depends on the nearby ISM density and local GCR flux, i.e. not very different from what occurs in the case of distributed sources.

Quantitatively, it is easy to extend our approach, as described in Sects. 2 and 3, to compute the distribution of secondary radioactive nuclei as well: one simply needs to add a term −Γ_radΨ in the right hand side of Eq. (2), where Γ_rad = ln(2)/(γ τ_1/2) is the decay constant, taking into account time dilation for radioactive nuclei with Lorentz factor γ = γ(K).

The results are shown in Fig. 10 for the two models Kr2.4 and Kol2.55 models and halo sizes L = 5, 10 and 20 kpc. As can be seen, ¹⁰Be/⁹Be ratios similar to the observed ones are obtained, for halo sizes L ≳ 4 kpc, in conformity with the results of Sect. 4.2.2. This confirms that a significant fraction of GCCRs could be present among the local cosmic-rays without disturbing significantly the familiar phenomenology.

Fig. 10. 10Be/9Be ratio at Earth as a function of kinetic energy per nucleon K/n for the Kr2.4 (left panel) and Kol2.55 (right panel) models. The different curves correspond to different halo size, L = 5, 10 and 20 kpc, as indicated (the other parameters are fixed as in Table 1).

5. Summary and discussion

In the absence of a definite model that would be supported by clear evidence, it appears important to revisit in a broader perspective what is often considered common knowledge regarding cosmic ray origin and propagation in the Galaxy. Most of the GCR phenomenological studies assume that cosmic ray sources are distributed throughout the Galatic disk, with a space and time granularity that can be neglected, at first order, with respect to the typical source seperation distances and repetition times. While this is a reasonable assumption, especially since such distributed sources are indeed known (notably SNRs and superbublles), it does not exclude that other types of sources may also contribute to the GCRs.

In the framework of GCR studies assuming distributed sources, reasonable estimates of the GCR diffusion coefficient and confinement time could be obtained, which in turn provided an estimate of the required GCR source power, comparable to a fair fraction of the total power released by SN explosions in the Galaxy. It turns out, however, that similar or even possibly larger average power is now known to be released around the Galactic centre through episodic events, and that these events lead to astrophysical conditions that appear to be potentially propitious for particle acceleration. According to Guo & Mathews (2012) and Miller & Bregman (2016), the event that led to the so-called Fermi bubbles released a total energy of 10⁵⁷ erg s⁻¹, for an average injection rate of 1–7  ×  10⁴² erg s⁻¹, which exceeds the total kinetic power of SN explosions in the interstellar medium by large amounts. It is thus natural to investigate whether such events could play a role in the global GCR phenomenon, and if so how the phenomenology would be affected.

The main result of this paper is the confirmation that, with the same type of parameters as those derived from distributed-sources GCR studies, the contribution of the GCCRs cannot be a priori neglected. We found, indeed, that save for models with very thin halos (e.g. MIN), the required injection power to match the observed GCR fluxes (𝒫_inj ≲ 10⁴¹ erg s⁻¹) is only a fraction of the inferred available power. For some transport parameters, the GCCRs could even dominate the GCR flux over a large range of energies. If this is the case, then they will have to satisfy the observational constraints directly. To see if this could be possible, we investigated the caracteristics of the GCCRs as if they were alone in the Galaxy, and found that, even with simple generic assumptions and first order modeling (see Sect. 2), it is indeed possible to reproduce the locally observed primary CR spectra as well as secondary-to-primary ratios (B/C) and secondary radioactive-to-stable ratio (¹⁰Be/⁹Be), although the latter is only true for a Kraichnan-like energy dependence of the diffusion coefficient.

However, even a sub-dominant (but non-negligible) contribution of the GCCRs may have important consequences on the general CR phenomenology. In particular, the source composition and source spectrum of the GCCRs may be expected to be somewhat different from those of the other sources. In regions of the spectrum where both contributions are roughly similar in magnitude, this may lead to gradual changes in the composition (hence different spectra for different nuclei), or to specific features in the elemental spectra. Such effects will be studied in a forthcoming paper.

If the GCCRs are able to contribute at a non-negligible level, some of the conclusions of the standard approach regarding the transport parameters of the energetic particles in the Galaxy may also have to be revised, which can result in the relaxation of some of the usual constraints. In particular, the modeling of the multi-wavelength emission of SNRs may not have to be done with the requirement that the associated energetic particles be also compliant with the entire GCR phenomenology, including the maximum energy problem or the elemental and isotopic ratios. Likewise, it may be possible to relax the constraints associated with the study of the so-called cosmic-ray clocks (radioactive secondary nuclei), if two main components with different timescales are mixed together among the GCRs.

The fact that GCCR typically need more time to reach the Solar System than the average cosmic ray from more evenly distributed sources has also some consequences on the link between CR nuclei and leptons. Since electrons and positrons rapidly lose energy as they propagate in the Galactic magnetic field, those observed at Earth must have been generated relatively nearby. GCCRs can only contribute secondary e^±, produced by collisions of high energy protons and helium nuclei on the atoms of the ISM (see e.g. Delahaye et al. 2009), while the bulk of the energetic leptons should still be due to local Galactic sources. Of course, such a “decoupling” between nuclei and lepton sources is not specific to our approach, and is also suggested for instance in the context of the observed positron anomaly (Adriani et al. 2009) where nearby pulsar sources can be invoked (Hooper et al. 2009; Profumo 2012; Linden & Profumo 2013). We note that this decoupling may also alleviate some of the problems that arise when trying to reconcile observations with models of dark matter (Boudaud et al. 2015a; Cirelli et al. 2009).

In this paper, we have estimated the average contribution of the GCCRs in the vicinity of the solar system, assuming steady state. Because of the very nature of the potential GCCR sources, the assumption of a constant particle injection is clearly wrong. However, for not too high energies, i.e. as long as the “confinement time” of the particles is much larger than the repetition time between acceleration events, the steady state solution provides an acceptable approximation.

To estimate the energy range where the steady state assumption can be expected to be valid, we can compare the relevant timescales. The most recent major event in the Galatic centre is suggested to have occurred ∼3 Myr ago (Miller & Bregman 2016) and to have lasted for ∼0.1 − 0.5 Myr (Guo & Mathews 2012), leading to the so-called Fermi Bubbles (Acero et al. 2016). We may thus assume a repetition time of the order of a few Myr. For each event, GCCR injection occurs on a much shorter timescale and can be approximated here as being instantaneous. This gives rise to a diffusion front propagating outwards from the Galactic centre. If diffusion were isotropic in a homogeneous medium with diffusion coefficient D, the density of particles at distance r from the source, at a time t after injection, would be simply given by the Green function:

(20)

where N₀ is the number of particles injected at r = 0 at t = 0.

At any given position r, this density rapidly increases to a maximum and decreases more slowly as the diffusion sphere expands. The typical duration of the event, as seen at radius r, can be estimated as the time during which the particle density is larger than half of its maximum value: Δt_1/2 ≃ 0.45r²/D.

Thus, at radius r, the steady-state solution provide a good approximation of the actual GCCR flux at any given time up to an energy K such that D(K)< 0.45r²/Δt_s, where Δt_s is the typical time interval between two source episodes. Numerically, at the solar radius, this gives:

(21)

With the parameterization of the diffusion coefficient adopted above, , with D₀ ∼ 0.07–0.1 kpc² Myr⁻¹, this corresponds to rigidities

(22)

In particular, for the Kr2.4 model, the above solution is roughly valid up to ∼40 TeV, and for the Kol2.55 model, up to 8 PeV.

At higher energy, however, the intermittent nature of the central source will affect the main characteristics of the GCCR. Qualitatively, one should expect a reduction of the high energy particles at any time when the diffusion front from the last event at these energies has already passed the solar radius, producing a knee-like feature. Ankle-like features can also be obtained at energies where the diffusion front from the last event is arriving at the solar radius at the time of observation, while the previous one has already left. These effects will be studied in more detail together with associated composition features in a forthcoming paper (Jaupart et al., in prep.).

The previous considerations can be extended to any position inside the Galactic disk. For any galactocentric radius r <  20 kpc, there is a critical rigidity R_crit (given by Eq. (22)) above which the steady state approximation is no longer valid. More specifically, in the inner region of the Galaxy, at r ≃ 1 kpc for example, this yields R_crit ∼ 2.8^1/δ GV, corresponding to a few GeV of kinetic energy in the Kr2.4 and Kol2.55 models. The diffusive front from the last event has thus gone through the inner regions of the Galaxy (r ≲ 1 kpc) for GCCRs with a kinetic energy above the GeV level. Thus, these regions are depleted in these GCCRs, which acts to flatten the distribution given by Fig. 4. This is also the energy range which dominates the production of gamma-rays from π⁰ decay. Therefore, in order to compute the gamma-ray background associated with the GCCRs and compare it to the observation, a more detailed, time-dependent treatment is needed. This will be addressed in a separate paper.

Acknowledgments

The authors are grateful to Pierre Salati for his useful advice and comments. They also acknowledge support from Agence Nationale de la Recherche (grant ANR-17-CE31-0014). This work is partially supported by an ENS scholarship to EJ from the University of Lyon.

References

Appendix A:

Numerical procedure to account for energy losses

We start from Eq. (16) and define τ_D, j and Q̃_j by

(A.1.)

(A.2.)

We then re-write Eq. (16) to obtain:

(A.3.)

where α_j = K/τ_D, j, P̄_j = τ_D,j × Q̃_j, x = log(K) and γ = χ(K)/K.

We then discretize the previous equation on a logarithmic scale with N + 1 values between K_min and K_max with step of . Hence the kth value on that grid is x_k = log(K_min) + k Δx. We obtain the discretized equation:

(A.4.)

that is:

(A.5.)

(A.6.)

(A.7.)

(A.8.)

namely:

(A.9.)

where [P̄_j] = M(a, b, c)[P̄_j] is a matrix to be inverted.

Regarding the boundary conditions, we need to specify them on our energy grid to perform the inversion. At high kinetic energy, we expect that energy losses and diffusive reacceleration are not significant so that:

(A.10.)

At low energy, we note that the primary cosmic-ray data can be fitted by power law spectra, and thus impose this as an ad hoc boundary condition, which can be expressed as:

(A.11.)

These boundary conditions are similar to those found in other generic studies (Donato et al. 2001) of GCRs phenomenology and allows to invert the relation [P̄_j] = M(a, b, c)[P̄_j] to obtain each Bessel order P_j.

All Tables

All Figures

	Fig. 1. Geometry of the diffusive volume.
In the text

	Fig. 2. H protons flux Φp at Earth, rescaled by K2, for the three standard models (at the obtained injection powers) against kinetic energy K. Data from Pamela and AMS02 are also displayed.
In the text

	Fig. 3. Left panel: H protons flux Φp at Earth, for the Kr2.4 model with and without energy losses and re-acceleration against kinetic energy (K). Solar modulation was turned off in both cases, such that ΦF = 0. Right panel: timescales ratios against kinetic energy for the Kr2.4 model.
In the text

	Fig. 4. Density profiles for a various set of halo size L for a purely diffusive equation. The profile for an infinite box (in both direction) is also displayed and is labeled 1/r as Ψ = N0/4πDr in this case.
In the text

	Fig. 5. H protons flux Φp at Earth, rescaled by K2, for the Kr2.4 model (see Table 1) against kinetic energy K.
In the text

	Fig. 6. Carbon flux ΦC at Earth, rescaled by (K/n)2, for the Kr2.4 model (see Table 1) against kinetic energy per nucleon K/n.
In the text

	Fig. 7. H protons flux Φp at Earth, rescaled by K2, for a Kolmogorov model, against kinetic energy K. Propagation parameters are displayed and we obtain an injection power 𝒫inj = 6 × 1041 erg s−1.
In the text

	Fig. 8. Boron to Carbon (B/C) ratio, at high energy (K ≥ 200 GeV), at Earth, vs. the kinetic energy per nucleon K/n, for L = 10 kpc, δ = 0.5 (Kraichnan turbulence) and two values of the parameter D0.
In the text

	Fig. 9. Boron to Carbon (B/C) ratio at Earth against kinetic energy per nucleon K/n for the Kr2.4 (left panel) and Kol2.55 (right panel) models (see Table 1).
In the text

	Fig. 10. 10Be/9Be ratio at Earth as a function of kinetic energy per nucleon K/n for the Kr2.4 (left panel) and Kol2.55 (right panel) models. The different curves correspond to different halo size, L = 5, 10 and 20 kpc, as indicated (the other parameters are fixed as in Table 1).
In the text