EGRB from starforming galaxies: do empirical scalings suffice?

# Extragalactic Gamma-ray Background from Star-forming Galaxies: Will Empirical Scalings Suffice?

## Abstract

Despite the influx of unprecedented-quality data from the Fermi Gamma-Ray Space Telescope that have been collected over nine years of operation, the contribution of normal star-forming galaxies to the extragalactic gamma-ray background (EGRB) remains poorly constrained. Different estimates are discrepant both their underlying physical assumptions and their results. With several detections and many upper limits for the gamma-ray fluxes of nearby starforming galaxies now available, estimates that rely on empirical scalings between gamma-ray and longer-wavelength luminosities have become possible and increasingly popular. In this paper we examine factors that can bias such estimates, including: a) possible sources of nontrivial redshift dependence; b) dependence on the choice of star-formation tracer; c) uncertainties in the slope and normalisation of empirical scalings. We find that such biases can be significant, pointing towards the need for more sophisticated models for the star-forming galaxy contribution to the gamma-ray background, implementing more, and more confident, physics in their buildup. Finally, we show that there are large regions of acceptable parameter space in observational inputs that significantly overproduce the gamma-ray background, implying that the observed level of the background can yield significant constraints on models of the average cosmic gamma-ray emissivity associated with star formation.

###### keywords:
gamma-rays: diffuse background – gamma-rays: galaxies – methods: analytical – galaxies: luminosity function, mass function –

## 1 Introduction

The gamma-ray sky consists of resolved point sources, Galactic diffuse emission, i.e. photons produced in interactions between cosmic rays and the interstellar medium and interstellar radiation field, and an isotropic, presumably extragalactic diffuse emission, the extragalactic gamma-ray background (EGRB, Clark et al. 1968; Fichtel et al. 1978; Sreekumar et al. 1998; Strong et al. 2004; Abdo et al. 2010; Ackermann et al. 2015). The EGRB is a superposition of individual unresolved point sources and truly diffuse emission, and it encodes valuable information about high energy processes in the universe. The dominant components of the EGRB are likely active galaxies (blazars), which are the brightest extragalactic sources and the most numerous population of resolved gamma-ray sources, and normal star-forming galaxies. These two classes of gamma-ray sources are guaranteed to contribute to the EGRB.

The relatively large number of resolved gamma-ray blazars has made possible the construction of empirical gamma-ray luminosity functions, which, when extrapolated to fainter fluxes, place the contribution of unresolved blazars around of the Fermi LAT measurement of the EGRB (Ajello et al. 2015). Star-forming galaxies on the other hand are individually faint in gamma rays, and only a handful of sources have been individually resolved (Ackermann et al. 2012; Tang et al. 2014; Peng et al. 2016). As a result, the construction of empirical gamma-ray luminosity functions remains infeasible, and indirect ways of estimating their contribution to the EGRB have to be employed.

In normal star-forming galaxies the dominant emission mechanism of gamma-rays is through hadronic interactions between cosmic rays (CR) and interstellar gas (IG) producing neutral pions, i.e., . Gamma-rays are produced through pion decay (Stecker 1971). A critical assumption in the case of normal star-forming galaxies is that the level of cosmic ray flux is set by the balance between cosmic-ray acceleration (likely in supernova remnants) and diffusive escape of cosmic rays to the intergalactic medium (e.g., Pavlidou & Fields 2001). In order for this assumption to hold, cosmic-ray losses due to escape have to dominate over losses due to pion production. In this case, the flux of gamma rays depends on: a) the flux of projectiles (cosmic rays) and b) the number of targets (IG particles). Since both these quantities are related to the star formation rate (SFR) of a galaxy, the gamma-ray emissivity of a galaxy is expected to depend on its star-formation properties.

In contrast to normal star-forming galaxies, starburst galaxies tend to have much higher gas fractions, and losses due to pion production may dominate. In this case, the flux of cosmic rays is set by the balance between acceleration and pion production, and the flux of gamma rays is simply proportional to the cosmic ray acceleration rate. The gamma-ray emissivity of such a calorimetric galaxy will still depend on its star formation rate, but following a different scaling than that of normal galaxies (e.g., Thompson et al. 2007, Wang & Fields 2016). For this reason, normal and starburst galaxies are treated as distinct classes with respect to their contribution to the EGRB.

Before the launch of Fermi, the only star-forming galaxies with confirmed gamma-ray emission were the Milky Way and the Large Magellanic Cloud. Of necessity, then, early studies of the star-forming galaxy contribution to the EGRB relied on theoretical arguments (e.g., Pavlidou & Fields 2002; Fields et al. 2010; Makiya et al. 2011; Stecker & Venters 2011). However, Fermi has now made possible the detection of several nearby star-forming galaxies, and provided upper limits for the gamma-ray flux of many more (Ackermann et al. 2012). As a result, it is now possible to obtain empirical scalings between gamma-ray luminosity and luminosities at lower, star-formation-tracing frequencies, and use such scalings together with lower-frequency luminosity functions to obtain estimates of the EGRB contribution from star-forming galaxies (e.g., Stecker & Venters 2011; Ackermann et al. 2012; Linden 2017).

However, there is still no consensus among these different studies on the fraction of the EGRB that can be attributed to star-forming galaxies. The difference is not just arithmetical: rather, there is tension between the theoretically expected and the empirically observed slope of the correlation between SFR and gamma-ray luminosity (). While Fields et al. (2010), using the Kennicutt-Schmidt (KS) relation, find a slope of , Ackermann et al. (2012) and Linden (2017) report a best-guess slope closer to . Given that star-forming galaxies are very faint in gamma rays, it is not expected that enough members of this class will be individually resolved in the foreseeable future to construct empirical luminosity functions. It is therefore important to understand the origin of any such tensions; to investigate whether empirical scalings and longer-wavelength luminosity functions can act as a reliable substitute of empirical luminosity functions; and identify any potential biases in the construction of models for the star-formation-related EGRB components from empirical scalings, as well as ways to address and remedy such biases. This is the scope of this paper.

