GeV OBSERVATIONS OF STAR-FORMING GALAXIES WITH THE FERMI LARGE AREA TELESCOPE

Table 3 summarizes results from the maximum likelihood analyses. Source detection significance is determined using the test-statistic (TS) value, TS ≡ −2(ln (L₀) − ln (L₁)), which compares the likelihoods of models including the galaxy under consideration (L₁) and the null hypothesis of no gamma-ray emission from the galaxy (L₀; Mattox et al. 1996). Significant (TS > 25) gamma-ray excesses above background were detected in directions coinciding with four galaxies in the sample: two prototypical starburst galaxies, M82 and NGC 253, and two starburst galaxies also containing Seyfert 2 nuclei, NGC 1068, and NGC 4945. Each of the four gamma-ray sources is associated with the corresponding galaxy from our sample according to the 2FGL catalog (Nolan et al. 2012).

We additionally obtain relatively large TS values in the directions of NGC 2146 (∼20) and M83 (∼15). These excesses do not pass the conventional threshold of TS > 25 required to claim a source detection. Lenain & Walter (2011) noted an excess in the vicinity of M83, but determined that the best-fit position was inconsistent with the galaxy location and instead proposed an association with the blazar 2E 3100. We will return to these two galaxies in Section 6.

CR-induced gamma-ray emission from galaxies is expected to be steady on the timescales of our observations. Variability is tested for the four significantly detected sources by partitioning the full observation period into 12 time intervals of ∼90 days each and performing a separate maximum likelihood fit for each of the shorter time periods. The resulting light curves are shown in Figure 1. None of the four sources show significant changes in flux over the 3 years of LAT observations. This finding is consistent with the results obtained by Lenain et al. (2010) for NGC 1068 and NGC 4945.

Zoom In Zoom Out Reset image size

Spectral energy distributions of M82, NGC 253, NGC 1068, and NGC 4945 are shown in Figure 2. Each flux measurement represents a separate maximum likelihood fit following the procedures described above, but for six logarithmically spaced energy bins. Each energy bin is further divided into five logarithmically spaced sub-bins of energy for binned likelihood analysis to maintain ∼10 bins for each power of 10 in energy. Using the power-law spectral models, the maximum likelihood photon index values for the four galaxies are in the range 2.1–2.4. For M82 and NGC 253, the extrapolated power-law fits from the LAT energy range naturally connect with spectral measurements obtained by imaging air-Cerenkov telescopes at higher energies (Acciari et al. 2009; Acero et al. 2009; Ohm et al. 2011).

Zoom In Zoom Out Reset image size

Integral flux upper limits in the range 0.1–100 GeV at the 95% confidence level are provided in Table 3 for galaxies not significantly detected by the LAT. Limits are computed using the profile likelihood technique, for which the flux of a target source is varied over a range while simultaneously maximizing the model likelihood with respect to all other parameters. One-sided 95% CL upper limits correspond to a decrease in profile likelihood, L_p, of 2Δln (L_p) = 2.71 relative to the maximum likelihood. The flux upper limits presented in Table 3 are derived using a power-law spectral model with photon index Γ = 2.2, corresponding to the typical index of the four LAT-detected starbursts. The upper limits increase by 10%–20% (less constraining) on average when assuming a photon index of 2.3.

NGC 1068 and NGC 4945 deserve special attention as galaxies with both circumnuclear starbursts and radio-quiet obscured AGNs. The gamma-ray spectra of NGC 1068 and NGC 4945 are similar to those of M82 and NGC 253, and none of the four galaxies show indications of gamma-ray variability. Lenain et al. (2010) suggested a model for NGC 1068 in which gamma rays are produced mainly through inverse Compton scattering of IR photons by electrons in the misaligned jet at distances ∼0.1 pc from the central engine. In this scenario, NGC 1068 is viewed as the first member of a new class of gamma-ray sources.

Studies of X-ray-selected Seyfert galaxies using LAT data demonstrate that Seyfert galaxies hosting radio-quiet AGNs are generally gamma-ray quiet as a population (Teng et al. 2011; Ackermann et al. 2012b). Excepting NGC 1068 and NGC 4945, no other radio-quiet Seyfert galaxies have been detected by the LAT. The upper limits in 0.1–100 GeV luminosity inferred for several other nearby Seyfert galaxies are below the gamma-ray luminosity of NGC 1068 (Ackermann et al. 2012b). The tightest constraint is found for NGC 4151, located at a distance of 11.2 Mpc, with gamma-ray luminosity <2 × 10⁴⁰ erg s⁻¹ (a factor ∼10 below that of NGC 1068).

The relative contributions of CR interactions versus AGN activity to the gamma-ray emission of NGC 1068 and NGC 4945 are not yet definitively established. In the multiwavelength analysis which follows, we will consider both the full sample of analyzed galaxies, and a subsample with galaxies hosting Swift-BAT-detected AGNs removed.

Many authors have proposed scaling relationships between galaxy SFRs and gamma-ray luminosities motivated by the connections between interstellar gas, star formation, SNe, and CRs (e.g., Pavlidou & Fields 2002; Torres et al. 2004; Thompson et al. 2007; Stecker 2007; Persic & Rephaeli 2010; Lacki et al. 2011). We now examine relationships between gamma-ray luminosity and several photometric tracers of star formation taking advantage of the 3 year accumulation of LAT data.

The radiative output of star-forming galaxies across much of the electromagnetic spectrum is related to the abundance of short-lived massive stars. Many photometric estimators of the recent star formation histories of galaxies have been used. Total IR luminosity 8–1000 μm is one well-established tracer of the SFR for late-type galaxies (reviewed by Kennicutt 1998a). The conversion proposed by Kennicutt (1998b),

$$\begin{equation} \frac{\rm SFR}{M_{\odot }\,{\rm yr}^{-1}}=\epsilon 1.7 \times 10^{-10} \frac{L_{8\hbox{--}1000 \,\mu {\rm m}}}{L_{\odot }}, \end{equation} \tag{ 1 }$$

assumes that thermal emission of interstellar dust approximates a calorimetric measure of radiation produced by young (10–100 Myr) stellar populations. All SFRs quoted in this work consider the stellar mass range 0.1–100 M_☉. The factor depends on the assumed initial mass function (IMF), with = 1 for the Salpeter (1955) IMF originally used by Kennicutt (1998b). We will use = 0.79 to convert⁵⁷ to the Chabrier (2003) IMF, proposed after further studies of star formation in the 0.1–1 M_☉ mass range.

Several other photometric SFR estimators have been calibrated using the Kennicutt (1998b) total IR luminosity relation. These include RC luminosity at 1.4 GHz produced by synchrotron-emitting CR electrons (Yun et al. 2001),

$$\begin{equation} {\rm SFR}(M_{\odot }\,{\rm yr}^{-1})= \epsilon (5.9\pm 1.8) \times 10^{-22} L_{1.4 \,{\rm GHz}} ({\rm W\:Hz}^{-1}), \end{equation} \tag{ 2 }$$

