Breaking the Radio - Gamma-Ray Connection in Arp 220

TL;DR

This paper investigates the discrepancy between gamma-ray and radio observations in Arp 220, proposing that increased cosmic ray energy injection and infrared optical depth can reconcile the data while maintaining consistency with infrared observations.

Contribution

It introduces a model that accounts for higher cosmic ray energy injection and optical depth to explain gamma-ray and radio emission discrepancies in Arp 220.

Findings

Higher cosmic ray injection rates are needed to match gamma-ray flux.

Infrared optical depth greater than one affects cosmic ray energy losses.

Optically thick medium explains gamma-ray and radio data simultaneously.

Abstract

Recent analyses of the gamma-ray spectrum from the ultra-luminous infrared galaxy Arp 220 have revealed a discrepancy in the cosmic ray energy injection rates derived from the gamma-rays versus the radio emission. While the observed radio emission is consistent with the star formation rate inferred from infrared observations, a significantly higher cosmic ray population is necessary to accurately model the measured gamma-ray flux. To resolve this discrepancy between the radio and gamma-ray observations, we find that we must increase the cosmic ray energy injection rate and account for an infrared optical depth greater than unity. Raising the energy injection rate naturally raises the total gamma-ray flux but also raises the radio flux unless there is also an increase in the energy loss rate for cosmic ray leptons. A optically thick medium results in an increase in energy losses via inverse Compton for cosmic ray leptons and preserves agreement with submillimeter, millimeter, and infrared wavelength observations.

Tables

Table 1. Table 1: Fiducial Inputs

Parameters	East	West	CND	Ref.
Radius	85	65	15	1,2
$R$ [pc]
Scale Height	30	30	10	1
$h$ [pc]
Column Density	${10}^{25}$	${10}^{25}$	${10}^{26}$	1,3
$N_{H_2}$ [cm-2]
Molecular Gas Massa	10	5.9	5.9	4
$M_{H_2}$ [${10}^8{M}_{\odot }$]
Average ISM Densitya	1.1	1.2	65
$n_{H_2}$ [${10}^4$ cm-3]
FIR Luminosity	3.0	3.0	6.0	1,3
$L_{\mathrm{IR}}$ [${10}^{11}{L}_{\odot }$]
FIR Energy Densitya	2.8	4.2	150
$U_{\mathrm{IR}}$ [${10}^4$ eV cm-3]
Dust Temperaturea	50	55	135
$T_{\mathrm{IR}}$ [K]
Notes. aDerived from above parameters; References:
(1) Scoville et al. (2017); (2) Barcos-Muñoz et al. (2015);
(3) Wilson et al. (2014); (4) Bolatto et al. (2013)

Table 2. Table 2: Cosmic Ray Energy Injection Rates

	East	West	CND
Scenario A (Protons & Electrons)
Range	$0.5-10$	$0.5-10$	$1.0-20$
Step Size	0.5	0.5	1.0
Scenario B (Protons Only)
Range	1.0	1.0	$2.0-40$
Step Size	0	0	2.0
Scenario C (Electrons Only)
Range	1.0	1.0	${10}^1-{10}^5$
Step Size	0	0	${10}^{0.2}$
Notes: In all cases, the units for the injection rates
are ${10}^{50}$ erg yr-1.

Table 3. Table 3: Data for Star-Forming Galaxies with Gamma-Ray Detections

Galaxy	$S_{60\mathrm{\mu m}}$	${\sigma }_{60\mathrm{\mu m}}$	$S_{100\mathrm{\mu m}}$	${\sigma }_{100\mathrm{\mu m}}$	$S_{1.4\mathrm{GHz}}$	${\sigma }_{1.4\mathrm{GHz}}$	$F_{0.1-100\mathrm{GeV}}$	${\sigma }_{0.1-100\mathrm{GeV}}$
	(Jy)	(mJy)	(Jy)	(mJy)	(mJy)	(mJy)	(erg cm-2 s-1)	(erg cm-2 s-1)
NGC 0224	595	$2.8\times {10}^4$	3073	$3.16\times {10}^5$	8600	1100	$5.88\times {10}^{-12}$	$0.863\times {10}^{-12}$
NGC 0253	968	65	1290	644	2990	114	$8.65\times {10}^{-12}$	$9.22\times {10}^{-13}$
NGC 1068	196	108	257	228	4850	170.4	$6.97\times {10}^{-12}$	$8.55\times {10}^{-13}$
NGC 2146	147	63	194	240	1070	39.9	$3.58\times {10}^{-12}$	$6.75\times {10}^{-13}$
NGC 3034	1480	74	1370	694	6200	209	$1.07\times {10}^{-11}$	$8.37\times {10}^{-13}$
NGC 4945	625	133	1330	797	6450	250	$1.38\times {10}^{-11}$	$1.24\times {10}^{-12}$
CIRCINUS	249	24.9	316	31.6	1370	–	$1.30\times {10}^{-11}$	$2.01\times {10}^{-12}$
ARP 220	104	112	115	138	327	9.8	$1.92\times {10}^{-12}$	$4.30\times {10}^{-13}$
SMC	$6.69\times {10}^3$	–	$1.50\times {10}^4$	–	–	–	$2.90\times {10}^{-11}$	$1.0\times {10}^{-12}$
LMC	$8.29\times {10}^4$	–	$1.85\times {10}^5$	–	–	–	$2.24\times {10}^{-10}$	–
Notes: NGC 0224 = M31, NGC 3034 = M82.

Table 4. Table 4: Data for Star-Forming Galaxies without Gamma-Ray Detections

Galaxy	$S_{60\mathrm{\mu m}}$	${\sigma }_{60\mathrm{\mu m}}$	$S_{100\mathrm{\mu m}}$	${\sigma }_{100\mathrm{\mu m}}$	$S_{1.4\mathrm{GHz}}$	${\sigma }_{1.4\mathrm{GHz}}$	$F_{0.1-100\mathrm{GeV}}$
	(Jy)	(mJy)	(Jy)	(mJy)	(mJy)	(mJy)	(${10}^{-12}$ erg cm-2 s-1)
IC 1623	22.9	62	31.6	113	249	9.8	0.481
NGC 0520	31.5	30	47.4	146	176	5.3	0.396
NGC 0660	65.5	88	115	134	373	13.9	1.06
NGC 0695	7.59	31	13.6	167	74.8	3.1	1.36
NGC 0828	11.5	53	25.3	113	104	3.7	0.715
MRK 1027	5.28	39	8.57	108	53.2	1.6	0.521
NGC 0891	66.5	71	172	353	239	7.9	2.04
NGC 1022	19.7	35	27.3	111	46.7	1.5	1.12
NGC 1055	23.4	51	65.3	82	199	6.7	1.56
NGC 1144	5.06	38	11.5	198	155	5.4	0.328
NGC 1365	94.3	33	166	126	376	13	1.33
IC 0342	181	70	392	434	191	7.3	1.83
NGC 1530	9.88	43	25.8	128	67.3	2.7	1.77
NGC 1614	32.1	83	34.3	430	137	4.9	1.14
VII Zw 31	5.51	31	10.1	126	41.4	1.3	0.151
NGC 2276	14.3	58	29.0	133	267	8.6	0.522
ARP 055	6.07	48	10.3	109	36.8	1.2	0.583
NGC 2903	60.5	52	130	208	445	13.9	1.08
UGC 05101	11.7	34	19.9	137	170	5.8	0.383
NGC 3079	50.7	64	105	141	769	27.1	1.40
NGC 3147	8.17	38	29.6	222	89.9	3.4	0.192
ARP 148	6.38	34	10.3	106	36.4	1.2	0.797
NGC 3556	32.6	68	76.9	231	217	7.1	1.11
NGC 3627	66.3	59	137	118	453	11.2	1.01
NGC 3628	54.8	77	106	240	291	8.7	0.866
ARP 299	113	52	111	133	677	25.4	1.08
NGC 3893	15.6	54	36.8	63	139	5	0.431
NGC 4030	18.5	46	50.9	117	154	5.3	0.399
NGC 4041	14.2	19	31.7	126	102	3.8	0.875
NGC 4414	29.6	80	70.7	72	240	7.9	0.578
NGC 4631	85.4	62	160	260	446	14.1	0.635
MRK 0231	30.8	42	29.7	108	309	12.1	0.651
NGC 4826	36.7	77	81.7	99	99.9	3.7	0.755
NGC 5005	22.2	21	63.4	131	181	6.2	1.11
NGC 5055	40.0	49	140	356	350	–	0.626
ARP 193	17.0	88	24.4	120	104	3.2	0.983
NGC 5135	16.9	46	31.0	177	200	7.1	1.28
NGC 5194	97.4	193	221	329	1570	–	0.723
NGC 5236	266	64	524	302	405	15	0.609
MRK 0273	22.5	42	22.5	70	145	5.1	0.449
NGC 5678	9.67	29	25.7	85	110	3.9	0.496
NGC 5713	22.1	65	37.3	88	158	5.7	0.203
NGC 5775	23.6	46	55.6	97	280	9.1	0.680
NGC 6240	22.9	54	26.5	174	426	15	0.269
NGC 6701	10.1	28	20.1	67	88.9	3.1	0.426
NGC 6946	130	71	291	458	187	6.7	0.551
IC 5179	19.4	62	37.3	90	169	6.2	0.293
NGC 7130	16.7	44	25.9	132	190	7.5	0.227
NGC 7331	45.0	91	110	468	328	8.1	1.19
NGC 7469	27.3	40	35.2	599	181	5.4	0.847
NGC 7479	14.9	49	26.7	247	99.0	3.7	1.95
NGC 7771	19.7	135	40.1	839	141	5.2	0.905
MRK 0331	18	46	22.7	194	70.7	2.2	0.303
Notes: NGC 5194 = M51, NGC 5236 = M83.