Here, we examine two factors as potential sources of bias in EGRB estimates of star-forming galaxies based on empirical scalings between gamma-ray luminosity and longer-wavelength luminosities. First, possible sources of nontrivial redshift dependence, such as the dependence of the supernova progenitor mass cutoff on metallicity, or a redshift-dependent performance of different star formation tracers. And, second, uncertainties in the slope and normalisation of empirical scalings, arising either due to the extremely-low-number statistics involved, or due to the evolution with cosmic time of the physical processes setting these scalings.

This paper is organised as follows. In §2 we describe the general formalism for a luminosity-function–based calculation of the contribution of a certain class to the EGRB. In §3 we discuss how the different factors discussed above enter the construction of a gamma-ray luminosity function. Our results are presented in §4, and they are discussed in the context of recent work in gamma-ray astrophysics and the cosmic history and star formation in §5, where a roadmap towards future improvements in the modeling of gamma-ray emission for star-forming galaxies is also suggested.

## 2 Formalism

The diffuse gamma-ray photon intensity produced by primary emission from unresolved star-forming galaxies with luminosity function , where is the comoving number density, is given by the the line-of-sight integral over the gamma-ray luminosity function:

 IE(E)=zmax∫0dzLγ,max∫Lγ,mindLγΦ(Lγ,z)dVdzdΩF1[Lγ,E(1+z),z] (1)

where is the observed photon energy, is the differential photon flux of an individual galaxy of redshift and gamma-ray luminosity , and is the comoving volume element per unit redshift and unit solid angle. In Eq. (1) we have assumed that all galaxies have identical photon spectral shapes, and that we only treat energies where absorption of gamma rays by pair production on the intergalactic background light (see, e.g., Venters 2010) is negligible. In this work, we will express the gamma-ray luminosity function in terms of the differential photon luminosity at a rest-frame energy of 0.6 GeV: . We use this same quantity to normalize the photon spectrum, , where . The photon flux from a single galaxy is then given by

 F1[Lγ,E(1+z),z]=Lγ4πξ2dNdE[E(1+z)] (2)

where is the coordinate distance that enters the Robertson-Walker metric,

 ds2=−c2dt2+a2(t)[dξ2/(1−kξ2)+ξ2dΩ2] (3)

so for a flat Universe the comoving volume element per unit solid angle and per unit redshift is so Eq. (1) simplifies to

 IE(E)=c4πH0zmax∫0dzLγ,max∫Lγ,mindLγΦ(Lγ,z)Lγ√Ωm(1+z)2+ΩΛdNdE[E(1+z)]. (4)

## 3 The Gamma-ray Luminosity Function

Since cannot be determined directly from Fermi data (because resolved star-forming galaxies are too few), it is common practice to rescale an infrared luminosity function assuming some relation between gamma-ray luminosity and infrared luminosity . A very significant fraction of the uncertainty in the calculation of the contribution of normal star-forming galaxies to the EGRB enters through the assumed relation between and . There are various approaches to deriving such a relation, which lead to different results.

### 3.1 Empirical scaling between Lγ and LIR

Ackermann et al. (2012) optimised simple power law scalings between gamma-ray luminosity (energy luminosity integrated between 100 MeV and 100 GeV, ) and total infrared luminosity, :

 Missing or unrecognized delimiter for \left (5)

In the equation above, the energy luminosity integrated between 0.1 and 100 GeV that enters the scaling is related to the differential photon luminosity at 0.6 GeV, , that enters Eq. (1) through .

Ackermann et al. (2012), depending on whether galaxies with AGN were included in (excluded from) their analysis, found a slope (), using the Expectation-Maximization (EM) algorithm [e.g., Dellaert (2002) and references therein], which is similar to the well known least-square fitting. A more sophisticated statistical analysis by Linden (2017) taking explicitly into account the possible spread in the optimised scaling, also resulted in a consistent result () for the slope. These analyses take into account a significant number of star-forming galaxies for which only upper limits, rather than statistically significant measurements, are available for their gamma-ray flux. However, most of these upper limits are weak, and we have confirmed that as a result the best-fit value for the scaling slope is naturally dominated by the effect of the (few) resolved galaxies. For this reason, in our own tests for sources of biases entering through an adopted empirical scaling, we limit ourselves to least-square fitting of power laws using the resolved sources only.

In using Eq. (5) and an infrared luminosity function to obtain the gamma-ray luminosity function, we adopt the implicit assumption that the scaling itself does not show a non-trivial redshift evolution. A few ways such a redshift dependence could enter, include, for example, the following.

• In escape-dominated galaxies, where the gamma-ray flux depends both on the gas mass and the SFR, the gas fraction of a galaxy of given mass increases with redshift, compounding the effect of increased star formation.

• The performance of infrared flux as a star-formation tracer is known to be redshift-dependent, so that infrared luminosity functions may not adequately represent the cosmic star formation history (which in turn sets the history of cosmic ray acceleration and thus star-formation-related gamma-ray emission).

• The supernova mass cutoff is metallicity-dependent, so that at high redshift / lower metallicity environments, a different fraction of stars act as supernova progenitors that will subsequently produce supernova remnants that will participate in cosmic-ray acceleration.

• The fraction of star formation that takes place in environments so gas-rich that gamma-ray production is calorimetric may increase with increasing redshift. This may impart a redshift-dependent change is the scaling slope between gamma-ray and infrared luminosities.

In this paper we quantitatively assess effects (a)-(c) above, while we qualitatively discuss the potential impact of (d) and ways it can be addressed in the future.

### 3.2 Analytically derived scaling between Lγ and LIR