and HCN J = 1–0 line luminosity indicating the quantity of dense molecular gas available to form new stars (Gao & Solomon 2004b),

$$\begin{equation} {\rm SFR}(M_{\odot }\,{\rm yr}^{-1})= \epsilon 1.8 \times 10^{-7} L_{{\rm HCN}} ({\rm K\,km\,s^{-1}\,pc^{2}}). \end{equation} \tag{ 3 }$$

Although these three SFR estimators are intrinsically linked, each explores a different stage of stellar evolution and is subject to different astrophysical and observational systematic uncertainties.

Figures 3 and 4 compare the gamma-ray luminosities of galaxies in our sample to their differential luminosities at 1.4 GHz, and total IR luminosities (8–1000 μm), respectively. A second abscissa axis has been drawn on each figure to indicate the estimated SFR corresponding to either RC or total IR luminosity using Equations (2) and (1). The upper panels of Figures 3 and 4 directly compare luminosities between wavebands, whereas the lower panels compare luminosity ratios. Taken at face value, the two figures show a clear positive correlation between gamma-ray luminosity and SFR, as has been reported previously in LAT data (see in this context Abdo et al. 2010b). However, sample selection effects, and galaxies not yet detected in gamma rays must be taken into account to properly determine the significance of the apparent correlations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We test the significances of multiwavelength correlations using the modified Kendall τ rank correlation test proposed by Akritas & Siebert (1996). This method is an example of "survival analysis" techniques, suitable for the analysis of partially censored data sets (i.e., containing a mixture of detections and upper limits). The Kendall τ coefficient, τ ∈ [ − 1, 1], indicates the degree of positive or negative correlation between two quantities by comparing the ordering of each pair of points in the data set. For example, τ = 1 describes a set of uncensored data points which are monotonically increasing. This non-parametric approach makes no assumption regarding the particular mathematical form for the relationship between compared quantities. See Appendix A for a detailed description of the method.

The multiwavelength relationships are tested in luminosity space because we are primarily interested in the intrinsic galaxy properties, and importantly, because multiwavelength comparisons in flux space can either falsely suggest a non-existent correlation or obscure a genuine physical correlation if the underlying relationship between intrinsic luminosities is nonlinear (Feigelson & Berg 1983; Kembhavi et al. 1986; and see Appendix B). We performed a series of Monte Carlo simulations to investigate potential biases resulting from selection effects for flux-limited samples of objects compared in luminosity space. The generalized Kendall τ correlation test is found to be robust against the selection effects expected for our application. A brief discussion of results from that study is presented in Appendix C.

Significances of the multiwavelength relationships are computed by comparing the τ correlation coefficient of the actual data to the distribution of τ correlation coefficients which could be obtained under the null hypothesis of independence between wavebands. Null hypothesis data sets are generated by scrambling derived gamma-ray luminosities among galaxies. Only "observable" permutations are retained to account for the truncation of measurable gamma-ray luminosities imposed by the LAT flux sensitivity (following the method of Efron & Petrosian 1999). Specifically, we require that the resultant gamma-ray flux (or upper limit) of each galaxy in a scrambled data set could have been measurable by the LAT in order for the scrambled data set to be included in the null hypothesis distribution. In practice, each permutation represents an exchange of gamma-ray luminosities between two randomly selected galaxies so that statistically meaningful null hypothesis distributions can be efficiently generated. Gamma-ray flux sensitivity thresholds are determined by randomly sampling from the distribution of gamma-ray flux upper limits presented in Table 3 for each galaxy in the pair to be exchanged. We perform 1000 such exchanges on the actual data before counting permutations toward the null hypothesis distribution.

Inaccuracies in the distance measurements to the galaxies are included in our analysis because these errors are propagated in the transformation of measured fluxes into luminosities. Scatter in distance measurements applied to a collection of objects tends to broaden the distribution of objects in luminosity space, and can therefore induce positive correlation. For the galaxies in our sample at distances 1 Mpc ≲ D ≲ 100 Mpc, the typical distance measurement precision for individual galaxies is 10%–20% (Freedman et al. 2001, and references therein). For our particular sample, we can compare the distances provided in the HCN survey of Gao & Solomon (2004a) to those reported in the IRAS Revised Bright Galaxies Sample (Sanders et al. 2003), and find that for galaxies beyond the Local Group, the average discrepancy between reported distances is 10%. We assume that measurement errors are normally distributed and centered on the true distances to the galaxies.

Figure 5 shows an example of the correlation significance test applied to the relationship between gamma-ray and IR luminosity using the full sample of star-forming galaxies. The distribution of correlation coefficients for 10⁶ "observable" permutations of the actual data is plotted for three levels of distance measurement scatter. The null hypothesis distributions tend toward positive values due to the truncation of measurable gamma-ray luminosities, effectively requiring a larger correlation coefficient for the actual data in order to claim a significant relationship between wavebands. As the level of scatter in distance measurements is increased, the tails of null hypothesis distributions extend to larger positive values. The probability that a null hypothesis data set would have a correlation coefficient larger than that of the actual data corresponds to the P-value measure of correlation significance.

Zoom In Zoom Out Reset image size

Given uncertainties for the measured gamma-ray fluxes of galaxies in our sample, the correlation coefficient of the actual data is more appropriately viewed as a probability density as opposed to a single value (the Kendall τ statistic does not explicitly include measurement errors). To account for this uncertainty, we draw a flux value for each LAT-detected galaxy from a normal distribution centered on the best-fit flux with standard deviation equal to the reported flux uncertainty. A correlation coefficient is then computed using the sampled gamma-ray fluxes. This process is repeated many times to generate a distribution of correlation coefficients representing the actual data, as shown by the gray bands in Figure 5. Finally, we integrate over the probability density of correlation coefficients representing the actual data to obtain significances for the multiwavelength relations.

The results of the correlation tests are summarized in Table 4. Using the full sample of 69 star-forming galaxies, the null hypothesis of independence between gamma-ray luminosity and either RC or total IR luminosity can be rejected at the ≳99% confidence level allowing for 20% uncertainty on the distance measurements to the galaxies. No evidence for a correlation between gamma-ray luminosity and HCN-line luminosity is found in this study, an apparent discrepancy considering the previously established correlations between RC, IR, and HCN-line luminosities (Liu & Gao 2010). However, global HCN-line luminosity estimates are not available for the SMC and LMC, M31, and M33, which effectively reduces the luminosity range over which the correlation between wavebands can be tested, and excludes three LAT-detected galaxies from the analysis.

We also considered the more conservative case in which galaxies hosting Swift-BAT-detected AGNs, including NGC 1068 and NGC 4945 along with seven others, were removed from the sample given the potential contributions of the AGN to the broadband emission of those galaxies. In that case, the P-values for the correlations between gamma-ray luminosity and either RC or IR luminosity are ∼0.05 allowing for 20% uncertainty on the distance measurements.

Even in the most conservative case described above, the data prefer a correlation between gamma-ray luminosity and both RC and IR luminosity. Although strong conclusions regarding the significance of the multiwavelength correlations cannot be made at the time, at minimum, scaling relations derived from the current sample have predictive value and offer a testable hypothesis for further observations.