Equations

$$q_p(E_p)=\frac{{\Gamma }_p-2}{E_{\min }^{-{\Gamma }_p}}\times \frac{\eta {\epsilon }_{\mathrm{CR}}}{V}\times {\left(\frac{E_p}{1\mathrm{GeV}}\right)}^{-{\Gamma }_p},$$

$$j_{\nu }^{\mathrm{Sy}}={\int }_{E_{\min }}^{\infty }d{E}_e{N}_e(E_e)P_{\nu }^{\mathrm{Sy}}({\gamma }_e),$$

$$P_{\nu }^{\mathrm{Sy}}({\gamma }_e,\psi )=\frac{\sqrt{3}{q}_e^{3}B}{m_e{c}^2}\sin (\psi )F(x).$$

$$P_{\nu }^{\mathrm{Sy}}({\gamma }_e)=\frac{32\pi c}{9\Gamma (4/3)}{\left(\frac{q_e^2}{m_e{c}^2}\right)}^2{U}_B{\beta }_e^2{\gamma }_e^2\frac{{\nu }^{1/3}}{{\nu }_c^{4/3}}{e}^{-\nu /{\nu }_c},$$

$$N(E)=q(E)\times t(E).$$

$$q_e^{\mathrm{tot}}(E_e)=q_e^{\mathrm{prim}}(E_e)+q_e^{\sec }(E_e),$$

$$q_e^{\mathrm{prim}}=\frac{{\Gamma }_e-2}{E_{\min }^{-{\Gamma }_e}}\times \frac{\xi \eta {\epsilon }_{\mathrm{CR}}}{V}\times {\left(\frac{E_e}{1\mathrm{GeV}}\right)}^{-{\Gamma }_e}.$$

$$q_e^{\sec }(E_e)\propto \int d{E}_p{N}_p(E_p){\sigma }_{\mathrm{p}-\mathrm{p}}(E_p).$$

$$t_e(E_e{)}^{-1}=t_{\mathrm{loss}}(E_e{)}^{-1}+t_{\mathrm{diff}}(E_e{)}^{-1}+t_{\mathrm{adv}}^{-1}.$$

$$t_e(E_e)=t_{\mathrm{loss}}(E_e)\simeq E_e/(-dE/dt).$$

$$-{\left(\frac{d{E}_e}{dt}\right)}_{\mathrm{Ion}}=\frac{9}{4}{m}_e{c}^3{\sigma }_T{n}_{\mathrm{mol}}\left[6.85+\ln ({\gamma }_e)\right],$$

$$-{\left(\frac{d{E}_e}{dt}\right)}_{\mathrm{Br}}^{\mathrm{ion}}=\frac{3\alpha }{2\pi }c{\sigma }_T{n}_{\mathrm{ion}}Z(Z+1)E_e$$

$$\times \left[\ln (2{\gamma }_e)-\frac{1}{3}\right],$$

$$-{\left(\frac{d{E}_e}{dt}\right)}_{\mathrm{Br}}^{\mathrm{mol}}=\frac{3.9\alpha }{8\pi }c{\sigma }_T{\varphi }_{\text{Hi}}^{s-s}{n}_{\mathrm{mol}}{E}_e,$$

$$-{\left(\frac{d{E}_e}{dt}\right)}_{\mathrm{IC}}=\frac{4}{3}c{\sigma }_T{U}_{\mathrm{rad}}\frac{{\gamma }_{\mathrm{K}}^2{\gamma }_e^2}{{\gamma }_{\mathrm{K}}^2+{\gamma }_e^2},$$

$$-{\left(\frac{d{E}_e}{dt}\right)}_{\mathrm{Sy}}=\frac{4}{3}c{\sigma }_T\frac{B^2}{8\pi }{\gamma }_e^2,$$

$$M_{{\mathrm{H}}_2}=\pi RH\times m_{{\mathrm{H}}_2}{N}_{{\mathrm{H}}_2}\frac{X_{\mathrm{CO}}^{\mathrm{ULIRG}}}{X_{\mathrm{CO}}^{\mathrm{MW}}}.$$

$$U_{\mathrm{IR}}=\frac{L_{\mathrm{IR}}}{4\pi R^2\times c}.$$

$$T_{\mathrm{IR}}={\left(\frac{L_{\mathrm{IR}}}{(2\pi R^2+2\pi RH)\times ca}\right)}^{1/4},$$

$$U_{\text{rad}}$$

$$T_{\mathrm{rad}}$$

$$-{\left(\frac{dE}{dt}\right)}_{\mathrm{IC}}$$

$$q_{\mathrm{IC}}(E_{\gamma })$$

$${\tau }_{\mathrm{IR}}\approx {\tau }_{\nu }{\left(\frac{\lambda (\nu )}{60\mu m}\right)}^{1.6}.$$

$$F_{\mathrm{FIR}}=1.26\times 10^{-11}(2.58S_{60}+S_{100}),$$

$$q_{\mathrm{Br}}=c{n}_{\mathrm{ISM}}{\int }_{E_{\min }}^{\infty }d{E}_e{N}_e(E_e)\frac{d\sigma }{d{E}_{\gamma }},$$

$$\frac{d\sigma }{d{E}_{\gamma }}=\frac{3\alpha {\sigma }_T}{8\pi E_{\gamma }}\left[\left[1+{\left(1-\frac{E_{\gamma }}{E_e}\right)}^2\right]{\varphi }_1+\left(1-\frac{E_{\gamma }}{E_e}\right){\varphi }_2\right].$$

$${\varphi }_u=4\left[\text{ln}\left(\frac{2E_e}{m_e{c}^2}\left(\frac{E_e-E_{\gamma }}{E_{\gamma }}\right)\right)-\frac{1}{2}\right].$$

$$n_{ph}(ϵ)=\frac{C_{\mathrm{dil}}}{{\pi }^2{\hslash }^3{c}^3}\frac{{ϵ}^2}{e^{ϵ/k{T}_{\mathrm{rad}}}-1}{\left(\frac{ϵ}{{ϵ}_0}\right)}^{\beta },$$

$$U_{\mathrm{rad}}$$

$$=\frac{C_{\mathrm{dil}}}{{\pi }^2{\hslash }^3{c}^3}{ϵ}_0^{-\beta }{\int }_0^{\infty }\frac{{ϵ}^{3+\beta }}{e^{ϵ/k{T}_{\mathrm{rad}}}-1}dϵ$$

$$=\frac{C_{\mathrm{dil}}{ϵ}_0^{-\beta }}{{\pi }^2{\hslash }^3{c}^3}(k{T}_{\mathrm{rad}}{)}^{4+\beta }\Gamma (4+\beta )\zeta (4+\beta ),$$