An alternative approach to relating and is to assume that is proportional to the SFR, and then use our understanding of the physics of gamma-ray production in star-forming galaxies to relate the SFR to . This is the approach used by Fields et al. (2010).

For galaxies with escape-dominated cosmic-ray losses, the gamma-ray production rate per interstellar H-atom, , is proportional to the galaxy’s volume averaged cosmic-ray proton flux, (e.g., Pohl 1994; Pavlidou & Fields 2001; Persic & Rephaeli 2010), so

 Lγ∝M\scriptsize{gas}Φ\scriptsize{cr} (6)

where is the gas mass of a galaxy and its cosmic-ray flux. If in addition we assume that the initial mass function (IMF), the supernova progenitor mass cutoff, and the confinement efficiency are comparable in all galaxies, while all supernova remnants accelerate, on average, the same number of cosmic rays1,

 Φcr∝RSN∝ψ (7)

where is the supernova rate (SNR) and is the SFR. Thus,

 Lγ∝M\scriptsize{gas}ψ (8)

can also be related to . For example, the Kennicutt-Schmidt scaling (Schmidt 1959; Kennicutt 1998) relates the SFR and gas surface densities,

 Σ\scriptsize{SFR}∝Σx\scriptsize{gas} (9)

which yields (Fields et al. 2010)

 M\scriptsize{gas}(ψ,z)∝(1+z)−βψω (10)

where and . The factor enters through the conversion of surface densities of gas and SFR in the KS law to volume densities. The gamma-ray luminosity of a galaxy then will be,

 Lγ(ψ,z)∝(1+z)−βψω+1. (11)

Note that this equation is valid only for normal escape-dominated galaxies. In starburst galaxies, which have very high cosmic-ray intensities and gas within small volumes, inelastic collisions compete with and sometimes even dominate over escape to regulate cosmic-rays losses (Paglione et al. 1996; Torres et al. 2004; Stecker 2007; Thompson et al. 2007; Lacki et al. 2011; Persic & Rephaeli 2010; Wang & Fields 2016). It is important then to keep in mind that unless we explicitly exclude starburst galaxies from our calculation, we will end up overestimating the gamma-ray signal.

The total infrared luminosity is a well-established tracer of the SFR for late type galaxies (Kennicutt 1998a). The conversion proposed by Kennicutt (1998b) is the following:

 ψ1M⊙{yr}−1=ϵ1.7×10−10L8−1000μmL⊙ (12)

This conversion assumes that thermal emission of interstellar dust approximates a calorimetric measure of radiation produced by young, i.e. Myr, stellar populations. The factor depends on the assumed initial mass function (IMF). Throughout this work we use Salpeter IMF (Salpeter 1955) and we consider it unchanging through space and time. Hence, the scaling relation between gamma-ray luminosity and TIR luminosity is

 Lγ(L8−1000μm,z)∝(1+z)−β(L8−1000μmL⊙)ω+1 (13)

In order then to calculate , which enters Eq. (1), we need to (1) adopt a luminosity function, (2) determine the slopes and (either from KS scaling or empirically) and (3) determine the normalisation of the scaling in Eq. (13). The latter can be derived using, for example, data from the Milky Way, where the local cosmic-ray flux, gas mass, and gamma-ray emission are well measured, or observations of nearby galaxies resolved in gamma rays.

In the model of Ackermann et al. (2012) the gamma-ray data used in order to derive the scaling relation between the gamma-ray luminosity and TIR luminosity of a galaxy include, in addition to upper limits, eight resolved galaxies. However, two of them (NGC 4945 and NGC 1068) are galaxies with active galactic nuclei (AGN) so, we are not going to consider them in this work. Our sample of galaxies consists of four normal star-forming galaxies (SMC, LMC, MW, M31) and four starburst (NGC 253, M82, NGC 2146, and Arp 220); effects of including the starburst galaxies in deriving the slope and normalisation of the scaling will be discussed in detail below. Detections of NGC2146 (Tang et al. 2014) and Arp 220 (Peng et al. 2016) were not available at the time of the Ackermann et al. (2012) analysis.

#### Infrared Luminosity Function

We begin by considering the adopted luminosity function. Fields et al. (2010) use an H luminosity function while Ackermann et al. (2012) use the luminosity function of Rodighiero et al. (2010):

 Φ(L)dlog10(L)=Φ∗(LL∗)1−αexp[−12σ2log210(1+LL∗)]dlog10(L) (14)

where, the parameter sets the slope at the faint end, is the characteristic luminosity and is the normalisation factor.

The choice of a luminosity function in some star-formation-tracing frequency is the first point of divergence between the different models: each choice sets an (explicit) redshift dependence of the gamma-ray emissivity, even if all other elements of a model remain the same. However, this effect can be controlled and corrected in a straight-forward way, described below.

A luminosity function in some star-formation tracer yields a history of cosmic star formation rate density, once we convert luminosity to star formation rate and integrate over all luminosities. Conversely, if we decide on a “preferred” cosmic star formation history, we can adjust any luminosity function using an overall redshift-dependent multiplicative factor so that, once integrated, it reproduces the desired cosmic star formation history. This approach has the advantage that it minimizes the sensitivity to biases relevant to each tracer. Computing a star-formation history using multiple tracers effectively averages over several different indicators thus generating a more robust result.

This is the approach we demonstrate and use in this work. As our “preferred” choice we use the star-formation history as surmised by Hopkins & Beacom (2006), considering all available SFR tracers. The luminosity function of Rodighiero et al. (2010) (which we choose to use for the remainder of this work) is not in agreement with this result (for example, Rodighiero et al. 2010 find no significant redshift dependence for , while Hopkins & Beacom 2006 report a strong peak at and a sharp decline at higher ). To bring the two in agreement, we introduce a redshift-dependent dimensionless normalisation correction factor in the scaling of Eq. (12):

 ψ1M⊙{yr}−1=ϵ1.7×10−10h(z)L8−1000μmL⊙ (15)