We fit scaling relationships between wavebands using simple power-law forms. For example,

$$\begin{equation} \log \left(\frac{L_{0.1\hbox{--}100 \,{\rm GeV}}}{{\rm erg \,s}^{-1}}\right) = \alpha \log \left(\frac{L_{8\hbox{--}1000 \,\mu {\rm m}}}{10^{10} L_\odot }\right) + \beta \end{equation} \tag{ 4 }$$

parameterizes the relationship between gamma-ray and total IR luminosity.

Two regression methods are employed: the expectation-maximization (EM) algorithm, and the Buckley–James algorithm. The EM algorithm is similar to the least-squares fitting method, and assumes that the intrinsic residuals of gamma-ray luminosity are normally distributed in log space about the regression line for fixed values of either RC or IR luminosity. Non-detections are incorporated in the regression analysis by determining the degree to which the upper limits are compatible with the assumed dispersion about the regression line. The variance of the intrinsic residuals becomes a third parameter in the fit. The Buckley–James algorithm is similar in approach, but generalizes to cases for which the intrinsic distribution of residuals relative to the regression line is not normally distributed. Instead, the intrinsic dispersion is estimated using the Kaplan–Meier distribution obtained from the scatter within the data set itself. Isobe et al. (1986) provide details regarding the theory and implementation several commonly used survival analysis methods including the EM and Buckley–James regression algorithms.

Table 5 reports the best-fits of gamma-ray luminosity versus either RC luminosity or total IR luminosity, which were evaluated using the ASURV Rev 1.2 code (LaValley et al. 1992).⁵⁸ The two regression methods yield consistent results. For both the RC and total IR tracers of SFR, a nearly linear power-law index of α = 1.0–1.2 is preferred. The scaling indices between gamma-ray luminosity and SFR tracers found for the larger sample of galaxies is consistent within uncertainties with that obtained considering the Local Group galaxies alone, α = 1.4 ± 0.3 (Abdo et al. 2010b). In particular, the gamma-ray luminosity upper limits found for LIRGs and ULIRGs disfavor a strongly nonlinear scaling relation. Table 5 also contains regression parameters fitted after removing galaxies hosting Swift-BAT-detected AGNs from the sample, which would be appropriate if the broadband emissions of those galaxies were heavily influenced by AGN activity. The fitted parameters in both cases are consistent within statistical uncertainties.

Table 5. Scaling Relationships between Global Luminosities: Power-law Fits with Full Sample

	Expectation-maximization Method			Buckley–James Method
	α	β	$\sqrt{{\rm Variance}}$	α	β	$\sqrt{{\rm Variance}}$
Full sample of galaxies
L1.4 GHz: L0.1–100 GeV	1.10 ± 0.05	38.82 ± 0.06	0.17	1.10 ± 0.06	38.81	0.20
(1021 W Hz−1): (erg s−1)
L8–1000 μm: L0.1–100 GeV	1.17 ± 0.07	39.28 ± 0.08	0.24	1.18 ± 0.10	39.31	0.31
(1010 L☉): (erg s−1)
Excluding galaxies hosting Swift-BAT-detected AGN
L1.4 GHz: L0.1–100 GeV	1.10 ± 0.07	38.81 ± 0.07	0.19	1.09 ± 0.11	38.80	0.24
(1021 W Hz−1): (erg s−1)
L8–1000 μm: L0.1–100 GeV	1.09 ± 0.10	39.19 ± 0.10	0.25	1.10 ± 0.14	39.22	0.33
(1010 L☉): (erg s−1)

Notes. Fitted parameters for relationships between gamma-ray luminosity and multiwavelength tracers of star formation. Using the RC case as an example, the scaling relations are of the form log L_{0.1–100 GeV} = αlog L_1.4 GHz + β, with luminosities expressed in the units provided in the leftmost column. The square root of the variance provides an estimate of the intrinsic dispersion of gamma-ray luminosity residuals in log space about the best-fit regression line. For the EM algorithm, the intrinsic residuals about the best-fit line are assumed to be normally distributed in log space. The Buckley–James algorithm uses the Kaplan–Meier method to estimate the distribution of residuals. See Section 4.3 for a description of the fitting methods. Both the complete sample of 69 galaxies (containing eight LAT sources) and a subsample of 60 galaxies excluding the AGN detected by the Swift BAT are analyzed (containing six LAT sources).

Download table as:  ASCIITypeset image

Note that the intrinsic dispersion values presented in Table 5 should be viewed as upper limits because we have not made an attempt to account for measurement uncertainties in the gamma-ray fluxes.

The best-fit power laws obtained using the EM algorithm are plotted in Figures 3 and 4. The darker shaded regions represent uncertainty in the fitted parameters from the regression and the lighter shaded regions indicate the addition of intrinsic residuals to one standard deviation. None of the limits from non-detected galaxies are in strong conflict with the best-fit relations.

Luminosity ratios of the form shown in the lower panels of Figures 3 and 4 have a direct physical meaning in terms of the energy radiated in different parts of the electromagnetic spectrum. For a galaxy with non-thermal emission powered mainly by CR interactions, and having an SFR of 1 M_☉ yr⁻¹ (similar to the MW), the corresponding luminosity ratios between wavebands are

$$\begin{equation} \log \left(\frac{L_{0.1\hbox{--}100 \,{\rm GeV}}}{L_{1.4 \,{\rm GHz}}}\right) = 1.7 \pm 0.1_{\rm (statistical)} \pm 0.2_{\rm (dispersion)}, \end{equation} \tag{ 5 }$$

and

$$\begin{equation} \log \left(\frac{L_{0.1\hbox{--}100 \,{\rm GeV}}}{L_{8\hbox{--}1000 \,\mu {\rm m}}}\right) = -4.3 \pm 0.1_{\rm (statistical)} \pm 0.2_{\rm (dispersion)}, \end{equation} \tag{ 6 }$$

found using the full sample of galaxies with the EM regression algorithm.

Further gamma-ray observations are required to conclusively establish the multiwavelength correlations examined in this section. Additional data may also show that scaling relations beyond simple power-law forms are required once farther and more active sources are either detected or are constrained by more stringent gamma-ray upper limits.

A discussion of scaling relations for diffuse gamma-ray radiation from star-forming galaxies invites comparison with the extensively studied RC–IR correlation. The RC–IR correlation extends over multiple galaxy types including irregulars, spirals, and ellipticals with ongoing star formation (de Jong et al. 1985; Wunderlich et al. 1987; Dressel 1988; Wrobel & Heeschen 1988, 1991), covers galaxies with magnetic energy densities likely ranging at least four orders of magnitude (Condon et al. 1991; Thompson et al. 2006), and applies to galaxies out to at least z ∼ 1 (Appleton et al. 2004). High-resolution imaging has also shown the correlation to hold on ∼100 pc scales within individual galaxies (Marsh & Helou 1995; Paladino et al. 2006; Murphy et al. 2006).