$$C_{\mathrm{dil}}=\frac{{\pi }^2{\hslash }^3{c}^3{ϵ}_0^{\beta }{U}_{\mathrm{rad}}}{(k{T}_{\mathrm{rad}}{)}^{4+\beta }\Gamma (4+\beta )\zeta (4+\beta )},$$

$$n_{ph}(ϵ)=\frac{U_{\mathrm{rad}}}{(k{T}_{\mathrm{rad}}{)}^{4+\beta }\Gamma (4+\beta )\zeta (4+\beta )}\times \frac{{ϵ}^{2+\beta }}{e^{ϵ/k{T}_{\mathrm{rad}}}-1}.$$

$$q_{\mathrm{IC}}(E_{\gamma })=\frac{3}{4}c{\sigma }_T{\int }_0^{\infty }dϵ\frac{n_{ph}(ϵ)}{ϵ}{\int }_{E_{\min }}^{\infty }d{E}_e\frac{N_e(E_e)}{{\gamma }_e^2}F(q,\Gamma ),$$

Full text

1 Introduction

Over the past decade, eight star-forming galaxies (mostly giant spirals) have been detected in $\gamma$-rays by Fermi, including M31 (NGC 0224), NGC 0253, NGC 1068, NGC 2146, M82 (NGC 3034), NGC 4945, Arp 220, and Circinus (Abdo et al., 2010b; Abdo et al., 2010a; Lenain et al., 2010; Hayashida et al., 2013; Tang et al., 2014; Peng et al., 2016; Griffin et al., 2016). Observations of these galaxies show a tentative correlation between their $\gamma$-ray and far-infrared (FIR) luminosities (Ackermann et al., 2012; Rojas-Bravo & Araya, 2016), indicating a physical relationship between these two quantities in star-forming galaxies as predicted by Thompson et al. (2007). Recent simulations show further evidence of a link between the FIR emission and $\gamma$-rays from galaxies via their star formation rate (SFR; e.g., Martin, 2014; Pfrommer et al., 2017).

The essential physical processes underlying this correlation are tied to massive stars, which dominate the diffuse emission in star-forming galaxies. Massive stars radiate predominately in the ultraviolet; much of this emission is absorbed by the surrounding interstellar medium (ISM) and re-radiated in the infrared. At the ends of their lives, massive stars explode as supernovae, and the resulting remnants accelerate electrons and protons to relativistic energies (e.g., de Jong et al., 1985; Helou et al., 1985). Cosmic ray electrons interacting with galactic magnetic fields produce non-thermal synchrotron emission, while cosmic ray protons collide with the ISM to produce neutral pions that subsequently decay into $\gamma$-rays.

In our current understanding, the correlation between radio and FIR emission seen in star-forming galaxies requires two additional conditions to hold. First, either star-forming galaxies must be cosmic ray electron calorimeters, meaning that cosmic ray electrons are confined to the galaxy and radiate the majority of their energy within their injection volume (Voelk, 1989), or both cosmic ray electrons and ultraviolet radiation escape from star-forming galaxies in proportional amounts (Lacki et al., 2010). Secondly, in most star forming galaxies, energy losses via synchrotron must be larger than losses via inverse Compton for cosmic ray electrons, implying that the energy density in magnetic fields must be slightly larger than the energy density in radiation for star-forming galaxies ($U_{B}\gtrsim U_{\rm rad}$; Lisenfeld et al., 1996; Thompson et al., 2006).

While star-forming galaxies are generally believed to be electron calorimeters, most normal star-forming galaxies are not proton calorimeters (e.g., Pfrommer et al., 2017). In quiescent systems, the diffusion time-scale for cosmic ray protons is significantly shorter than their energy loss time-scale, and thus, cosmic ray protons preferentially escape into the halo. Additionally, these are systems in which primary cosmic ray electrons are likely dominant over secondary cosmic ray electrons and positrons produced in interactions of cosmic ray protons with the ISM.

In contrast, intensely star-forming systems (such as Arp 220) are at least partial proton calorimeters (e.g., Lacki et al., 2011; Yoast-Hull et al., 2013; Eichmann & Becker Tjus, 2016) and secondary cosmic rays dominate over primary cosmic ray electrons (e.g., Lacki & Thompson, 2013; Peretti et al., 2018). As a result, in these systems, we would expect both the non-thermal radio and the $\gamma$-ray emission to be proportional to the population of cosmic ray protons. Thus, for starburst systems in which secondary cosmic rays dominate, there may be a constant ratio between the $\gamma$-ray and radio fluxes, regardless of whether the galaxy is a complete proton calorimeter.

Observations clearly show that strong correlations exist between the FIR, radio, and $\gamma$-ray emission in star-forming galaxies (e.g., Ackermann et al., 2012; Schöneberg et al., 2013; Rojas-Bravo & Araya, 2016), see Fig. 11 in Appendix A. Examining the ratio of $\gamma$-ray flux to radio fluxes in star-forming galaxies detected with Fermi, we find tentative evidence for a constant $\gamma$-ray to radio flux ratio in starbursts galaxies, see Fig. 1. The ratio for normal star-forming galaxies is lower than that for starbursts, as expected due to the dominance of primary cosmic ray electrons and the lack of proton calorimetry in quiescent galaxies. In contrast, the ratio of $\gamma$-ray to radio flux for the Seyfert galaxy Circinus is significantly higher than the starburst ratio, likely due to the dominance of the central active galactic nucleus (AGN; Wojaczyński et al., 2015).

Fig. 1 also shows Arp 220, the closest ultra-luminous infrared galaxy (ULIRG), as slightly more than a standard deviation away from the mean $\gamma$-ray to radio luminosity ratio for starbursts. In investigating Arp 220 in more detail, modeling of the $\gamma$-ray and radio spectra reveal a larger cosmic ray energy injection rate is necessary to reproduce the observed $\gamma$-ray flux than is necessary for the radio flux, see Fig. 2. Of course, the magnitude of this discrepancy between the inferred energy injection rates depends on assumptions about the cosmic ray populations in the nuclei of Arp 220 (Yoast-Hull et al., 2017).

To examine these assumptions in detail, we begin by defining the injection spectrum for cosmic ray protons such that

[TABLE]

where $\Gamma_{p}$ is the cosmic ray (injection) spectral index, $\varepsilon_{\rm CR}$ is the cosmic ray energy injection rate, and $\eta=0.1$ is the acceleration efficiency of cosmic ray protons. The critical factor here, besides the energy injection rate (see Fig. 3), is the cosmic ray proton spectral index: as the spectral index steepens, the energy injection rate required to supply a set $\gamma$-ray flux increases (because the factor $1/E_{\rm min}^{\Gamma_{p}}$ decreases).

Thus, for flatter spectral indices ($\Gamma_{p}\sim 2.1-2.2$), there exist models where the radio and $\gamma$-ray fluxes imply cosmic ray energy injection rates which are consistent with each other (e.g., Torres, 2004; Lacki et al., 2011; Lacki & Thompson, 2013), and for steeper spectral indices ($\Gamma_{p}\sim 2.3-2.5$), the energy injection rates derived from the $\gamma$-ray and radio fluxes are inconsistent with each other (e.g., Yoast-Hull et al., 2017). As the observed $\gamma$-ray spectrum is best fit by spectral indices of $\Gamma_{\gamma}=2.35\pm 0.16$ (Peng et al., 2016), both sets of models are consistent with the observations within the statistical uncertainty on the $\gamma$-ray spectral index111It should be noted here that $\Gamma_{p}\sim\Gamma_{\gamma}$, as the cross section for neutral pion production depends on energy such that $\sigma_{p-p}(E_{p})\propto lnE_{p}$ (e.g., Kelner et al., 2006).. Throughout the rest of this paper, we will be concentrating on solutions to the tension between the radio and $\gamma$-ray observations for steeper spectral indices.

In this paper, we examine possible mechanisms which can resolve the discrepancy between the observed radio and $\gamma$-ray fluxes from Arp 220. Using our previously developed models for starburst galaxies (Yoast-Hull et al., 2013, 2014, 2015), we explore how changes in energy injection rate, infrared optical depth, and magnetic field strength affect the cosmic ray electron populations and the resulting non-thermal emission, see Fig. 3. We find that accounting for an infrared optical depth greater than unity in the models, as required by submillimeter and millimeter wavelength observations, easily resolves the tension between the energy injection rates derived from the radio and $\gamma$-ray spectra, respectively.

