Galactic Diffuse Emission from Radio to Ultra-high-energy γ-Rays in Light of Up-to-date Cosmic-Ray Measurements

Cosmic rays (CRs) travel throughout the Galaxy, leaving traces from radio to ultra-high-energy γ-rays due to interactions with the interstellar gas, radiation field, and magnetic field. Therefore, it is necessary to utilize multiwavelength investigations on the Galactic diffuse emission to shed light on the physics of CR production and propagation. In this work, we present a spatially dependent propagation scenario, taking account of a local source contribution, while making allowances for an additional CR component freshly accelerated near their sources. In this picture, after reproducing the particle measurements at the solar system, we calculated the intensity and compared the spectral energy distribution to observations from Fermi-LAT and LHAASO-KM2A in the γ-ray band, and from WMAP and Planck among other radio surveys at lower energies. Multiband data considered in conjunction, the former comparison exhibits sufficiently good consistency in favor of our model, while the latter calls for improvement in data subtraction and processing. From this standpoint, there remains potential for advanced observations at energies from milli-eVs to MeVs toward the Galactic plane, in order to evaluate our model further and more comprehensively in the future.


Introduction
Interactive and ubiquitous, cosmic rays (CRs) play an impactful role in varieties of celestial events.Nonetheless, the origin of CRs remains a century-long mystery.CRs below the knee (around a few PeVs) are expected to be of Galactic origin.Away from their energetic accelerators, particles propagate in the magnetic halo of the Milky Way and interact with the interstellar gas, radiation field (ISRF) and magnetic field (IMF), generating secondary particles and photons.Unlike the particle measurements, which lose most of their initial directional information, the secondary photon emission preserves the spatial distribution of the progenitor CRs and thus turns out to be a unique and irreplaceable probe of CR propagation.
Wide-band diffuse emission from radio frequencies to ultra-high-energy (UHE) γ rays is yielded via interactions between CRs and the interstellar gas, the ISRF, and the IMF.In general, for a typical random magnetic field of a few µGs, the synchrotron radiation of CR electrons and positrons (CREs) results in radio emission from MHz to THz.The bremsstrahlung of CREs in the ISM generates high-energy emission from X rays to γ rays.The inverse Compton scattering (ICS) between CREs and the ISRF, together with the inelastic hadronic interactions between CR nuclei and the ISM give diffuse γ rays in a wide energy range.The multi-wavelength diffuse emission observations can therefore be used to constrain the source distribution and propagation of Galactic CRs.
With forefront space-borne and ground-based instruments, γ-ray observations have advanced into higher energy domains, enabling the possibility to study the GDE in multiwavelength approaches.Radio (Haslam et al. 1981a;Remazeilles et al. 2015) and microwave (Planck Collaboration et al. 2016a,b) surveys of early years, in conjunction with recent measurements at above hundreds of MeVs (Ackermann et al. 2012) and at even higher energies of multi-TeVs (Lemiere et al. 2015;Smith & VERITAS Collaboration 2015) can be comprehensively investigated to renovate and reconstruct the current theoretical framework of CRs.While some previous studies have taken data-driven, phenomenological methodologies, others have proposed variant refined models beyond our traditional understanding of CRs, postulating exotic origins (Bringmann & Salati 2007) and/or novel interaction mechanisms (Calore et al. 2015), but only a few of them show the capability to interpret the observed GDE spectrum across different wavelengths and other unforeseen anomalous phenomena simultaneously (Strong et al. 2010;Carlson et al. 2016;Orlando 2016).Straightforward comparison of high-level models each derived from within a certain energy range corresponding to a certain series of astronomical entities and processes could sometimes be misleading, due to the interdependence among the involved physical quantities and uncertainties introduced by different assumptions.Thus, it might be a safer and more refreshing perspective that we aim at fitting all available data at different energy ranges inside a unified configuration when trying to construct or assess our models, which partially motivates this work.Meanwhile, multi-wavelength observations have been used in many recent studies on both point and extended sources, which rely heavily on the accurate modeling of the GDE.
Featuring a spatially dependent propagation (SDP) scenario with extra CR origins beyond the standard paradigm, our model, which has been developed based on principally the up-to-date measurements of secondary-to-primary ratios (Zhang et al. 2023a), could be finer tested on this wise.In recent years, it has been argued extensively that the excesses in the observed CR spectra could be naturally explained by a rather young supernova remnant (SNR) located near us (Guo et al. 2016;Liu et al. 2018), whose contribution is also regarded in this work.The suggested CR confinement around the accelerators also influences the GDE intensity across different frequencies and will be depicted more specifically hereafter.
The rest of this paper is organized as follows.Section 2 describes the model setting and obtains the model parameters.Section 3 presents and discusses the results of wide-band diffuse emission.Finally, Section 4 concludes this work.

