© The Authors 2024

1. Introduction

Building a consistent description of the cosmic-ray (CR) spectrum has been a challenging theoretical task for decades. While it is well established that these relativistic charged particles are accelerated by magnetic shocks (Fermi 1949; Bell 1978; Hillas 1984; Blandford & Eichler 1987), a comprehensive understanding of the relative contributions of all these sources is still lacking. It is thought that the Galactic cosmic-rays (GCRs) dominate the overall CR population for energies below ∼10¹⁸ eV, while substructures in the CR spectrum, the so-called knees and ankle, are considered as indicators of a transition between Galactic and extragalactic CRs (Aloisio et al. 2012).

In common scenarios, GCRs find their origin in magnetic shocks produced by supernovae (SNe) remnants (Strong et al. 2007). The total thermal energy released by an such event, of order 10⁵¹ erg, is believed to be partially converted into CRs with a yield of ∼10%, and with energies of up to a few 10–100 TeV (Lagage & Cesarsky 1983).

The case for the existence of other significant CR sources within the Galaxy has recently been bolstered by the discovery of large hemispherical structures above and below the Galactic plane by the ROSAT, Fermi and eRosita space telescopes (Bland-Hawthorn & Cohen 2003; Su et al. 2010). The soft-X-ray component of these so-called Fermi bubbles (FBs) consists of spherical bright shells extending up to 14 kpc away from both sides the galactic plane. The edge of these shells seems to coincide with that of a gamma-ray component, filling its inner volume with a relatively uniform surface brightness. Two main scenarios have been put forward to explain these features: starbust-driven disk outflows (Sofue 2000) or a past activity of the central supermassive black hole (SMBH), Sgr A^⋆ (Guo & Mathews 2012). Both scenarios involve a total energy injection E_F of at least 10⁵⁵ erg, with the energy linked to the latter scenario ranging between 10⁵⁶ and 10⁵⁷ erg (Guo & Mathews 2012; Yang et al. 2022). Also, hydrodynamical simulations allow constraints to be placed on their age, which is estimated to be ∼1–6 Myr for the SMBH jet-driven scenarios (Guo & Mathews 2012; Zhang & Guo 2020), and up to 50 Myr for nuclear star formation-driven bubbles (Miller & Bregman 2016). Whether or not the expanding X-ray shells and the volume-filling gamma-ray component correspond to a single event is still unclear, because accretion from the SMBH is a stochastic phenomenon that can leads to repetitive emissions of outbursts separated by time periods of as little as ∼10⁵ yr (King & Nixon 2015).

Theoretical studies suggest that two main scenarios can be proposed to explain the observations of the FBs, where the gamma-ray emission is the result of either (i) CR protons interacting with the surrounding medium (“hadronic” model) or (ii) inverse Compton scattering by CR electrons (“leptonic” model) (Yang et al. 2012). In both cases, the event leading to the formation of the FBs is expected to have accelerated a large quantity of CRs. Whether or not this specific population of proton and electron CRs could make a significant contribution to the CR spectrum measured on Earth remains unclear. Theoretical work following the isotropic propagation of CRs in a simplified Milky Way (MW) model led the authors to conclude that only a small fraction of the assumed injected energy E_F is found to be sufficient to significantly affect the CR spectrum measured in the solar neighborhood (Jaupart et al. 2018).

In the present paper, we revisit these conclusions using a realistic magnetohydrodynamical numerical model of the MW. In particular, our model includes magnetic fields and CR diffusion along the magnetic field lines. More specifically, we study the diffusion of CRs accelerated by SN remnants or by a central source (CS). The paper is organized as follows. In Sect. 2 we present our numerical setup. In Sect. 3 we present how CRs of different origin diffuse into the Galaxy for two different energies, and discuss the implications for the CR spectrum observed in the solar neighborhood. Conclusions and perspectives are presented in Sect. 4.

2. Simulations and method

2.1. Magnetohydrodynamics with anisotropic cosmic-ray diffusion