This paper is structured as follows: In Section 2, we lay out the equations relating to the cosmic ray spectra and radio emission. In Section 3, we discuss observational constraints on the physical quantities that are input into the models. In Section 4, we describe how our models are constructed and tested, and we report our results. In Section 5, we compare our results with observed and other quantities. Our conclusions are summarized in Section 6.

2 Theoretical Framework

To resolve the difference between the cosmic ray energy injection rates inferred from $\gamma$-ray and radio spectra, we must lower the radio flux while maintaining or raising the $\gamma$-ray flux for a given cosmic ray energy injection rate. Unfortunately, we cannot simply attribute the radio and $\gamma$-ray fluxes to different cosmic ray populations.

The primary hadronic $\gamma$-ray emission mechanism is neutral pion decay. In the process of producing neutral pions, cosmic ray protons also create charged pions which subsequently decay into secondary electrons and positrons, along with various neutrinos. Hence, the initial interactions leading to hadronic $\gamma$-rays also lead to radiation from secondary leptons and thus radio emission. Therefore, whether the $\gamma$-rays are hadronic or leptonic in origin, there should also be a substantial amount of corresponding radio emission.

To explore all the variables affecting radio emission, we break down the various equations governing the radio synchrotron emissivity in Fig 3. Inputs affecting the total synchrotron flux include: the magnetic field strength ($B$), the average ISM gas density ($n_{\rm ISM}$), the radiation field energy density ($U_{\rm rad}$), the electron-to-proton (number) ratio ($\xi$), the cosmic ray injection spectral index ($\Gamma_{e}$), and the cosmic ray energy injection rate ($\varepsilon_{\rm CR}$).

Both the radio and $\gamma$-ray spectra provide constraints on the cosmic ray spectral index, and the $\gamma$-ray flux, in particular, provides a lower limit on the the cosmic ray energy injection rate as discussed in Section 1, limiting our avenues of altering the radio flux. Additionally, lowering the electron-to-proton (number) ratio has limited effect due to the presence of secondary electrons and positrons. Further observational constraints on the magnetic field strength, ISM density, and radiation field energy density are discussed in Section 3.

2.1 Cosmic Ray Leptons

Non-thermal emission dominates the radio spectrum at GHz frequencies, and so, we concentrate on the synchrotron emissivity,

[TABLE]

where $N_{e}(E_{e})$ is the total leptonic energy spectrum. $P_{\nu}^{\rm Sy}$, the full synchrotron radiative power for a single electron, is given by (Ginzburg, 1969; Blumenthal & Gould, 1970; Rybicki & Lightman, 1979; Schlickeiser, 2002; Dermer & Menon, 2009; Boettcher et al., 2012; Ghisellini, 2013)

[TABLE]

Here, $\psi$ is the pitch angle, and $F(x)=x\int_{x}^{\infty}d\xi K_{5/3}(\xi)$. $K_{5/3}(\xi)$ is a modified Bessel function of the second kind, and $x=\nu/\nu_{\psi}$ where $\nu_{\psi}=3q_{e}B/4\pi m_{e}c\times\gamma_{e}^{2}\sin(\psi)$. Assuming an isotropic pitch angle distribution and adopting the asymptotic approximation by Boettcher et al. (2012), the synchrotron power is given by

[TABLE]

where $\Gamma(4/3)$ is the Gamma function and the critical frequency $\nu_{c}=3q_{e}\gamma_{e}^{2}B/4\pi m_{e}c$.

The energy spectrum for cosmic ray leptons is given by the product of the cosmic ray injection spectrum, $q(E)$, and the cosmic ray confinement time-scale, $t(E)$, see Fig. 4.

[TABLE]

The cosmic ray injection spectrum, $q(E)$, is assumed to be a power-law, proportional to the energy injection rate and inversely proportional to the volume of the injection region, as seen in the previous section.

The cosmic ray lepton spectrum includes both primary electrons, accelerated alongside the primary protons, and secondary electrons and positrons, the decay products of charged pions produced during collision of cosmic ray protons with the surrounding ISM. Thus, the total cosmic ray lepton injection spectrum is given by

[TABLE]

where the injection spectrum of the primary electrons is given by

[TABLE]

Here, $\xi$ is the electron-to-proton (number) ratio, which is usually taken to be $\xi=0.01$ for $E_{e}=E_{p}=1$ GeV. For simplicity, we will assume that the spectral indices for the injection spectra for the cosmic ray protons and electrons are the same, $\Gamma_{p}=\Gamma_{e}$. Then, the injection spectrum for secondary cosmic rays depends directly on the cosmic ray proton energy spectrum such that (e.g., Schlickeiser, 2002; Kelner et al., 2006)

[TABLE]

The cosmic ray confinement time is comprised of the energy-loss time-scale ($t_{\rm loss}(E_{e})$), the diffusion time-scale ($t_{\rm diff}(E_{e})$), and the energy-independent advective time-scale ($t_{\rm adv}$) such that

[TABLE]

Because of the extreme environment of Arp 220, the energy-loss time-scale is significantly shorter than either the diffusion or advective time-scales, see Fig. 4. Further, an advective time-scale competitive with the energy-loss time-scale would only result in the need for a larger cosmic ray energy injection rate for the observed $\gamma$-ray spectrum, increasing the discrepancy between the injection rate derived from the $\gamma$-ray spectrum and the observed supernova rate. Thus, to obtain a conservative lower limit on the cosmic ray injection rate, we will assume that the total lifetime for the cosmic ray electrons is simply

[TABLE]

Based on the above equations, lowering the radio flux requires lowering either the synchrotron power ($P_{\nu}^{\rm Sy}$) or the lepton energy spectrum ($N_{e}(E_{e})$), see Fig. 3. From equation 4, we can see that the only means by which we can decrease the synchrotron power is by decreasing the magnetic field strength.

The situation in regard to the energy spectrum for cosmic ray leptons is more complicated. To decrease the leptonic energy spectrum, we could decrease the energy injection spectrum. However, the observed $\gamma$-ray spectrum provides stringent constraints on the cosmic ray spectral index ($\Gamma_{p}$) and a lower limit on the cosmic ray energy injection rate ($\varepsilon_{\rm CR}$). Thus, lowering the population of secondary electrons and positrons is not possible within the confines of the observed $\gamma$-ray flux. Hence, the primary avenue for lowering the leptonic energy spectrum is via the confinement time-scale, see Fig. 3. As the energy-loss time-scale dominates the total confinement time-scale, see Fig. 4, we will concentrate our efforts on the energy loss rate for cosmic ray electrons.

2.2 Energy Loss Rates

As in previous papers, we include energy losses via ionization, bremsstrahlung, inverse Compton, and synchrotron in our models. The energy loss rates for these mechanisms are directly proportional to the gas density ($n$), the radiation field energy density ($U_{\rm rad}$), and the magnetic field strength ($B$) such that (Ginzburg, 1969; Schlickeiser, 2002; Dermer & Menon, 2009; Schlickeiser & Ruppel, 2010)

[TABLE]

where $\phi_{\rm\text{H\,{i}}}^{s-s}\approx 45$ is the scattering function for strongly shielded gas.

The loss rates for ionization and bremsstrahlung both scale with the gas density (note that $n_{\rm mol}\gg n_{\rm ion}$) and because of the relative energy dependencies, bremsstrahlung losses dominate over ionization losses starting at $E_{e}\gtrsim 0.5$ GeV, see Fig. 5. Similarly, the loss rates for inverse Compton and synchrotron have the same energy dependence, and thus, inverse Compton losses dominate for $U_{\rm rad}>U_{B}$ and synchrotron losses dominate for $U_{\rm rad}<U_{B}$, where $U_{B}=B^{2}/8\pi$.