Model
During the active phases of varieties of astrophysical objects, particles are accelerated up to very high energies (VHE) and then go through a long voyage across the Galaxy before some of them enter the solar system.The whole process can be described from three major aspects: the injection from sources, the propagation in the Milky Way, and the production of secondaries.

Source injection
In this work, the SNRs are considered as the dominant sources of CRs and a continuous distribution (Case & Bhattacharya 1996) is adopted.The spatial distribution can be parameterized as where r ⊙ = 8.5 kpc is the distance of the Solar system to the Galactic center.Other parameters are adopted as a = 1.69, b = 3.33, and z s = 0.2 kpc (Case & Bhattacharya 1996).Besides this constituent defined as the background component, the contribution from individual nearby sources is also considerable, as evidenced by the curiosities in the energy spectra (Yoon et al. 2017;Atkin et al. 2018;Adriani et al. 2019;An et al. 2019) and large-scale anisotropies (Bartoli et al. 2015a;Aartsen et al. 2016;Amenomori et al. 2017;Abeysekara et al. 2019).To reproduce these anomalies, a local source is incorporated to our model.Though many astrophysical objects exhibit the ability to accelerate particles up to VHE, the SNRs are generally considered as the dominant ones.The injection spectra are assumed to follow a broken power-law energy distribution: where Q 0 is the normalization factor, ν 1 (ν 2 ) is the spectral index at rigidities lower (higher) than the break rigidity R br , and R c is the cut-off rigidity.The detailed information of the injection spectra is listed in Table 1 and Table 2 for the background and the local source respectively.

CR propagation
As mentioned above, we adopt an SDP model to describe the propagation of CRs in the Milky Way.The SDP model is primarily motivated1 by γ-ray observations of pulsar halos, which suggest very slow diffusion for regions surrounding those pulsars (e.g.Abeysekara et al. 2017;Aharonian et al. 2021) compared with the average diffusion coefficient inferred from CR secondary-to-primary ratios.It is a natural expectation that the slow diffusion regions in the Galactic plane are abundant due to the population of such middle-aged pulsars.Therefore, a two-zone diffusion (slow disk plus fast halo) scenario is reasonable to describe the propagation of CRs.It was also shown that this two-zone diffusion model can suppress the amplitudes of the dipole anisotropies, and the spatial variations of CR intensities and spectra from Fermi observations (Guo & Yuan 2018).
Following the previous work (Guo & Yuan 2018), we assume an anti-correlation between the diffusion coefficient with the source distribution.The diffusive halo is divided into two regions: the slow diffusion inner halo (IH) and the fast diffusion outer halo (OH).At the outer halo border (r = r h , z = ±z h ), the free-escape condition, ψ(r h , z, p) = ψ(r, ±z h , p) = 0, is imposed.With z h , ξz h , and (1 − ξ)z h being the half-thickness of the entire halo, the IH, and the OH respectively, the diffusion coefficient is given by where β is the velocity of the CR particle in unit of light speed c, and where g(r, z) = Nm 1+f (r,z) .The propagation parameters are summarized in Table 3.For a more complete and detailed description of the SDP model, we refer the readers to Zhang et al. (2023a).
As CRs enter the heliosphere, their energy spectra will be modified by the solar magnetic field.This so-called solar modulation effect is accounted for using the prevalent force field approximation (Gleeson & Axford 1968).In this work, without considering the chargesign dependent modulation effect, a constant modulation potential of ∼ 550 ± 150 MV is Table 3 0.56 0.55 0.1 4.0 0.05 4.0 5.0 6.0 assumed, which, together with other parameters, suffices to fit the observed CR spectra properly.However, it should be noted that this simplified treatment results in uncertainties of the calculated diffuse emission, particularly for the low-frequency (radio) and low-energy (γ-ray) parts.

Secondary production
As CRs travel through the Milky Way, they leave substantial imprints of secondary nuclei, leptons (positrons and electrons), and photon emission of radio to γ rays, which provide important probes to study their propagation.Relatively heavy secondary nuclei, such as lithium, beryllium and boron, are produced through fragmentation of the primaries chancing upon interstellar gas molecules, in which case each nucleon is considered to take after the energy of its parent.The production of these secondary particles can be described as follows: where ψ i is the density of the primaries, v i is the velocity of the parent particle, n H,He is the number density of hydrogen/helium in the ISM, and σ i+H/He→j denotes the cross section of the fragmentation process i → j.
Besides the fragmentation of heavy nuclei, inelastic collisions of light nuclei with the ISM will also produce secondary particles, such as antiprotons, electrons, positrons, and γ rays.The source term is the convolution of the primary spectra and the relevant differential cross sections: where ψ i (p i ) is the solution of the propagation equation of the primaries in the Milky Way.
The interactions near the acceleration sites are also considered, in which case, we assume that these particles have not experienced adequate propagation before they escape from the source regions, and ψ i = Q pri,i × τ , where τ is the effective confinement time of the particles around the source, which is estimated to be 0.2 Myr in our model through fitting the observed CR spectra.Approximately, the average grammage accumulated by all escaped primary particles from t = 0 to t = τ is X = ρcτ ≈ m p n H cτ = 0.32 g cm −2 , assuming a constant density n H = 1 cm −3 in the proximity of the sources for such injections.The charged secondary particles experience the same propagation procedure as primary CRs in the Milky Way.
The γ-ray emissivity from pp collisions can also be calculated with Eq. ( 6).For the production cross section, we use a recently developed interpolation routine based on Monte Carlo simulations, AAfrag, which employs the updated QGSJET-II-04m model tuned with the LHC data on hadronic interactions (Kachelrieß et al. 2019).For proton energies below 4 GeV and helium energies below 5 GeV, where AAfrag does not cover, an old cross section model of Dermer (1986) is used.
The bremsstralung emissivity of CREs is ) (7) where α is the fine structure constant, r 0 is the classical electron radius, c is the speed of light, n s is the number density of gas of species s, E e is the energy of CREs, ψ e is the differential density of CREs, ϵ is the energy of emitted photons, ϕ 1,s and ϕ 2,s are shielding factors which can be found in Blumenthal & Gould (1970).
The gas distributions used in this work are the ones embedded in GALPROP.They are HI distribution from Gordon & Burton (1976) renormalized to the results given by Dickey & Lockman (1990), the molecular gas distribution from Bronfman et al. (1988) with a conversion factor from CO to H 2 of ∼ 1.9 × 10 20 mols cm −2 /(K km s −1 ) (Strong & Mattox 1996), and the ionized hydrogen distribution from Cordes et al. (1991).As for the calculation of the diffuse emission, the gas column densities for each line of sight are further renormalized to match the data from surveys of HI (Kalberla et al. 2005) and molecular gas (Dame et al. 2001).
In addition to the interactions depicted above, CREs also emit γ rays via the ICS process in the ISRF.The ICS γ-ray emissivity is (Moskalenko & Strong 2000) where γ is the Lorentz factor of parent electrons, ϵ 1 is the energy of seed photons, ϵ 2 is the energy of scattered photons, (n γ ϵ 2 1 f γ ) and (n e γ 2 f e ) are the differential spectra of seed photons and parent electrons, and is the yield function of ICS photons for electron-photon interactions of fixed energies, for an isotropic distribution of photons, where q ′ = ϵ 2 /[4ϵ 1 γ 2 (1 − ϵ 2 /γ)] and 1/4γ 2 < q ′ ≤ 1.
VHE γ rays would experience attenuation when traveling through the ISRF due to the pair production effect (Zhang et al. 2006;Moskalenko et al. 2006).We calculate the absorption fraction of the three dimensional emissivity of γ rays, and then integrate along the line-of-sight to get the fluxes (Zhang et al. 2023b).
The Milky Way is filled with randomly oriented magnetic field, and the electrons will lose energy and produce wide-band emission through the synchrotron process.After averaging over the pitch angle for an isotropic electron distribution, the synchrotron emissivity of a single electron integrated over all directions is (Ghisellini et al. 1988) where ν is the radiation frequency, γ is the electron Lorentz factor, ν B = eB/(2m e c), x = ν/(3γ 2 ν B ), and K 4/3 , K 1/3 are modified Bessel functions.The Galactic magnetic field strength is modeled as B = B 0 e −R/30 kpc−|z|/4 kpc .The observable synchrotron intensity is then obtained by integrating the emissivity along the line of sight: To be conclusive, the major ingredients of the model include an SDP propagation scheme, a nearby source component which mainly accounts for the spectral bumps of primary CRs, and the secondary production around acceleration sources which is to explain the B/C hardening and diffuse γ-ray excesses.Note that the spectral breaks of primary and secondary particles are likely a coincidence.With higher energy measurements of the primary spectra by DAMPE (An et al. 2019;Alemanno et al. 2021) and CALET (Adriani et al. 2019(Adriani et al. , 2020(Adriani et al. , 2023)), it was shown that the break rigidity is about 500 GV for protons, about 650 GV for helium, 400 -500 GV for carbon and oxygen nuclei.The break rigidity of B/C and B/O measured by DAMPE is about 200 GV (Alemanno et al. 2022).Therefore, the break energies for primary and secondary nuclei may indeed have different origins.

Galactic diffuse emission
Starting from the model setting described in Sec. 2, we adjust the model parameters to match the measured spectra of CR protons, positrons, CREs, and boron-to-carbon ratio, as shown in Figure A1 in Appendix A. The diffuse emission from neV to PeV energies is then calculated, which will be discussed in detail as follows.2018), based on original measurements: radio surveys (green stars) (Guzmán et al. 2011;Landecker & Wielebinski 1970;Haslam et al. 1981bHaslam et al. , 1982;;Reich et al. 1999Reich et al. , 2001)), Planck synchrotron map (teal circle) (Planck Collaboration et al. 2016a) and WMAP (magenta arrows) (Bennett et al. 2013).In the left panel, the black solid line represents the overall intensity with modulation potential of 550 MV.The colour orange and green mark the radiation from background and fresh interactions respectively; the dotted, dashed and long-dashed lines correspond to contributions by CREs at various energy ranges: [0.1 GeV, 1 GeV], [1 GeV, 10 GeV] and [10 GeV, 100 GeV].In the right panel, the total synchrotron fluxes for different solar modulation potentials are shown.Note that the magnetic field strength for each case is slightly re-scaled to better fit the wide-band data.

