Unstable cosmic ray nuclei constrain low-diffusion zones in the Galactic disc

ABSTRACT

1 INTRODUCTION

One of the cornerstones of the standard picture of cosmic ray (CR) transport (see e.g. Strong, Moskalenko & Ptuskin 2007; Grenier, Black & Strong 2015; Gabici et al. 2019 for a review) is related to the overabundances of certain, so-called secondary species with respect to Solar system values. Such a difference is commonly interpreted as resulting from spallation reactions of primary CRs with interstellar gas. The measurement of secondary-to-primary ratios like B/C therefore allows us to constrain the matter traversed during propagation, called the grammage. As the amount of grammage accumulated is about three orders of magnitude larger than the column density of the disc, CRs must propagate diffusively. In fact, flux ratios involving unstable secondaries such as ¹⁰Be allow us to infer that CRs are confined in a larger region called the Galactic halo of CRs, or halo for short, extending to kiloparsec distances above and below the formalized. All these interpretations can be formalized in a standard framework where CRs are produced in the Galactic disc and diffuse in the halo by scattering upon turbulent Galactic magnetic field, first proposed by Ginzburg & Syrovatskii (1964). The combination of B/C and ¹⁰Be then allows us to determine the two most important parameters of this model which are the diffusion coefficient κ and the halo height H (Donato, Maurin & Taillet 2002; Lipari 2022).

As magnetic turbulence controls both κ and H, it is essential for a comprehensive picture of Galactic CR propagation. Yet, the question of how magnetized turbulence is generated and distributed on large scales remains open. Turbulence can be generated both by supernova explosions and by inhomogeneities of CR density via the CR streaming instability (Kulsrud & Pearce 1969; Skilling 1975; Marcowith, van Marle & Plotnikov 2021). For Galactic propagation, the former is believed to be particularly important for CR particles with high rigidity above a few hundred GV while the latter might be more relevant for CRs below a few hundred GV (Blasi, Amato & Serpico 2012). This has been employed to explain the change in the rigidity dependence of κ at |${\mathcal {R}}\simeq 300$| GV as inferred from B/C data (Génolini et al. 2017). We note, however, that most of the analyses adopt the simplified assumption that κ is homogeneous, specifically that it is the same in the disc and in the halo (see, however, Tomassetti 2012; Evoli & Yan 2014). More recent studies suggest that the spatial dependence of κ might be more complicated if turbulence self-generated via the CR streaming instability is taken into account (Recchia, Blasi & Morlino 2016; Evoli et al. 2018b).

In fact, observations seem to indicate that the diffusion coefficient could be different between the halo and the disc, at least for certain regions in the disc. Gamma-ray observations of potential CR sources such as supernova remnants or pulsars have been interpreted to have regions with a suppressed isotropic diffusion coefficient in their surroundings. For example, HAWC (Abeysekara et al. 2017), Fermi-LAT (Di Mauro, Manconi & Donato 2019), and H.E.S.S. (Mitchell et al. 2022) observed haloes of gamma rays (referred to as gamma-ray haloes to be distinguished from the Galactic halo of CRs) at GeV and TeV energies around the pulsars Monogem and Geminga. These haloes originate from high-energy e⁺/e⁻ pairs accelerated by the pulsars up to energies of |$100\, \mathrm{TeV}$| (Abeysekara et al. 2017) and produce gamma-rays via inverse Compton scattering upon the cosmic microwave background and interstellar radiation fields. These observations provide some support for the hypothesis that pulsars are sources of primary positrons, capable of explaining the positron excess observed by AMS-02 (Aguilar et al. 2013) and PAMELA (Adriani et al. 2010). Interestingly, the radial profile of gamma-ray haloes cannot be explained by models where electrons diffuse isotropically from pulsars with a diffusion coefficient as inferred by fits to local CR data, but they require a diffusion coefficient lower by about two to three orders of magnitude (Abeysekara et al. 2017). Under the initial assumption that the diffusion coefficient is that low in the entire disc, electrons at TeV energies could not reach Earth, since the diffusion time would be larger than the energy loss time. However, many experiments observe these electrons (e.g. DAMPE Collaboration 2017; Adriani et al. 2018), hence suppression of the diffusion coefficient by that order can be excluded (Hooper & Linden 2018). In order to simultaneously explain the positron excess and the gamma-ray haloes a two-zone diffusion model was introduced (Fang et al. 2018; Tang & Piran 2019) with inhibited diffusion in a spherical region at least the size of the gamma-ray halo and the usual interstellar medium value outside. Within this two-zone approach, it is possible to explain the positron excess with a distribution of pulsars (Manconi, Di Mauro & Donato 2020). However, the origin of these low-diffusion regions remains uncertain (see, however, Evoli, Linden & Morlino 2018a).