For ‘low’ densities ($n\sim 200$ cm*-3*) and ‘high’ magnetic or radiation fields ($U>10^{5}$ eV cm*-3*), synchrotron or inverse Compton losses dominate for all relevant energies (see Fig. 5). However, for ‘moderate’ densities ($n\sim 1000$ cm*-3*) and more ‘moderate’ magnetic or radiation fields ($U>1000$ eV cm*-3*), bremsstrahlung losses are dominate at low energies ($E_{e}<1.0$ GeV) and synchrotron or inverse Compton losses dominate at higher energies, see Fig. 5.

In equation (14), $\gamma_{\rm K}$ is the critical Klein-Nishina Lorentz factor which is given by $\gamma_{\rm K}=3\sqrt{5}m_{e}c^{2}/8\pi kT_{\rm rad}$ (Schlickeiser & Ruppel, 2010). For $\gamma\ll\gamma_{\rm K}$, the inverse Compton cross section is well approximated by the Thompson cross section, while for $\gamma\geq\gamma_{\rm K}$ the full Klein-Nishina cross section must be used and the energy loss rate from inverse Compton becomes energy independent. As shown in Fig. 5, this change in the cross section becomes important for either high energies ($E_{e}\gtrsim 10$ TeV), high dust temperatures ($T_{\rm rad}\gtrsim 250$ K), or high optical depths.

3 Observational Constraints

As in previous work (Yoast-Hull et al., 2015, 2017), we model the eastern nucleus with a one-zone model and the western nucleus as two zones: an inner circumnuclear disc (CND) and an outer torus. Here, we will discuss the relevant properties for both nuclei, and we have listed these properties in Table LABEL:tab:obs.

Based on observations of the molecular media and high-resolution radio continuum measurements, we assume radii of $R=85$ pc, $R=65$ pc, and $R=15$ pc for the eastern nucleus, the outer western nucleus, and the CND, respectively (Barcos-Muñoz et al., 2015; Scoville et al., 2017). We assume scale-heights of $h=30$ pc in the eastern and outer western nuclei and $h=10$ pc for the CND (Scoville et al., 2017; Sakamoto et al., 2017).

3.1 Molecular Gas Masses

Millimeter wavelength observations provide probes of the H2 column density ($N_{H_{2}}$) ranging from $10^{25}$ cm*-2* to greater than $10^{26}$ cm*-2* (Wilson et al., 2014; Scoville et al., 2017), assuming a Galactic CO-to-H2 conversion factor of $X_{\rm CO,20}=2$. Modeling and observations suggest that the conversion factor should be lower in ULIRGs, with potential values ranging from $X_{\rm CO,20}\simeq 0.2-1.0$ (Bolatto et al., 2013, and references therein). We require that the molecular gas mass be less than the dynamical mass. So, for a disc of radius $R$ and total height $H$, we calculate the molecular gas mass such that

[TABLE]

For the eastern and outer western nuclei, we adopt moderate column densities of $10^{25}$ cm*-2*, and for the CND, we assume a density of $10^{26}$ cm*-2* based on Scoville et al. (2017). By choosing a moderate conversion factor of $X_{\rm CO,20}^{\rm ULIRG}=0.6$, we settle on molecular gas masses of $\sim 10^{9}$ M*⊙* for each of the nuclei, see Table LABEL:tab:obs, which is consistent with other estimates of the gas mass (e.g., Barcos-Muñoz et al., 2015).

3.2 Radiation Fields

Measurements of the infrared luminosity give us a means to calculate a lower limit on the radiation field energy density ($U_{\rm IR}$) such that

[TABLE]

From the Stefan-Boltzmann law ($U_{\rm IR}=aT_{\rm IR}^{4}$), we can also calculate an effective dust temperature from the infrared luminosity such that

[TABLE]

where $a$ is the radiation constant $a=4\sigma/c$.

For the CND, estimates from millimeter wavelength observations put the infrared luminosity at $L_{\rm IR}\simeq 6\times 10^{11}$ L*⊙* (Wilson et al., 2014; Scoville et al., 2017). This corresponds to an effective dust temperature of 135 K with an energy density of $1.5\times 10^{6}$ eV cm*-3*. Adopting slightly less intense luminosities for the eastern nucleus and the outer western nucleus gives us effective dust temperatures of $\sim 50$ K and energy densities of $\sim 10^{4}$ eV cm*-3*, see Table LABEL:tab:obs.

3.3 Optical Depths

To lower the total radio emission while simultaneously preserving a high energy injection rate, one could increase the radiation field energy density (or the infrared optical depth), increase the average ISM density, or decrease the magnetic field strength. From Zeeman splitting in OH masers in Arp 220, observations confirm that the magnetic field strengths in the Arp 220 nuclei are of mG strength (McBride et al., 2015). Magnetic fields of similar strengths are found even in extremely dense molecular gas clouds in our own Milky Way (e.g., Crutcher, 2012, and references therein). Thus, we cannot decrease the magnetic field strength significantly from previous models (Yoast-Hull et al., 2016).

Similarly, estimates of the dynamical masses for the nuclei are $\sim 2\times 10^{9}~{}\rm M_{\odot}$ which provides a rough upper limit on the molecular gas mass (Scoville et al., 2017). Therefore, we cannot increase the average ISM density significantly above previous models (Yoast-Hull et al., 2017).

Based on these observational constraints, we will be focusing on the optical depth and thus the radiation field energy density $U_{\rm rad}$. For an optically thick disc, the radiation field energy density in the disc is higher than the emergent radiation field energy density, inferred from observations using equation (17). The energy density and dust temperature in the disc are related to the emergent flux and radiation temperature such that

[TABLE]

The effective optical depth ($\tau_{\rm IR}$) used here is the optical depth of the integrated infrared spectrum, which peaks near 60 $\mu$m (González-Alfonso et al., 2004). We will discuss the equivalent optical depth for millimeter wavelengths and compare with observational estimates in Section 5.1.

Using these relations, we find that for the optically thick case ($\tau_{\rm IR}\geq 1$)

[TABLE]

Additionally, from equations (34) & (36), we know that the inverse Compton emissivity ($q_{\rm IC}$) depends on equation (35), $I(u\times\tau_{\rm IR}^{1/4}kT_{\rm IR},E_{\gamma})$, which depends on the electron loss time-scale and thus the electron energy loss rates, see Appendix B for details.

4 Results

4.1 Model Setup

The goal of this study is to find a solution or set of solutions in parameter space via $\chi^{2}$ tests that simultaneously fit the $\gamma$-ray and radio spectra. For the $\gamma$-ray spectra, the most critical parameters are the cosmic ray spectral index and the cosmic ray energy injection rate. Because the $\gamma$-ray observations provide no information about regarding spatial distribution, we test three different scenarios for where and how energy is injected into cosmic rays.

Previously, we began with injection rates of $0.7\times 10^{50}$ erg yr*-1*, $0.7\times 10^{50}$ erg yr*-1*, and $1.3\times 10^{51}$ erg yr*-1* for the eastern nucleus, western nucleus, and CND, respectively (Yoast-Hull et al., 2015; Varenius et al., 2016). Thus, in Scenario A, we begin by assuming similar initial injection rates of $0.5\times 10^{50}$ erg yr*-1*, $0.5\times 10^{50}$ erg yr*-1*, and $10^{50}$ erg yr*-1* and increase the injection rates in equal proportion up to a total injection rate of $4\times 10^{51}$ erg yr*-1*, see Table 2. This model includes primary protons along with primary and secondary electrons, with the assumption that the primary electron-to-proton (number) ratio is $\xi=0.02$ (Yoast-Hull et al., 2013).

In Scenario B, we fix the injection rates in the eastern and outer western nuclei to $10^{50}$ erg yr*-1* in each region to be consistent with the observed supernova rates, see Varenius et al. (2016). In the CND, we test a range of injection rates from $2\times 10^{50}$ to $4\times 10^{51}$ erg yr*-1*, see Table 2, and inject only primary protons, including secondary electrons and positrons but not primary electrons. In Scenario C, we similarly fix the injection rates in the eastern and western nuclei and test a range of rates from $10^{51}-10^{55}$ erg yr*-1* in steps of $10^{0.2}$ yr*-1* in the CND, see Table 2. For Scenario C, we only include primary electrons and no primary protons or secondary electrons and positrons.