Diffuse radio emission
The calculated synchrotron fluxes are shown in the left panel of Figure 1, compared with the data taken from Orlando (2018), for a sky region slightly above the Galactic plane (10 • < |b| < 20 • ).The model prediction is roughly consistent with the observations.At high frequencies the model flux is slightly lower than the Planck data.This may be solved by a slight re-scale of the magnetic field strength and lowering CRE fluxes below ∼ 1 GeV (see discussion on the uncertainty of the solar modulation in the next paragraph).Note that the WMAP results seem to be higher than the Planck flux at similar frequency (23 GHz).As discussed in Orlando (2018), there might be degeneracies among different components in such frequency bands, e.g., synchrotron, free-free, thermal dust, and anomalous microwave emissions.The WMAP synchrotron intensities may be over-estimated.We therefore use WMAP results as upper limits (shown by arrows).
To better see the mapping relation between CREs and synchrotron emission, we show the contributions from CREs in three energy bands, 0.1 − 1 GeV, 1 − 10 GeV, and 10 − 100 GeV, respectively.The emission below 1 GHz mainly comes from CREs with energies smaller than 1 GeV, in which range the fluxes are uncertain due to solar modulation.In the right panel of Figure 1, we show the synchrotron fluxes for three different values of the modulation potentials, 550 MV as the benchmark, 400 MV, and 700 MV, respectively.For each modulation parameter, the source parameters are tuned to reproduce the measured CR spectra (see Table 1).The magnetic field strength is slightly re-scaled to better fit the radio data.For Φ = 400 MV, the low energy CRE fluxes are lower (Figure A1), and the resulting synchrotron spectrum is harder, which can better match the data in a wide frequency band.