Additionally, observations of gamma-rays from dense molecular clouds near supernova remnants can be interpreted with a suppressed diffusion coefficient in this area (Gabici et al. 2010). In this scenario, hadronic CRs accelerated by the supernova remnant escape into the interstellar medium and propagate to nearby dense molecular clouds, where they produce gamma-rays via pion production. The diffusion coefficient required to explain the intensity of gamma-rays is suppressed by around two orders of magnitude with respect to the Galactic average. This suppression can be self-generated via CR streaming instability (Malkov et al. 2013; Nava et al. 2016, 2019; Brahimi, Marcowith & Ptuskin 2020; Recchia et al. 2021) and could last for several |$100\, \mathrm{kyr}$| at GeV energies (Jacobs, Mertsch & Phan 2022). Furthermore, suppression of the diffusion coefficient at GeV energies has been inferred in dense molecular clouds and in stellar bubbles (Aharonian, Yang & de Oña Wilhelmi 2019; Yang et al. 2023). Hence, there is plenty of evidence of inhomogeneous diffusion and the impact of these low diffusion zones on Galactic CR propagation might be non-negligible.

Here, we will devise a simple semi-analytical model of CR propagation in a scenario where we have a suppressed averaged diffusion in the Galactic disc and investigate constraints of upcoming AMS-02 (Jiahui Wei 2022) and HELIX (Park et al. 2019) ¹⁰Be/⁹Be measurements. We then show how to calculate a large-scale averaged diffusion coefficient from smaller-scale features using the numerical solution of a stochastic differential equation (SDE). Finally, we explain the implications of our findings and conclude. Detailed calculations and additional results are presented in a number of appendices.

2 MODEL

We use a 1D transport equation to describe CR propagation in the Galaxy following: Ginzburg & Syrovatskii (1964) and Parker (1965). CRs are assumed to propagate within a magnetized halo of size 2H extending symmetrically above and below the Galactic disc. The halo is assumed to have a negligible gas density and a free escape boundary condition is imposed at the halo’s boundaries. The sources of CRs are considered to be located only within the disc of height h ≪ H. Similarly, the interstellar gas, which is the target of interactions inducing energy loss for CRs, is confined to this region. Since all observed and potential sources of suppressed diffusion are confined to the disc, we distinguish between two zones, the disc, and the halo in the following.

Analytically solving the transport equation for a diffusivity with randomly distributed regions of suppression, as shown in Fig. 1, is intractable. A numerical solution with the finite-difference technique that is commonly used in CR transport codes is also unfeasible: due to the large dynamical range in diffusivity, unrealistically small times-steps would be required for the algorithm to remain accurate. Here, we therefore resort to a coarse-graining approach where we adopt a spatially homogeneous diffusion coefficient |$\kappa _{\mathrm{disc}}({\mathcal {R}})$| in the disc which is related to an unsuppressed diffusion coefficient |$\kappa _{\mathrm{high}}({\mathcal {R}})$| by a scaling factor α,

Here, we have assumed the same rigidity-dependence for |$\kappa _{\mathrm{disc}}({\mathcal {R}})$| and |$\kappa _{\mathrm{high}}({\mathcal {R}})$|⁠. This can be motivated, e.g. for gamma-ray haloes around PWNe, by the fact that the rigidity-dependence in the low-diffusion zone has been found to be compatible with a Kolmogorov phenomenology that is commonly assumed for κ_high (Di Mauro et al. 2019). We note that this approach is also flexible enough to cover other potential sources of suppressed diffusion, like CR self-confinement, as long as the suppressed and unsuppressed diffusion coefficients have the same rigidity-dependence. We posit that this scenario deserves investigation on its own merits. For the scenario where α results from the presence of bubbles of suppressed diffusion, we provide some relations between the ratio (κ_low/κ_high), the filling fraction f of bubbles in the disc and the resulting α in Section 4.

In the halo, the absence of low diffusion zones implies κ_halo = κ_high. The rigidity dependence of the diffusion coefficient is modelled as a broken power law to account for spectral breaks observed in CR data at low (Vittino et al. 2019) and high rigidities (Evoli, Aloisio & Blasi 2019),