To simulate the magnetohydrodynamical evolution of an isolated MW Galaxy, we used the adaptive mesh refinement code RAMSES (Fromang et al. 2006). We used the Harten–Lax–van Leer Discontinuity (HLLD) (Miyoshi & Kusano 2005) Riemann solver to calculate the MHD flux. Computation of the CR energy density flux was performed using the implicit numerical solver introduced in Dubois & Commerçon (2016) and Dubois et al. (2019), which we modified to support several CR groups. The full system is then described by the following set of equations:

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·[\rho \mathbf{u}]=0,\hfill \\ \hfill & \frac{\partial \rho \mathbf{u}}{\partial t}+\mathrm{\nabla }·[\rho \mathbf{u}\otimes \mathbf{u}+P_{\text{tot}}\mathbb{I}-\frac{\mathbf{B}\otimes \mathbf{B}}{4\pi }]=0,\hfill \\ \hfill & \frac{\partial E_{\text{tot}}}{\partial t}+\mathrm{\nabla }·[\mathbf{u}(E_{\text{tot}}+P_{\text{tot}})-\frac{(\mathbf{u}·\mathbf{B})\mathbf{B}}{4\pi }]\hfill \\ \hfill & =-P_{\text{CR}}\mathrm{\nabla }·\mathbf{u}+\mathrm{\nabla }·(\mathbb{D}\mathrm{\nabla }{E}_{\text{CR}})-\mathcal{L}(\rho ,T),\hfill \\ \hfill & \frac{\partial \mathbf{B}}{\partial t}-\mathrm{\nabla }\times [\mathbf{u}\times \mathbf{B}]=0,\hfill \\ \hfill & \frac{\partial E_{\text{CR}}}{\partial t}+\mathrm{\nabla }·(\mathbf{u}{\mathbf{E}}_{\mathbf{\text{CR}}})=-P_{\text{CR}}\mathrm{\nabla }·\mathbf{u}+\mathrm{\nabla }·(\mathbb{D}\mathrm{\nabla }{E}_{\text{CR}}).\hfill \end{array}$$(1)

Here, ρ is the gas density, u its velocity field, B is the magnetic field vector, P_tot is the total pressure including the gas pressure P = (1 − γ)E, the magnetic field pressure P_mag = B²/(8π), and the CR pressure P_CR = (γ_CR − 1)E_CR, where γ and γ_CR are the adiabatic index of the gas and the CR, respectively. As the CRs are considered here as a relativistic fluid, we take γ_CR = 4/3, while the gas is nonrelativistic and thus characterized by the index γ = 5/3. Here, E is the gas internal energy density, E_CR is the CR energy density and E_tot designates the total energy density, and E_tot = E + 0.5ρu² + B²/(8π)+E_CR. Also, ℒ(ρ, T) is a cooling function modeling the heating and cooling processes taking place inside the ISM, and 𝔻 describes the diffusion tensor of the CR, and can be decomposed following:

$$\begin{array}{c}\hfill \mathbb{D}=D_{\perp }{\delta }_{\mathit{ij}}-(D_{\perp }-D_{\Vert })\mathbf{b}\times \mathbf{b},\end{array}$$(2)

where b = B /∥B∥ is the normalized magnetic field vector and D_∥ is the parallel diffusion coefficient, accounting for the diffusion along magnetic field lines. D_⊥ is the perpendicular diffusion coefficient, effectively taken as 10⁻² × D_∥, which accounts for the diffusion of CRs perpendicularly to the magnetic field lines.

2.2. Initial conditions

The initial conditions are similar to those in Renaud et al. (2013). We initially set a box of width L_box = 120 kpc filled with an axisymmetric distribution of gas and stellar particles. Particles position and velocity distributions of the particles are calculated using the Multi-Gaussian Expansion (MGE) method (Emsellem et al. 1994), which uses a set of Gaussian functions to model the various component of the MW based on a given input, here chosen as the Besançon model (Robin et al. 2003). The magnetic field is initially seeded with a toroidal component of amplitude 1 μG.

A spherically symmetric dark matter halo is modeled by the following density profile:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \rho (r)\propto {\{1+{\left(\frac{r}{r_c}\right)}^2\}}^{-1},\end{array}\end{array}$$(3)

with a characteristic radius r_c of 2.7 kpc and up to a truncation radius at 50 kpc, using 4 × 10⁶ particles of 1.13 × 10⁵ M_⊙. A pre-existing star population of 4.6 × 10¹⁰ M_⊙ is modeled with 4 × 10⁶ particles of 1.15 × 10⁴ M_⊙. This population is decomposed into a bulge, a spheroid and a thick and a thin disk following Renaud et al. (2013). In addition, a particle of mass 4 × 10⁶ M_⊙ is located at the center of the box to model the SMBH. The gas component is modeled by an exponential profile with a radial (respectively vertical) scale length of 10 kpc (0.25 kpc) and truncation at 30 kpc (2.5 kpc), with a total mass of 5.94 × 10⁹ M_⊙.