A simple model proposed to explain the RC–IR correlation posits that galaxies are effective "calorimeters" of both UV photons and CR electrons (Völk 1989). In this picture, if the ISM is optically thick to UV photons, reprocessed starlight emitted in the IR becomes a proxy for the abundance of massive stars, and hence SFR. The luminosity of CR electrons is assumed to be proportional to the SFR in the basic calorimeter model. If either synchrotron losses dominate for CR electrons, or if the ratio of synchrotron losses to other loss mechanisms is uniform for many galaxy types, then the RC luminosity should also increase together with the SFR. The latter possibility requires some degree of fine-tuning (but see Lacki et al. 2010). Numerical models of CR propagation in the MW suggest that our Galaxy is an effective CR electron calorimeter (Strong et al. 2010).

Gamma-ray observations extend our study of CRs to include hadronic populations beyond the solar system, and notably, to measure the large-scale propagated spectra of CR nuclei in the MW and other galaxies. At energies above the threshold for pion production, hadronic gamma rays have the same spectral index as the underlying CR proton spectrum in the thin-target regime relevant for proton–proton interactions in the ISM. Diffuse pionic gamma rays presumably offer direct access to the majority of CR energy content, although leptonic diffuse emission and individual sources within galaxies must also contribute to the total emission. The gamma-ray spectra of the SMC and LMC are both consistent with models in which pionic gamma rays dominate the high-energy emission (Abdo et al. 2010e, 2010h). Many authors have argued that the majority of starburst emission at energies >0.1 GeV is produced by interactions of CR nuclei in the ISM, assuming primary CR proton to electron ratios similar to values found in the MW (e.g., Paglione et al. 1996; Torres 2004; Persic et al. 2008; Lacki et al. 2010).

For the MW in particular, global CR-induced gamma-ray emission above 0.1 GeV stems mostly from neutral pion decay (Bloemen 1985; Bertsch et al. 1993), with leptonic processes contributing 30%–40% in the energy range 0.1–100 GeV (Strong et al. 2010). The diffuse neutral pion decay component for the MW has spectral index ∼2.75 above ∼1 GeV, matching the locally measured CR proton spectrum (Simpson 1983; Sanuki et al. 2000; Adriani et al. 2011). CR escape from the MW plays a major role in the formation of the measured CR spectrum, but other galaxies with higher concentrations of ambient gas may work differently, behaving as calorimeters of CR nuclei, as suggested by Thompson et al. (2007). Approximately, 10% (possibly up to 20%) of the global gamma-ray emission of the MW is thought to originate from individual sources (Strong 2007).

We now review models for the interactions of CR nuclei to examine physical quantities involved in the scaling of hadronic gamma-ray emission, following the formalism of recent works by Persic & Rephaeli (2010) and Lacki et al. (2011). Neutral pion decay emission is related both to the energy density of CR nuclei, U_p, and to the ambient gas density, n, which serves as target material for inelastic collisions. The total hadronic gamma-ray luminosity (erg s⁻¹) can be expressed as

$$\begin{equation} L_{\gamma,\pi ^{0}} = \int _V\int E_{\gamma }\frac{dq}{dE_{\gamma }}\;n(\textbf {\em r})\;U_{\rm p}(\textbf {\em r})dE_{\gamma }dV, \end{equation} \tag{ 7 }$$

integrating over gamma-ray energies, E_γ, and over the whole galactic volume, V, which comprises the disk and surrounding halo occupied by CRs. The factor q(⩾ E_γ) denotes gamma-ray emissivity (s⁻¹ H-atom⁻¹ eV⁻¹ cm³) normalized to the CR energy density (Drury et al. 1994). Next, we define an average effective density of ambient gas encountered by CR nuclei:

$$\begin{equation} \langle n_{\rm eff}\rangle = \frac{\int _V n(\textbf {\em r})U_{\rm p}(\textbf {\em r}) dV}{\int _V U_{\rm p}(\textbf {\em r}) dV}. \end{equation} \tag{ 8 }$$

Note that 〈n_eff〉 can be substantially lower than the average gas density found in the galaxy disk. CR nuclei of the MW are thought to spend most of their time in the low-gas-density halo, based upon combined measurements of unstable isotopes and spallation products (e.g., the B/C ratio) in the locally observed CR composition (reviewed by Strong et al. 2007). Using the above definition for the average effective gas density, Equation (7) separates to

$$\begin{equation} L_{\gamma,\pi ^{0}} = \int _V U_{\rm p} dV \times \langle n_{\rm eff}\rangle \times \int E_{\gamma }\frac{dq}{dE_{\gamma }} dE_{\gamma }. \end{equation} \tag{ 9 }$$

The leading factor of Equation (9) represents the total energy of CR nuclei, and can be re-written as the product of the total luminosity of CR nuclei, L_p, and their characteristic residence time in the galactic volume, τ_res:

$$\begin{equation} \int _V U_{\rm p} dV = L_{\rm p}\;\tau _{\rm res}. \end{equation} \tag{ 10 }$$

The trailing factor of Equation (9) depends upon the parent spectrum of CR nuclei, for which we assume a power-law spectral distribution, $dN_{\rm p}/dP_{\rm p} \propto (P_{\rm p} \,{\rm GeV}^{-1}c)^{-\Gamma_{\rm p}}$, where P_p is the momentum. Gamma-ray emissivities (s⁻¹ H-atom⁻¹ cm³) in the energy range 0.1–100 GeV,

$$\begin{equation} \mathcal {Q}_{0.1\hbox{--}100 \,{\rm GeV}} = \int _{0.1 \,{\rm GeV}}^{100 \,{\rm GeV}} E_{\gamma }\frac{dq}{dE_{\gamma }} dE_{\gamma }, \end{equation} \tag{ 11 }$$

are computed using the proton–proton interaction cross section parameterizations of Kamae et al. (2006). Note that the emissivities are normalized to the average energy density of CR nuclei (kinetic energies 1–10⁶ GeV) within the entire galactic volume. A nuclear enhancement factor of 1.85 is included to account for heavy nuclei in both CRs and target matter (Mori 2009). The emissivities for Γ_p = 2–3 range from $\mathcal {Q}_{0.1\hbox{--}100 \,{\rm GeV}} = (5.7\hbox{--}7.9) \times 10^{-17}$ s⁻¹ H-atom⁻¹ cm³.

Substituting Equations (10) and (11) into Equation (9) yields

$$\begin{equation} L_{0.1\hbox{--}100 \,{\rm GeV},\pi ^{0}} = L_{\rm p}\;\tau _{\rm res}\;\langle n_{\rm eff}\rangle \;\mathcal {Q}_{0.1\hbox{--}100 \,{\rm GeV}}. \end{equation} \tag{ 12 }$$

In the paradigm that SNRs are the primary sources of galactic CRs (Ginzburg & Syrovatskii 1964), CR luminosity is the product of the SN rate, Γ_SN, the kinetic energy released per SN, E_SN, and the fraction of kinetic energy going into CR nuclei with kinetic energies >1 GeV (i.e., acceleration efficiency), η:

$$\begin{equation} L_{\rm p} = \Gamma _{\rm SN}\;E_{\rm SN}\;\eta. \end{equation} \tag{ 13 }$$

Core-collapse SN rates can be estimated from SFRs given an IMF and a minimum stellar mass required to produce a core-collapse SN, which we take to be 8 M_☉. For a Chabrier (2003) IMF, the transformation between SFR in the mass range 0.1–100 M_☉ and core-collapse SN rate is Γ_SNΨ = SFR, with Ψ = 83 M_☉. The average kinetic energy released per core-collapse SN is E_SN ∼ 10⁵¹ erg (Woosley & Weaver 1995). The gamma-ray luminosity becomes

$$\begin{equation} L_{0.1\hbox{--}100 \,{\rm GeV},\pi ^{0}} = \Psi ^{-1}\,{\rm SFR}\;E_{\rm SN}\;\eta \;\tau _{\rm res}\;\langle n_{\rm eff}\rangle \;\mathcal {Q}_{0.1\hbox{--}100 \,{\rm GeV}}. \end{equation} \tag{ 14 }$$

Greater physical intuition can be realized by inserting values for the physical quantities similar to those of the MW,

$$\begin{eqnarray} L_{0.1\hbox{--}100 \,{\rm GeV},\pi ^{0}} &=& 9.4 \times 10^{38}\,{\rm erg\,s}^{-1} \left(\frac{\rm SFR}{M_{\odot }\,{\rm yr}^{-1}}\right) \left(\frac{\mathit {E}_{\rm SN}}{10^{51}\,{\rm erg}}\right) \nonumber\\ &&\times\left(\frac{\eta }{0.1}\right) \left(\frac{\tau _{\rm res}}{10^7\,{\rm yr}}\right) \left(\frac{\langle n_{\rm eff}\rangle }{{\rm cm}^{-3}}\right), \end{eqnarray} \tag{ 15 }$$

assuming a power-law index for CR nuclei of Γ_p = 2.75 corresponding to $\mathcal {Q}_{0.1\hbox{--}100 \,{\rm GeV}} = 7.8 \times 10^{-17}$ s⁻¹ H-atom⁻¹ cm³.

A consistency check for this simple approach can be performed for the specific case of the MW, for which the column density of material traversed by CR nuclei, or "grammage," is well known (e.g., Dogiel et al. 2002). The mean escape length for CR nuclei in the MW is $\bar{x} = \tau _{\rm res}\;m_{\rm p}\;\langle n_{\rm eff}\rangle \;c \approx 12$ g cm⁻² (Jones et al. 2001; cf. the interaction length is ∼55 g cm⁻²). We may equivalently express Equation (15) in terms of the escape length as

$$\begin{eqnarray} L_{0.1\hbox{--}100 \,{\rm GeV},\pi ^{0}} &=& 6.0 \times 10^{38}\,{\rm erg\,s}^{-1} \left(\frac{\rm SFR}{M_{\odot }\,{\rm yr}^{-1}}\right) \left(\frac{\mathit {E}_{\rm SN}}{10^{51}\,{\rm erg}}\right) \nonumber\\ &&\times \left(\frac{\eta }{0.1}\right) \left(\frac{\bar{x}}{10\,{\rm g\;cm}^{-2}}\right). \end{eqnarray} \tag{ 16 }$$

If we take SFR = 1.6 M_☉ yr⁻¹ corresponding to an SN rate of 1.9 (± 1.1) per century (Diehl et al. 2006), then our estimate for the pionic gamma-ray luminosity of the MW becomes

$$\begin{equation} L_{0.1\hbox{--}100 \,{\rm GeV},\pi ^{0}} = 1.1 \times 10^{39}\,{\rm erg\,s}^{-1} \left(\frac{\mathit {E}_{\rm SN}}{10^{51}\,{\rm erg}}\right) \left(\frac{\eta }{0.1}\right). \end{equation} \tag{ 17 }$$

For comparison, the luminosity estimated from a full numerical treatment (GALPROP code) incorporating an ensemble of multiwavelength and CR measurements is $L_{0.1\hbox{--}100 \,{\rm GeV},\pi ^{0}} = (5\hbox{--}7) \times 10^{38}$ erg s⁻¹, with E_SNη = (0.3–1) × 10⁵⁰ erg (Strong et al. 2010).

It is commonly assumed that the kinetic energy released per SN, and acceleration efficiency are universal from galaxy to galaxy. In that case, the model expressed in Equations (15) and (16) implies that the hadronic gamma-ray luminosity of the ISM scales linearly with the SFR, and with the product of residence time and effective ambient gas density (or equivalently, grammage). These quantities are likely correlated in real galaxies. For example, several authors have argued that SFRs and gas densities are related (e.g., via the Schmidt–Kennicutt law) so that a nonlinear scaling of gamma-ray luminosities with SFRs should be expected: $L_{0.1\hbox{--}100 \,{\rm GeV},\pi ^{0}} \propto {\rm SFR}^{1.4\hbox{--}1.7}$ (Persic & Rephaeli 2010; Fields et al. 2010).

In general, the residence time of CR nuclei in a galaxy, τ_res, can be expressed as

$$\begin{equation} \tau _{\rm res}^{-1}=\tau _{\rm esc}^{-1}+\tau _{\rm pp}^{-1}. \end{equation} \tag{ 18 }$$

A galaxy becomes a "calorimeter" of CR nuclei when the residence time approximately equals the collisional energy loss timescale, τ_res ≈ τ_pp, i.e., energy losses are dominated by inelastic collisions with interstellar matter rather than by escape of energetic particles. The proton–proton collisional energy loss timescale depends primarily on the average gas density encountered by CR nuclei, and is nearly independent of energy because the inelastic cross section, σ_pp, is only weakly energy dependent for CRs with kinetic energies greater than ∼1 GeV. The collisional energy loss timescale is (Mannheim & Schlickeiser 1994)

$$\begin{equation} \tau _{\rm pp} \approx (0.65\langle n_{\rm eff}\rangle \;\sigma _{\rm pp}\;c )^{-1} \approx 5 \times 10^7\,{\rm yr}\;\left(\frac{\langle n_{\rm eff}\rangle }{{\rm cm}^{-3}}\right)^{-1}. \end{equation} \tag{ 19 }$$

By contrast, the escape energy loss timescale,

$$\begin{equation} \tau _{\rm esc}^{-1}=\tau _{\rm dif}^{-1}+\tau _{\rm adv}^{-1}, \end{equation} \tag{ 20 }$$

may be energy dependent since the diffusion time, τ_dif, depends upon particle rigidity in interstellar magnetic fields. For the MW, the main contribution to the gamma-ray luminosity comes from photons with energies ∼1 GeV, corresponding to CR proton energies of ∼10 GeV. The diffusive escape timescale can be estimated as τ_dif ∼ z²_h/6D ∼ (1–2) × 10⁷ yr for a halo height of z_h = 4 kpc, and taking diffusion coefficients in the range D = (4–8) × 10²⁸ cm² s⁻¹ for a particle rigidity of ∼10 GV (Ptuskin et al. 2006). Given a grammage of 12 g cm⁻³, the average gas density corresponding to an escape time of 2 × 10⁷ yr is 〈n_eff〉 ≈ 0.3 cm⁻³ (implying that τ_pp ∼ 1.5 × 10⁸ yr). In spite of modeling differences, it is generally understood that the diffusive escape timescale is much shorter than the collisional loss timescale for the MW (Strong et al. 2007, and references therein).