such that, , where is the Hopkins & Beacom (2006) cosmic star-formation history, and is the star-formation rate distribution function obtained by the luminosity function of Eq. (14) and the scaling of Eq. (12). Fitting an analytic expression to , we obtain:

 Unknown environment '% (16)

This approach allows for a fair comparison between models using different luminosity functions, but also between models that are luminosity-function–based (including, e.g., the Fields et al. 2010, Stecker & Venters 2011, and Ackermann et al. 2012 models) and models that are based on the overall cosmic star formation history (such as the Pavlidou & Fields 2002 and the Ando & Pavlidou 2009 models).

The correction factor is the first redshift dependence we identify (cosmic evolution of the performance of a single star formation tracer) that is unaccounted for when using empirical scalings alone.

#### Scaling slopes ω & β

Perhaps the most puzzling discrepancy between the theoretical approach of Fields et al. (2010) and the empirical scalings of Ackermann et al. (2012) and Linden (2017) is the discrepancy in the scaling slope between and .

From a physics perspective, for escape-domintated galaxies, if is indeed proportional to the scaling slope should deviate significantly from unity to reflect the compounded effect of both star-formation ( cosmic ray accelerators flux of projectiles) and gas ( availability of targets). On first inspection, this could be the effect of including starburst galaxies in the dataset from which the slope is calculated: starbursts are expected to be calorimetric and hence to exhibit a scaling (see e.g. Wang & Fields 2016); since Ackermann et al. (2012) and Linden (2017) consider both normal and starburst galaxies when determining a best-guess slope, it is reasonable to expect a slope closer to 1 than when considering normal galaxies alone. However, this is not the case: when fitting only the normal starforming galaxies detected by Fermi (Milky Way, LMC, SMC, and M31), we find a slope of 0.9, even flatter than the slope derived for all galaxies!

If the empirical scaling for normal galaxies alone does indeed exclude a steeper slope, then this would have important implications: it could imply, for example, that confinement of cosmic rays in normal galaxies is not only variable, but SFR-dependent (with higher SFR galaxies exhibiting poorer confinement properties); or that the IMF is SFR-dependent; or that any scatter in the scaling is dependent on gas content; or finally, that the primary contribution in the flux from star-forming galaxies is leptonic rather than hadronic and thus dependent on cosmic-ray flux but not on gas.

Before such scenarios are entertained, however, we need to test the extent to which the scaling slope is well-constrained. There are two important systematic effects that could bias the best-fit slope more than what statistical uncertainties allow for.

The first concerns the choice of star formation tracer. For example, it has been suggested that the sum of near-ultraviolet (NUV) plus luminosity of each galaxy () is a better star-formation tracer than the total infrared luminosity , as the luminosity term corrects for possible existence of dust (Kennicutt & Evans 2012). Indeed, we find that for normal galaxies, the best-guess slope in the scaling is , significantly steeper than the slope obtained when using . However, because the NUV luminosity of the Milky Way cannot be measured, this is not a fair comparison, as we have to omit the Milky Way from our fit.

To remedy this, and improve the statistics of the fit (as, if we exclude the Milky Way, there are only three normal galaxies detected by Fermi), we use the following sample to compare best-fit slopes: (SMC, LMC, M31, NGC253, M82). These are all galaxies resolved by Fermi and used by Ackermann et al. (2012), regardless of their starburst status, but excluding the Milky Way as well as the ones that appear to host an active galactic nucleus. We compare best-fit slopes of the scaling between and three star-formation tracers: , , and (see Hao et al. 2011). The best-fit slopes are , , and respectively, while the statistical error on the fitted slope is in all cases . These slopes diverge more than expected from statistical fluctuations alone, and it is therefore clear that the choice of star-formation tracer affects the empirically determined scaling slope.

The second effect involves the fact that the scaling itself will have finite scatter. We are sampling this scaling with very few detected objects, while the weak upper limits in other nearby galaxies do not add significant constraints either on the slope or the scatter of the relation. In this way, we suffer from a form of “cosmic variance”: if some of the very few galaxies we did happen to be able to “draw” in the nearby universe lie in the “outskirts” of the scaling, the fitted slope may appear significantly steeper or shallower than the true one.

To test this effect, we performed Monte Carlo calculations to examine how the slope derived from a power-law fit is affected if we only fit a handful of sources drawn from a scaling with significant scatter. We assume that the scaling now takes the form [compare with Eq. (5)]

 Missing or unrecognized delimiter for \left (17)

where is the true slope (which we assume to be ), is the normalisation of the scaling, and is a random number, which we assume to be Gaussianly distributed around zero. The standard deviation of the Gaussian quantifies the amount of scatter in the scaling. We draw eight points from this scaling (equal to the number of Fermi-resolved starforming galaxies), and we fit a power law using least-square fitting. We repeat the process times.

Figure 1 shows the mean and standard deviation of fitted slopes as a function of assumed scatter (standard deviation in the distribution of ). If the scaling has a scatter of (as estimated, e.g., by Linden 2017), then the fitted slope carries a “cosmic variance” uncertainty of , in addition to any statistical uncertainty. For scatter closer to 1 dex, the spread of slopes would be .

Due to the two effects discussed above, at this stage a difference between ( dependent on ) and ( dependent on alone) cannot be confidently claimed. The scaling slope between and remains poorly constrained.