Some of the parameters have little effect on our results and so we consider values which have been found in other studies. Specifically, we fix s_l = 0.04, which has minimal impact on the results, as shown by Weinrich et al. (2020). Furthermore, the high energy break parameters can be determined by primary nuclei (Evoli et al. 2019): |$s_h=0.5, \ \delta _h=0.34, \ {\mathcal {R}}_h=312\, \mathrm{\mathrm{GV}}$|⁠. All other parameters remain as free parameters.

In the steady state, the transport equation for Galactic CRs reads (Berezinskii et al. 1990):

where ψ_j is the density differential in momentum of CRs of species j, connected to the phase-space density f_j by ψ_j = 4πp²f_j. The terms on the right-hand side describe advection, diffusion, (adiabatic) energy losses, decay with boosted decay time γτ_j and spallation of CRs, respectively. While we distinguish between diffusion at z < h and z > h, we neglect the spatial distribution of spallation, energy losses and the sources in the disc, since we expect the effect to be small and approximate them with a delta distribution δ(z). We make use of the cross-sections provided by Evoli et al. (2019, 2020). The advection speed v(z) = sgn(z)v_c is directed vertically out of the Galactic disc and left as a free parameter, giving dv/dz = 2v_cδ(z). The dominating energy loss process is ionization on neutral gas in the disc, which can again be assumed to be infinitesimal (Mannheim & Schlickeiser 1994; Evoli et al. 2019):

where we use the formula given by Mannheim & Schlickeiser (1994) and correct a typo as described in Jacobs et al. (2022).

There are three ways to produce CRs of species j, accounted for by the right-hand side of equation (3):

The first term describes injection by a population of sources, where the dominating injection mechanism is believed to be diffusive shock acceleration, which in the linear strong shock case results in a |${\mathcal {R}}^{-2}$| spectrum (Malkov & Drury 2001). To account for non-linear modification effects (e.g. Diesing & Caprioli 2021), here we assume |${\mathcal {R}}^{-2.3}$| for all primaries with the abundances given by Evoli et al. (2020). Unstable heavier elements can decay into species j, resulting in the second term of (5). Furthermore, heavier CRs can spallate into species j by interacting with the interstellar gas in the disc (the third term of equation 5). The uncertainty on those interaction cross-sections will be a source of systematic errors, which is discussed in Section 4. The disc height is |$100\, \mathrm{pc}$| and the hydrogen number density is assumed to be |$n_{\mathrm{gas}}=1\, \mathrm{cm^{-3}}$|⁠, similar to what was adopted in Ferrière (2001) and Phan et al. (2021).

We solve the diffusion equation separately in the disc and in the halo, demanding continuity of the CR density and flux at the interface and in the disc, in addition to enforcing free escape at the halo’s boundaries. We take into account the entire nuclear chain up to iron. At low energies, the fluxes in the Solar system are altered by solar modulation. To account for this, we use a force field model by Gleeson & Axford (1967) with a solar modulation potential of |$\psi =700\, \mathrm{MeV}$| that we obtained by a fit to AMS-02 Oxygen data (Aguilar et al. 2021). We use these fluxes to calculate the B/C ratio, where we make the same assumptions as the AMS-02 collaboration on isotopic composition, namely pure ¹²C as well as a combination of ¹⁰B and ¹¹B (Aguilar et al. 2016). These results are fitted to AMS-02 data (Aguilar et al. 2016) and recent Voyager data (Cummings et al. 2016). The latter have been measured outside the heliosphere and are therefore not subject to solar modulation.

It can be shown that for high energies the B/C ratio depends on an averaged diffusion ratio:

Since the averaged diffusion coefficient in the Galaxy cannot change significantly, as shown by Hooper & Linden (2018) and Profumo et al. (2018), the impact of a suppressed diffusion in the disc on B/C is rather small and is degenerate from a change of H/κ_halo. To break the degeneracy between H and κ_halo, |$\strut ^{10} \text{Be}$| data are needed, since they allow to determine the propagation time. Because a suppressed diffusion coefficient in the disc inhibits the escape of freshly produced |$\strut ^{10} \text{Be}$| into the halo, these measurements can also potentially allow to determine α.

For kinetic energies per nucleon below a few GeV, we use available data presented in Garcia-Munoz, Mason & Simpson (1977), Wiedenbeck & Greiner (1980), Garcia-Munoz, Simpson & Wefel (1981), Lukasiak et al. (1994), Connell (1998), Lukasiak (1999), Yanasak et al. (2001), and Nozzoli & Cernetti (2021), but there are no published data at energies above |$E_{k/n}\approx \, 1\, {\rm GeV/n}$|⁠. At this point, in addition to the data above – which can be found in the CRDB (Maurin et al. 2023) we will consider preliminary AMS-02 data (Jiahui Wei 2022). In addition, we will forecast the constraining power of upcoming HELIX data (Park et al. 2019) by fitting to mock data. For generating these, we have adopted the model parameters as determined by fitting to all other data sets. The energy binning is chosen in accordance with Park et al. (2019) and for each energy bin, we have drawn a sample from a normal distribution centred on the model and with errors as given in Park et al. (2019).

