Interpreting the relation between the gamma-ray and infrared luminosities of star-forming galaxies

TL;DR

This paper investigates the relationship between gamma-ray and infrared luminosities in star-forming galaxies, revealing a steepening of the correlation in low-SFR galaxies due to cosmic-ray escape effects.

Contribution

It introduces a model considering galaxy properties to explain the non-linear gamma-ray to infrared luminosity relation, especially in low-SFR galaxies.

Findings

The gamma-ray to infrared luminosity relation steepens in low-luminosity galaxies.

Cosmic-ray escape significantly affects gamma-ray production in low-SFR galaxies.

The model's predictions align with observed data.

Abstract

It has been found that there is a quasi-linear scaling relationship between the gamma-ray luminosity in GeV energies and the total infrared luminosity of star-forming galaxies, i.e. $L_{\gamma}\propto L_{\rm IR}^{\alpha}$ with $\alpha\simeq 1$. However, the origin of this linear slope is not well understood. Although extreme starburst galaxies can be regarded as calorimeters for hadronic cosmic ray interaction and thus a quasi-linear scaling may hold, it may not be the case for low star-formation-rate (SFR) galaxies, as the majority of cosmic rays in these galaxies are expected to escape. We calculate the gamma-ray production efficiency in star-forming galaxies by considering realistic galaxy properties, such as the gas density and galactic wind velocity in star-forming galaxies. We find that the slope for the relation between gamma-ray luminosity and the infrared luminosity gets steeper for low infrared luminosity galaxies, i.e. $\alpha\rightarrow 1.6$, due to increasingly lower efficiency for the production of gamma-ray emission. We further find that the measured data of the gamma-ray luminosity is compatible with such a steepening. The steepening in the slope suggests that cosmic-ray escape is very important in low-SFR galaxies.

Tables

Table 1. Table 1: Parameters of GeV-detected star-forming galaxies

Name	$L_{1-500\mathrm{GeV}}$	$L_{\mathrm{IR}}$	SFR
	$\mathrm{erg}{s}^{-1}$	${10}^9{L}_{\odot }$	$M_{\odot }{\mathrm{yr}}^{-1}$
SMC	$\left(7.57\pm 0.34\right)\times {10}^{36}$	$0.07\pm 0.01$	0.04-0.08aafootnotemark:
LMC	$\left(2.94\pm 0.05\right)\times {10}^{37}$	$0.7\pm 0.1$	0.20-0.25bbfootnotemark:
M31	$\left(6.31\pm 0.21\right)\times {10}^{37}$	$2.4\pm 0.4$	0.35-1ccfootnotemark:
NGC 253	$\left(3.51\pm 0.40\right)\times {10}^{39}$	21	3.5-10.4ddfootnotemark:
M82	$\left(7.74\pm 0.28\right)\times {10}^{39}$	46	13-33eefootnotemark:
NGC 2146	$\left(3.39\pm 1.20\right)\times {10}^{40}$	100	26.6-79.7fffootnotemark:
Arp 220	$\left(6.24\pm 2.21\right)\times {10}^{41}$	1400	254.8-764.3fffootnotemark:

Equations

$$\log \left(\frac{L_{\gamma }}{\mathrm{erg}{\mathrm{s}}^{-1}}\right)=\alpha \log \left(\frac{L_{\mathrm{IR}}}{10^{10}{L}_{\odot }}\right)+\beta ,$$

$$f_{\pi }=1-\exp (-\frac{t_{\mathrm{esc}}}{t_{\mathrm{loss}}}),$$

$$t_{\mathrm{loss}}=3.36\times 10^5\frac{H}{1\mathrm{kpc}}{\left(\frac{{\Sigma }_g}{1\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}\right)}^{-1}\mathrm{yr}.$$

$$t_{\mathrm{diff}}=1.97\times 10^6{\left(\frac{H}{\mathrm{kpc}}\right)}^2{\left(\frac{E_p}{\mathrm{GeV}}\right)}^{-\frac{1}{3}}\mathrm{yr}.$$

$$t_{\mathrm{adv}}=\frac{H}{v_w}=1.96\times 10^6\frac{H}{1\mathrm{kpc}}{\left(\frac{v_w}{500{\mathrm{km}\mathrm{s}}^{-1}}\right)}^{-1}\mathrm{yr},$$