An additional constraint to the scaling slope may arise from the level of the EGRB itself. If an assumed value of the slope results in an overproduction of the background compared to its level established by Fermi (Ackermann et al. 2015), then such slopes would have to be excluded. To test whether such constraints are possible, we will examine both best-fit values of using the latest data on Fermi-resolved star-forming galaxies, as well as other possible values of between 1 (i.e., no effect of gas) and 2 (i.e., , SFRgas, maximum possible effect of gas). In each case, the value of is taken always to be consistent with that of , i.e., . Physically, this is equivalent to changing the slope of the Kennicutt-Schmidt scaling between gas and star formation (and thus this change also propagates to ). Galaxies are still assumed to be escape-dominated. A possible convergence of the slope towards 1 is because and become uncorrelated (any is possible for a given gas mass), rather calorimetry setting in.

The factor is the second, explicit, redshift dependence we identify (cosmic evolution of the average size of galactic disks) that is unaccounted for empirical scalings between and . Note that it is possible that the Kennicutt-Schmidt scaling itself evolves with redshift (see, e.g., Gnedin & Kravtsov 2010) or with gas density (e.g. Gao & Solomon 2004; Tassis 2007 and references therein); however we do not treat these effects here.

#### Normalisation of Lγ−L\scriptsize{TIR} scaling

For each slope that we will examine, we have to determine the normalisation of the scaling of Eq. (13). We do so by performing least-squares fitting on the sample of resolved star-forming galaxies that are relevant in each case (normal or normal + starburst), while assuming that the scaling slope is fixed at the desired value. We call the normalisation resulting in this way .

### 3.3 Effect of Metallicity

One of the assumptions in Fields et al. (2010) is that the ratio of cosmic-ray flux to star-formation rate is constant for all normal star-forming galaxies. In general, this will be only true on average, as the proportionality between the two depends on various properties that may be different between galaxies of the same redshift (such as the stellar IMF, the supernova progenitor mass cutoff, the acceleration efficiency in supernova remnants, or the confinement efficiency). Variations of these properties between galaxies will result in scatter in the final scaling. If in addition the average value among many galaxies any of these properties evolves with redshift, then this will introduce an extra, unaccounted-for redshift dependence in the gamma-ray luminosity function.

Here we examine the possible redshift dependence of one of these factors: the supernova progenitor mass cutoff, which depends on the metallicity, the average value of which in turn depends on redshift. We assume that the IMF retains a Salpeter functional form (Salpeter 1955), and we explore the effect of an evolving metallicity, which can alter the minimum mass of a star that undergoes a supernovae explosion. The relation between and thus becomes

 RSN∝f(Z)ψ (18)

where encodes the effects of metallicity, . In order to specify we follow Ibeling Heger (2013), who calculate the dependence of the low mass limit for making core-collapse supernovae on the initial stellar metallicity. Their main conclusion is that for a fixed IMF may be higher at than at , where . We are interested in the minimum mass required for a star to undergo a classical core-collapse event; following the notation of Ibeling Heger (2013), this mass is denoted by . The relation they suggest for the metallicity-dependent transition mass is

 Unknown environment '%' (19)

where the coefficients are the best-fit parameters of a sixth-order polynomial they use to fit their numerical results. Following standard practice, we take the lowest-mass star to be , and highest-mass star to be . Hence, the fraction of stars undergoing a core-collapse supernova will be the integral of the IMF from to , divided by the integral of the IMF for the entire range. Since in general and the Salpeter IMF scales as , will scale as

 f(Z)=⎛⎝M\scriptsize{up}′(Z)Mup′(Ztoday)⎞⎠−1.35 (20)

Following Kistler et al.(2013), we quantify the metallicity evolution with redshift using

 Z(z)=0.03×10−0.15z (21)

Hence, using Eqs. (19) and (21), we obtain .

### 3.4 Scaling including explicit redshift dependences

Combining all the effects discussed above we finally obtain the following expression for the scaling relation between the gamma-ray luminosity and the total infrared luminosity of a galaxy:

 Lγ=Lγ,0(1+z)−β[h(z)]ω+1f(z)−1.35(L8−1000μmL⊙)ω+1, (22)

where may either be an assumed value or a best-fit value from the resolved galaxy dataset, while is always obtained from the resolved galaxy dataset given a value of , and is consistent with the value of [i.e., ].

## 4 Results

In this section we calculate the collective contribution of normal starforming galaxies to the EGRB that results from Eqs. (1) and (22). We discuss the effect on this calculation of different scaling slopes and explicit dependences on redshift entering Eq.(22) and we constrain these effects by comparing our results with Fermi EGRB data from Ackermann et al. (2015).

For the photon spectrum we have used a broken power law with the break, in the source rest frame, at GeV:

 dNdE=⎧⎪ ⎪⎨⎪ ⎪⎩(E0.6GeV)−1.9,E<0.6GeV(E0.6GeV)−2.3,E≥0.6GeV. (23)

The high-energy spectral index is consistent with the average spectral index of detected star-forming galaxies (whether starburst galaxies are included in the sample or not). The existence and location of the break is consistent with observations of diffuse emission from the Milky Way (Abdo et al. 2009) and encodes the pionic origin of the signal (see, e.g., Prodanović & Fields 2004).

### 4.1 Effect of scaling slope ω+1

In Fig. 2 we show the estimated normal star-forming galaxies contribution to the EGRB, for different values of . The normalisation constant for each case is obtained using only the four normal star-forming galaxies.

The orange line corresponds to (, i.e., linear correlation between the surface density of SFR and the surface density of gas). The solid cyan line corresponds to (, i.e., the surface density of SFR and the surface density of gas are completely uncorrelated, or, equivalently, the gas mass does not enhance the gamma-ray emission of a galaxy). These two extreme scenarios set the bounds to the possible contribution of normal star-forming galaxies to the EGRB if all other assumptions in our calculation hold.