In order to find the best-fitting parameters we adopt a Gaussian log-likelihood,

where d runs over the observables B/C and ¹⁰Be/⁹Be and i runs over rigidities |${\mathcal {R}}_i$| for B/C AMS-02 data or kinetic energies per nucleon E_k/n,i else. The errors are the statistical and systematic uncertainties summed in quadrature. We have refrained from taking into consideration the possibility of correlated errors (see e.g. Derome et al. 2019) as no official information on this has been provided by the collaborations; we also consider this to be of lesser importance than, for instance, in searches for excesses well-localized in rigidity (e.g. Boudaud et al. 2020).

Apart from the best-fitting values, we want to quantify the allowed parameter intervals. Therefore, we have performed a Monte Carlo Markov Chain (MCMC) study which also provides robustness in the identification of the best-fitting parameter values. We make use of the emcee package (Foreman-Mackey et al. 2013) and use uniform priors as listed in Table 1. For model comparison, we use a likelihood ratio test, modified to take into account the fact that the parameter α is bounded (Cowan et al. 2011).

In the following we will investigate the null hypothesis, i.e. that there is no suppression, and the full model for different data sets. Each of them contains the AMS-02 and Voyager B/C data, together with all published ¹⁰Be/⁹Be data. The first setup (‘w/o prelim’.) reflects our best knowledge to date and only accounts for existing data. The second setup (‘w/ prelim’.) additionally takes into account preliminary AMS-02 data on ¹⁰Be/⁹Be (Jiahui Wei 2022). We expect the latter to be able to break the degeneracy of κ and H. Additionally, we will predict the effects of upcoming HELIX data as described above in the setup called ‘w/o forecast’.

3 RESULTS

The corner plot of the MCMC for the full model with preliminary data (‘full w/ prelim’.) is shown in Fig. 2. The lower left triangle displays the 2D marginalized posteriors of all parameter combinations. The dark blue (bright blue) regions represent the |$68\, ~{{\ \rm per\ cent}}$| (⁠|$95\, {{\ \rm per\ cent}}$|⁠) quantiles. On the diagonal the 1D marginalized posteriors are shown with the |$16\, {{\ \rm per\ cent}}$| and |$84\, {{\ \rm per\ cent}}$| quantiles indicated by blue-dashed lines. Both in the 1D and 2D posteriors the best-fitting values are marked with solid black lines, see also Table 2, where we give the χ² of the setups as well. The first column, which shows the posterior as a function of α clearly indicates that the null hypothesis is firmly ruled out. An anticorrelation between κ_halo and α exists, see equation (6), where a smaller α partially offsets the effects of a larger κ_halo. The degeneracy of κ_halo and halo height H is largely broken and the halo height can be constrained to |$H=8.2^{+1.1}_{-0.9}\, \mathrm{kpc}$|⁠, however, some correlation can still be observed. Another degeneracy exists between δ and v_c. This is a well-documented effect observed in advection–diffusion models, see e.g. fig. 2 of Putze, Derome & Maurin (2010). It is due to the fact that in order to match the secondary-to-primary ratio over a limited rigidity range, a larger advection speed v_c can be compensated by a larger δ, that is a stronger rigidity dependence of κ. The best-fitting values are generally in good agreement with the median of the distribution, but for the low energy break rigidity volume effects in the case δ = δ_l lead to an extended non-Gaussian tail. The corner plots of the other setups and the null hypothesis are given and explained in appendix  B and the spectra of B/C and ¹⁰Be/⁹Be for the ‘full w/ prelim’ model are shown in Fig. 10.

Figure 1.

Schematic sketch of low diffusion zones in the disc and the effect of course-graining. Zones of low diffusivity κ_low in the Galactic disc of height 2h are surrounded by a zone with high diffusivity κ_high, which extends into the Galactic halo of height 2H. Coarse-graining leads to an averaged diffusion zone in the disc κ_disc while in the halo κ_halo = κ_high.

Open in new tabDownload slide

Figure 2.

