High-energy gamma-ray and neutrino production in star-forming galaxies across cosmic time: Difficulties in explaining the IceCube data

Abstract

We present new theoretical modeling to predict the luminosity and spectrum of gamma-ray and neutrino emission of a star-forming galaxy, from the star formation rate (ψ), gas mass (M_gas), stellar mass, and disk size, taking into account production, propagation, and interactions of cosmic rays. The model reproduces the observed gamma-ray luminosities of nearby galaxies detected by Fermi better than the simple power-law models as a function of ψ or ψM_gas. This model is then used to predict the cosmic background flux of gamma-rays and neutrinos from star-forming galaxies, by using a semi-analytical model of cosmological galaxy formation that reproduces many observed quantities of local and high-redshift galaxies. Calibration of the model using gamma-ray luminosities of nearby galaxies allows us to make a more reliable prediction than previous studies. In our baseline model, star-forming galaxies produce about 20% of the isotropic gamma-ray background unresolved by Fermi, and only 0.5% of IceCube neutrinos. Even with an extreme model assuming a hard injection cosmic-ray spectral index of 2.0 for all galaxies, at most 22% of IceCube neutrinos can be accounted for. These results indicate that it is difficult to explain most of the IceCube neutrinos by star-forming galaxies, without violating the gamma-ray constraints from nearby galaxies.

1 Introduction

The IceCube Collaboration has revealed the existence of extraterrestrial high-energy neutrinos, though their origin still remains a mystery (Aartsen et al. 2013). The arrival distribution is consistent with being isotropic and the flavor ratio is consistent with ν_e: ν_μ: ν_τ = 1 : 1 : 1, suggesting an extragalactic, astrophysical origin. A variety of source models have been proposed so far, such as gamma-ray bursts, active galactic nuclei, star-forming galaxies, and galaxy clusters, though no individual source has been identified yet (for recent reviews, see, e.g., Mészáros 2017; Halzen 2017).

Star-forming galaxies are one possible origin of the IceCube neutrinos, in which cosmic rays (CRs) are produced by supernovae and they produce pion-decay neutrinos via inelastic collisions with the interstellar medium (ISM) (Loeb & Waxman 2006; Thompson et al. 2006; Stecker 2007; Lacki et al. 2011; Murase et al. 2013; He et al. 2013; Tamborra et al. 2014; Anchordoqui et al. 2014; Liu et al. 2014; Emig et al. 2015; Chang et al. 2015; Giacinti et al. 2015; Senno et al. 2015; Moharana & Razzaque 2016; Chakraborty & Izaguirre 2016; Xiao et al. 2016; Bechtol et al. 2017). Starburst galaxies are especially important in this context, because CRs are expected to produce pions efficiently, due to their high gas densities and confinement by strong magnetic fields. However, a definite conclusion has not yet been reached about whether star-forming galaxies are the dominant source of IceCube neutrinos. Tamborra, Ando, and Murase (2014) concluded that starburst sources could possibly explain a significant fraction of the IceCube flux based on an empirical relation between gamma-ray and infrared luminosities, while Bechtol et al. (2017) argued that galaxies cannot produce more than 30% of the IceCube data if the upper limit on the non-blazar fraction of the extragalactic gamma-ray background is considered. Xiao et al. (2016) found that star-forming galaxies might explain about 50% of the IceCube flux, but assumed a dominant contribution from hypernovae, which is rather uncertain. In these studies, however, little attention was paid to the consistency of the model with observed gamma-ray fluxes from nearby star-forming galaxies (Ackermann et al. 2012a; Tang et al. 2014; Peng at al. 2016), which should provide a useful constraint because gamma-rays are inevitably emitted if neutrinos are produced in a star-forming galaxy.

In this work, we present a new calculation of the contribution from star-forming galaxies to the diffuse gamma-ray and neutrino background, which takes into account cosmic-ray physical processes in galaxies and is consistent with the gamma-ray observations of nearby galaxies. For this purpose, we construct a new theoretical framework to predict both flux and spectrum of gamma-ray and neutrino emissions from a galaxy based on star formation rate (SFR), stellar mass, gas mass, and effective radius.

Gamma-ray flux from a galaxy has been modeled in a number of studies to predict the cosmic gamma-ray background from star-forming galaxies (Strong et al. 1976; Lichti et al. 1978; Dar & Shaviv 1995; Pavlidou & Fields 2002; Thompson et al. 2007; Ando & Pavlidou 2009; Fields et al. 2010; Makiya et al. 2011; Stecker & Venters 2011; Ackermann et al. 2012a; Chakraborty & Fields 2013; Lacki et al. 2014; Komis et al. 2017; Lamastra et al. 2017), and most studies estimated gamma-ray luminosity only from one or two physical quantities (e.g., SFR, gas mass, or infrared luminosity) assuming empirical relations. However, such an approach may induce some bias in the predictions of the cosmic background (Komis et al. 2017). The gamma-ray spectrum was also often treated in a simple way, using the Milky Way and M 82 spectra as templates for normal and starburst galaxies. Our model includes a larger number of physical quantities of a galaxy to take into account the production, propagation, and interactions of cosmic rays, and hence flux and spectrum can be calculated for diverse individual galaxies across cosmic time.

This modeling is compared with the observed gamma-ray luminosities of six nearby galaxies, and we will show that our modeling reproduces the observed gamma-ray luminosities better than the simple scaling relations with SFR and gas mass. The cosmic background flux and spectrum of gamma-rays and neutrinos are then calculated by using a semi-analytical model of cosmological galaxy formation (the Mitaka model, Nagashima & Yoshii 2004), which reproduces many observational data of local and high-z galaxies. The Mitaka model provides us with the physical quantities of galaxies, enabling us to predict gamma-ray and neutrino emissions based on a self-consistent calculation of formation and evolution of galaxies in a cosmological framework, taking into account key baryonic processes such as gas cooling, star-formation, galaxy merger, and subsequent starbursts.

We organize this paper as follows. The new model of gamma-ray and neutrino emission from a star-forming galaxy is described in section 2. The model prediction is compared with the gamma-ray observations of nearby galaxies in section 3. Before discussing the cosmic neutrino background, we examine the expected contribution to IceCube neutrinos from the Galactic disk in section 4. The results on the cosmic gamma-ray and neutrino background are then presented in section 5. Discussions on model dependence and uncertainties are given in section 6, followed by summary in section 7. We adopt a flat ΛCDM cosmology with Ω_M = 0.3, Ω_B = 0.05, and h = 0.7 throughout this work.

2 Methods

2.1 Production, propagation and interaction of cosmic rays

Suppose that the four quantities of SFR (denoted as ψ), stellar mass M_*, gas mass M_gas, and the effective disk radius R_eff (the radius including half of the total galactic light, related to the exponential scale R_d as R_eff = 1.68 R_d) are given for a galaxy. We first need to determine the production rate of pions by CR interactions as a function of the CR proton energy E_p.

We assume that the production rate of CRs is proportional to the SFR, and the CR spectrum at the time of injection into the ISM is described by a single power-law. Then the number of CRs produced in a galaxy per unit time and CR energy is expressed as

\begin{equation} \frac{dN_{\rm p}}{dtdE_{\rm p}} = C \left(\frac{\psi }{{M}_\odot \:{\rm yr}^{-1}} \right) \left(\frac{E_{\rm p}}{{\rm GeV}}\right)^{-\Gamma _{\rm inj}}, \end{equation}

(1)

where the normalization factor C will be fixed by fitting to the observed gamma-ray luminosities of nearby galaxies later. Observations of the cosmic-ray spectrum on Earth, gamma-ray spectrum of supernova remnants, and the Galactic diffuse gamma-ray emission favor the injection spectral index in a range of Γ_inj = 2.2–2.4 (see Ackermann et al. 2012b; Caprioli 2012 and references therein). In this work we adopt Γ_inj = 2.3 as the baseline value, but we will also test the case of the strong shock limit, Γ_inj = 2.

A fraction f_π(E_p) of CRs interact with the ISM before they escape into the intergalactic space, and it can be expressed as f_π(E_p) = 1 − exp (−t_esc/t_pp), where t_esc(E_p) is the escape time of a CR particle from the galaxy and t_pp(E_p) is the mean time taken to interact with the ISM by the proton–proton (pp) collisions. This can be written as t_pp = (n_gasσ_ppc)⁻¹, where n_gas is the number density of ISM and σ_pp(E_p) is the inelastic part of the pp cross-section, for which we use the formula given in Kelner, Aharonian, and Bugayov (2006). We calculate n_gas as |$n_{\rm gas} = M_{\rm gas}/(2 \pi R_{\rm eff}^2H_{\rm g} m_{\rm p})$|⁠, where m_p is the proton mass and H_g is the scale height of the gas disk.

Our galaxy formation model does not compute the scale height H_g. The mechanism to determine disk heights of galaxies is complex, depending on many physical processes including magnetic field and cosmic rays, gas pressure, and turbulent motion. In this work, we take a simple empirical approach to estimate H_g by assuming H_g ∝ R_eff, and the coefficient is determined by the observations for our Galaxy: H_g = 150 pc (Mo et al. 2010) and R_eff = 6.0 kpc (Sofue et al. 2009). Observations of nearby galaxies show that the stellar scale height of disks (H_*) is roughly proportional to R_eff, and the above H_g/R_eff ratio is consistent with observations if the difference between stellar and gas scale heights is taken into account (H_* ∼ 2H_g in our Galaxy). In section 6 we will also examine other modelings of H_g and the effects on our results. In the Mitaka galaxy formation model, starbursts occur at the time of a major merger of two disk galaxies, resulting in the formation of a spheroidal galaxy. We assume that such star-bursting galaxies are nearly spherical, and hence H_g = R_eff.

The escape time, t_esc, is determined as follows: |$t_{\rm esc}$| = min [|$t_{\rm diff}$|⁠, |$t_{\rm adv}$|], where |$t_{\rm diff}$| and |$t_{\rm adv}$| are the CR diffusion time and the advection time by galactic outflow, respectively. These are estimated from galactic properties as |$t_{\rm diff}=H_{\rm g}^2/[2D(E_{\rm p})]$| and t_adv = H_g/σ, where D(E_p) is the diffusion coefficient of CRs and σ is the escape velocity from the gravitational potential of the galactic disk. We estimate σ from H_g and the column density of total mass |$\Sigma = (M_* + M_{\rm gas})/(\pi R_{\rm eff}^2)$| assuming the relation for an isothermal sheet: GΣ = σ²/(2πH_g) (Mo et al. 2010).

The diffusion coefficient D is not well constrained by observations. Theoretically, D is expected to depend on the CR energy E_p and magnetic field strength B. Furthermore, fluctuation pattern of magnetic fields is important for CR diffusion. We consider two regimes regarding the proton Larmor radius R_L = 2.0 × 10⁻⁷(E_p/GeV)(B/6 μG)⁻¹ pc following Aloisio and Berezinsky (2004). Suppose that there is a coherent length of turbulence l₀, and magnetic fields are random and uncorrelated beyond l₀. This length should be smaller than the region size, i.e., l₀ ≤ H_g. Then the diffusion of high energy CRs with l₀ < R_L < (H_gl₀)^1/2 is described by the small-angle random scattering approximation with a mean free path |$l_{\rm mfp} \sim (R_{\rm L}^2/l_0)$|⁠, and hence |$D \sim c l_{\rm mfp}/3 \sim c R_{\rm L}^2 / (3l_0)$|⁠. It should be noted that when R_L > (H_gl₀)^1/2, the mean free path becomes larger than H_g, and hence in this regime we set l_mfp ∼ H_g, and hence D ∼ cH_g/3 and t_diff = 3H_g/(2c). Diffusion of lower energy CRs with R_L < l₀ is determined by resonant scattering in the turbulent magnetic field fluctuations, resulting in D ∼ cl₀(R_L/l₀)^δ/3, where δ depends on the spectrum of interstellar turbulence and we adopt δ = 1/3 for the Kolmogorov-type turbulence. Then the diffusion coefficient can be described as follows:

\begin{eqnarray} D(E_{\rm p}) = \left\lbrace \begin{array}{ll}\frac{c l_0}{3} \left[ \left( \frac{R_{\rm L}}{l_0}\right)^\frac{1}{3} + \left(\frac{R_{\rm L}}{l_0} \right)^2 \right] & \left( R_{\rm L} \le \sqrt{H_{\rm g} l_0} \right) \\ \frac{cH_{\rm g}}{3} & \left( R_{\rm L} > \sqrt{H_{\rm g} l_0} \right) \end{array} \right. \end{eqnarray}

(2)

(see also Parizot 2004).

The coherent length of turbulence l₀ is poorly understood theoretically, but l_outer ∼  100–200 pc has been suggested for turbulence generated by supernova remnants in the Milky Way (Haverkorn et al. 2008; Chepurnov & Lazarian 2010; Iacobelli et al. 2013), where l_outer is the outer scale of turbulence and related as l₀ = l_outer/5 for the Kolmogorov turbulence. If this is a result of the physics of supernova remnants interacting with the ISM, it may not strongly depend on the global properties of galaxies. Therefore, as a baseline model we use a common value of l_{0, max} = 30 pc for all galaxies whose H_g is larger than l_{0, max} but l₀ = H_g for smaller galaxies, i.e., l₀ = min (l_{0, max}, H_g). Uncertainties and other possible modelings about l₀ will be discussed in section 6. The Larmor radius becomes larger than l₀ for a very high energy of E_p > 1.5 × 10¹⁷(l₀/30 pc)(B/6 μG) eV, and hence the regime of R_L < l₀ is relevant for GeV gamma-rays and PeV neutrinos in most galaxies.

In order to estimate magnetic field strength, we assume that the energy density of the magnetic field is close to that of supernova explosions injected into ISM on the advection time scale, t_adv = H_g/σ. Then, using a dimensionless parameter η, the magnetic field is given as B²/(8π) = ηE_SNr_SNt_adv/V, where E_SN is the kinetic energy of a supernova, r_SN the supernova rate in a galaxy, and |$V = 2\pi R_{\rm eff}^2 H_{\rm g}$| the volume of the gas disk. We set η = 0.03 because it reproduces the observed magnetic field strength of our Galaxy (6 μG, Beck 2008; Haverkorn 2015), using ψ, R_eff, M_gas, and M_* in table 1. Here we made standard assumptions that stars heavier than 8 M_⊙ end their life by supernovae, E_SN = 10⁵¹ erg, and the Salpeter IMF (Woosley & Weaver 1995; Salpeter 1955) for all galaxies.

Table 1.

Properties of gamma-ray galaxies.

Objects	L  γ	Ref.‡	ψ	Ref.‡	M  gas	Ref.‡	M*	Ref.‡	R  eff	Ref.‡
	[1039 erg s−1]†		[M⊙ yr−1]		[109 M⊙]		[109 M⊙]		[kpc]
MW	0.82 ± 0.24	(1)	2.6	(3), (4)	4.9	(8)	50	(14)	6.0	(19)
LMC	0.047 ±  0.005	(1)	0.24	(4), (5)	0.53	(9)	1.5	(15)	2.2	(20)
SMC	0.011 ±  0.003	(1)	0.037	(4), (5)	0.45	(10)	0.46	(15)	0.7	(21)
NGC 253	6 ±  2	(1)	7.9	(4), (6)	4.3	(11)	21	(16)	3.7	(22)
M 82	15 ±  3	(1)	16.3	(4), (6)	1.3	(12)	8.7§	(17)	1.2	(22)
NGC 2146	40 ±  21	(2)	17.5‖	(7)	4.1	(13)	20	(18)	1.8	(23)

^†The photon energy range is 0.1–100 GeV except for NGC 2146, for which 0.2–100 GeV is adopted according to ref. (2).

^‡References: (1) Ackermann et al. (2012a), (2) Tang, Wang, and Tam (2014), (3) Diehl et al. (2006), (4) Makiya, Totani, and Kobayashi (2011), (5) Kennicutt et al. (2008), (6) Sanders et al. (2003), (7) Yun, Reddy, and Condon (2001), (8) Paladini et al. (2007), (9) Staveley-Smith et al. (2003) and Fukui et al. (2008), (10) Stanimirović et al. (1999) and Leroy et al. (2007), (11) Knudsen et al. (2007) and Springob et al. (2005), (12) Chynoweth et al. (2008) and Mao et al. (2000), (13) Tsai et al. (2009), (14) Bland-Hawthorn and Gerhard (2016), (15) McConnachie et al. (2012), (16) Lucero et al. (2015), (17) Sofue et al. (1992), (18) Skibba et al. (2011), (19) Sofue, Honma, and Omodaka (2009), (20) van der Marel (2006), (21) Gonidakis et al. (2009), (22) J-band half-light radius in Jarrett et al. (2003), converted into kiloparsec using distances taken from the NASA/IPAC Extragalactic Database, (23) Hα half-light radius in Marcum et al. (2001).

^§From a dynamical mass estimate (see text).

^‖Estimated from the radio continuum luminosity at 1.4 GHz and equation (13) of Ref. (7).

Now we can calculate the CR production rate dN_p/dtdE_p and the fraction f_π of these interacting with the ISM as a function of E_p, if M_*, M_gas, ψ, and R_eff of a galaxy are given.

2.2 Neutrinos and gamma-rays from galaxies

Following from the previous subsection, the number luminosity of neutrinos or gamma-rays, i.e., the number of particles produced per unit time and particle energy in a galaxy, is calculated as

\begin{equation} \frac{dL_{N,i}}{dE_i} = \int _{E_i}^\infty f_\pi \frac{dN_{\rm p}}{dtdE_{\rm p}} \frac{dn_i}{dE_i} dE_{\rm p}, \end{equation}

(3)

where the index i denotes neutrino or gamma-ray (i = ν, γ), dn_i/dE_i is the number of neutrinos or gamma-ray photons per unit particle energy E_i generated from one pp interaction, which is a function of E_i and E_p and is calculated using the formulae in Kelner, Aharonian, and Bugayov (2006). It should be noted that in equation (3) we only consider the initial pp collision of each cosmic ray with the ISM, because a CR particle loses a significant fraction of its energy in the first collision and subsequent interactions have minor energy contributions to the gamma-ray or neutrino spectrum. Furthermore, the contribution from the second and later interactions is effectively taken into account because we fix the normalization factor C to fit the observed gamma-ray luminosity of nearby galaxies.

The cosmic background radiation flux of gamma-rays or neutrinos per steradian and per unit E_i is then calculated as

\begin{eqnarray} \Phi _i(E_i) = \frac{c}{4\pi } \int dz \left| \frac{dt}{dz} \right| (1+z) \frac{d{\cal L}_{N,i} [(1+z)E_i; z]}{dE_i^{\rm rest}}, \end{eqnarray}

(4)

where |$E^{\rm rest}_i = (1+z) E_i$| is the rest-frame particle energy, and |${d{\cal L}_{N,i}[E_i^{\rm rest};z]}/{dE_i^{\rm rest}}$| is the neutrino/gamma-ray number luminosity density at redshift z per unit of comoving volume. Our galaxy formation model generates mock galaxy catalogues at various redshifts, in which the jth galaxy is given a weight |$n_j^{\rm gal}$| (comoving number density of such a galaxy). We use morphology, ψ, M_gas, M_*, and R_eff of these catalog galaxies to calculate the number luminosity density by summing up the catalog galaxies as:

\begin{equation} \frac{d{\cal L}_{N,i}[E_i;z]}{dE_i} = \sum _j n^{\rm gal}_j \left[\frac{dL_{N,i}(E_i)}{dE_i}\right]_j \ . \end{equation}

(5)

For the gamma-ray background, additional corrections to equation (4) are necessary to take into account the attenuation of the gamma-ray flux due to e^± creation by interaction with background optical and infrared photons, and the subsequent cascade emission produced by e^±s. We use the baseline model calculation of τ_γγ[E_γ; z] in Inoue et al. (2013) for the absorption optical depth, and calculate the cascade emissivity following the treatment of Inoue and Ioka (2012).

2.3 Semi-analytical modelling of galaxy formation

In this work we use a semi-analytical model of cosmological galaxy formation presented by Nagashima and Yoshii (2004), which is called the Mitaka model. This model first computes the merger history of dark matter haloes based on the extended Press–Schechter theory. Then baryonic processes in these haloes, such as gas cooling, star formation, supernova feedback, and galaxy mergers, are calculated with phenomenological prescriptions. Free parameters included in describing baryons are determined to match several observations of local galaxies. This model reproduces various observed properties of local and high-z galaxies, such as luminosity functions, size-luminosity relation, luminosity densities, and stellar mass densities (Nagashima & Yoshii 2004; Kashikawa et al. 2006; Kobayashi et al. 2007, 2010).

This model includes two modes of star formation: starburst and quiescent. Starburst is assumed to occur after major mergers of galaxies leading to spheroidal galaxies; otherwise, stars are formed in a disk (the quiescent mode). Contributions of starbursts to the total cosmic SFR in this model is about 5% at z ∼  0, which increases to 15% at z ∼  1 and 30% at z ∼  2.

3 Comparison with nearby galaxies

To fix the value of the normalization factor C in the equation (1), we use the data of six galaxies detected by the Fermi-LAT, including three normal galaxies (SMC, LMC, MW) and three starbursts (NGC 253, M 82, NGC 2146). It should be noted that the latter three galaxies are generally called “starbursts” by their intensive star formation activity, but their morphology is disk-like and without evidence of major mergers in recent past. Therefore we treat these as disk galaxies in the theoretical modeling, and hence the definition of ”starburst” is different from that (major merger) used in the galaxy formation model. Their physical quantities are shown in table 1.

There are several notes for this sample. A different photon energy range of gamma-ray luminosity L_γ is used for NGC 2146 in the literature, and hence we also change the range inthe theoretical calculation accordingly. Stellar mass M_* for M 82 is not yet observationally well constrained, and we used a dynamical mass M_dyn ∼  10¹⁰ M_⊙ (Sofue et al. 1992) and estimated M_* from M_* = M_dyn − M_gas. However, even if we set M_* = 0 for M 82, our results are hardly changed. In addition to the galaxies shown in table 1, M 31, NGC 4945, NGC 1068, and Arp 220 are also detected in gamma-rays. However, we do not use NGC 4945, NGC 1068, or Arp 220 in our calculation, because gamma-ray emissions from these galaxies are likely affected by active galactic nuclei activities (Ackermann et al. 2012a; Yoast-Hull et al. 2017). M 31 is also removed from the sample, because a recent analysis of gamma-rays from M 31 has shown that the emission is not correlated with regions where most of the atomic and molecular gases reside, suggesting that they did not originate from the CR–ISM interactions (Ackermann et al. 2017).

The parameter C is then determined by a standard χ²-fitting to the six observed values of L_γ. Figure 1 shows a comparison of gamma-ray luminosities between model predictions and observations. We also show the best-fitting proportional and power-law relationships to the observed data as a function of ψ or ψM_gas, because fits to these quantities were often made in previous studies, motivated by the expectation that CR production rate is proportional to ψ and the target gas mass for pp reactions scales with M_gas. Though our model is based on a simple modeling of physical processes on the whole galactic scale, our simple model nicely reproduces the observed luminosities compared with the power-law fits to ψ or ψM_gas. It may be rather surprising that this simple model predicts the correct luminosities for various types of galaxies, ranging widely from dwarf galaxies like LMC and SMC to intense starburst galaxies. Note that C can be converted to the energy injected into CRs from a supernova, as can be seen from equation (1). Setting the energy range of CRs to be 10⁹–10¹⁵ eV, the best-fitting C corresponds to the CR energy of 0.24E_SN with the same E_SN, IMF, and the threshold stellar mass for supernovae assumed in the model. This is comparable to the standard assumption that CRs carry ∼ 10% of supernova explosion energy. This slightly larger value is partly a result of considering only the first time pp interaction to produce gamma-rays, as mentioned above.

Fig. 1.

Predicted (crosses) and observed (circles) gamma-ray luminosities of nearby galaxies are compared in L_γ–ψM_gas (left-hand panel) and L_γ–ψ (right-hand panel) plots. The baseline model parameters described in the text are used in the model calculation. The best-fitting proportional (dotted) and power-law (dashed) relations to the observed data are also shown. (Color online)

4 Gamma-rays and neutrinos from the galactic disc

The gamma-ray sky measured by Fermi is dominated by the diffuse Galactic emission (DGE), which is mostly from the pion decay, and hence the Milky Way could also be a source of high-energy neutrinos. It is important to check the predicted neutrino flux from the Milky Way does not contradict the IceCube data. In figure 2 we show our model prediction of the Galactic diffuse gamma-ray and neutrino emission with our baseline model parameters, but changing Γ_inj. For comparison, we also show the pion-decay component of the ^SS^Z4^R20^T150^C5 model of DGE in Ackermann et al. (2012b) (the A12 model, hereafter) within the low-latitude regions of |b| < 8°, which was calculated by the more detailed GALPROP code and is in good agreement with the observed data. The flux of our model is normalized so that it matches the A12 model at 1 GeV. We also show the IceCube neutrino flux data, which is scaled into the |b| < 8° region from the whole IceCube data, assuming isotropy.

Fig. 2.

Predicted gamma-ray (left-hand panel) and neutrino (right-hand panel) spectra from the Galactic disk region are shown. The spectral shapes of gamma-rays and neutrinos predicted by our model are shown with some different values of Γ_inj. The black dashed curve shows the spectrum of the diffuse Galactic gamma-ray emission (DGE) in the disk regions of |b| <  8°, extracted from the GALPROP-based model of Ackermann et al. (2012b). The red squares are neutrino spectrum per flavor observed in all sky by the IceCube experiment, but multiplied by the solid angle fraction of the Galactic disk (|b| <  8°) as an estimate for the observed neutrino flux in this region. (Color online)

Figure 2 indicates that our model matches the A12 model for Γ_inj = 2.3, in which case the contribution from the Milky Way can possibly explain 13% of the IceCube energy flux. It should be noted here that we do not introduce any cutoff in the accelerated cosmic ray spectrum [equation (1)], which may be rather unrealistic. Indeed, if we introduce a cutoff at PeV, for example, the contribution decreases to 3%. Interestingly, some previous studies have also suggested a considerable contribution to IceCube neutrinos from the Milky Way (Gaggero et al. 2015; Neronov & Semikoz 2016; Palladino & Vissani 2016). Increased statistics of the IceCube neutrino events in the future may reveal an excess from the Galactic disk region compared with the isotropic component, which would give important information about the maximum cut-off energy of proton acceleration in the Milky Way. On the other hand, a hard injection spectrum of Γ_inj ≲ 2.1 is disfavored because the Galactic disk component of neutrinos would be too strong, unless there is a cut-off below E_p ∼  10¹⁴ eV. It should be noted that this constraint is from the average CR spectrum in our Galaxy, and we cannot rule out the possibility that the spectrum at production is different for other galaxies, especially for starburst galaxies that on average have dense ISM and strong magnetic fields. Hence, we also test optimistic cases where Γ_inj ≲  2.1 is realized in other galaxies, even though it is unlikely for our Galaxy.

5 Cosmic gamma-ray and neutrino background

Figure 3 presents the cosmic gamma-ray and neutrino background spectra predicted by our baseline model, in comparison with the Fermi and IceCube data. Our calculations show that the gamma-ray energy flux from star-forming galaxies is (5.1–7.0) × 10⁻⁴ MeV cm⁻²s⁻¹str⁻¹ above 100 MeV, which is 18%–25% of the IGRB flux observed by Fermi, in reasonable agreement with previous studies (Fields et al. 2010; Makiya et al. 2011; Stecker & Venters 2011; Ackermann et al. 2012a; Lacki et al. 2014).

Fig. 3.

Cosmic diffuse background of gamma-rays (left-hand panel) and neutrinos (right-hand panel) from star-forming galaxies predicted by our baseline model are shown for different values of Γ_inj. Solid and dashed lines correspond to the contributions from all galaxies and starburst galaxies respectively. Data points represent the gamma-ray spectrum of unresolved isotropic gamma-ray background observed by Fermi-LAT [blue, Ackermann et al. (2015)] and the astrophysical neutrino spectrum per flavor observed in the IceCube experiment [red, Aartsen et al. (2015a)]. For the purpose of presentation, the scale of the vertical axis is different between the left- and right-hand panels. (Color online)

The neutrino flux predicted by our baseline model (Γ_inj = 2.3) is only a 0.5% contribution to the IceCube data. If we assume a harder spectrum at injection Γ_inj = 2.1 and 2.2, the contributions increase to 8.4% and 2.1%, respectively. Even in the most optimistic (and extreme) case of Γ_inj = 2 in all galaxies, star-forming galaxies can account for only 22% of the IceCube data. It should be noted that such a hard injection spectrum in our Galaxy is in contradiction with the Galactic diffuse gamma-ray spectrum and the isotropy of the IceCube data. Therefore, we conclude that star-forming galaxies cannot be the major source of the IceCube neutrinos, and a reasonable estimate of their contribution is about 1%–8% or less.

6 Dependence on model parameters

Dependence of our results on different modelings of l_{0, max} and H_g is shown as the change of the gamma-ray luminosities of nearby galaxies (figure 4) and the cosmic gamma-ray/neutrino background flux (figure 5) from our baseline model. For l_{0, max} we test different values, 10 and 100 pc, from the baseline model (30 pc), and also test a model with l_{0, max} = H_g, as the case where the coherent scale of turbulence is determined by the system size. The model with l_{0, max} = 10 pc also agrees reasonably well with the L_γ data of nearby galaxies, and the models with l_{0, max} = 100 pc and l_{0, max} = H_g show larger deviations, especially for MW and NGC 253, which deviate by up to a factor of 3. The changes in background fluxes are small (less than 20%) for l_{0, max} = 10 and 100 pc, but the neutrino flux is reduced by a factor of 2–3 in the model of l_{0, max} = H_g. This is because l₀ becomes larger in galaxies of H_g > 30 pc, and the diffusion coefficient becomes larger with larger l₀ if R_L < l₀, which is valid in most galaxies.

Fig. 4.

Comparisons between predicted and observed gamma-ray luminosities of nearby galaxies for different modellings of l_{0, max} (upper panel) and H_g (lower panel). In both panels the red cross points correspond to our baseline model. For the purpose of presentation, gamma-ray luminosities of SMC, LMC, and MW are multiplied by the numbers indicated in parentheses in the figure. (Color online)

Fig. 5.

Same as figure 3, but showing dependence on the modelling of l_{0, max} (upper panels) and H_g (lower panels), using the baseline model value of Γ_inj = 2.3.

To check the dependence on the modelling of H_g, we test the following two cases for disk galaxies: (1) assuming that the height of gas disk in all disk galaxies are similar to that of the Milky Way, i.e., H_g = 150 pc, and (2) assuming energy equipartition between gas motion and energy produced by supernovae, as ρ_gasσ² = αE_SNr_SNt_adv/V, where ρ_gas = M_gas/V is the ISM gas density. In the latter case, both σ and H_g are determined by this relation combined with GΣ = σ²/(2πH_g), without using H_g ∝ R_eff. The dimensionless factor α = 4 is determined as reproducing H_g = 150 pc for physical quantities of the Milky Way.

The lower panel of figure 4 shows that different modelling of H_g results in clear discrepancies between predicted and observed L_γ for local galaxies. As shown in the lower panels of figure 5, the gamma-ray and neutrino background flux is reduced in the constant H_g model by a factor of about 2. This is likely because the gas density is underestimated in small galaxies whose density is higher if we assume H_g ∝ R_eff. The background fluxes are reduced by a larger factor of about 3–4, for the model assuming equipartition between gas and supernovae. This implies that the gas density is higher than that expected from equipartition. Such a situation may be expected, especially in starburst galaxies, if gas density becomes high enough before star formation starts and supernova energy is injected to the ISM. Figure 4 indeed indicates that the discrepancy in the gas–SN equipartition model is larger for starburst galaxies. In any case, different modelings of l_{0, max} and H_g lead to lower neutrino background flux, and hence it does not affect our conclusions that the majority of IceCube neutrinos cannot be explained by star-forming galaxies.

7 Summary

In this work we constructed a new theoretical model to predict flux and spectrum of gamma-rays and neutrinos by interactions of cosmic rays produced by supernovae in a star-forming galaxy, and applied it in order to predict the flux and spectrum of the cosmic gamma-ray and neutrino background radiation.

Our model calculates gamma-ray and neutrino spectra of a star-forming galaxy from four physical quantities, i.e., SFR, size, gas mass, and stellar mass, taking into account the production, propagation, and interactions of cosmic rays in the ISM. This model is tested against the gamma-ray luminosities measured by Fermi of the six nearby galaxies listed in table 1. In spite of this sample including a wide variety of star-forming galaxies, from dwarf galaxies like LMC and SMC to starburst galaxies such as M 82, our model reproduces the observed gamma-ray luminosities fairly well, within a factor of ∼2. The agreement is significantly better than the simple phenomenological modelings assuming a power-law relation between L_γ and SFR or SFR × M_gas. This model can be tested with a larger sample of nearby galaxies by high-sensitivity gamma-ray observations in future projects such as the Cherenkov Telescope Array (CTA).

We have examined the neutrino emission from the disk of our Galaxy predicted by our model, and found that the observed isotropic distribution of IceCube neutrinos constrains the injection cosmic-ray energy spectrum in our Galaxy to be softer than Γ_inj ∼  2.2, unless there is a cut-off at ≲  10¹⁴ eV in the CR energy spectrum.

The contribution from star-forming galaxies to the extragalactic gamma-ray and neutrino background is calculated by using a semi-analytical cosmological galaxy formation model. This model is quantitatively in agreement with many observations of galaxies at local and high redshifts. We have found that star-forming galaxies make about 20% of the isotropic gamma-ray background flux unresolved by Fermi, and only 1%–8% or less of the IceCube flux, with reasonable values of Γ_inj = 2.1–2.3. Even with the most optimistic case, where Γ_inj = 2.0 in all galaxies, the contribution is no more than 22%. Therefore, we conclude that the majority of IceCube neutrinos cannot be explained by star-forming galaxies. We also examined the dependence of these results on modelings about l_{0, max} (the maximum length of coherence in ISM turbulence) and H_g (gas scale height of a galactic disk), and found that alternate prescriptions give an even lower neutrino flux.

It may be worth mentioning that the Pierre Auger Collaboration has recently reported a possible correlation of ultra-high-energy cosmic rays (UHECRs) above E ∼  10¹⁹ eV with nearby starburst galaxies (Aab et al. 2018). In this work we considered cosmic-ray production from normal supernova remnants, and it is generally considered that such cosmic rays cannot reach the energy scale of UHECRs. The reacceleration of such low-energy cosmic rays by e.g., large-scale terminal shocks generated by galactic wind (Anchordoqui 2018) would be required to produce UHECRs originating from normal supernovae. Another possibility is that UHECRs are produced by other sources related to young stellar populations.

Our results demonstrate that star-forming galaxies make only a minor contribution to the IceCube flux, implying that other sources are required to explain most of the observed data. Since correlation analyses have shown that gamma-ray blazars and gamma-ray bursts are not the dominant source (Aartsen et al. 2015b, 2017), there is no strong candidate for the origin of IceCube neutrinos. Future detectors such as the IceCube-Gen2, as well as multi-messenger approaches with next-generation telescopes like CTA, will be the key to understanding the nature of high-energy neutrinos.

Acknowledgments

TT was supported by JSPS KAKENHI Grant Numbers JP15K05018 and JP17H06362. NK is supported by the Hakubi project at Kyoto University. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.