The green line in Fig. 2 is obtained using as a star formation tracer to determine empirically the slope of the scaling. From least-square fitting we obtain , i.e., . We examine this case because the NUV MIR (Mid-Infrared) luminosity is a better estimator of recent SFR than the TIR luminosity (Kennicutt & Evans 2012). In this case, however, we do not consider the Milky Way in our sample, since we cannot measure its NUV luminosity. Hence, to obtain the slope, we perform least-square fitting using the three other normal star-forming galaxies (SMC, LMC, M31). Then, requiring that the scaling has the same slope, we determine its normalisation as described in §3.2.3. The indigo line is based on the formalism of Fields et al. (2010), where it is assumed that and . This result is different from the result of Fields et al. (2010) since we have used the Rodighiero et al. (2010) luminosity function, and we have included the additional redshift dependence to ensure that the luminosity function yields a cosmic SFR history consistent with Hopkins & Beacom (2006).

For comparison, we overplot the results of the Ackermann et al. (2012) calculation (with the red dashed line) and of that of Ando & Pavlidou (2009) (with the blue dashed line), who used the SFR density as a function of redshift instead of a luminosity function. With the exception of the spectral slopes, these models are very close to our “fiducial” model (green line).

The indigo and orange lines are inconsistent with the EGRB Fermi LAT data (they over-predict the observed background). This is additional evidence that the scaling is shallower than the theoretically predicted one based on the Kennicutt-Schmidt law (i.e., that ). However, before we can conclude that these steep slopes are excluded, we have to test the sensitivity of the normalisation of the scaling, , to the number of resolved galaxies used to empirically determine its value. We do so in the next section.

We note here that the spectral shape of the unresolved emission will not be a power law for energies a few tens of GeV. The EGRB spectrum at higher energies is modified by three effects: i) because not all sources have the same spectral index, the hardest sources will dominate at the highest energies, giving the resulting spectrum upwards curvature (Stecker & Salamon 1996; Pavlidou et al. 2008; Pavlidou & Venters 2008); ii) absorption by the extragalactic background light (EBL) will eventually become important, giving the spectrum downwards curvature (e.g., Salamon & Stecker 1998; Venters, Pavlidou & Reyes 2009 and references therein); iii) electromagnetic cascades from the highest energy photons will alter the spectrum (e.g., Venters 2010). These effects are not treated in our calculation.

### 4.2 Sensitivity to number of resolved galaxies

The significant deviation (over three orders of magnitude) between predicted intensities corresponding to the extreme values of the scaling slopes in Fig. 2 has its origin in three factors: (a) the value of itself (i.e., the power to which an increase of the SFR at early times is raised to produce the corresponding rise in gamma-ray emissivity); (b) the suppression factor due to the changing galactic disc size (higher values of result to lower and thus smaller suppression at higher redshifts); (c) the change in normalisation factor , as the same dataset is fit with power laws of different slope.

It is expected that the effect of factor (c) above is disproportionally significant, due to the small number (four) of resolved normal galaxies that enter the normalisation factor calculation. It is therefore interesting to test whether (and to what extent) the deviation between extreme scenarios would decrease if more resolved galaxies were included in the fitting sample that is used to obtain .

To implement this test, we repeat the calculations of §4.1, but now we include the four starburst galaxies NGC253, M82, NGC 2146 and Arp220 in the dataset we use to determine for each value of . It should be stressed that this calculation is not representative of the contribution of starburst galaxies to the EGRB, as for starburst galaxies other important assumptions in our model will not hold (for example, the scaling slope between infrared and gamma-ray luminosities should be exactly 1, and the suppression factor would not enter the calculation).

Our results are shown in Fig. 3. Our “fiducial” case is only marginally affected, however the curves corresponding to the two limiting slope values now differ in predicted intensity by “only” 1.5 orders of magnitude. Qualitatively our results remain unchanged (the Kennicutt-Scmidt value of still significantly overpredicts the background), however the discrepancy is considerably smaller. We conclude that while our results clearly indicate that the level of the EGRB itself is quite constraining of the scaling slope between luminosities that trace star formation and gamma-ray luminosity, due to the extremely limited statistics of resolved, normal star-forming galaxies, we should be careful so as not to overestimate the significance of such constraints. An additional factor that could be alleviating these constraints is that at higher redshifts the transition to calorimetry may occur at lower star formation rates than today (see discussion in §5).

### 4.3 Effect of Metallicity

The effect of including the effect of a cosmically evolving average metallicity on the calculation of the collective intensity of unresolved star-forming galaxies is shown in Fig. 4. The blue line shows our “fiducial” calculation with the normalisation factor calculated including both normal and starburst galaxies. The red line shows the same calculation but omitting the metallicity-related factor from Eq. (22). The metallicity does not affect our result appreciably, thus the ratio of SNR-to-SFR can be indeed assumed to be constant for all redshifts, if the IMF is constant. It is plausible that an evolving IMF may affect this calculation mode significantly. However this case is more complicated to treat, since a redshift dependence of the IMF would have to be modeled, and then one would have to factor the changing IMF also in the scaling between infrared luminosity and SFR.

Note that our correction is made on an average sense and does not reflect the distribution of galaxy metallicities at a given z, nor the distribution of metallicities within a single galaxy.

### 4.4 Contributed intensity as a function of redshift - effect of the cosmic star formation history

The contribution to the collective intensity from higher redshift galaxies is shown in Fig. 5, where we plot the fraction of intensity contributed at observer-frame energy (which corresponds to the high-energy power-law part of the spectrum) by galaxies at redshifts between 0 and over total intensity contributed by galaxies at redshifts between 0 and (the maximum redshift of our integration). Solid line colours are as in Fig. 2.