Corner plot of the MCMC scan of the full model including preliminary data, ‘w/ prelim’. The 2D marginalized posteriors are shown on the lower triangle. The dark/bright blue areas corresponds to the |$68\, {{\ \rm per\ cent}}$|/|$95\, {{\ \rm per\ cent}}$| quantiles. Similarly, the principal diagonal displays the 1D marginalized posterior with the |$16\, {{\ \rm per\ cent}}$| and |$84\, {{\ \rm per\ cent}}$| quantiles marked as blue dashed lines. The black lines indicate the best-fitting values as listed in Table 2.

Open in new tabDownload slide

Fig. 3 shows the best-fitting ¹⁰Be/⁹Be spectra as a function of kinetic energy per nucleon for the full model with and without preliminary data (blue-dashed and red-solid lines) and the null hypothesis with preliminary data (green dotted–dashed lines). At low energies the ratio is determined by an interplay between diffusion and advection. Since Be is purely secondary, it is exclusively produced in the disc. A faster transport in the disc by advection or diffusion leads to freshly spallated ¹⁰Be being transported out of the disc into the halo, where it will decay. Since the Solar system is located in the disc, this results in a lower observed ratio. Therefore, the Be ratio at low energies is especially sensitive to the suppression of the diffusion coefficient, which has little effect in the propagation afterwards. The high energy limit, where ¹⁰Be can be considered as stable due to time dilatation, is given by the ratio of the production cross-sections (Maurin, Ferronato Bueno & Derome 2022) and is dominated by their uncertainties. The energy at which the transition happens is determined by the energy as which the decay time of particles becomes larger than the escape time. Furthermore, the steepness of the rise is determined by the escape time. Together the above effects explain, why our setup can constrain α, H, and κ_halo, as shown in Fig. 2. The decrease in α together with the increase in H of the full model ‘w/ preliminary’ compared to the null hypothesis results in a flatter transition to the high energy limit as preferred by the AMS-02 data. Without the AMS-02 data the transition is unconstrained, resulting in both α, κ, and H being essentially unconstrained, as visible in Fig. B1.

Figure 3.

Best-fitting spectra for the Beryllium ratio to the setups with only existing data (red solid), preliminary AMS-02 data (blue dashed), and the corresponding null hypothesis α = 1 (green dotted-dashed). The null hypothesis predicts a lower Be ratio at low energies, and a transition to the high energy limit earlier, since the best-fitting halo height is smaller (for more details see text). The lower panel is the pull plot.

Open in new tabDownload slide

The 1D marginal posterior as a function of α for all three data sets is shown in Fig. 4. Here, we have used a Gaussian kernel density estimation with a width of 0.1. The posterior for existing data is indicated by the red solid line, the one including preliminary data by the blue-dashed line and for the forecast by the dotted purple line.

Figure 4.

Marginalized posterior of the scaling factor α. The setup with existing data only (red solid) shows no preference for a suppression. The setup with preliminary AMS-02 data, ‘w/ prelim’., allows constraining the scaling factor to around |$0.20^{+0.10}_{-0.06}$|⁠, with a model preference of |$3.5\, \mathrm{\sigma }$|⁠. Upcoming HELIX data will be able to constrain the ‘w/ prelim’. best-fitting model to |$6.8\, \mathrm{\sigma }$|⁠. Due different effects the increase in the significance difficult to estimate from the posterior (see text for details).

Open in new tabDownload slide

We give the median and estimate its error by its distance to the 16 and |$84~{{\ \rm per\ cent}}$| quantiles in Table 2. For the setup ‘w/o prelim’. the posterior is flat and extends to α = 1. The observed decline towards α = 1 is a result of the chosen prior α < 1 and the convolution with the Gaussian kernel. Values larger than 1 are unconstrained by our model, since the residence time in the disc would become negligible. Thus without preliminary AMS-02 data, we do not find a preference for a suppressed diffusion coefficient. It comes as no surprise that H and κ are also still degenerate. Including preliminary AMS-02 data, allows to constrain the median to |$\alpha =0.20^{+0.10}_{-0.06}$|⁠. In order to compare models we use a likelihood ratio test and account for the null hypothesis being at the boundary of our parameter space (Cowan et al. 2011). This results in a |$3.5\, \mathrm{\sigma }$| preference for a model with suppressed diffusion. Should the trend continue with upcoming HELIX data, we forecast a detection by |$6.8\, \mathrm{\sigma }$|⁠. This increase of significance is in contrast to the naive expectation gleaned from looking at the posteriors. However, there are several reasons why we deem the modified test statistics to be more reliable. First, the MCMC scan lacks statistics so far in the tails of the posterior and secondly, we expect volume effects to contribute in the tail as well. Hence, the median and its errors are well estimated by the posterior, but the model significance has to be determined using test statistics.