$$t_{\mathrm{esc}}=\frac{1}{t_{\mathrm{diff}}^{-1}+t_{\mathrm{adv}}^{-1}}.$$

$$\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}=ϵ1.7\times 10^{-10}\frac{L_{\mathrm{IR}}}{L_{\odot }}$$

$$\begin{array}{r}{\Sigma }_{\mathrm{SFR}}=(2.5\pm 0.7)\times 10^{-4}{\left(\frac{{\Sigma }_g}{M_{\odot }{\mathrm{pc}}^{-2}}\right)}^{1.4\pm 0.15}\\ M_{\odot }{\mathrm{y}{\mathrm{r}}^{-1}\mathrm{kpc}}^{-2}\end{array}$$

$$\begin{array}{r}{\Sigma }_g=(10.43\pm 2.41)\times 10^{-5}{\left(\frac{L_{\mathrm{IR}}}{L_{\odot }}\right)}^{\frac{1}{1.4\pm 0.15}}\\ {\left(\frac{R}{\mathrm{pc}}\right)}^{-\frac{2}{1.4\pm 0.15}}{\mathrm{g}\mathrm{cm}}^{-2},\end{array}$$

$$\begin{array}{rl}{v}_w& =10^{1.56\pm 0.13}{\left(\frac{\mathrm{SFR}}{M_{\odot }{\mathrm{yr}}^{-1}}\right)}^{0.35\pm 0.06}\\ & =10^{1.56\pm 0.13}{\left(\frac{L_{\mathrm{IR}}}{5.8\times 10^9{L}_{\odot }}\right)}^{0.35\pm 0.06}{\mathrm{km}\mathrm{s}}^{-1}.\end{array}$$

$$\bar{R}=0.1{\left(\frac{1.26M_{\ast }}{M_{\odot }}\right)}^{0.14}{\left(1+\frac{1.26M_{\ast }}{3.98\times 10^{10}{M}_{\odot }}\right)}^{0.25}\mathrm{kpc},$$

$${\sigma }_{{\ln }_R}=0.34+\frac{0.13}{1+{\left(\frac{1.26M_{\ast }}{M_{\odot }}\right)}^2}\mathrm{kpc}.$$

$$M_{\ast }=10^{11}{\left(\frac{L_{\mathrm{IR}}}{5.38_{-2.69}^{+5.38}\times 10^{10}{L}_{\odot }}\right)}^{1.11}.$$

$$L_{\gamma }=L_{{\pi }_0}=C{f}_{\pi }{L}_{\mathrm{IR}},$$

$$f(x)=A{\left(\frac{x}{x_b}\right)}^{{\alpha }_1}{\left[\frac{1}{2}+\frac{1}{2}{\left(\frac{x}{x_b}\right)}^{1/\delta }\right]}^{({\alpha }_2-{\alpha }_1)\delta },$$

$$\delta ={\log }_{10}\frac{x_2}{x_b}={\log }_{10}\frac{x_b}{x_1}.$$

$${\log }_{10}{L}_{1-500\mathrm{GeV}}=A+\alpha ({\log }_{10}{L}_{\mathrm{IR}}-{\log }_{10}{L}_b)-\delta ({\alpha }_1-{\alpha }_2)\times {\log }_{10}\left[\frac{1}{2}+\frac{1}{2}{\left(\frac{L_{\mathrm{IR}}}{L_b}\right)}^{1/\delta }\right].$$

$${\log }_{10}{L}_{1-500\mathrm{GeV}}=A+\alpha ({\log }_{10}\mathrm{SFR}-{\log }_{10}{\mathrm{SFR}}_b)-\delta ({\alpha }_1-{\alpha }_2)\times {\log }_{10}\left[\frac{1}{2}+\frac{1}{2}{\left(\frac{\mathrm{SFR}}{{\mathrm{SFR}}_b}\right)}^{1/\delta }\right].$$

Full text

cosmic rays – gamma-rays: ISM – galaxies: star formation

1 Introduction

