Star-forming galaxies dominate the diffuse, isotropic 1 γ -ray background 2

8 The Fermi Gamma Ray Space Telescope has revealed a diffuse, isotropic γ -ray background at energies ranging from 0.1 GeV to 1 TeV [1] whose astrophysical sources remain uncertain. Previous efforts to understand the origin of this background have been hampered by the lack of physical models capable of predicting the γ -ray emission produced by the many candidate sources, which include star-forming galaxies (SFGs) [2–6], active galactic nuclei (particularly blazars [7–9]), millisecond pulsars [7], and dark matter annihilation [10]. In the absence of predictive models, estimates of the contribution from potential sources have relied on a highly-uncertain process of empirically scaling the emission from a small sample of local, resolved sources by their estimated cosmological abundances. Here we present the ﬁrst calculation of the contribution of SFGs to the γ -ray background that is based on a physical model for the γ -ray emission produced when cosmic ray ions accelerated in supernova remnants interact with the interstellar medium [11]. We validate the model by showing that it reproduces the γ -ray spectra, source count distribution and far infrared- γ -ray correlation observed for nearby, resolved SFGs. When we apply the model to the observed cosmological SFG population, we recover an excellent match to the γ -ray background from 1 GeV to 1 TeV. Our result shows that SFGs alone can explain the full diffuse γ -ray background over this energy range, and strongly suggests that emission in excess of our model at energies < 1 GeV originates from cosmic ray electrons produced in the same galaxies.


Introduction
supernova remnants [13]. This process transfers ∼10% of the supernova mechanical energy to relativistic particles, yielding on average ∼ 10 50 erg in CR ions per supernova [14,15], with a substantially smaller portion deposited in CR electrons. The 19 resulting CRs follow a power law distribution in particle momentum, well approximated for relativistic CRs as a power law in 20 total particle energy E [16,17], dN/dE ∝ E −p , with an index p that observations of individual supernova remnants, analytical 21 models, and numerical simulations all indicate to be in the range p ≈ 2.0 − 2.6, with a mean value p ≈ 2.2 − 2.3 [18,19]. Some of these CR ions collide inelastically with ISM nuclei, producing roughly equal numbers of π 0 , π + and π − mesons that rapidly 23 decay via the channels π 0 → 2γ, π − → µ − +ν µ , and π + → µ + + ν µ . The decay of π 0 particles is responsible for most of the 24 observed Galactic γ-ray foreground, which displays a characteristic spectrum that rises sharply from ∼ 0.1 − 1 GeV as a result 25 of the 135 MeV rest mass of the π 0 particle. 26 2 Galactic emission model 27 The diffuse γ-ray emission from a SFG depends on three factors: its total star formation rate (which determines its supernova 28 rate and thus the rate at which CRs are injected), the distribution of γ-ray energies produced when individual CRs collide with 29 ISM nuclei (which depends on the parent CR energy E) and the fraction of CRs (again as a function of E) that undergo inelastic 30 collisions before escaping the galaxy. The latter is known as the calorimetry fraction f cal , and is the most poorly-understood of 31 the three. 32 For the purposes of the calculations we present here, we compute the conversion from star formation rate to CR injection by 33 assuming that stars form with a Chabrier initial mass function [20], which gives the distribution of mass for newly-formed 34 stars, and that stars with initial masses of 8 − 50 M , where M is the mass of the Sun, end their lives as supernovae 35 [21]. Each supernova injects 10 50 erg of energy in CR ions, distributed in energy for CR energies E > m p c 2 as dṄ/dE = 36 φṀ * (E/E 0 ) −p exp (−E/E cut ), where E 0 = 1 GeV, the cutoff energy E cut = 10 8 GeV 1 ,Ṁ * is the galaxy's star formation rate, 37 and the spectral index p = 2.2 [18,19]. The normalisation factor φ that corresponds to our choice of initial mass function, 38 supernova mass range, and CR energy per supernova is φ = 6.52 × 10 42 s −1 GeV −1 M −1 yr. 39 Similarly, we calculate the γ-ray spectrum produced by CRs of energy E using the parameterised model of Ref. [22], which gives the differential cross section dσ γ /dE γ E γ , E for the production of γ-rays of energy E γ by CRs of energy E in inelastic collisions. Combining these factors, the rate at which a SFG produces γ-ray photons is where σ pp = 40 mbarn is the mean proton-proton inelastic cross-section. The calorimetry fraction f cal (E) depends on the 40 properties of the galaxy, and the lack of a model for this dependence has previously precluded proper evaluation of Equation 1.
[11] recently introduced a model for f cal (E), based on the idea that the rate at which CRs diffuse through the The γ-rays produced inside a galaxy do not propagate unhindered, but are subject to losses due to pair production in 48 collisions with far-infrared (FIR) photons inside the source galaxy and with photons that are part of the extragalactic background 49 light (EBL) that pervades intergalactic space; the latter is generally the more important effect. 50 We compute the optical depth τ γγ due to galactic FIR photons using the model of Ref.
[23], and the optical depth τ EBL due to the EBL using the model of Ref.
[24]. Taking these into account, Equation 1 can be used to compute the specific photon flux dF γ /dE γ (i.e., the number of photons per unit area, time and energy) received at the Earth from a galaxy at redshift z as where d L (z) is the luminosity distance of the source. The radiation that is absorbed by the galactic and extragalactic photon 51 fields is reprocessed to lower energies in the pair-production cascade. In this process, both the electron and positron, from the 52 initial high energy pair produced, inverse Compton scatter lower energy photons up to γ-ray energies and these, in turn, produce 53 further pairs, and so on. The photon spectrum from this cascade can be parameterised according to the method developed in 54 Ref. [25]. For the purposes of our calculation here, we include the effect of the cascade by adding a component to dF γ /dE γ 55 with a spectral shape as computed by Ref.
[25], and with a normalisation such that its energy is equal to the integrated energy 56 lost to photon-photon scattering. In our model the cascade contributes 5 -20% of the total emission between 1-100 GeV. 57 3 γ-rays from star-forming galaxies 58 We now have a model that predicts the γ-ray emission of a SFG. The next step in our analysis is to apply this model to a galaxy Ref.
[27]. The full sample contains 34,930 galaxies, but we exclude those whose parameters are uncertain because they contain 63 bright active galactic nuclei, have unreliable redshifts, or lack a good fit to the surface brightness profile. This leaves a sample 64 of 22,279 galaxies.
total masses, surface densities, and volume densities, we assume that the galactic disc has an area of 2πR 2 e and a thickness 70 of 2h g . We then derive the gas surface density from the stellar and star formation rate surface densities, Σ * = M * /2πR 2 e and 71Σ * =Ṁ * /2πR 2 e , using the empirical Kennicutt-Schmidt relation obtained by Ref. [28]. Likewise, we derive the gas velocity 72 dispersion from the star formation rate using the empirical relation obtained by Ref. [29]. Finally, we derive the scale height 73 assuming the galactic disc is in hydrostatic equilibrium [30]. We provide full details of these procedures in the Methods.