In each scenario, we test the models against both the $\gamma$-ray and radio spectra while varying the injection spectral index (2.3, 2.4, 2.5), magnetic field strength ($0.1-100$ mG), and optical depth ($1-10^{4}$). For the radio spectra, we also vary the absorption fraction ($0.2-1.0$) and the ionized gas density ($50-800$ cm*-3*), both of which affect thermal radio emission. For our $\chi^{2}$ tests against the $\gamma$-ray data, we assume values for the spectral index and optical depth are equal across each of the regions (East, West, CND) in order to reduce the computation time, and we assume that the magnetic field strength in the CND is $\sim 7$ times larger than in the East and the West based on the difference in ISM densities.

Regarding radio data, there exist separate spectra for the eastern and western nuclei. Thus, in Scenario A, we assume separate values for spectral index, optical depth, magnetic field, ionized gas density, and absorption fraction in the eastern nucleus versus the western nucleus, and we also assume that the spectral index, optical depth, and absorption fraction are constant across the outer western nucleus and the inner CND. We also assume that the ionized gas density increases by a factor of $\sim 50$ from the outer western nucleus to the inner CND.

In Scenarios B & C, we also assume that the bulk of the $\gamma$-ray emission and the minority of the radio emission comes from the CND. Thus, for the CND, we require that the radio emission be less than 1% of the total emission from the western nucleus. As such, we can test separate spectral indices, magnetic field strengths, and optical depths for the West versus the CND.

4.2 Gamma-Ray Spectra

The best-fitting results from the $\chi^{2}$ tests, performed against the $\gamma$-ray data from Peng et al. (2016), require energy injection rates of $2\times 10^{51}$ erg yr*-1* and $10^{51}$ erg yr*-1* for Scenarios A and B, respectively, see Fig. 6. These results for Scenarios A and B are independent of optical depth (see Fig. 6), as hadronic $\gamma$-rays dominate the emission in these models. Additionally, these models are best fit with spectral indices of $\Gamma_{p}\simeq 2.4$, closely matching the spectral index ($\Gamma_{\gamma}=2.35$) derived from fitting a pure power-law to the $\gamma$-ray spectrum (Peng et al., 2016).

For Scenario C, we assume injection of cosmic ray leptons, and thus the energy injection rate required to fit the $\gamma$-ray data changes with both optical depth and magnetic field strength, see Fig. 7. Overall, preferred models have energy injection rates of $\leq 10^{52}$ erg yr*-1*, and the best-fitting models require both high optical depths and high magnetic field strengths.

Examining Fig. 7, we can see that at low magnetic field strengths, the energy injection rate is nearly independent of optical depth. However, as magnetic field strength increases, the relationship between optical depth and injection rate nears a pure power-law. Further, these models show a clear requirement for an optical depth significantly greater than unity, see Fig. 7.

4.3 Radio Spectra

For the radio spectra, we have collected data from Downes & Solomon (1998); Rodríguez-Rico et al. (2005); Barcos-Muñoz et al. (2015), covering a frequency range from 1.7 to 43 GHz. Here, we limit our analysis to models with only the best-fitting energy injection rates from the $\gamma$-ray analysis. For Scenario A, our best-fitting models for both the eastern and western nuclei are restricted to magnetic field strengths greater than $10$ mG and optical depths greater than $10^{3}$, see Fig. 8. For Scenarios B and C, where we have assumed constant injection rates of $10^{50}$ erg yr*-1* for the East and West, we find solutions with magnetic field strengths as low as a few mG and optical depths ranging from unity to $\sim 10^{3}$.

In the case of the CND, we consider any model producing radio emission that is less than 1% of the total radio emission from the western nucleus to be acceptably negligible in regards to the total radio emission. These models are marked as having a $\chi^{2}$ value of 1 (corresponding to dark purple) in Fig. 8. If we allow for magnetic field strengths below 1 mG, then the minimum requirement for the optical depth for a joint solution is simply unity. If we require that the magnetic field strength be greater than or equal to 1 mG, then the optical depth must be larger than 10 for Scenario B and larger than 100 for Scenario C, see Fig. 8.

In all cases, it is clear that as magnetic field strength increases, optical depth must also increase, see Fig. 8. This is expected because as magnetic field strength increases, the total synchrotron emissivity increases. Thus, to keep the radio emission small, the optical depth must also increase. Above, we also found that for Scenario C, as the energy injection rate increases and assuming high ($\geq 10$ mG) magnetic field strengths, the optical depth must decrease, see Fig. 7. For the $\gamma$-rays, the $\gamma$-ray flux increases as the energy injection rate increases. Thus, to avoid over-predicting the $\gamma$-ray flux, the optical depth must decrease as the energy injection rate increases. Regardless, there is clearly a minimum optical depth necessary to provide a joint solution to the radio and $\gamma$-ray spectra.

5 Discussion

5.1 Observed Optical Depths

As stated above, the best-fitting optical depths derived from our $\chi^{2}$ analysis are effective optical depths over the whole of the infrared range. Thus, a full and complete comparison between our $\tau_{\rm IR}$ and optical depth values derived at specific wavelengths from observational measurements is complex. However, a rough comparison is simple enough if we assume that our effective optical depth is approximately representative of the optical depth at the peak of the infrared spectrum, $\sim 60$ $\mu$m in Arp 220 (González-Alfonso et al., 2004; Rangwala et al., 2011). Assuming optical depth depends on frequency to the 1.6 power in ULIRGs (Casey, 2012), we have

[TABLE]

In Wilson et al. (2014), the authors calculated optical depths of 1.7 and 5.3 for the eastern and western nuclei, respectively. The optical depths were derived from a combination of submillimeter observations with the Atacama Large Millimeter/submillimeter Array (ALMA) at 691 GHz (434 $\mu$m) and at 345 GHz (869 $\mu$m). Converting these optical depths to optical depths in the infrared range, we get $\tau_{60}=40.3$ and $\tau_{60}=125.7$ for the eastern and western nuclei, respectively. These values would suggest that Scenarios B and C are more likely than Scenario A, as optical depths larger than $10^{3}$ were necessary for joint solutions between the radio and $\gamma$-ray data in Scenario A.

5.2 Energy Densities

From our results, it is clear that there is significantly higher energy density in infrared radiation in the nuclei of Arp 220 than in the free-streaming emergent (observed) radiation field. This is consistent with submillimeter and millimeter wavelength estimates of the optical depth at those wavelengths (e.g., Wilson et al., 2014; Scoville et al., 2017). To evaluate the how the energy density in the trapped radiation field compares with the other components of the ISM, we evaluate the energy densities for turbulence, magnetic fields, radiation fields, and cosmic rays below.

The turbulent energy density is given as $U_{\rm turb}=\rho\sigma_{\rm turb}^{2}/2$, where $\rho$ is the ISM gas density ($\rho_{\rm H_{2}}=m_{\rm H_{2}}n_{\rm H_{2}}$). Measurements by Scoville et al. (2017) show the turbulent velocities to be 120 km s*-1* and 250 km s*-1* in the eastern and western nuclei, respectively. Combined with gas densities of $\sim 10^{4}$ and $6.5\times 10^{5}$ cm*-3*, these velocities give us turbulent energy densities of $1.7\times 10^{6}$ and $4.2\times 10^{8}$ eV cm*-3*.

For magnetic field strengths of $\sim 4-20$ mG, the corresponding energy densities ($U_{B}=B^{2}/8\pi$) range from $4\times 10^{5}$ eV cm*-3* to $10^{7}$ eV cm*-3*. Thus, the magnetic field energy density lags slightly behind the turbulent energy density in the eastern and outer western nuclei. However, in the CND, the magnetic field energy density clearly lags behind the turbulent energy density by more than an order of magnitude.

For the radiation field, the total energy density is given by $U_{\rm rad}=\tau_{\rm IR}U_{\rm IR}$. Thus, for the East and the West, the total radiation field energy densities range from $\sim 10^{4}-10^{7}$ eV cm*-3* based on the results of our $\chi^{2}$ analysis in Scenarios B and C. In the CND, the total radiation field energy density must be greater than $\sim 10^{7}$ eV cm*-3* and greater than $\sim 10^{8}$ eV cm*-3* in Scenarios B and C, respectively.

If we take our infrared optical depths approximated from the Wilson et al. (2014) measurements, the corresponding energy densities are $1.1\times 10^{6}$ eV cm*-3* for the eastern / outer western nucleus versus $1.9\times 10^{8}$ eV cm*-3* for the CND. These specific values are within the range of energy densities from our acceptable models and they agree well with the turbulent energy densities.