Nearby star-forming and starburst galaxies have been identified to be GeV-TeV gamma-ray sources (Acero et al., 2009; VERITAS Collaboration et al., 2009; Abdo et al., 2010; Ackermann et al., 2012). Cosmic rays (CRs) accelerated by supernova remnants or stellar winds interact with the interstellar medium (ISM) and produce neutral pions (schematically written as $p+p\rightarrow\pi^{0}$+other products), which in turn decay into high-energy gamma-rays ($\pi^{0}\rightarrow\gamma+\gamma$). Using the two-year observations obtained with the Large Area Telescope (LAT) aboard the Fermi mission, Abdo et al. (2010) found a correlation between the gamma-ray luminosity and star-formation-rate (SFR) for the Local Group galaxies. With a larger galaxy sample, Ackermann et al. (2012) found a quasi-linear scaling relation between gamma-ray luminosity and infrared luminosity which applies to both quiescent galaxies of the Local Group and low-redshift starburst galaxies. Since young stars in star-forming galaxies emit ultraviolet emission which is absorbed by dust in the ISM and reprocessed into infrared (IR) emission, the IR luminosity is a good indicator of SFR. Recently, the scaling relation was extended to an even higher IR luminosity regime, with detection of GeV emission from nearby luminous and ultraluminous infrared galaxies (Tang et al., 2014; Peng et al., 2016; Griffin et al., 2016).

The quasi-linear correlation between the gamma-ray luminosity $L_{\gamma}$ in the energy range $0.1-100\ \rm GeV$ and the total infrared luminosity $L_{\rm IR}$ ($8-1000\ \mu m$) has the form (Ackermann et al., 2012)

[TABLE]

where a nearly linear power-law index of $\alpha=1.0-1.2$ is obtained (Ackermann et al., 2012). Since the infrared luminosity reflects the cosmic ray injection rate in the galaxy, the quasi-linear relationship implies that a constant fraction of cosmic ray energy is converted into gamma-rays across all these galaxies. While, theoretically, extreme starburst galaxies, such as Arp 220, can be considered as cosmic ray calorimeter (e.g., (Torres, 2004; Lacki et al., 2011; Yoast-Hull et al., 2015; Wang & Fields, 2018)), i.e., almost all of the cosmic ray energy is converted into secondary particles due to high gas density, CRs in low-SFR galaxies are expected to transfer only a small fraction of energy into secondary particles due to a much lower gas density in these galaxies. For starburst galaxies (SBGs) such as M82, the optical depth for $pp$ interactions over typical SBG size is low and CRs escape mostly unscathed, e.g., for a high average gas density, $n_{\rm H}=10^{2}\ \rm cm^{-3}$, the mean free path of a few TeV proton is $65\ \rm kpc$, much larger than the typical disk radius (e.g., Persic et al. (2008)). A consequence is that these SBGs are inefficient calorimeters (at $\sim 30\%$ level, see Lacki et al. (2011)), and normal star-forming galaxies have even lower fractions of cosmic-ray energy converted into gamma-rays. Thus, the origin of this linear slope is puzzling.

The aim of this paper is to interpret this apparent quasi-linear relationship. We first calculate the gamma-ray production efficiency in star-forming galaxies by considering realistic galaxy properties, such as the gas density and galactic wind velocity in star-forming galaxies. It is shown that the efficiency of producing gamma-rays becomes increasingly low for low infrared luminosity galaxies. We further find that this theoretic model can reproduce the observed relation between the gamma-ray luminosity and infrared luminosity. Motivated by such a gradually changing efficiency, we also use an empirical function, i.e., a smoothly broken power-law function, to fit the observation data, and find that the scaling index at the low luminosity end deviates significantly from the linear relation, which is consistent with our theoretical expectation.

The rest of this paper is structured as follows. In Section 2, we describe the formulas of calculating the gamma-ray emission related to IR luminosity. In section 3, we present the results. We give our discussions in Section 4 and conclusions in Section 5.

2 Gamma-ray production in star-forming galaxies

Massive stars in galaxies end their lives as core-collapse supernovae, whose remnants can accelerate CRs and inject them into the ISM. CRs produce high-energy gamma rays through inelastic collisions with ISM ($pp$ collisions), electron and positron bremsstrahlung, and inverse Compton (IC) scattering of the primary and secondary electrons. Detailed calculations have shown that the pionic decay gamma-rays dominate the emission above $100\ \rm MeV$ for star-froming galaxies (e.g., Domingo-Santamaría & Torres (2005); Rephaeli et al. (2010)), although leptonic emission is expected to become increasingly important at lower energies.