74
To verify that this approach predicts reasonably accurate γ-ray spectra, we apply it to four local, resolved galaxies, chosen 75 to span a wide range of gas and star formation surface densities: Arp 220, NGC 253, M31, and NGC 4945. The input data we 76 use for these calculations are summarised in Extended Data Table 1, and we show the results of the computation in Figure 1,

77
where the solid lines show the spectra derived using only stellar data and, for comparison, the dotted lines show the results we 78 obtain if we add directly-measured gas properties. We see that the fits are slightly improved if we make direct use of gas data 79 but, even for the stellar data only, our model reproduces the observed γ-ray spectra to better than a factor of 2 for all galaxies at 80 energies > 1 GeV, and within a factor of ≈ 1.5 for the two more rapidly star-forming galaxies, which, as we show below, are 81 more akin to the population that dominates the γ-ray background. We do not expect agreement at energies 1 GeV because 82 we have not included the leptonic component of the γ-ray emission, which dominates below the ∼ 0.1 − 1 GeV cutoff in the 83 hadronic component set by the π 0 rest mass.

84
Having verified that we can obtain accurate predictions of γ-ray spectra from stellar data alone, we carry out two additional 85 validation steps. First, we examine the correlation between galaxies' FIR and γ-ray luminosities. We calculate the total γ-ray Clouds. This is not surprising, since these are not part of the field galaxy population we are simulating, and are much closer to 106 the Milky Way than we expect any randomly-chosen field galaxy to be. Having validated our model, we are now ready to compute the contribution of SFGs to the diffuse, isotropic γ-ray background.
We do this by calculating an observed spectrum for each galaxy in the CANDELS sample using Equation 2, and summing over all galaxies in bins of redshift: Here Ω S = 173 arcmin 2 is the solid angle surveyed by CANDELS, n S,j the number of surveyed galaxies in the jth redshift bin, 109 and n zbin the number of redshift bins. The factor f corr,j is the ratio of the expected total star formation rate in each redshift bin 110 (based on the measured cosmic star formation history [36]) to the sum of the star formation rates of CANDELS galaxies in that 111 bin; its purpose is to correct for the fact that, due to its limited field of view and various observational biases, the distribution of 112 star formation with respect to redshift in CANDELS does not precisely match the total star formation history of the Universe. 113 We use redshift bins of size ∆z = 0.1, chosen to ensure that the number of sample galaxies in each bin is large enough that the 114 uncertainty in the mean spectrum due to Poisson sampling of the galaxy population is small. 115 We present the results of this calculation in Figure 4. The figure shows that the expected contribution of SFGs to the diffuse 116 isotropic γ-ray background fully reproduces both the intensity and the spectral shape of the observations from ≈ 1 GeV to ≈ 1 117 TeV. We emphasise that we obtain this agreement from our model with no free parameters: our only inputs are the CR injection obtain simply by setting f cal = 1 for all galaxies at all energies, which clearly both overestimates the intensity and yields a 122 spectral slope that is flatter than observed. By contrast, the EBL is important for the spectral shape only for energies > 100 GeV, 123 as illustrated by the green dashed line in Figure 4, which shows the result of a calculation where we omit γγ opacity effects. 124 We show the relative contributions to the background made by galaxies with differing star formation rates and redshifts in for their redshift. contributor. It is nonetheless instructive to examine the precise reasons why our conclusions differ from some earlier work.