The coarsest cells have a size of 937.5 pc. We allow refinement down to 14 pc following two criteria. We refine a cell if (i) there are more than 20 particles from initial conditions or more than 2 × 10⁴ M_⊙ total mass in the cell, or (ii) the local thermal Jeans length is smaller than four times the cell size (Truelove et al. 1997).

As gas and stars do not initially share the same gravitational potential, we first need to simulate a phase of “warm relaxation” for about 1 Gyr at the beginning of our simulation with no star formation and an isothermal equation of state at T = 5 × 10⁴ K so that any structure formed during the relaxation of the gas+particles system can dissipate. During this phase, instabilities lead to the formation of the asymmetric structures of the Galaxy, namely the bar and the spiral arms.

Once these structures are formed, we activate gas cooling and heating, and star formation. Gas cooling is computed via cooling tables tabulated at solar metallicity available in the public version of RAMSES, with a temperature floor of 50 K. Star formation is modeled by a density based criterion: each cell with gas density of above 1 cm⁻³ and colder than 2 × 10⁴ K forms stars following the Schmidt (1959) law with 10% efficiency per free-fall time. New stars have a mass of 10⁴ M_⊙.

After an additional ∼50 Myr, the star formation rate (SFR) reaches a realistic order of magnitude (i.e. around 1 M_⊙/yr). The corresponding gas density and magnetic field maps are shown in Fig. 1. CRs injection from SNe is then turned on and the simulation is continued for an extra 5 Myr so that the CRs injected from the SNe have enough time to diffuse significantly. Time at the end of this phase is hereafter referred to as t = 0.

Fig. 1. Initial conditions for CR injection. The left panel shows the hydrogen density, while the right panel shows the magnetic field amplitude overlaid with its associated field lines.

2.3. Stellar feedback

We include feedback from type II SNe. Injection is performed for stellar particles with an age of more than 10 Myr and a mass superior to 10 M_⊙. At each injection event, mass and specific thermal energy are injected into the parent gas cell of the stellar particle following:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill m_{\mathrm{SNe}}& =f_{\mathrm{SNe}}\times m_{\star },\hfill \\ \hfill e_{\text{SNe}}& ={10}^{51}\text{erg}/(10{\text{M}}_{\odot }),\hfill \end{array}\end{array}$$(4)

where m_⋆ is the mass of the stellar particle, and f_SNe = 0.2.

2.4. Cosmic-ray injection and diffusion

To assess the impact of a CS of CRs on the overall CR spectrum measured at the solar neighborhood, we followed the propagation of two distinct populations of CRs injected by either the SNe or the CS, each of them being split into two energy groups.

In the following, we define the solar neighborhood as the ensemble of cells of the simulation whose position fulfills −0.5 kpc < z < 0.5 kpc and 7.5 kpc < r < 8.5 kpc, where z is the axis orthogonal to the galactic disk and r is the radial position of the cell.

2.4.1. Energy groups and injection spectrum

We follow two energy groups of CRs, defined as intervals of kinetic energies following:

$$\begin{array}{c}\hfill \begin{array}{cc}& {10}^{1.5}\le E_{\mathrm{CR},1}\le {10}^{2.5}\mathrm{GeV}\hfill \\ \hfill & {10}^{3.5}\le E_{\mathrm{CR},2}\le {10}^{4.5}\mathrm{GeV}.\hfill \end{array}\end{array}$$(5)

These groups are referred to below as either E₁ or E₂, or as their central value (100 GeV or 10 TeV respectively). We assume that the injection spectrum from the sources –that is, the SNe and the CS– follows a power law given by

$$\begin{array}{c}\hfill \frac{\text{d}n(K_{\mathrm{CR}})}{\text{d}{K}_{\mathrm{CR}}}\propto K_{\mathrm{CR}}^{-\alpha },\end{array}$$(6)

where n(K_CR) is the volumetric density of CRs, K_CR is the kinetic energy of the CRs, and α is the intrinsic source injection spectrum.

For each SNe event, the amount of CR energy theoretically injected into the whole CR energy spectrum is η_SNe × e_SNe, with η_SNe = 0.1 (Dermer & Powale 2013). The fractions of this energy to be injected into our two energy bins E₁ and E₂ is then calculated from Eq. (6), taking a slope of α = 2.4 (Lagutin et al. 2014; Blasi & Amato 2012).