Therefore, we expect the injection power of CRs is proportional to SFR, i.e., $L_{\rm CR}\propto{\rm SFR}$. The gamma-ray luminosity produced by $pp$ collisions between CRs and ISM can be parameterized by $L_{\gamma}\propto f_{\pi}\,L_{\rm CR}$ , where $f_{\pi}$ is the efficiency of CR energy transferred to secondary pions. On the other hand, CRs can escape from galaxy through diffusion and galactic wind advection. These two processes compete to generate the efficiency

[TABLE]

where $t_{\rm esc}$ is the escape time of CRs and $t_{\rm loss}$ is the energy loss time of CRs via $pp$ collisions.

The CR energy loss timescale depends on the average gas density and the inelastic collision cross section, given by $t_{\rm loss}=(0.65n\sigma_{pp}c)^{-1}$ (Mannheim & Schlickeiser, 1994), where $\sigma_{pp}$ is the inelastic cross section. For CR energy larger than $1\ \rm GeV$, $\sigma_{pp}$ is nearly a constant with energy, and we take the value $2.5\times 10^{-26}\ \rm cm^{2}$. The ISM number density $n$ relates with the gas surface density $\Sigma_{g}$ by $\Sigma_{g}=nm_{p}H$, where $m_{p}$ is mass of proton, $H$ is the scale height of galaxy disk. Then the CR energy loss time is given by

[TABLE]

CRs escape out of the galaxy in the form of diffusion or galactic wind advection. In the case of diffuse process, CRs are scattered by small-scale inhomogeneous magnetic fields. Diffusive time scales can be approximated as $t_{\rm diff}=H^{2}/{4D}$, where $D=D_{0}\,(E/E_{0})^{\delta}$ is the diffusion coefficient, $D_{0}$ and $E_{0}$ are normalization factors. We take the standard diffusion coefficient for ISM as $D=3.86\times 10^{28}(E_{p}/{\rm GeV})^{1/3}\ {\rm cm^{2}\ s^{-1}}$, where the value of $\delta=1/3$ is for the Kolmogorov type turbulence. Then the diffusion time is given by

[TABLE]

In the case of advection process, CRs are transported outward with galactic wind on a characteristic timescale

[TABLE]

where $v_{w}$ is the velocity of the galactic wind. The timescale for CRs escaping out of the galaxy is parameterized as

[TABLE]

Next, we calculate the pion production efficiency by considering realistic galaxy properties, such as gas surface density $\Sigma_{g}$, galactic wind velocity $v_{w}$ and galaxy scale height $H$. The total IR luminosity in 8-1000$\mu m$ is one well-established tracer of the SFR for late-type galaxies, so we take (Kennicutt, 1998)

[TABLE]

where the factor $\epsilon=0.79$ is for initial mass function (IMF) derived from Chabrier (2003), and $\epsilon=1$ is for initial mass function used by Kennicutt (1998). Below we use the IR luminosity instead of SFR in the subsequent calculation.

To determine the gas surface density $\Sigma_{g}$, we use the classical Kennicutt-Schmidt law (K-S law) that extends over several orders of magnitude in SFR and gas density (Kennicutt, 1998):

[TABLE]

where $\Sigma_{\rm SFR}$ is disk-averaged SFR density. $\Sigma_{\rm SFR}$ can be obtained from $\Sigma_{\rm SFR}={\rm SFR}/\pi R^{2}$. Replacing SFR with the IR luminosity, K-S law can be written as

[TABLE]

where $\epsilon=1$ for the IMF used by Kennicutt (1998) is adopted. Galactic winds or gaseous outflows from starburst galaxies can accelerate CRs to high speed. For luminous infrared galaxies at low redshift, winds from more luminous starbursts have higher speeds roughly as $v_{w}\propto{\rm SFR}^{0.35}$ (Martin, 2005). A similar relation is found for star-forming galaxies at $z=1.4$ (Weiner, 2009). Thus, we have

[TABLE]

Considering that the relation between galaxy scale height and SFR is not straightforward, we assume galaxy height relates with galaxy radius as $H=(0.21\pm 0.02)\ R$ (Padilla & Strauss, 2008), where $R$ is the disk radius of a galaxy. For late-type galaxies, there is a relation between the radius of galaxy and the total stellar mass, given by (Shen et al., 2003)

[TABLE]

with the dispersion of

[TABLE]

For star-forming galaxies in local universe, based on the tight relationship between total stellar mass and SFR of galaxy (Peng et al., 2010), we get ${\rm SFR}=7.94(\rm M/{10^{11}\rm M_{\odot}})^{0.9}\ \rm M_{\odot}\ {\rm yr}^{-1}$, where the dispersion is 0.3 dex. The SFR used in Peng et al. (2010) is computed for the Kroupa IMF, so we convert it to the case for the Chabrier IMF by using ${\rm log}({\rm SFR}_{\rm Chabrier})={\rm log}({\rm SFR}_{\rm Kroupa})-0.04$. Then the total stellar mass relates with the IR luminosity by

[TABLE]

Using the above relations, we take galaxy parameters, such as gas surface density $\Sigma_{g}$ , galactic wind velocity $v_{w}$ and galaxy scale height $H$, as input parameters and generate the pion production efficiency $f_{\pi}$ as an output, which is in turn used to produce the $L_{\gamma}-L_{\rm IR}$ relation

[TABLE]

where $L_{\gamma}$ is the total energy of gamma-rays as integrating over corresponding energy range, $L_{\pi_{0}}$ is the total energy of $\pi_{0}$ produced by $pp$ collisions and $C$ is a normalization factor. The normalization factor $C$ can be obtained using the data of M82.

3 Results

The pion production efficiency $f_{\pi}$ is shown in Fig. 1. The uncertainty given in Fig. 1 takes into account of all the uncertainties in the relation between the specific galaxy parameter and IR luminosity. One can see that the $pp$ interaction is quite inefficient in low IR luminosity galaxies, with efficiency being less than $5\%$ for $L_{\rm IR}<10^{8}L_{\odot}$. For galaxies with higher IR luminosity $>10^{11}L_{\odot}$, the efficiency increases to about $25\%$, albeit with a somewhat larger uncertainty.

Now we compare our theoretical prediction with the observed gamma-ray luminosity of star-forming galaxies (Fig. 2). Considering that gamma-rays at $100\ \rm MeV$ may be contaminated by leptonic emission or point sources like pulsars, we use a higher threshold energy of $1\ \rm GeV$ for our study. It is believed that $1-500\ \rm GeV$ gamma-ray emission of star-forming galaxies and star-burst galaxies has a hadronic origin. The parameters of the GeV-detected star-forming galaxies, including $L_{1-500\ \rm GeV}$, $L_{\rm IR}$ and SFR, are listed in Table 1. The data of $L_{1-500\ \rm GeV}$ are taken from Peng et al. (2019). From Fig. 2, one can see that the theoretical model for the gamma-ray luminosity ($L_{1-500\ \rm GeV}$) (the blue lines) is in compatible (within the low available statistics) with observations, especially for galaxies with IR luminosity $L_{\rm IR}$ above $10^{9}L_{\odot}$. The slope for the relation between the gamma-ray luminosity and IR luminosity is not a constant. For galaxies with higher IR luminosity, the slope of the curve approaches to $1$, but it steepens to $\sim 1.6$ for low infrared luminosity galaxies, which reflects an increasingly lower efficiency for the production of gamma-ray emission (see Fig. 1). This means that high luminosity galaxies are closer to (but still short of) being CR calorimeters (shown by the gray line), while low luminosity galaxies deviate from the calorimetric limit significantly. The model can explain the data obtained in Ackermann et al. (2012) and Peng et al. (2016). We will discuss this in more detail in the next section.

4 Discussions

As shown in Fig 2, the Small Magellanic Cloud (SMC) is an apparent outlier of the correlation. We suggest that this could be due to the underestimate of SFR for SMC. The IR luminosity can be regarded as a well-established tracer of SFR only when the IR emission of interstellar dust is nearly calorimetric measure of radiation produced by young stellar populations (Bell, 2003). That is, for galaxies with IR luminosity larger than $10^{9}L_{\odot}$, the IR luminosity is a robust SFR indicator, but for low IR luminosity galaxies such as SMC, the substantially low metallicity and dust content lead to a low optical depth for IR photons. As a result, the observed IR emission only reflects a fraction of the star formation activity, and thus the SFR of these galaxies are underestimated. After considering the combined observations of $H_{\alpha}$ and IR emission, we estimate that the SFR of SMC is $0.04-0.08\ M_{\odot}\ {\rm yr}^{-1}$ (Wilke et al., 2004). Using Eq.(7) to convert the SFR into IR luminosity, we find that the IR luminosity of SMC is corrected to $\sim 4.5\times 10^{8}L_{\odot}$. The new data for SMC is shown by a red dot in Fig. 2 and one can see that it agrees well with our theoretical model.

The finding of a changing slope for the relation between the gamma-ray luminosity and the total IR luminosity suggests that the mechanism dominating the CR energy loss process is different for different IR-luminosity galaxies. We now discuss the physical process accounting for the scaling slope change. The CR energy loss timescale through $pp$ collision and escape timescales through diffusion and advection are shown in Fig. 3. One can see that the diffusion escape is the fastest when the IR luminosity is less than $10^{11}L_{\odot}$, leading to a low efficiency of pion production and a steep slope. This situation changes when the IR luminosity is larger than $10^{11}L_{\odot}$. The diffusion timescale increases sharply due to the increase of the galaxy disk scale height. The advection escape becomes increasingly important in more luminous star-forming galaxies due to higher galactic wind velocities. For the timescale of $pp$ collision, it declines first due to increases of the gas density and then rise due to the relatively larger galaxy disk scale height. The timescales of the advection and the energy loss through $pp$ collisions are on the same order of magnitude, thus CRs collide with ISM effectively before escaping out of the galaxy.

Since the derived relation between $L_{1-500\ \rm GeV}$ and $L_{\rm IR}$ indicates a gradually changing slope, we try to use a smoothly broken power law (SBPL) function to fit the data. For more details about the SBPL fit, please refer to the Appendix. As it is difficult to constrain all the parameters of the SBPL model for a small number of the sample galaxies, we fix the slope at the high IR-luminosity end as $\alpha=1$, which is motivated by the calorimetric limit for these galaxies. We further remove SMC from the sample as its IR luminosity is not an accurate SFR tracer. We use the maximum likelihood approach as illustrated in our previous work (Peng et al., 2019). We find that the slope $\alpha$ for low IR-luminosity galaxies approaches to $1.56\pm 0.15$ as the IR luminosity decreases. The result is shown in Fig. 2 by the red lines. The slope found for this empirical fit function is well consistent with our theoretical model.

As mentioned above, the IR luminosity underestimates SFR for low IR-luminosity galaxies. To eliminate this issue, we study the relation using SFR directly. The result is shown in Fig. 4. The slope changes from $1$ at high SFR end to $1.41\pm 0.15$ at low SFR end for a fit with a smoothly broken power-law function. We compare the maximum likelihood of the two fits with a smoothly broken power-law function and a single power-law function, and find that the difference in the maximum likelihood between the two models is not significant ($<2$). From this point of view, we think that the smoothly broken power law function is equally acceptable by the data.

5 Conclusions

We calculate the gamma-ray luminosity of star-forming galaxies by considering realistic galaxy properties, such as gas surface density, galaxy scale height and wind velocity. The derived relation between the gamma-ray luminosity and IR luminosity shows a gradually changing slope, as expected from increasingly lower efficiency of gamma-ray production for low luminosity galaxy. We further find that the measured data is well consistent with such a changing slope. As an comparison, we use a smoothly broken power law function to fit the data and find that the slope for low IR luminosity galaxies deviates from the linear relation significantly, but agrees well with our theoretical calculation. Our result suggests that most CRs escape before significant collisional losses in low-SFR galaxies, mainly through the diffusion process. This is consistent with the findings in some recent numerical simulations (Pfrommer et al., 2017; Chan et al., 2018).

Acknowledgments

X.Y.W. is supported by the National Key R & D program of China under the grant 2018YFA0404203 and the NSFC grants 11625312 and 11851304. F.K.P acknowledges support from the Doctoral Starting up Foundation of Guizhou Normal University 2017 (GZNUD[2017] 33).