Finally, for the cosmic rays, the energy density involves an integral over the energy spectrum such that $U_{\rm CR}=\int EN(E)dE$. For Scenarios B and C, the cosmic ray energy densities for the East and West are $1.8\times 10^{3}$ and $3\times 10^{3}$ eV cm*-3*. Because Scenarios B and C have low assumed energy injection rates, along with high energy loss rates for cosmic rays, the corresponding energy densities are low and consistent with our previous study in Yoast-Hull et al. (2015). For the CND, in Scenarios B and C, the total cosmic ray energy densities are $\sim 4\times 10^{4}$ eV cm*-3*. Thus, even with the additional cosmic ray population inferred from the $\gamma$-ray flux, the cosmic ray energy density is well below the other energy densities.

5.3 Further Interactions with Infrared Photons

Having established that the infrared optical depth in the nuclei of Arp 220 is greater than unity, we must assess how this increase in the radiation energy density affects other interaction rates. Specifically, we examine how the optical depth for photon-photon interactions and the source function for photopion production increase with increasing infrared optical depth.

5.3.1 Photon-Photon Interactions

From previous studies, we know that even considering only the emergent infrared energy density, the radiation fields in Arp 220 are large enough for photon-photon ($\gamma-\gamma$) interactions to create significant absorption in the TeV energy $\gamma$-ray spectrum (Torres, 2004; Lacki & Thompson, 2013; Yoast-Hull et al., 2015). The optical depth for $\gamma-\gamma$ interactions depends on the integral over the photon number density over the square of the photon energy and is thus directly proportional to the radiation field energy density and inversely proportional to the temperature to the third power, see Appendix C for details. Hence, $\tau_{\gamma-\gamma}$ depends on the effective infrared optical depth such that $\tau_{\gamma-\gamma}\propto\tau_{\rm IR}^{1/4}$.

Calculating the photon-photon optical depth for various infrared optical depths, we naturally find that as $\tau_{\rm IR}$ increases, the peak optical depth for $\gamma-\gamma$ interactions shifts up in value and down in energy, see Fig. 9. However, we find that for even very large values for the infrared optical depth, $\tau_{\gamma-\gamma}$ does not rise above unity for energies below 100 GeV. The $\gamma$-ray energy spectrum measured in Peng et al. (2016) extends from $0.1-10$ GeV. Thus, while there is a significant increase in the $\gamma$-ray absorption due to photon-photon interactions above 100 GeV, the spectral region covered by the Fermi $\gamma$-ray observations is optically thin.

At present, there are no $\gamma$-ray detections at energies above 10 GeV, only an upper limit at 100 GeV from VERITAS (Fleischhack & VERITAS Collaboration, 2015). However, a future detection by VERITAS or the Cherenkov Telescope Array (CTA) could provide valuable additional constraints on the infrared optical depth and the cosmic ray spectral index.

5.3.2 Proton-Photon Interactions

In an environment with extremely intense radiation fields, pions produced in proton-photon ($p-\gamma$) interactions can be the dominant source of $\gamma$-ray and neutrino production. Here, we compare the source function for pions from $p-\gamma$ processes against that from $p-p$ processes, see Fig. 10. Using the simplified models for photo-hadronic interactions from Hümmer et al. (2010), we plot the source function for photopions for a variety of effective infrared optical depths and for varying maximum cosmic ray energies, see Appendix C for details about the source function.

Examining Fig. 10, we find that even at large optical depths, the pion flux from $p-p$ interactions is orders of magnitude larger than the pion flux from $p-\gamma$ interactions. Additionally, the maximum cosmic ray energy must be greater than $10^{18}$ eV for the cosmic ray protons in order for the photopion source function to even approach the source function of pions from proton-proton interactions.

6 Summary

As established above and in our previous analysis (Yoast-Hull et al., 2017), there is a significant discrepancy between the cosmic ray energy injection rates derived from the radio spectra versus the $\gamma$-ray spectrum in Arp 220, see Fig. 2. Observational constraints on the average ISM density and the magnetic field strength do not allow for adjustment of the radio spectrum through decreases in the synchrotron power spectrum or through increases in the bremsstrahlung energy loss rates, see Fig. 3.

However, observational diagnostics show that the nuclei of Arp 220 are optically thick over the extended range of the infrared through to millimeter wavelengths. It follows then that the intrinsic infrared energy densities inside the nuclei is significantly larger than the emergent energy densities (e.g., Wilson et al., 2014; Scoville et al., 2017). By testing different energy injection scenarios with varying values for an effective infrared optical depth, we have shown above that an optical depth greater than unity can be combined with a large cosmic ray energy injection rate to simultaneously fit the radio and $\gamma$-ray spectra, see Fig. 6 – 8.

Infrared optical depths calculated from submillimeter wavelength measurements (Wilson et al., 2014) are within the range of effective optical depths considered acceptable from our $\chi^{2}$ tests. Further, using these optical depths to calculate the total radiation field energy density, we find equipartition exists in both nuclei between the turbulent energy densities and the radiation field energy densities.

These results are, of course, dependent on several assumptions about the total dense gas mass, dust-to-gas ratio, and additional properties of the dust (e.g., grain size) in Arp 220. Additionally, as the image of $\gamma$-ray emission from Arp 220 is consistent with that of a point source, we lack any spatial information regarding the $\gamma$-ray flux or the composition of the additional cosmic ray population needed to produce such a $\gamma$-ray flux. Because of these factors and the complexities related to radiative transfer, we cannot use our results to pinpoint an exact value for the optical depths.

Regardless of the limitations of these models, we have clearly established that the infrared energy density in the nuclei is significantly larger than that in the emergent radiation field. Thus, models for cosmic ray interactions must properly account for optically thick radiation in calculating radio and $\gamma$-ray emission. This may have major implications for methods for deriving SFRs from radio emission and calculations of $\gamma$-ray and neutrino emission at high redshift, provided those regions or systems host very compact regions of star formation (i.e., large star-formation surface densities).

Acknowledgements

We thank Jay Gallagher, Todd Thompson, and Almog Yalinewich for useful discussions on Arp 220 and on greybodies and radiative transfer. We also thank Mubdi Rahman for discussions on statistics. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611. We acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC). This research was undertaken, in part, thanks to funding from the Canada Research Chairs program.

Appendix A Nearby Galaxies

A.1 Regression Analysis

The correlation between the diffuse FIR and radio emission in star-forming galaxies is well-established. Recent detections of star-forming galaxies in $\gamma$-rays have shown a second correlation between the diffuse FIR and $\gamma$-ray emission, see Fig. 11. Additionally, we show here that there is a third correlation between the diffuse $\gamma$-ray and radio emission in star-forming galaxies.

For the galaxies listed in Tables LABEL:tab:data1 & LABEL:tab:data2, we perform a Bayesian linear regression analysis on the radio and FIR luminosities. Assuming normal priors, the mean slope and intercept are $m_{A}=1.148\pm 0.047$ and $b_{A}=-7.324\pm 0.501$, respectively, with a standard deviation of $\sigma_{A}=0.258$. The slope derived here is slightly different from the best-fitting slope in Yun et al. (2001); this is expected due to the difference in samples. Here, we are concentrating on a small sample ($N\sim 60$) of nearby IR-bright galaxies, while general studies of the FIR–radio correlation have sample sizes of $\gtrsim 1800$.

For the galaxies listed in only Table LABEL:tab:data1, we also perform a Bayesian linear regression analysis on the FIR and $\gamma$-ray luminosities and the radio and $\gamma$-ray luminosities. Again assuming normal priors, we obtain mean slopes and intercepts of $m_{B}=1.235\pm 0.093$ and $b_{B}=-6.125\pm 0.920$ with a standard deviation of $\sigma_{B}=0.297$ for the FIR and $\gamma$-ray luminosities and $m_{C}=1.112\pm 0.223$ and $b_{C}=1.739\pm 1.061$ with a standard deviation of $\sigma_{C}=0.481$ for the radio and $\gamma$-ray luminosities.