For our “fiducial” model, over of the intensity at comes from and from . As expected, the contribution of higher redshifts is larger for higher values of , which reflects the effect of the increase of gas mass at higher redshifts. However the overall qualitative behaviour remains the same: the cosmic star-formation history is the factor that primarily dictates the contribution of different redshifts to the total intensity. This is further emphasised by the fact that the Ando & Pavlidou 2009 calculation (blue dashed line) also follows the general trend, although it uses very different model assumptions. The reason is that it encodes the cosmic star formation history of Hopkins & Beacom (2006), which has been enforced in all of the models represented by solid lines through the function .

As an extreme counter-example, we plot, with the red dashed line, the results for (as in Ackermann et al. 2012) but omitting the function and, unlike Ackermann et al. 2012 (who take ), taking an extreme value of . Because the luminosity function of Rodighiero et al. (2010) does not decline above and because we do not have information of the behaviour of this luminosity function at very high redshift, without the function the result is sensitively dependent on the assumed value of . As a result, depending on the choice of , there can be significant contribution to the overall estimated intensity from unphysically high redshifts (in our extreme example, over of the intensity at comes from ). Conversely, a significant contribution from higher redshifts may be missed if is taken to be too low.

That the observed spectral break occurs at GeV while it is hard-coded to be at GeV in the source rest-frame indicates that most photons come from redshifts around , consistent with the location of the peak of the Hopkins & Beacom (2006) star formation history.

## 5 Summary and Discussion

The detection of several nearby normal star-forming and starburst galaxies by Fermi LAT (Ackermann et al. 2012) has made for the first time possible to establish empirical scalings between the galaxies’ gamma-ray luminosity and their luminosity at star-formation–tracing wavelengths. This development represents an unprecedented breakthrough in our ability to estimate the contribution to the EGRB from other galaxies where cosmic rays accelerated by star-formation products (supernova remnants) interact with these galaxies’ interstellar gas and radiation fields. Such estimates are naturally becoming the standard in the field (e.g., Stecker & Venters 2011; Ackermann 2012; Linden 2017). It is therefore important to study and understand any biases that may enter these calculations, as well as identify ways that these biases can be eliminated or alleviated. This has been the scope of this work.

We have tested whether the scaling slope between gamma-ray and infrared luminosity is well constrained, assuming that the underlying correlation has finite scatter, and given that we currently sample the underlying correlation with very few resolved galaxies. We found that even in the case where all (normal + starburst) currently resolved galaxies are considered, there is a significant “cosmic variance” effect: depending on whether the resolved galaxies happen to be drawn directly from the “core” or the “edges” of the correlation, the resulting fitted slopes show an appreciable spread (for scatter of 1 dec, the spread of slopes is ).

Our estimate of the strength of this effect is a lower limit: we have not considered errors of measurement in the luminosities; and we have not treated in any way the expected transition from a steeper slope at lower luminosities to a shallower slope at higher luminosities. This should occur because lower luminosities correspond to normal, escape-dominated star-forming galaxies, where an increased star formation rate results in a multiplicative enhancement of gamma-ray luminosity through both higher cosmic-ray acceleration and more targets (gas). On the other hand, higher luminosities correspond to starburst, calorimetric galaxies, where the gamma-ray luminosity scales as the cosmic ray flux (and hence the SFR) only. The “true” correlation is therefore theoretically expected to be a broken power law. In addition, the infrared luminosity where the break occurs may be redshift-dependent: due to the higher gas abundance of higher-redshift galaxies (e.g. Magdis et al. 2012), calorimetry may set in at lower star-formation rates.

Fermi observations of the EGRB itself constrain , as, for example, the theoretically predicted value of 1.7 (derived assuming escape-dominated losses and a correlation between star formation rate and gas mass given by the Kennicutt-Schmidt relation) results to an overprediction of the observed background. However, the following factors may affect these result. First, the scaling normalisation has to be derived using normal galaxies alone, as all the Fermi-observed starbursts are calorimetric (Wang & Fields 2016) and will obey a different scaling. There are only four Fermi-detected normal galaxies, and hence the normalisation is very poorly constrained. Second, to obtain the contribution from normal galaxies alone, the luminosity functions have to only be integrated up to the luminosity where calorimetry sets in. In our calculations we have assumed that the result of Fields et al. 2010 that the contribution from calorimetric galaxies to this calculation is very small and can be neglected holds at all redshifts. However, that result assumed that calorimetry sets in at a fixed value of the star formation rate at all redshifts, which, as discussed above, is not obvious.

We have found that the cosmic evolution of the supernova mass cutoff due to an increasing average metallcity does not appreciably affect this calculation, and can safely be neglected.

On the other hand, we have found that the different and redshift-dependent performance of various star formation tracers can considerably affect the calculation. For example, although the Kennicutt (1998b) scaling between total infrared luminosity and star-formation rate holds only for optically thick galaxies, it is uniformly applied to all sources, and infrared luminosity functions have been commonly used due to their availability rather than the robustness with which they trace star formation. One possible solution to this problem, which we have applied here, is to apply a redshift-dependent correction factor to the infrared luminosity function that brings the resulting cosmic star formation history in agreement with combination estimates that use all available star formation tracers, such as the one by Hopkins & Beacom (2006).

An evolving stellar IMF is another possible bias that may affect the ERGB calculation, however we have not treated such an evolution in this work. Instead, we used a Salpeter IMF (Salpeter 1955) for all redshifts. A proper treatment of an evolving IMF should model its impact on the scaling between star-formation-tracing luminosity and SFR, as well as on the SNR - SFR scaling.

We note that while our results indicate a rather robust spectral break (a peak in a - plot) tracing the redshift of the peak in cosmic star formation history, such a feature is not visible in Fermi data. This may be an indication that gamma rays from star-forming galaxies, even if they are a significant contribution to the EGRB (Linden 2017), are not the dominant one. A contribution of normal galaxies to the EGRB low enough that the spectral peak would be “hidden” in the overall observed spectrum is clearly well within the parameter space allowed by other observational inputs in our model.