142
One contributing factor is that earlier models were forced to adopt single power laws for the emitted γ-ray spectrum 143 in different classes of galaxies [39, 41]. By contrast, Figure 1 demonstrates that none of the four nearby resolved galaxies 144 shown have spectra that are well described by a γ-ray spectrum in the form of a pure power law over the energy range from 145 E γ = 1 − 1000 GeV; our model correctly captures this behaviour, but earlier pure power law models did not. Similarly, we 146 calculate f cal as a function of energy directly, rather than relying on an empirical FIR-γ correlation, and our calorimetry 147 fractions are on average larger than those implicitly assumed in earlier works. This is because many of the lower estimates 148 for the contribution from SFGs to the γ-ray background rely on a FIR-γ relation derived from early Fermi detections of < 10 149 individually-resolved SFGs [3] that yields somewhat lower γ-ray luminosities than more recent fits using a larger (but still 150 small) sample of SFGs [33], and with which our model agrees (Figure 2). Thus the reason we find that SFGs can produce the 151 full background, whereas earlier models could not, is that our model predicts γ-ray emission that is both somewhat brighter and 152 has a more complex spectral shape than the values adopted in earlier work. active galactic nuclei using a radio-γ relation derived from a sample of 16 resolved objects, coupled with a radio luminosity 156 function extrapolated to redshifts considerably higher than those well-sampled by observations [43]. By contrast, our assignment explain the background below 1 GeV as arising from the leptonic counterpart to the hadronic emission we have measured.

162
Comparing the total power in the γ-ray background from 0.1 − 1 GeV to the total available power in hadronic cosmic rays (the 163 integral under the orange dashed line in Figure 4) indicates that a small fraction of the hadronic power would be sufficient to 164 explain the observed signal, so this explanation is plausible on energetic grounds. However, we leave detailed modelling of the 165 leptonic emission for future work. The second question is whether SFGs also dominate the observed neutrino background 166 at energies > 1 TeV. We address this issue in the Supplementary Information,   ; in the left panel we show the local starburst galaxies Arp 220 and NGC 253, and in the right panel we show the local quiescent galaxies M31 and NGC 4945. The solid lines show model predictions using only stellar data of the type we have available for the CANDELS sample, while the dotted lines shows results predicted if we supplement this with observed gas data. We list the full set of observed quantities used in computing these models in Extended Data Table 1.     [1]; the blue line shows our prediction for the background due to SFGs; the red dot-dashed line shows the contribution to this spectrum due to the γγ scattering cascade, and the green dashed line shows the spectrum we would obtain in the absence of γγ opacity. For comparison, the orange line is the model prediction for full calorimetry, i.e., with f cal (E) = 1 for all galaxies at all CR energies.