In Fig. 11, we show FIR, $\gamma$-ray, and radio luminosities for 60 nearby galaxies. Galaxies with $\gamma$-ray detections are shown in filled symbols and those with only $\gamma$-ray upper limits are shown in open symbols. Pure star-forming galaxies are denoted with circular symbols and active galaxies are denoted with triangular symbols.

For each analysis, we ran 10,000 samples with a burn sample of 1000. Modeling for the radio and FIR luminosities converged. However, there were $\sim 2000$ and $\sim 200$ divergences for the FIR and $\gamma$-ray luminosities and the radio and $\gamma$-ray luminosities, respectively. With an extremely small sample size of 8 galaxies and probable outliers (e.g., Circinus), there is little to be done to create convergence. Thus, we see a significantly larger spread in the posterior models for the bottom two figures in Fig. 11.

A.2 Multi-Wavelength Observations

The total far-infrared fluxes for each galaxy are calculated such that

[TABLE]

where $F_{\rm FIR}$ is given in units of erg cm*-2* s*-1* (Sanders & Mirabel, 1996, and references therein). Infrared data are taken from the IRAS Revised Bright Galaxy Sample (Sanders et al., 2003) and radio data are taken from the NRAO VLA Sky Survey (NVSS; Condon et al., 1998). $\gamma$-ray data in Table LABEL:tab:data1 are taken from the Third Fermi Source Catalog (Acero et al., 2015) and from Lopez et al. (2018) and Ackermann et al. (2016) for the SMC and the LMC, respectively. Upper-limits for $\gamma$-ray fluxes in Table LABEL:tab:data2 are taken from Rojas-Bravo & Araya (2016). Exceptions:

M31 (NGC 0224): The infrared fluxes are taken from Xu & Helou (1996), and the radio flux is taken from Dennison et al. (1975).

NGC 1144: The infrared fluxes are taken from the IRAS Bright Galaxy Sample (Soifer et al., 1989).

NGC 3627: The radio flux used is the sum of both NVSS components associated with the galaxy, following Condon et al. (2002).

NGC 4945: The radio flux is taken from Elmouttie et al. (1997).

NGC 5055: The radio flux is taken from Condon et al. (2002) and is inconsistent with the flux listed in the NVSS catalog.

M51 (NGC 5194): The radio flux is taken from Condon et al. (2002) and references therein.

Circinus: The infrared fluxes are taken from the IRAS Point Source Catalog (Helou & Walker, 1988), and the radio flux is taken from Elmouttie et al. (1998).

NGC 7331: The radio flux used is the sum of both NVSS components associated with the galaxy, following Condon et al. (2002).

Appendix B Leptonic Emissivities

B.1 Bremsstrahlung

The emissivity for bremsstrahlung is given as

[TABLE]

where $E_{\rm min}=\sqrt{E_{\gamma}(2m_{e}c^{2}+E_{\gamma})}$ and the differential cross section is given by (Ginzburg, 1969; Blumenthal & Gould, 1970; Stecker, 1971; Schlickeiser, 2002)

[TABLE]

The scattering functions $\phi_{1}$ and $\phi_{2}$ are related such that $\phi_{2}=-2/3\phi_{1}$ and are equivalent to to $Z^{2}\phi_{u}$ where $\phi_{u}$ is given by (Ginzburg, 1969; Blumenthal & Gould, 1970; Schlickeiser, 2002)

[TABLE]

B.2 Photon Number Density

Throughout the paper, we assume a greybody distribution for optical depths of unity or less. The photon number density for such a distribution is given by (Rybicki & Lightman, 1979)

[TABLE]

where $\beta=1.6$ and $\lambda_{0}=200$ µm, parameters which are adopted from Casey (2012). The number density is normalized such that (e.g., Ginzburg, 1969; Rybicki & Lightman, 1979; Ghisellini, 2013)

[TABLE]

where $\Gamma(x)$ is the gamma function and $\zeta(x)$ is the Riemann zeta function. Solving for the dilution factor, $C_{\text{dil}}$, we have

[TABLE]

and substituting back into the number density gives us,

[TABLE]

In the optically thin case, $\beta=1.6$ and thus $\Gamma(5.6)\zeta(5.6)\approx 63$. In the optically thick case, $\beta=0$ and thus $\Gamma(4)\zeta(4)=\pi^{2}/15\approx 6.5$.

B.3 Inverse Compton Emissivity

The emissivity for inverse Compton emission is given as (Blumenthal & Gould, 1970; Schlickeiser, 2002)

[TABLE]

with units of GeV*-1* cm*-3* s*-1*, where the minimum cosmic ray electron energy is given by

[TABLE]

Here, $E_{\gamma}$ is the energy of the resulting $\gamma$-ray and $\epsilon$ is the energy of the incident photon. The function $F(q,\Gamma)$ is part of the Klein-Nishina cross section and is given by (Ginzburg, 1969; Blumenthal & Gould, 1970; Stecker, 1971; Rybicki & Lightman, 1979; Schlickeiser, 2002; Dermer & Menon, 2009; Boettcher et al., 2012; Ghisellini, 2013)

[TABLE]

where

[TABLE]

Here, $\Gamma$ is a product of the photon and electron energies and is not to be confused with the gamma function, $\Gamma(x)$, or the cosmic ray spectral index, $\Gamma_{p}$.

To better understand the effects of the dust temperature ($T_{\rm rad}$) and to simplify our expression for the inverse Compton emissivity, we start by designating the second integral as $I(\epsilon,E_{\gamma})$ such that

[TABLE]

where

[TABLE]

Then, substituting our expression for the greybody distribution from above into our expression for emissivity, we have

[TABLE]

where $K(\beta)=\Gamma(4+\beta)\zeta(4+\beta)$. Changing variables from $\epsilon$ to $u=\epsilon/kT_{\rm rad}$ gives

[TABLE]

Thus, the inverse Compton emissivity decreases as the square of the dust temperature and also changes as $I(\epsilon,E_{\gamma})$ changes. As $I(\epsilon,E_{\gamma})$ depends on $N_{e}(E_{e})=q_{e}(E_{e})t_{e}(E_{e})$, the expression will change as the electron lifetime changes with dust temperature.

Appendix C Photon Interactions

For photon-photon ($\gamma-\gamma$) interactions, the cross section is given by

[TABLE]

where $\beta_{\rm cm}=(1-\gamma_{\rm cm}^{-2})^{1/2}$ and $\gamma_{\rm cm}$ is the center-of-momentum frame Lorentz factor of the produced electron and positron (Dermer & Menon, 2009; Boettcher et al., 2012). Then, the optical depth for $\gamma-\gamma$ interactions is

[TABLE]

where $s_{0}=\epsilon\epsilon_{1}$ is the maximum center-of-momentum energy for head-on collisions. Dermer & Menon (2009) define $d\tau_{\gamma-\gamma}/dx$ as the absorption probability per unit path length, so the approximate optical depth is merely $\tau_{\gamma-\gamma}\simeq l_{\rm path}\times d\tau_{\gamma-\gamma}/dx$. For an inclined ellipse, we take the path length through the disc to be $l_{\rm path}=\pi/2\times R~{}\sec(90^{\circ}-i)$.

If we substitute into equation 38 with our expression for the greybody distribution from Section B.2 and change variables from $\epsilon$ to $u=\epsilon/kT_{\rm rad}$, we have

[TABLE]

Thus, for $\gamma-\gamma$ interactions, the optical depth is directly proportional to $U_{\rm rad}/(kT_{\rm rad})^{3}$.

For proton-photon ($p-\gamma$) interactions, we utilize the parameterizations presented in Hümmer et al. (2010), which is in turn based on the SOFIA code (Mücke et al., 2000). Their parameterizations cover photopion production through both direct and multi-photon production and via the $\Delta$-resonance and higher order $\Delta$- and $N$- resonances. The parameterized source functions (equation 28) take the form

[TABLE]

where $y=E_{p}\epsilon/m_{p}$ and $\chi,M$ are parameterization constants given in Tables 4-6 in Hümmer et al. (2010). The functions $f(y)$ are further parameterizations given in their equations 31-35 and 38-39. The parameterization functions are given for each of the types of production (resonances, direct, multi-pion) and the parameterization constants are further broken down for each type of pion ($\pi^{+},\pi^{-},\pi^{0}$).