Large uncertainties exist for the escape time due to uncertainty in diffusion parameters, and because CRs might also be advected outward by bulk motion "superwinds" which are observed for many starbursts (Lehnert & Heckman 1996). The advective timescale, τ_adv, like the collisional energy loss timescale, is nearly independent of energy.

The spectrum of pionic gamma rays will be affected by a transition between dominant energy loss mechanisms, becoming softer for galaxies in which energy-dependent diffusive losses are important. In particular, the harder gamma-ray spectra found for starburst galaxies may be understood if the energy loss rate due to either proton–proton interactions or advection is considerably faster than the diffusion timescale. A scaling relation for gamma-ray luminosity should include the possibility that residence times for CR nuclei in different galaxies may vary substantially.

In the calorimetric limit, the residence time and effective density of target material are inversely proportional, and accordingly, the anticipated scaling of pionic gamma-ray luminosity becomes a linear scaling of the SFR:

$$\begin{eqnarray} L_{0.1\hbox{--}100 \,{\rm GeV},\pi ^{0}}\vert _{{\tau _{\rm res}} \approx \tau _{\rm pp}} &=& 5 \times 10^{39}\,{\rm erg\,s}^{-1} \left(\frac{\rm SFR}{M_{\odot }\,{\rm yr}^{-1}}\right)\nonumber\\ &&\times \left(\frac{\mathit {E}_{\rm SN}}{10^{51}\,{\rm erg}}\right) \left(\frac{\eta }{0.1}\right). \end{eqnarray} \tag{ 21 }$$

For the calorimeter case, we assume an underlying spectrum of CR nuclei described by a power law with Γ_p = 2.2 ($\mathcal {Q}_{0.1\hbox{--}100 \,{\rm GeV}} = 7.9 \times 10^{-17}$ s⁻¹ H-atom⁻¹ cm³), representative of the LAT-detected starbursts. The gamma-ray luminosity expected in the calorimetric limit for CR nuclei is plotted in Figures 3 and 4 assuming an average CR luminosity per SN of E_SN η = 10⁵⁰ erg.⁵⁹

To determine how well the theoretical calorimetric limit above describes nearby starburst galaxies, we can transform the observed scaling relation between gamma-ray luminosity and total IR luminosity to a scaling relation between gamma-ray luminosity and SFR using Equation (1). This produces

$$\begin{equation} L_{0.1\hbox{--}100 \,{\rm GeV}}\vert _{\rm scaling\;relation} = N \left(\frac{\rm SFR}{M_{\odot }\,{\rm yr}^{-1}}\frac{1}{1.7 \epsilon }\right)^{\alpha }, \end{equation} \tag{ 22 }$$

where N = 1.85^+0.37_{− 0.31} × 10³⁹ erg s⁻¹ and α = 1.16 ± 0.07, fit using the EM algorithm for the complete sample of 69 galaxies. Substituting best-fit values for the normalization and index, and assuming a Chabrier IMF ( = 0.79) yields

$$\begin{eqnarray} L_{0.1\hbox{--}100 \,{\rm GeV}}\vert _{\rm scaling\;relation} &\approx& (1.3 \pm 0.3) \times 10^{39}\,{\rm erg\,s}^{-1} \nonumber\\ &&\times\left(\frac{\rm SFR}{M_{\odot }\,{\rm yr}^{-1}}\right)^{1.16 \pm 0.07}. \end{eqnarray} \tag{ 23 }$$

The ratio between the gamma-ray luminosity estimated via the observed scaling relation with total IR luminosity (Equation (23)) and predicted in the calorimetric limit for CR nuclei (Equation (21)) provides an estimate for the calorimetric efficiency of CR nuclei:

$$\begin{eqnarray} && F_{\rm cal} \lesssim \frac{L_{0.1\hbox{--}100 \,{\rm GeV}}\vert _{\rm scaling\;relation}}{L_{0.1\hbox{--}100 \,{\rm GeV},\pi ^{0}}\vert _{\tau _{\rm res}\approx \tau _{\rm pp}}} \approx (0.3 \pm 0.1) \nonumber\\ && \quad\times \left(\frac{\rm SFR}{M_{\odot }\,{\rm yr}^{-1}}\right)^{0.16 \pm 0.07} \left(\frac{E_{\rm SN}}{10^{51}\,{\rm erg}}\right)^{-1} \left(\frac{\eta }{0.1}\right)^{-1}. \end{eqnarray} \tag{ 24 }$$

This calorimetric efficiency estimate should be viewed as an upper limit because leptonic processes and individual sources must also contribute to the observed global gamma-ray emission of galaxies. We may infer from Equation (24) that starburst galaxies with SFR ∼10 M_☉ yr⁻¹ have maximum calorimetric efficiencies of 30%–50% whereas dwarf galaxies with SFR ∼0.1 M_☉ yr⁻¹ have lower calorimetric efficiencies of 10%–20%, when assuming an average CR luminosity per SN of E_SN η = 10⁵⁰ erg, using a Chabrier IMF, and attributing all observed gamma rays to hadronic processes. Real galaxies are expected to exhibit some dispersion about the trend. The estimates presented here are in good agreement with the calculations of Lacki et al. (2011) for M82 and NGC 253. Our own MW is the only galaxy for which global gamma-ray luminosity can be compared to direct CR measurements, and the CR energetics are in fair agreement with our estimates above considering that the estimated gamma-ray luminosity is below (but consistent with) the best-fit scaling relation.

Energy-independent cooling mechanisms for CR nuclei capable of preserving the injected spectral shape are preferred for the LAT-detected starburst galaxies given their relatively hard gamma-ray spectra. An enhancement of interaction efficiency is suggested by the slightly nonlinear scaling relations found in Section 4.3. However, the overall normalizations of these scaling relations imply that most starburst galaxies are not fully calorimetric for E_SN η = 10⁵⁰ erg. Our present analysis does not tightly constrain the CR calorimetry scenario for ULIRGs. Energy-independent CR transport should also be examined as a potentially important energy loss mechanism for CR nuclei in starburst galaxies (see in this context Persic & Rephaeli 2012). Our observations suggest that a substantial fraction of CR energy escapes into intergalactic space for most galaxies since CR nuclei are expected to dominate the CR energy budget, but that some starburst systems such as M82, NGC 253, NGC 1068, and NGC 4945 have substantially higher calorimetric efficiencies relative to the MW. A detailed discussion of particular starburst galaxies as potential CR calorimeters is presented by Lacki et al. (2011).