Diffuse γ-ray emission
The comparisons between model calculations and observed γ-ray data in different sky regions are shown in Figure 2. The Fermi-LAT measurements in two latitude belts, 10 • < |b| < 20 • and 8 • < |b| < 90 • , for the whole longitude range (Ackermann et al. 2012), the ARGO-YBJ (Bartoli et al. 2015b) and Tibet-ASγ (Amenomori et al. 2021) measurements in the inner Galactic plane region (25 • < l < 100 • , |b| < 5 • ), the Milagro measurement in a smaller region in the inner Galactic plane (30 (Abdo et al. 2008), and the CASA-MIA upper limits a region covering mostly the outer plane (50 (Borione et al. 1998), are employed for comparison.It is shown that in general the model can reproduce the observations relatively well.Only for ∼GeV energies the model fluxes are slightly higher.We keep in mind that the uncertainty of the solar modulation may take effect in such an energy range.
The agreement is satisfactory to an extent: at intermediate latitudes, the contributions from GCRs and extragalactic backgrounds are comparable; at the Galactic disk, the predominant emission is from the Galactic background CR nuclei, with a small portion of electronic origins.The good consistency between the predictions and the observations of γ rays through the sky further demonstrates the reliability of our model.In addition, the Galactic background CRs and the fresh ones make major contributions below and above tens of GeVs respectively.This feature has already been reflected in the hardening of the CR spectra and secondary-to-primary ratio spectra as shown in Figure A1 in Appendix A. It suffices to say that our model has been tested in the γ-ray band from tens of MeVs to hundreds of TeVs.
-Diffuse emissions from radio to PeV γ rays.The black solid line is the total radiation.For the rest of the lines, the colour orange, pink, green and teal denotes comfrom synchrotron, Bremmstrahlung, ICS and π 0 decay; the dots and dashes tell whether they are generated from background or fresh interactions; the magenta and gray solid lines are contributions by the IGRB and resolved sources detected by Fermi-LAT.Data are from: Fermi-LAT (Ackermann et al. 2012), ARGO-YBJ (Bartoli et al. 2015b), Tibet ASγ (Amenomori et al. 2021), CASA-MIA (Borione et al. 1998), Milagro (Abdo et al. 2008).
Fig. 3.-Diffuse γ-ray emission of the inner region (left) and the outer region (right) of the Galaxy.Areas containing resolved sources detected by LHAASO-KM2A is masked, the masked region is referred to FIG. 1 in Cao et al. (2023).The red, black and blue solid lines represent results predicted by models with modulation potentials of 400 MV, 550 MV and 700 MV.

Spatial distributions
The energy spectra from radio to PeV γ-ray emission have been calculated and compared with the measurements.There exists satisfactory consistency between the calculations and the observations in the γ-ray band over the sky.Moreover, with the aim of conducting a more comprehensive and detailed study the CR propagation, we should delve deeper into their spatial distribution, especially those at different energy bands, to study the evolution of CR composition.Figure 4 shows the diffuse emission skymaps for radiation at the energy of 1.5×10 −9 MeV, 1.4×10 3 MeV and 1.4×10 6 MeV, where two distinct features are noticeable.Firstly, the distribution appear rather smooth in the radio band and yet somewhat uneven in the range of γ rays.Secondly, only at above TeVs does the fresh component dominate.To explore these features a step further, the 1-d spatial distributions along the Galactic latitudes and longitudes are explored and presented in Figure 5, showing similar attributes to those exhibited already in Figure 4.In addition, as the injected source distribution is assumed to be identical for the CREs and CR nuclei during our calculation, the heterogeneity of the γ rays must have originated from the distribution of the ISM.Hopefully the Fermi-LAT and the LHAASO experiments will provide cleanly subtracted diffuse skymaps in the γ-ray range and the space-borne experiments will give similar measurements in the radio and X-ray bands in the future.

Comparison with other models
Similar to ours, modifications on the diffusion coefficient have been performed in other studies, e.g., a linear dependence on the distance to the Galactic center (Gaggero et al. 2015).Alternatively, despite minimized by our model, the influence of processes such as reacceleration and convection can also be highlighted (e.g.Orlando 2018;Qiao et al. 2022).The unresolved γ-ray sources have also been brought up to explain the high-energy GDE (Vecchiotti et al. 2022;Schwefer et al. 2023), but the correspondent radiation mechanism remain to be further understood.Some models take account of variants less recognized in this context, e.g., a spatially dependent X CO (Gaggero et al. 2015), with uncertain factors such as Galactic chemical evolution.Predictions from some of these models raised above are drawn together with those from this work, as shown in Figure 6.At radio frequencies, we compare with the baseline DRE model of Orlando (2018) instead of the best-fit DRElowV model which employed some arbitrary tuning of the Alfven velocity for particular particle species.Our model gives relatively lower fluxes at low frequencies, and can better match the data.This is primarily due to that the reacceleration in our case is not as strong as the DRE model of Orlando (2018).Note that the simplified magnetic field model we adopt in this work may affect the results as explored in (Orlando & Strong 2013).As for the diffuse γ-ray emission, previous works tend to give higher fluxes above 100 TeV, to fit the ASγ data (Amenomori et al. 2021).Updated fitting to the LHAASO measurements may give slightly lower fluxes at the UHE band (Cao et al. 2023).Our model prediction differs from those works mainly in the TeV region, featuring by a bump-like structure of our model prediction, which could be tested with future measurements by LHAASO (Li et al. 2023).The one-dimensional distributions along Galactic latitudes are compared with the results of a homogeneous diffusion model of Zhang et al. (2023b).The results show that the SDP model in this work gives a faster decrease from the disk to the pole regions, due to a slower diffusion in the Galactic disk.(Orlando 2018;Lipari & Vernetto 2018;De La Torre Luque et al. 2023;Schwefer et al. 2023;Zhang et al. 2023b).Left: of synchrotron radiation; middle: fluxes of diffuse γ-ray emission; right: Galactic latitude distribution of diffuse γ-ray emission at 1.5 × 10 7 MeV energy.(Cummings et al. 2016), AMS-02 (Aguilar et al. 2019b(Aguilar et al. ,a, 2015(Aguilar et al. , 2018)), CREAM (Yoon et al. 2017) and DAMPE (Ambrosi et al. 2017;Alemanno et al. 2022;An et al. 2019), and CALET (Adriani et al. 2022).
Fig.1.-Themodel predicted synchrotron emission, compared with data compiled or processed inOrlando (2018), based on original measurements: radio surveys (green stars)(Guzmán et al. 2011;Landecker & Wielebinski 1970;Haslam et al. 1981bHaslam et al. , 1982;;Reich et al. 1999Reich et al. , 2001)), Planck synchrotron map (teal circle)(Planck Collaboration et al. 2016a) and WMAP (magenta arrows)(Bennett et al. 2013).In the left panel, the black solid line represents the overall intensity with modulation potential of 550 MV.The colour orange and green mark the radiation from background and fresh interactions respectively; the dotted, dashed and long-dashed lines correspond to contributions by CREs at various energy ranges: [0.1 GeV, 1 GeV], [1 GeV, 10 GeV] and [10 GeV, 100 GeV].In the right panel, the total synchrotron fluxes for different solar modulation potentials are shown.Note that the magnetic field strength for each case is slightly re-scaled to better fit the wide-band data.