4 DISCUSSION

We have seen in the previous section that preliminary data on Beryllium isotopes from the AMS-02 experiment favour a smaller diffusion coefficient in the disc than in the halo. As far as the interpretation is concerned, our analysis does not allow distinguishing two scenarios for the origin of this suppressed diffusion: either the diffusion coefficient is suppressed (more or less) everywhere in the disc, for instance due to different turbulence generation or damping than in the halo; or the diffusion coefficient is strongly suppressed in localized regions of a certain filling fraction.

For the scenario where the suppression is due to bubbles of low diffusivity, we have estimated α for a given ratio (κ_low/κ_high) by numerically simulating the transport of a large number of CRs with SDEs. To this end, we have modelled the disc as a periodic medium with a diffusion coefficient κ_high and spherical regions of suppressed diffusion coefficient κ_low, arranged on a cubic grid. For the ratio of (κ_low/κ_high), we have adopted the values 10⁻³ or 10⁻². The sources of CRs have been assumed to be uniformly distributed in the medium. We have computed the coarse-grained diffusion coefficient as κ_disc from the mean-square displacement of particles,

We have checked that the exact configuration (e.g. simple cubic, body-centred cubic or face-centred cubic) does not matter for the asymptotic value of κ_disc, but only the filling fraction.

In Fig. 5, we show the suppression α achieved as a function of the volume-filling fraction f of bubbles of low diffusion in the disc. The value |$\alpha = 0.2^{+0.10}_{-0.06}$|⁠, found in the CR parameter study can be realised for a filling fraction of |$\sim 66~{{\ \rm per\ cent}}$|⁠. We also compare to the arithmetic mean achieved for a given filling fraction.

Figure 5.

Effective diffusion coefficient κ_disc from the SDE computation as a function of the volume filling fraction f. The solid blue (dashed orange) line shows κ_disc in unit of κ_high, assuming a ratio κ_low/κ_high of 10⁻³ (10⁻²). The sources of CRs are assumed to be distributed uniformly in the disc. Note that the lines end at the maximum possible filling fraction of |$\sim 74~{{\ \rm per\ cent}}$| for spherical inclusions. For comparison, the green dotted (red dotted–dashed) line shows the arithmetic mean (κ_low/κ_high)f + (1 − f) assuming κ_low/κ_high of 10⁻² (10⁻³). The grey line and band highlight the suppression found from the CR fit, |$\alpha = 0.2^{+0.10}_{-0.06}$|⁠.

Open in new tabDownload slide

In the literature, the value typically estimated for the filling fraction f is lower,

where N_bubbles is the number of bubbles, |$\mathcal {R}_{\text{bubble}}$| their rate, R_bubble, V_bubble, and τ_bubble are the radius, volume, and lifetime of one bubble, respecticely, and R_disc, V_disc denote the radius and volume of the Galactic disc, respectively. In fact, a filling fraction of |$\sim 60~{{\ \rm per\ cent}}$|⁠, which would be in agreement with the determined value of α, would require a bubble radius of |$150 \, \text{pc}$|⁠.

We note that such a high filling fraction does not contradict H.E.S.S. observations (Hooper et al. 2017), as long as the sources of electrons are distributed in both high and low diffusion zones. Hence, zones of low diffusivity in our Galaxy might be more common and extended than previously assumed (Profumo et al. 2018; Jóhannesson, Porter & Moskalenko 2019; Manconi et al. 2020). Interestingly, recent studies of TeV gamma ray haloes around pulsars allow for similar sizes of |$100\, \mathrm{pc}$| of the low diffusion region (Schroer, Evoli & Blasi 2023). Additionally, simulations of galaxy evolution require high filling fractions in order to suppress clump formation in star-forming regions of the interstellar medium (Semenov, Kravtsov & Caprioli 2021). The setup we used here can be seen as a simplification of a more complex model, including magnetic field simulations, which in the future will supersede the simplistic approach used here (Ruszkowski & Pfrommer 2023).

As discussed in Section 3, the reason that the data prefer a suppressed diffusion coefficient in the disc is that the transition from the low to the high energy limit in the observed |$\strut ^{10} \mathrm{Be}/\strut ^{9} \mathrm{Be}$| is relatively flat. Since the Evoli et al. (2018b) cross-sections used here result in a relatively high |$\strut ^{10} \mathrm{Be}/\strut ^{9} \mathrm{Be}$|⁠, other cross-section parametrizations will likely lead to a smaller filling fraction, resulting also in a smaller bubble size.