Celestial gamma-ray radiation not resolved into individual sources is often treated as the sum two components: a spatially structured component originating mainly from CR interactions in the MW, and a nearly uniform all-sky component often called the IGRB. Shortly after the discovery of an isotropic component using the OSO-3 satellite (Clark et al. 1968) and early spectral characterization with SAS-2 (Fichtel et al. 1975), several authors proposed that multitudinous extragalactic sources too faint to be individually resolved must contribute to the observed flux (Strong et al. 1976; Lichti et al. 1978). The spectrum of the isotropic component has been subsequently measured with EGRET (Sreekumar et al. 1998; Strong et al. 2004), and most recently with the Fermi-LAT (Abdo et al. 2010f). Note that the intensity attributed to the isotropic diffuse component is instrument dependent in the sense that more sensitive instruments are capable of extracting fainter individual sources. In this work, the term IGRB will refer specifically to the most recent measurement reported by the Fermi-LAT Collaboration, consistent with a featureless power law of spectral index of 2.41 ± 0.05 between 0.2–100 GeV with integral intensity above 0.1 GeV of 1.03 ± 0.17 × 10⁻⁵ ph cm⁻² s⁻¹ sr⁻¹.

The origin of the IGRB flux is not yet fully understood, in part because the contribution from the most prominent extragalactic gamma-ray source class, namely, blazars, is well constrained by population synthesis techniques. More than one thousand sources at high Galactic latitudes, predominantly blazars, have now been discovered at GeV energies (Ackermann et al. 2011). The total intensity of these resolved sources averaged over the full sky is 0.44 × 10⁻⁵ ph cm⁻² s⁻¹ sr⁻¹. Based on the empirically determined flux distribution of high-latitude sources, blazars with fluxes below the LAT detection threshold are expected to contribute less than 30% of the IGRB intensity (Abdo et al. 2010c).

Star-forming galaxies far outnumber AGNs in number density, but are more challenging to detect in high-energy gamma rays due to their comparatively modest luminosities and lack of beamed emission. Even with the improved sensitivities of contemporary gamma-ray telescopes, a limited number of exclusively low-redshift star-forming galaxies are expected to be individually detected. Estimates for the collective intensity of unresolved galaxies consequently rely upon scaling relations between gamma-ray luminosity and galaxy properties (such as SFR and gas content), and studies of the evolving cosmological population of galaxies. Calculations by Pavlidou & Fields (2002) and Thompson et al. (2007) during the EGRET era, and by Fields et al. (2010), Makiya et al. (2011), and Stecker & Venters (2011) incorporating additional constraints from the first year of the Fermi mission demonstrated that star-forming galaxies could make a substantial contribution to the IGRB, comparable to that of the blazars.

Fields et al. (2010) pointed out that uncertainties in the estimated contribution can arise from a degeneracy between density and luminosity evolution of star-forming galaxies over cosmic time. In other words, a change in either the density of galaxies having a given SFR (density evolution), or a change in the SFR of individual galaxies (luminosity evolution) could account for variations in total SFR density with redshift. If the relationship between gamma-ray luminosity and SFR is nonlinear, such distinctions are important.

Stecker & Venters (2011) proposed an approach combining global gamma-ray luminosity scaling relations with multiwavelength luminosity functions which describe the abundances of individual galaxies within different luminosity classes over redshift. As an example of this approach, we now estimate the contribution of non-AGN-dominated star-forming galaxies to the IGRB using the relationship between gamma-ray luminosity (0.1–100 GeV) and total IR (8–1000 μm) luminosity identified in Section 4.3, together with an IR luminosity function provided by Spitzer observations of the VIMOS VLT Deep Survey and GOODS fields (Rodighiero et al. 2010).