We followed the same procedure for the central source. A fraction f_CS = 0.1 of E_F = 10⁵⁶ erg, the minimum total energy required to form the FBs is assumed to be redistributed as CR energy. Note that this amount of injected CR energy can be rescaled in post-processing to explore the effect of the exact value of E_F; see Sect. 3.2. The corresponding amount of energy to be injected into each energy group is obtained from the power law assuming α = 2.4, in accordance with measurements a hadronic model for the gamma-ray emission of the FBs (Ackermann et al. 2014). Several energy pulses can be injected with a fixed time interval Δt.

2.4.2. Energy losses

In this study, we only evaluate the propagation of the CRs in the Galactic disk, and not their effect on the ISM. Thus, we intentionally neglect any feedback process on the gas and any possible energy losses. In particular, we do not consider the cooling of proton CRs through pion production from proton–proton collision. The typical cooling timescale for this process is given by:

$$\begin{array}{c}\hfill t_{\mathrm{pp}}=\frac{1}{{\kappa }_{{\pi }^0}{n}_{\mathrm{H}}{\sigma }_{\mathrm{pp}}c}\simeq 5·{10}^{15}\times {\left(\frac{n_{\mathrm{H}}}{1{\mathrm{cm}}^{-3}}\right)}^{-1}\mathrm{s},\end{array}$$(7)

where κ_π⁰ ≃ 0.2 is the fraction of the protons energy transferred to the pions, σ_pp ≃ 40 mb is the collision cross section, and n_H is the hydrogen number density (Rieger 2023). Taking n_H = 10 cm⁻³ as an upper limit for the gas density in our simulated MW, we obtain a typical cooling timescale of t_pp ∼ 10⁷ yr, which is also the typical duration of our simulations. The results presented in this paper should thus be considered as an upper limit on the contribution of the CS to the local CR population.

2.5. Simulations parameters

In total, we ran six simulations whose various parameters are summarized in Table 1. The diffusion coefficient of the CRs in the ISM for each simulation is fixed and taken assuming a power law of slope β = 0.3 (Strong et al. 2007):

$$\begin{array}{c}\hfill D(E)={10}^{28}{\text{cm}}^2{\text{s}}^{-1}{\left(\frac{E}{1\text{GeV}}\right)}^{\beta }.\end{array}$$(8)

The obtained values for the 100 GeV and 10 TeV energy groups are 10²⁹ cm² s⁻¹ and 7 × 10²⁹ cm² s⁻¹, respectively. As solving the anisotropic diffusion equation of the CRs is particularly computationally expensive, our five main simulations are run with a maximum resolution of 29 pc. Two additional simulations with a maximum resolution of 15 pc and 59 pc were also run to evaluate numerical convergence. Finally, we performed one additional run with isotropic CR diffusion for comparison.

3. Results

3.1. Disk propagation and outflows

Fig. 2 presents projections of the CR energy density in both SNe and CS populations for our s15pc run and at four different times. Variable t here refers to the time since the injection of the unique pulse of 10⁵⁷ erg. The 100 GeV CRs are found to diffuse much more slowly than the 10 TeV group, which is expected due to their lower diffusion coefficient. After ∼18 Myr, most 10 TeV CRs have diffused out of the Galactic center, while a large fraction of the 100 GeV CRs are still confined within the bar.

Fig. 2. Evolution of the several CR groups at four different times after injection of a pulse from the CS. The first and third columns present the energy density maps of the CR groups from SNe, at 100 GeV and 10 TeV respectively, while the second and fourth rows presents the same quantity for the CR groups associated with the CS. Each maps are normalized by the mean energy density in the SNe group at t = 0. The gray annulus in the second column represents the solar neighborhood as defined in the main text. (Movie available online).

As we want to quantify the effect of CRs within the disk, we first evaluated how much of the injected energy stays within the disk and conversely what fraction of the CR energy outflows from the Galactic disk. As visible in Fig. 3, a significant amount of CR energy is propagates out of the disk, and forms large bubble-shaped structures in both sides of the Galactic plane. Remarkably, the spatial extent of these outflowing structures is consistent with the estimated size of the FBs (Su et al. 2010), which we evaluate to roughly 10 kpc at t = 3 Myr.

Fig. 3. Edge-on view of the 10 TeV CR group energy density injected by the CS in our s15pc run, 12 Myr after injection. As visible, Bubble-like structures can be seen to naturally emerge in our model due to the topology of the magnetic field.

To evaluate the impact of the outflow on the overall CR energy budget in the disk, we introduced f_i, outflow, the fraction of the energy injected by the CS escaping the disk for the energy group E_i, defined by:

$$f_{i,\text{outflow}}(t)=1-\frac{{\sum }_{j\in \text{disk}}{\rho }_{i,j}(t)·V_j}{{\sum }_{j\in \text{box}}{\rho }_{i,j}(t)·V_j},$$

with ρ_i, j(t) the CR energy density in cell j, for the E_i energy group and from the CS population. V_j designates the volume of the cell j. Here, the disk is defined as the set of gas cells contained in a slice 1 kpc in thickness vertically centered on the midplane of the box.

Results are shown in Fig. 4. The choice of the propagation mode is found to have minor influence on the fraction of energy staying in the disk as more than 80% of the injected CR energy escapes from the disk after at most 8 (respectively 2.5) Myr for the E₁ (E₂) energy bin for both anisotropic and isotropic runs. Intermittence is found to increase the proportion of energy contained in the disk, as recurrent injection from the CS constantly refuels the disk with additional energy that has not yet had time to outflow from the disk. One can verify this assumption by computing the diffusion length L_diff after 1 Myr, the time between two consecutive pulses:

$$\begin{array}{c}\hfill L_{\text{diff}}=\sqrt{\mathrm{\Delta }{t}_{\mathrm{p}}·D_{\Vert }}\sim 0.5\mathrm{kpc},\end{array}$$

Fig. 4. Evolution of foutflow as a function of time and for various parameters. The left panel presents the influence of the propagation model (i.e. either anisotropic or isotropic) for our 29 pc resolution, single-pulse run (s29pc). The right panel shows the influence of the intermittence of the source.

which is of the same order of magnitude as the half-height of the galactic disk. In any case, the CR energy escaping the disk and filling the bubble-like structures shown in Fig. 3 is expected to largely dominate the total CR energy injected by the CS. We also observe that the magnetic field topology might locally redirect outflowing streams of CRs back into the disk. However, due to the Eulerian nature of RAMSES, this effect remains difficult to quantify.

3.2. How much does the central source contribute to the local cosmic-ray population?

In Fig. 5, we evaluate the relative contribution of a single CS pulse to the total amount of CR energy measured in the solar neighborhood. Because the total energy E_F required to form the FBs is not well constrained and is thought to range from 10⁵⁶ to 10⁵⁷ erg, we rescaled the amount of energy injected from the CS in post-processing to explore this whole interval. The resulting range of values is represented as the colored surfaces. The dotted (dashed) line corresponds to the case where E_F = 10⁵⁶ erg (10⁵⁷ erg) is injected from the CS. The orange surface shows the estimated age τ_F of the FBs, which spans from 1 to 6 Myr (Guo & Mathews 2012; Zhang & Guo 2020). For each of the two energy bins, the values obtained in the case where the CRs propagate isotropically is shown by the shaded surface.

Fig. 5. Proportion of the total CR energy measured in the solar neighborhood attributed to the CS source as a function of time and for each energy bin. As current estimates of the energy required to form the FBs span between 1056 and 1057 erg, this quantity is left as a free parameter, the obtained range of values being represented by the colored surfaces. For each of them, the dotted (dashed) line represents the scenario in which 1056 erg (1057 erg) is required to form the FBs. The estimated age of the FBs (1–6 Myr) (Guo & Mathews 2012; Zhang & Guo 2020) is represented by the orange surface.

In the lower-energy bin, the contribution of the CS is strictly negligible at any time in the τ_F estimate range. This contribution might eventually increase up to a few percent in the isotropic case after 15 Myr, much beyond the larger estimate of τ_F. For the higher energy bin, the contribution of the CS reaches up to a few percent within the first 6 Myr in the anisotropic case, and up to 70% in the isotropic case. We conclude that the current estimates of the age of the FBs and of the total CR energy injected from the CS imply that the CS is not able contribute significantly to the observed total CR energy budget. We note that this results from the anisotropic CR diffusion, and that considering an isotropic diffusion significantly changes the result. Taking the source intermittence into account does not change these conclusions.

We emphasize that our simulation with isotropic diffusion provides results consistent with previous theoretical work (Jaupart et al. 2018), which are thus likely to overestimate the CS contribution in the solar neighborhood by up to several orders of magnitude.

3.3. Impact on the slope of the spectrum

To assess the effect of the central source on the CR flux in the solar neighborhood, we computed the slope of the spectrum using the energy density populating our two energy bins within all the cells contained in the solar neighborhood. The configuration is illustrated in Fig. 6.

Fig. 6. Method used to measure the spectrum slope. The central SMBH (•) particle and the SNe (⋆) injected CRs with an initial spectrum of slope α0. While propagating along the magnetic field lines, the slope of this spectrum can exhibit some modifications due to the energy dependence of the diffusion coefficient, and to a fraction of the CRs escaping the disk. The slope of the spectrum in the solar neighborhood αS(t) (later simply denoted α) is then measured as a function of time for both the SNe and the CS populations by comparing the relative amounts of energy contained in each energy bin, averaged over all solar neighborhood gas cells.