It is worth mentioning that the filling fraction of low diffusion zones fitted here differs from previous analyses (Jóhannesson et al. 2019). The main reason is the inclusion of Be data in our analysis, which allows us to constrain the diffusion coefficient in the disc and estimate the filling fraction of low diffusion zones. Furthermore, we adopt a different coarse-graining procedure to relate the diffusion coefficient in the disc κ_disc with the filling fraction of low diffusion zones f. Previous analyses (e.g. Jóhannesson et al. 2019) assume that particles start in the centre of a low diffusion zone, which is surrounded by a high-diffusion zone and then use the characteristic diffusion times to calculate an effective diffusion coefficient in the disc used for CR propagation afterwards. This procedure could lead to κ_disc being a few times suppressed with respect to κ_halo even for a filling fraction of a few per cent only. However, this approach ignores the fact that particles, after leaving the high diffusion zone and propagating in the halo, could come back and cross the disc many times before finally escaping the Galaxy. When particles reenter the disc, they have a lower than one probability to travel in other low diffusion zones and this could lead to κ_disc larger than previously expected for the same filling fraction (see also Fig. 5 and the discussion above).

As a remark, we would like to stress also that our results have been obtained using a linear model where the diffusion coefficients are determined by a fit to CR data. There exists also non-linear Galactic CR propagation models in which the diffusion coefficient is determined by turbulence generated via CRs themselves. Many of these models, however, predict κ_disc > κ_halo at least for certain regions in the halo which seems to contrast with the results in our linear model. For example, Recchia et al. (2016) found that balancing the growth rate of CR-induced turbulence and the non-linear Landau damping rate (see e.g. Ptuskin & Zirakashvili 2003 for more details) results in the diffusion coefficient being highest in the disc and decreasing linearly with |z| up to |z| = H. Recently, Evoli et al. (2018b) revisited this type of non-linear model without imposing the free-escape boundary condition at |z| = H to study the origin of the halo. The authors took into account also other types of turbulence-damping mechanisms and, more importantly, the injection of turbulence due to supernova explosions. In this case, the diffusion coefficient is also relatively high in the disc and reaches a minimum around |z| ≃ 1 kpc but increases again at larger |z|. We note, however, that these models do not consider the geometry of the Galactic magnetic field which could be important in modelling CR-induced turbulence.¹ Given that a suppressed diffusion coefficient in the disc is preferred by preliminary AMS-02 data, it would be interesting to revisit these non-linear models with a more precise geometry of the Galactic magnetic field and confront them with unstable secondary CR data.

5 CONCLUSION

We have presented an extended model of Galactic CR transport, taking into account the presence of regions of suppressed diffusivity in the Galactic disc. To this end, we have solved the CR transport equation assuming two diffusion coefficients, one in the disc (κ_disc) and one in the halo (κ_halo). Fitting this model to existing B/C and ¹⁰Be/⁹Be data by AMS-02, Voyager and several other experiments, we have found no preference for a suppressed diffusion coefficient in the disc. This has to be expected since these data sets do not even break the degeneracy between halo height and diffusion coefficient. Adding preliminary AMS-02 data on ¹⁰Be/⁹Be reveals a |$3.5\, \mathrm{\sigma }$| preference for a suppressed diffusion coefficient in the disc with a best-fitting value of |$\alpha = \kappa _{\rm disc}/\kappa _{\rm halo} = 0.20^{+0.10}_{-0.06}$|⁠. Additionally, these data constrain the halo height to |$H=8.9^{+1.3}_{-1.0}\, \mathrm{kpc}$| and the diffusion coefficient in the halo to |$\kappa _{\rm halo} = 0.28^{+0.06}_{-0.05}\, \mathrm{pc^2\,yr^{-1}}$|⁠. We predict that, if the preliminary findings by AMS-02 are confirmed, upcoming HELIX data will be able to detect a suppressed diffusion coefficient in the disc with |$6.8\, \mathrm{\sigma }$| significance. The inclusion of additional elemental ratios might additionally increase the significance of our model.

The suppressed diffusion coefficient in the disc found in this work can be interpreted as the collective effect of low diffusion zones present inside the Galactic disc. In fact, these low diffusion zones have been inferred from gamma-ray observations of regions surrounding supernova remnants or PWNe (e.g. Gabici et al. 2010; Abeysekara et al. 2017) and also inside dense molecular clouds (e.g. Yang et al. 2023). Having found the scaling factor α = κ_disc/κ_halo, we proceed to constrain the filling fraction of low diffusion zones f. The relation between κ_disc and f is found using simulations of CR transport in an inhomogeneous diffusion medium. Within a given medium with low diffusion zones of filling fraction f, the effective diffusion coefficient of CRs is as presented in Fig. 5. This allows us to constrain the filling fraction to be around 66 per cent, i.e. about 2/3 of the disc should be filled with low diffusion zones. Our results suggest that low diffusion zones surrounding potential CR sources such as supernova remnants and pulsars might be more ubiquitous than previously expected.