11/19
Methods Model for calorimetry 192 We employ the model of Ref. [11] for CR transport and calorimetry in SFGs. The basic premise of the model is that, in the neutral phase that dominates the mass of the ISM and thus the set of available targets for γ-ray production, CR transport is primarily by streaming along magnetic field lines. However, this yields approximately diffusive transport when averaged over scales comparable to or larger than the coherence length of the magnetic field, with a diffusion coefficient D V st is the CR streaming speed, h g is the gas scale height, and M A is the Alfvén Mach number of the turbulence. For diffusive transport with losses in a disc geometry, the calorimetry fraction is given by (using the favoured parameters of Ref. [11]) where 0 F 1 is the generalised hypergeometric function and τ eff is the dimensionless effective optical depth of the ISM, given by Here σ pp = 40 mbarn is the mean proton-proton inelastic cross section, η pp = 0.5 is the elasticity of pp collisions, Σ g is the gas 193 surface density of the galactic disc, c is the speed of light, D 0 is the diffusion coefficient at the galactic midplane, µ p = 1.17 is 194 the number density of nucleons per proton, and m H = 1.67 × 10 −24 g is the mass of a hydrogen atom.

195
To evaluate the calorimetry fraction for a CR of energy E, we must therefore determine the midplane diffusion coefficient D 0 for CRs of that energy, which in turn depends on the streaming speed V st . This speed is dictated by the balance between excitation of the streaming instability and dissipation of the instability by ion-neutral damping, the dominant dissipation mechanism in the weakly-ionised neutral ISM. Balancing these two effects yields a CR proton streaming velocity where V Ai is the ion Alfvén speed, γ d = 4.9 × 10 13 cm 3 g −1 is the ion-neutral drag coefficient, χ is the ionised mass fraction, 196 ρ = Σ g /2h g is the midplane mass density of the ISM, C is the midplane number density of CRs, e is the elementary charge, 197 u LA is the velocity dispersion of Alfvén modes in the ISM at the outer scale of the turbulence, µ i is the atomic mass of the 198 dominant ion species, p is the index of the CR energy distribution, and γ = E/m p c 2 is the Lorentz factor of the CR. Since 199 γ-ray production in our model is dominated by galaxies with high star formation rates and gas surface densities within which i) 200 the ISM is molecule-dominated, ii) the main ionised species is C + , and iii) the magnetic field is set by a turbulent dynamo, 201 we adopt the fiducial parameters of Ref.

13/19
At this point we have specified all the ingredients required to compute f cal (E) for a galaxy of known Σ g , σ g , and h g , save one: C, the CR number density. We estimate this as follows: consistent with our discussion in the main text, for a galaxy with star formation rate per unit areaΣ * , the relativistic cosmic ray injection rate per unit area is The CR number density at the midplane is then given by where t loss is the CR loss time. This is given by t loss = 1/ t −1 Our calorimetry model requires, as input, the gas surface density Σ g , scale height h g , and velocity dispersion σ g , along with the 212 surface density of star formationΣ * . However, the CANDELS data set that we use provides only the cosmological redshift z, 213 stellar mass M * , half-light or effective radius R e (corrected to 5000 Å according to Ref. [48]), and total star formation rateṀ * 214 for our sample galaxies. We must therefore estimate the gas properties from observed correlations between gas and stellar 215 properties. We do so as follows.

216
The half light radius R e at 5000 Å serves as a first order estimate of how the star formation and matter are distributed throughout the galactic disc. We therefore estimate the star formation rate surface density asΣ * =Ṁ * /2πR 2 e and the stellar surface density as Σ * = M * /2πR 2 e . We estimate the gas surface density from the observed correlation between gas, stellar, and star formation surface densities given by Ref.
[28]: Similarly, there is a strong correlation between galaxy star formation rates and velocity dispersions, which we use to derive σ g .
For this purpose we fit the relationship using the MaNGA galaxy sample [29]. A powerlaw fit to the data obtained in this survey Finally, we derive the gas scale height under the assumption that the gas is in vertical hydrostatic equilibrium, in which case the scale height is [30] h g = σ 2 g πG (Σ g + Σ * ) When applying this model to the galaxies in the CANDELS sample, we take the stellar properties from Ref.
[27] and correct 217 the half-light radii to 5000 Å following the procedure outlined in Ref. [48]. This gives us a sufficient set of input parameters to 218 apply our model.