This procedure requires an IR luminosity function for exclusively non-AGN-dominated star-forming galaxies. In obtaining luminosity functions from Spitzer data, Rodighiero et al. utilized a catalog of spectral templates to simultaneously determine the redshift and spectral classification of each galaxy in their data set. Objects classified as type-I AGNs were subsequently removed from the sample. Although some AGNs likely remain, the fractional contamination of IR-bight obscured AGNs is expected to be small compared to non-AGN galaxies (Ballantyne & Papovich 2007; Petric 2010), such that the final IR luminosity functions accurately represent the cosmic star formation history. The IR luminosity function for non-AGN galaxies published by Rodighiero et al. (given in the rest frame of the galaxies) is consistent with those found from other IR surveys (Sanders et al. 2003; Chapman et al. 2005; Le Floc'h et al. 2005; Huynh et al. 2007; Magnelli et al. 2009).

The collective intensity (ph cm⁻² s⁻¹ sr⁻¹ GeV⁻¹) of unresolved star-forming galaxies at observed photon energy E₀ can be computed using the line-of-sight integral,

$$\begin{eqnarray} I(E_0)&=&\int _0^{z_{\rm max}} \int _{L_{{\gamma },{\rm min}}}^{L_{{\gamma },{\rm max}}} \Phi (L_{\gamma },z) \frac{d^2V}{dz d\Omega } \frac{dN}{dE} \nonumber\\ \ms \ms &&\times(L_{\gamma },E_0(1+z)) dL_{\gamma }dz, \end{eqnarray} \tag{ 25 }$$

where Φ(L_γ, z) = d²N/(dVdL_γ) expresses the number density of galaxies per unit luminosity interval and dN/dE(L_γ, E₀(1 + z)) is the differential photon flux of an individual galaxy with integral gamma-ray luminosity L_γ at redshift z. The factor d²V/(dzdΩ) represents the comoving volume element per unit redshift and unit solid angle. We transform the integral over gamma-ray luminosity into an integral over IR luminosity using the scaling relation Equation (4) and fitted parameters from Table 3.

Given the limited number of star-forming galaxies detected in high-energy gamma rays, the spectral transition from quiescent to starburst galaxies has not yet been mapped in detail. Accordingly, two spectral model choices for star-forming galaxies in the gamma-ray energy band are used in our calculation. First, we adopt a power-law spectral model with photon index 2.2, characteristic of the LAT-detected starbursts. The second spectral model is based on a model of the global emission from the MW (Strong et al. 2010) scaled to the appropriate corresponding IR luminosity. These two spectral models should be viewed as bracketing the expected contribution since multiple galaxy types with different gamma-ray spectral characteristics contribute to the IGRB, e.g., dwarfs, quiescent spirals, and starbursts.

Only the redshift range 0 < z < 2.5 is considered in this work because IR luminosity functions are not yet well constrained beyond z ∼ 2.5. We explicitly assume that the relationship between gamma-ray and IR luminosity found for galaxies at z ⩽ 0.05 holds for galaxies at higher redshifts. Although we are unable to validate this assertion directly, the closely related RC–IR correlation does not show signs of evolution up to z ∼ 2 (Ivison et al. 2010), and possibly to redshifts z ≳ 4 (Sargent et al. 2010). Those observations are indicative of a consistent relationship between star formation and CRs in galaxies over the past 10 Gyr.

We account for the attenuation of gamma rays by interactions with the extragalactic background light using the model of Franceschini et al. (2008). The contribution of resulting cascade emission to the IGRB is expected to be faint compared to the primary component (Makiya et al. 2011).

Estimates for the collective intensities of unresolved star-forming galaxies, including both quiescent galaxies and starbursts, in the energy range above 0.1 GeV are shown in Figure 7 relative to the first-year Fermi-LAT IGRB measurement. Shaded regions denote the range of statistical and systematic uncertainties. We use the scaling relation between gamma-ray luminosity and total IR luminosity obtained using the EM algorithm with the full sample of 69 galaxies, allowing for one standard deviation of variation in the fitted parameters, and include log-normal intrinsic scatter of gamma-ray luminosities. An additional 30% uncertainty is allocated for the normalization of the IR luminosity function. Using the methodology described above with the scaled MW spectral model, the estimated integral photon intensity of unresolved star-forming galaxies with redshifts 0 < z < 2.5 above 0.1 GeV is 1.2^+1.2_{− 0.6} × 10⁻⁶ ph cm⁻² s⁻¹ sr⁻¹. If the power-law model is assumed, the estimated integral photon intensity is 0.8^+0.7_{− 0.4} × 10⁻⁶ ph cm⁻² s⁻¹ sr⁻¹. These estimates cover a range of 4%–23% of the integral IGRB intensity above 0.1 GeV measured with the LAT, 1.03 ± 0.17 × 10⁻⁵ ph cm⁻² s⁻¹ sr⁻¹ (Abdo et al. 2010f). Differential intensities for the two spectral models are presented in Table 6. The range of intensities predicted in this work is consistent with previous estimates for the intensity of unresolved star-forming galaxies shown in Figure 7 for comparison.

Zoom In Zoom Out Reset image size

The relative contributions of star-forming galaxies to the IGRB according to their redshift and IR luminosity class are shown in Figure 8. In both panels, fractional contributions are normalized to the estimated total intensity of star-forming galaxies with redshifts 0 < z < 2.5. These panels provide important checks to the reliability of our total estimate given the uncertainty in the relationship between gamma-ray and IR luminosity for all galaxy types (e.g., no ULIRGs are detected at GeV energies), and our incomplete knowledge of the IR luminosity function at intermediate redshifts. The black dashed curve indicates the total IR luminosity above which Rodighiero et al. (2010) report completeness in the sample of galaxies used to derive their published IR luminosity functions. The luminosity function in the region of phase space with IR luminosity above that threshold is well constrained by Spitzer data. Although the IR luminosity function for galaxies below this threshold is not directly determined by the detections of individual galaxies, the total density of IR emission is constrained by the observed extragalactic background light. The IR luminosity functions presented by Rodighiero et al. (2010) integrated over IR luminosity and redshift are consistent with the cosmic SFR history estimated from analyses of IR observations by Pérez-González et al. (2005), Le Floc'h et al. (2005), and Caputi et al. (2007), and that estimated with far-ultraviolet data (Tresse et al. 2007).

Zoom In Zoom Out Reset image size

Although the contribution from star-forming galaxies is only considered in the redshift range from 0 < z < 2.5 in this work, Figure 8 makes apparent that the contribution from star-forming galaxies at redshifts z > 1.5 is diminishing. The particular shape of the cumulative contribution curve in the lower panel of Figure 8 is determined mainly by the IR luminosity function.

Estimates for the combined intensities of unresolved blazars and star-forming galaxies are plotted in Figure 9. The summed contribution falls short of explaining the IGRB, suggesting that important aspects of these populations are not included in current models and/or that other high-energy source classes or diffuse processes contribute a significant fraction of the observed IGRB intensity. Dermer (2007) provides an extensive discussion of possible contributions to the IGRB.

Zoom In Zoom Out Reset image size

Besides the beamed emission from blazars, in which the relativistic jet direction coincides with our line of sight, the gamma-ray emission from the cores of misaligned AGNs must also constitute part of the IGRB. More than 10 such radio galaxies have now been detected in Fermi-LAT data at redshifts up to z ∼ 0.7 (Abdo et al. 2010i). A relationship between the gamma ray and radio emission of LAT-detected misaligned AGNs has been used to estimate that gamma-ray-loud radio galaxies account for ∼25% of the unresolved IGRB above 0.1 GeV (Inoue 2011). However, estimates for this contribution are not yet tightly constrained due to the small sample of LAT-detected objects and uncertainties regarding the scaling relation between radio and gamma-ray luminosities.

Large-scale structure formation shocks leading to the assembly of galaxy clusters represent another interesting candidate population to explain the remaining IGRB intensity (Colafrancesco & Blasi 1998). The tenuous intergalactic medium of galaxy clusters is thought to be a reservoir for CR nuclei which accumulate over cosmological timescales. Although non-thermal leptonic populations are well established by observation of Mpc-scale radio features, no galaxy cluster has been detected in the GeV energy range (e.g., Ackermann et al. 2010). Other extragalactic source classes have been considered, such as gamma-ray bursts (Le & Dermer 2007), and extended emission from the lobes of radio galaxies (Stawarz et al. 2006; Massaro & Ajello 2011), but these contributions are expected to be less than 1% in the GeV energy range. Unresolved radio-quiet AGNs are also not expected to make a significant contribution to the IGRB (Teng et al. 2011).

An example of a truly diffuse component is the cumulative emission resulting from electromagnetic cascades of ultra-high energy CRs interacting with cosmic microwave background photons (Berezinskii & Smirnov 1975) and very high energy photons interacting with the extragalactic background light (Coppi & Aharonian 1997). The gamma rays produced by those cascades are expected to have a relatively hard spectrum. Already, the LAT IGRB measurement has been used to constrain this component and to predict the cosmogenic ultra-high energy neutrino flux originating from charged pion decays of the ultra-high energy CR interactions (Ahlers et al. 2010; Berezinsky et al. 2011; Wang et al. 2011).

Galactic sources, such as a population of unresolved millisecond pulsars at high Galactic latitudes, could become confused with isotropic diffuse emission as argued by Faucher-Giguère & Loeb (2010). Part of the IGRB may also come from our Solar System as a result of CR interactions with debris of the Oort Cloud (Moskalenko & Porter 2009).

Finally, a portion of the IGRB may originate from "new physics" processes involving, for instance, the annihilation or decay of dark matter particles (Bergström et al. 2001; Ullio et al. 2002; Taylor & Silk 2003).

Studies of anisotropies in the IGRB intensity on small angular scales provide another approach to identify IGRB constituent source populations (Siegal-Gaskins 2008). The fluctuation angular power contributed by unresolved star-forming galaxies is expected to be small compared to other source classes because star-forming galaxies have the highest spatial density among confirmed extragalactic gamma-ray emitters, but are individually faint (Ando & Pavlidou 2009). Unresolved star-forming galaxies could in principle explain the entire IGRB intensity without exceeding the measured anisotropy (Ackermann et al. 2012a). By contrast, the fractional contributions of unresolved blazars and millisecond pulsars to the IGRB intensity are constrained to be less than ∼20% and ∼2%, respectively, due to larger angular power expected for those source classes.