The slope α of the CR spectrum in the solar neighborhood (⊙) is then obtained following:

$$\begin{array}{c}\hfill \alpha (t)=\frac{1}{log{E}_2/E_1}log\left(\frac{{\sum }_{j\in \odot }{\rho }_{E_1,j}(t)·V_j}{{\sum }_{j\in \odot }{\rho }_{E_2,j}(t)·V_j}\right).\end{array}$$(9)

Here, the logarithmic term is the ratio between the CR energy contribution from both SNe to that from the CS contained in the solar neighborhood for the two energy bins. This fraction is obtained by summing the CR energy contained in each cell of the solar neighborhood, here labeled with the index j. The evolution of the slope is presented in Fig. 7. Here, the exact quantity of energy injected by the CS is left as a free parameter, which can vary between 10⁵⁶ and 10⁵⁷ erg. The colored surfaces represent the range of values covering this whole interval, while the colored lines represent the median values. The left panel shows the influence of the propagation mode on the evolution of the slope for the emission of a single pulse at t₀. The gray dotted line presents the injection spectrum α₀ = 2.4. The purple curve presents the slope of the CR spectrum when only accounting for the two SN groups. This slope is relatively constant, with a mean value of $\overline{\alpha }=2.6$, which is consistent with measured values of the average spectrum in the local Galaxy; that is, in the ISM within 1 kpc distance from the Sun (Neronov et al. 2017). The green line shows the slope of the total CR spectrum, i.e. which contains both the contribution of the SNe and the CS, for the isotropic propagation mode, while the blue line indicates the same quantity but the anisotropic propagation mode.

Fig. 7. Evolution of the spectrum slope under various conditions. Because the exact quantity of energy injected at each pulse is only weakly constrained, this quantity is left as a free parameter, which can vary between 1056 and 1057 erg. The corresponding range of values obtained are represented by the colored surfaces, while the dotted (dashed) line represents the case where a total CR energy of 1056 erg (1057 erg) is injected from the CS at t = 0. The left panel shows the influence of the propagation mode (i.e. isotropic or anisotropic). The right panel shows the influence of the source intermittence. The results of our two simulations with intermittent source, i.e. with multiple injections separated in time by period Δt are presented for Δt = 1 Myr and Δt = 3 Myr. The simulation with single anisotropic injections is also represented by the solid line for comparison.

The propagation mode is found to have a critical impact on the ulterior shape of the spectrum. In the isotropic case, the spectrum becomes significantly harder after ∼10 Myr. This hardening effect is due to the fact that the CRs belonging to the higher energy group E₂ diffuse faster than the E₁ group because of their larger diffusion coefficient. They thus populate the solar neighborhood more rapidly, leading to a flattening of the spectrum. In the anisotropic case, the propagation is more sensitive to the complex morphology of the Galaxy, because the CRs propagate preferentially along magnetic field lines, which follows its spiral shape (see Fig. 1). Direct effect observed here is weaker impact on the slope of the spectrum, which is likely to be undetectable for the first 6–8 Myr after the energy injection.

We also model the eventuality that the central source of CRs is intermittent, and present the result of our two associated simulations, with injection periods of 1 Myr and 3 Myr. Results are presented in the right panel of Fig. 7. Remarkably, the latency period between the energy injection and the time at which the spectrum in the solar neighborhood starts being modified by the contribution of the CS CRs is found to be independent of the source period, and roughly equal to 6 Myr. At later times, the effect on the slope increases with the frequency of injection from the central source. About 15 Myr after the first injection, modification of the slope is found to span from ∼7% for the single-pulse simulation, to 21% for the intermittent run with Δt_p = 1 Myr. Considering that current measurements of the CR spectrum slope have a precision of a few percent (Neronov et al. 2017), the contribution of the CS in the run with intermittent injection is unlikely to be measurable less than 10 Myr after the first injection. As most estimations of the age of the FBs ranges between 1 and 6 Myr, we can nevertheless conclude that a significant contribution of CRs emitted from the Galactic center to the overall CR spectrum measured in the solar neighborhood is unlikely to be observed by current surveys. These conclusions could however be modified observational evidence of an older outburst were recovered; that is, with an estimated age of more than 10 Myr.

4. Conclusion

We investigated whether or not CR outbursts from the Galactic center are likely to contribute significantly to the CR spectrum observed at the solar neighborhood. By comparing simulations run with isotropic and anisotropic diffusion, we conclude that the diffusion mode strongly influences the extent to which the CS will contribute to the local CR flux. When taking into account the preferential diffusion of CRs along magnetic field lines, CRs are found to require a much longer timescale to reach the solar neighborhood. By comparing the typical timescale required for the CRs to significantly impact the local CR flux slope with the estimated age of the Fermi-bubble in the jet-driven scenario, we conclude that the period of strong activity from the Galactic center suspected to explain the observations made by Fermi and eRosita telescopes is unlikely to affect the CR flux up to the TeV regime measured at present. Nonetheless, these conclusions do not exclude the possibility that hypothetical active phases of the Galactic center at older times contributed to the local CR spectrum.

We emphasize that both our CS injection scheme, and our propagation model remain simplified. We do not include some internal properties of the source, such as its possible anisotropy, stochasticity, or continuity. We also do not include production and acceleration of CRs within and at the edge of the outflowing bubbles, or CR reacceleration, streaming, energy loss, or secondary CR production. Further research work should focus on evaluating whether acceleration of CRs within the shock waves observed by eRosita is likely to challenge our conclusions. Treatment of CRs with a larger number of energy bins and energy losses might also help to constrain the possible effects of such an event on the local CRs flux more precisely. Finally, more observational constraints might help improving the modeling of the CS, especially regarding its intrinsic injection spectrum.

Data availability

Movie associated to Fig. 2 is available at https://www.aanda.org

Movie

Movie 1 associated with Fig. 2 (CR_total) Access here

Acknowledgments

We thank Alexandre Marcowith, Joakim Rosdahl, Yohan Dubois, Etienne Jaupart, Volker Heesen, Marcus Brüggen, Franco Vazza and Maria Werhahn for insightful discussions. We thank Eric Emsellem who provided us with the initial conditions for the Milky Way. This work was performed thanks to HPC resources provided by the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon. We also thank Florent Renaud for providing the open-source visualization tool RDRAMSES.

References

Appendix A: Energy conservation and numerical convergence

While running our first simulations, we found nonphysical production and destruction of CR energy in our simulated box, which we associated with incorrect estimates of CR flux at interfaces between cells. Visual artifacts known as chessboard patterns are usually visible under these conditions and are shown in Fig. A.1. The origin of these patterns is discussed in Sharma & Hammett (2007).

Fig. A.1. CR energy density projection for one of our preliminary, low-resolution runs. Several artifacts of diffusion produced by injection from SNe are visible. A typical feature of such events is a spike in CRs energy propagating perpendicularly to the magnetic field lines, and endowed with a chessboard pattern.

These sources of nonconservation, known to dramatically impact the reproducibility of our results, can be efficiently reduced by the use of flux limiters. Energy conservation with and without the use of slope limiter is presented in Fig. A.2. For these tests, we theoretically injected a pulse of 10⁵⁶ erg in both energy bins. Maroon and orange curves represent the energy effectively injected into the 100 GeV and the 10 TeV groups, respectively. As visible on the left panel, numerical artifacts lead to dramatic energy production, which multiplies the amount of energy theoretically injected by a factor of nearly five. After activating the slope limiter, we obtained the results presented in the right panel, which fulfill energy conservation over our tests run by nearly 90 % in both energy groups for the anisotropic runs. We verified that the subsequent energy loss, likely associated with nonconservation due to Dirichlet boundary conditions at coarse-to-fine level interface (Commerçon et al. 2014) only weakly affect our results.

Fig. A.2. Comparison of energy conservation between test runs without (top panel) and with (bottom panel) flux limiter. Test runs were performed by injecting a pulse of energy 1056 erg in both energy bins at t = 0 Myr. The total CR energy contained in the simulated box is measured at various times for each of the two energy bins.

We verified numerical convergence by running several simulations with the same physical parameters and different numerical resolutions, namely 59 pc, 29 pc, and 15 pc. In Fig. A.3, we present the total amount of energy measured in the solar neighborhood for our three single pulse simulations, with the above resolutions. While not perfectly matching, runs with resolutions of 29 pc and 15 pc seem to come to reasonable agreement. For this reason, we derived most of our results in this paper from our 29 pc resolution runs, which allowed us to run more simulations in the limited time we had for this project. We emphasize that the physical properties of our simulated Galaxy cannot be identical from one run to another. The stochasticity of star formation, as well as the resolution-dependent effect of SNe on the gas dynamics lead to star formation rates that can differ by up to ± 0.5 M_⊙ yr⁻¹ from one resolution to another. Consequently, the morphology of the Galaxy, as well as the topology of its magnetic field vary from one run to another, and can affect CR kinematics. This variability is likely to justify the lack of perfect convergence between our 29 pc and 15 pc resolution runs. Also, the amount of energy measured in the solar neighborhood in the lowest energy bin (100 GeV) is only a very small fraction of the injected energy. The dependence of the simulated Galaxy and its magnetic field on morphological properties is thus likely to be much larger than for the highest energy bin.

Fig. A.3. Numerical convergence test. A single pulse of ECS = 1056 erg is injected from the central source and the total amount of CR energy contained in the solar neighborhood cells is followed over time.

All Tables

All Figures

	Fig. 1. Initial conditions for CR injection. The left panel shows the hydrogen density, while the right panel shows the magnetic field amplitude overlaid with its associated field lines.
In the text

Fig. 2. Evolution of the several CR groups at four different times after injection of a pulse from the CS. The first and third columns present the energy density maps of the CR groups from SNe, at 100 GeV and 10 TeV respectively, while the second and fourth rows presents the same quantity for the CR groups associated with the CS. Each maps are normalized by the mean energy density in the SNe group at t = 0. The gray annulus in the second column represents the solar neighborhood as defined in the main text. (Movie available online).

In the text

	Fig. 3. Edge-on view of the 10 TeV CR group energy density injected by the CS in our s15pc run, 12 Myr after injection. As visible, Bubble-like structures can be seen to naturally emerge in our model due to the topology of the magnetic field.
In the text

	Fig. 4. Evolution of foutflow as a function of time and for various parameters. The left panel presents the influence of the propagation model (i.e. either anisotropic or isotropic) for our 29 pc resolution, single-pulse run (s29pc). The right panel shows the influence of the intermittence of the source.
In the text

Fig. 5. Proportion of the total CR energy measured in the solar neighborhood attributed to the CS source as a function of time and for each energy bin. As current estimates of the energy required to form the FBs span between 1056 and 1057 erg, this quantity is left as a free parameter, the obtained range of values being represented by the colored surfaces. For each of them, the dotted (dashed) line represents the scenario in which 1056 erg (1057 erg) is required to form the FBs. The estimated age of the FBs (1–6 Myr) (Guo & Mathews 2012; Zhang & Guo 2020) is represented by the orange surface.

In the text

Fig. 6. Method used to measure the spectrum slope. The central SMBH (•) particle and the SNe (⋆) injected CRs with an initial spectrum of slope α0. While propagating along the magnetic field lines, the slope of this spectrum can exhibit some modifications due to the energy dependence of the diffusion coefficient, and to a fraction of the CRs escaping the disk. The slope of the spectrum in the solar neighborhood αS(t) (later simply denoted α) is then measured as a function of time for both the SNe and the CS populations by comparing the relative amounts of energy contained in each energy bin, averaged over all solar neighborhood gas cells.

In the text

Fig. 7. Evolution of the spectrum slope under various conditions. Because the exact quantity of energy injected at each pulse is only weakly constrained, this quantity is left as a free parameter, which can vary between 1056 and 1057 erg. The corresponding range of values obtained are represented by the colored surfaces, while the dotted (dashed) line represents the case where a total CR energy of 1056 erg (1057 erg) is injected from the CS at t = 0. The left panel shows the influence of the propagation mode (i.e. isotropic or anisotropic). The right panel shows the influence of the source intermittence. The results of our two simulations with intermittent source, i.e. with multiple injections separated in time by period Δt are presented for Δt = 1 Myr and Δt = 3 Myr. The simulation with single anisotropic injections is also represented by the solid line for comparison.

In the text

	Fig. A.1. CR energy density projection for one of our preliminary, low-resolution runs. Several artifacts of diffusion produced by injection from SNe are visible. A typical feature of such events is a spike in CRs energy propagating perpendicularly to the magnetic field lines, and endowed with a chessboard pattern.
In the text

	Fig. A.2. Comparison of energy conservation between test runs without (top panel) and with (bottom panel) flux limiter. Test runs were performed by injecting a pulse of energy 1056 erg in both energy bins at t = 0 Myr. The total CR energy contained in the simulated box is measured at various times for each of the two energy bins.
In the text

	Fig. A.3. Numerical convergence test. A single pulse of ECS = 1056 erg is injected from the central source and the total amount of CR energy contained in the solar neighborhood cells is followed over time.
In the text