ACKNOWLEDGEMENTS

The authors would like to thank Jakob Böttcher for advice on likelihood ratio tests in the presence of parameter bounds and acknowledge helpful conversations with Dieter Breitschwerdt, Stefano Gabici, Yoann Génolini and Pierre Salati. This project was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 490751943.

DATA AVAILABILITY

The best-fitting spectra and mock data will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: TWO ZONE DIFFUSION MODEL

The transport equation for CRs in steady state reads:

In a two zone model, this equation is solved using the following steps:

• The transport equation is solved in the halo (h < |z| < H) with the boundary condition ψ(H) = 0.

• We solve the equation in the disc and demand that both the density ψ and its spatial flux be conserved at |z| = h.

• We integrate the transport equation over an infinite volume around z = 0.

In the following we will denote all properties in the halo with a o for outer and in the disc with i for inner. Starting at (i), the transport equation at z > 0 reads:

where we have neglected decay into species j. Now we make the Ansatz ψ_j = Cexp (λz) and introduce the decay rate Γ_j = (γτ_j)⁻¹. Plugged into equation (A2) this results in:

Using the boundary condition ψ(z = H) = 0 gives:

Now we can move to step (ii). The solution for equation (A3) can be written in a manner analogous to equation (A2) as

Now we match the solutions at z = h:

resulting in:

If we use the conservation of the flux we get:

which can be solved for B_o to give:

If we now eliminate B_o and express A_i in terms of B_i:

In order to find B_i, we need to integrate equation (A1) over an infinite volume around z = 0. Denoting ψ_j(z = 0) ≡ ψ_{j, 0}, this gives:

where the derivative can be computed for the solution for z > 0 and the advection velocity is changing sign at the disc. Then the above equation reads:

Now we use this equation to calculate B_i. The differential equation that has to be solved is:

This is done numerically using Scipy’s solve|$\_$|ivp function (Virtanen et al. 2020).

APPENDIX B: ADDITIONAL PLOTS

The corner plots of the MCMC scans for the ‘Full w/o prelim’. setup is shown in Fig. B1. Due to the lack of Be data in the GeV region, the degeneracy between H and κ_o coming from B/C data cannot be broken. This is visible in the second column, third row. As a consequence, the suppression factor α is unconstrained as well. Due to the posterior being non-Gaussian the marginalized posteriors are subject to volume effects. Including both preliminary AMS-02 and a forecast for HELIX data, ‘w/ forecast’, results in the corner plot shown in fig. B2. Here, the degeneracy between H and κ_o is broken and α can be constrained, as explained in Section 4.

Figure B1.

Same as Fig. 2, but for the full model excluding preliminary data, called ‘w/o prelim’. The extension of the posterior to α = 1 indicates that the data cannot constrain α. Furthermore, the ratio of H/κ_o is largely unconstrained.

Open in new tabDownload slide

Figure B2.

Same as Fig. 2, but for the full model including a forecast of upcoming HELIX data.

Open in new tabDownload slide

For the null hypothesis with preliminary data, ‘Null w/ prelim’., the corner plot is given in Fig. B3. Here, the degeneracy between H and κ_o is broken again. Compared to the full model, smaller values for both parameters are preferred. In the case δ ≈ δ_l the lower break rigidity R_l becomes unconstrained.

Figure B3.

Same as Fig. 2, but for the null model including preliminary AMS-02 data.

Open in new tabDownload slide

Removing preliminary AMS-02 data results in the corner plots shown in Fig. B4. The degeneracy between H and κ_o still exist, but is less pronounced than in the full model. The halo height is similar to the value found by Wei (2022). Again the same volume effect as for the model with preliminary data is visible.

Figure B4.

Same as Fig. 2, but for null model with existing data only.

Open in new tabDownload slide

Figure B5.

Best-fitting spectra of B/C (top panel) and |$\strut ^{10}\text{Be}/\strut ^{9}\text{Be}$| (bottom panel) including statistical |$68\, {{\ \rm per\ cent}}/95\, {{\ \rm per\ cent}}$| quantiles for the full model ‘w/ prelim’. The red line indicates the best fit. The lower panels are pull plots.

Open in new tabDownload slide