We conclude that empirical scalings between star-formation tracing luminosities and gamma-ray luminosities of star-forming galaxies, despite being a significant breakthrough that allows the construction of gamma-ray luminosity functions, do not constrain the latter enough to produce precision estimates of the normal-galaxy contribution to the EGRB. Once all possible biases are considered, these estimates are at best accurate within an order of magnitude, with the main limitation being that the starforming galaxies that can be individually resolved in gamma rays are (a) very few and (b) all part of the local Universe. Further progress towards a more accurate determination of the normal galaxy contribution to the EGRB necessarily passes through more sophisticated models incorporating all recent progress in our understanding of a cosmically-evolving star-formation process, and a unified treatment of normal and starburst galaxies, including the gradual transition to calorimetry.

## Acknowledgements

I.K. would like to thank A. Tritsis for useful suggestions and discussions. Support for this work was provided by the European Research Council under the European Unions Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 617001 and by a Marie Sklodowska-Curie RISE grant ”ASTROSTAT” (project number 691164).

### Footnotes

1. In reality, all of these quantities may show variations, but as long as they are not dependent on star formation rate or gas content of a galaxy, they will simply contribute to the scatter of the relation, or, if they happen to evolve with cosmic time, introduce a redshift dependence.

### References

1. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, Physical Review Letters, 103, 251101
2. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, Physical Review Letters, 104, 101101
3. Ackermann, M., Ajello, M., Allafort, A., et al. 2012, ApJ, 755, 164
4. Ackermann, M., Ajello, M., Albert, A., et al. 2015, ApJ, 799, 86
5. Ajello, M., Gasparrini, D., Sánchez-Conde, M., et al. 2015, ApJL, 800, L27
6. Ando, S., & Pavlidou, V. 2009, MNRAS, 400, 2122
7. Clark, G. W., Garmire, G. P., & Kraushaar, W. L. 1968, ApJL, 153, L203
8. Dellaert, F. 2002, Tech.Rep., Georgia Institute of Technology
9. Fichtel, C. E., Simpson, G. A., & Thompson, D. J. 1978, ApJ, 222, 833
10. Fields, B. D., Pavlidou, V., & Prodanović, T. 2010, ApJL, 722, L199
11. Gao, Y., & Solomon, P. M. 2004, ApJ, 606, 271
12. Gnedin, N. Y., & Kravtsov, A. V. 2010, ApJ, 714, 287
13. Hao, C.-N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124
14. Hopkins, A. M., & Beacom, J. F. 2006, ApJ, 651, 142
15. Ibeling, D., & Heger, A. 2013, ApJL, 765, L43
16. Kennicutt, R. C., Jr. 1998a, ARA&A, 36, 189
17. Kennicutt, R. C., Jr. 1998b, ApJ, 498, 541
18. Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531
19. Kistler, M. D., Stanek, K. Z., Kochanek, C. S., Prieto, J. L., & Thompson, T. A. 2013, ApJ, 770, 88
20. Lacki, B. C., Thompson, T. A., Quataert, E., Loeb, A., & Waxman, E. 2011, ApJ, 734, 107
21. Linden, T. 2017, Phys. Rev. D, 96, 083001
22. Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6
23. Makiya, R., Totani, T., & Kobayashi, M. A. R. 2011, ApJ, 728, 158
24. Paglione, T. A. D., Marscher, A. P., Jackson, J. M., & Bertsch, D. L. 1996, ApJ, 460, 295
25. Pavlidou, V., & Fields, B. D. 2001, ApJ, 558, 63
26. Pavlidou, V., & Fields, B. D. 2002, ApJL, 575, L5
27. Pavlidou, V., Siegal-Gaskins, J. M., Fields, B. D., Olinto, A. V., & Brown, C. 2008, ApJ, 677, 27-36
28. Pavlidou, V., & Venters, T. M. 2008, ApJ, 673, 114-118
29. Peng, F.-K., Wang, X.-Y., Liu, R.-Y., Tang, Q.-W., & Wang, J.-F. 2016, arXiv:1603.06355
30. Persic, M., & Rephaeli, Y. 2010, MNRAS, 403, 1569
31. Pohl, M. 1994, A&A, 287, 453
32. Prodanović, T., & Fields, B. D. 2004, Astroparticle Physics, 21, 627
33. Rodighiero, G., Vaccari, M., Franceschini, A., et al. 2010, A&A, 515, A8
34. Salamon, M. H., & Stecker, F. W. 1998, ApJ, 493, 547
35. Salpeter, E. E. 1955, ApJ, 121, 161
36. Schmidt, M. 1959, ApJ, 129, 243
37. Sreekumar, P., Bertsch, D. L., Dingus, B. L., et al. 1998, ApJ, 494, 523
38. Stecker, F. W. 1971, NASA Special Publication, 249
39. Stecker, F. W., & Salamon, M. H. 1996, ApJ, 464, 600
40. Stecker, F. W. 2007, Astroparticle Physics, 26, 398
41. Stecker, F. W., & Venters, T. M. 2011, ApJ, 736, 40
42. Strong, A. W., Moskalenko, I. V., & Reimer, O. 2004, ApJ, 613, 956
43. Tang, Q.-W., Wang, X.-Y., & Tam, P.-H. T. 2014, ApJ, 794, 26
44. Tassis, K. 2007, MNRAS, 382, 1317
45. Thompson, T. A., Quataert, E., & Waxman, E. 2007, ApJ, 654, 219
46. Torres, D. F., Reimer, O., Domingo-Santamaría, E., & Digel, S. W. 2004, ApJL, 607, L99
47. Venters, T. M., Pavlidou, V., & Reyes, L. C. 2009, ApJ, 703, 1939
48. Venters, T. M. 2010, ApJ, 710, 1530
49. Wang, X., & Fields, B. D. 2016, arXiv:1612.07290
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters