A Herschel/PACS Far Infrared Line Emission Survey of Local Luminous Infrared Galaxies

TL;DR

This study analyzes far-infrared line emissions in about 240 local luminous infrared galaxies, revealing how line deficits relate to galaxy properties and identifying a transition point in radiation field intensity linked to starburst activity.

Contribution

It provides a comprehensive analysis of FIR line deficits and PDR conditions in LIRGs, establishing a relation between radiation field intensity, gas density, and infrared surface brightness.

Findings

Line deficits increase with FIR color and surface density.

Fractional PDR contribution to [CII] increases in warmer LIRGs.

G0/nH ratio correlates with infrared surface brightness, with a critical transition at ~5x10^10 Lsun/kpc^2.

Abstract

We present an analysis of [OI]63, [OIII]88, [NII]122 and [CII]158 far-infrared (FIR) fine-structure line observations obtained with Herschel/PACS, for ~240 local luminous infrared galaxies (LIRGs) in the Great Observatories All-sky LIRG Survey (GOALS). We find pronounced declines -deficits- of line-to-FIR-continuum emission for [NII]122, [OI]63 and [CII]158 as a function of FIR color and infrared luminosity surface density, $\Sigma_{\rm IR}$. The median electron density of the ionized gas in LIRGs, based on the [NII]122/[NII]205 ratio, is $n_{\rm e}$ = 41 cm$^{-3}$. We find that the dispersion in the [CII]158 deficit of LIRGs is attributed to a varying fractional contribution of photo-dissociation-regions (PDRs) to the observed [CII]158 emission, f([CII]PDR) = [CII]PDR/[CII], which increases from ~60% to ~95% in the warmest LIRGs. The [OI]63/[CII]158PDR ratio is tightly correlated with the PDR gas kinetic temperature in sources where [OI]63 is not optically-thick or self-absorbed. For each galaxy, we derive the average PDR hydrogen density, $n_{\rm H}$, and intensity of the interstellar radiation field, in units of G$_0$, and find G$_0$/$n_{\rm H}$ ratios ~0.1-50 cm$^3$, with ULIRGs populating the upper end of the distribution. There is a relation between G$_0$/$n_{\rm H}$ and $\Sigma_{\rm IR}$, showing a critical break at $\Sigma_{\rm IR}^{\star}$ ~ 5 x 10$^{10}$ Lsun/kpc$^2$. Below $\Sigma_{\rm IR}^{\star}$, G$_0$/$n_{\rm H}$ remains constant, ~0.32 cm$^3$, and variations in $\Sigma_{\rm IR}$ are driven by the number density of star-forming regions within a galaxy, with no change in their PDR properties. Above $\Sigma_{\rm IR}^{\star}$, G$_0$/$n_{\rm H}$ increases rapidly with $\Sigma_{\rm IR}$, signaling a departure from the typical PDR conditions found in normal star-forming galaxies towards more intense/harder radiation fields and compact geometries typical of starbursting sources.

Tables

Table 1. Table 1 Table Content of Emission Line and Continuum Data

Column #	Quantity	Units
(1)	Id.	…
(2)	Galaxy name	…
(3,4)	Raster (x,y)	[pixel]
(5,6)	Central Spaxel (x,y)	[pixel]
(7,8)	R.A. and Dec of (4,5)	…
(9)	Dist to [O i]63  cont.	[″]
(10)	Central line flux	$\times$ 10-17 [W m-2]
(11)	Central cont. flux dens.	[Jy]
(12)	Companion Galaxy?	…
(13,14,15)	R-C, R-L, C-L	…
(16)	Best line flux	$\times$ 10-17 [W m-2]
(17)	Best cont. flux dens.	[Jy]
(18)	Best measurement	…
(19)	AOR	…
(20)	Program	…

Table 2. Table 3 Best Parameters From The Line Deficit vs. S 63 / S 158 subscript 𝑆 63 subscript 𝑆 158 \hbox{$S_{\rm 63}$}/\hbox{$S_{\rm 158}$} Fits

Line	${ϵ}_0$	$\delta$	1$\sigma$ dispersion
	($\times$ ${10}^{-3}$)	($e$-fold)	($\times$ ${10}^{-3}$)
[O i]63	4.37 $\pm$ 0.93	1.15 $\pm$ 0.30	1.32
[O iii]88	0.94 $\pm$ 0.17	N/A	0.94
[N ii]122	1.27 $\pm$ 0.16	0.67 $\pm$ 0.08	0.11
[C ii]158	14.0 $\pm$ 0.9	0.68 $\pm$ 0.04	1.69
[N ii]205	1.26 $\pm$ 0.16	0.50 $\pm$ 0.04	0.07

Table 3. Table 4 Best Parameters From The Line Deficit vs. Σ IR subscript Σ IR \Sigma_{\rm IR} Fits

Line	${\alpha }_0$	${\alpha }_1$	${\alpha }_2$	1$\sigma$ dispersion
				[dex]
[O i]63	–15.83	2.711	–0.1412	0.24
[O iii]88	–8.78	1.272	–0.0721	0.32
[N ii]122	–7.66	1.023	–0.0613	0.22
[C ii]158	–11.47	2.043	–0.1145	0.17
[N ii]205	–8.24	1.165	–0.0719	0.20

Equations

$$\begin{array}{r}{L}_{\mathrm{line}}/L_{\mathrm{FIR}}={ϵ}_0{e}^{-(S_{63}/S_{158})/\delta }\end{array}$$

$$\begin{array}{r}{T}_{\mathrm{dust}}=20.24+14.54(S_{63}/S_{158})-3.75(S_{63}/S_{158}{)}^2\\ +0.46(S_{63}/S_{158}{)}^3\end{array}$$

$$\begin{array}{r}log(L_{\mathrm{line}}/L_{\mathrm{IR}})={\alpha }_0+{\alpha }_{1}log{\Sigma }_{\mathrm{IR}}+{\alpha }_2(log{\Sigma }_{\mathrm{IR}}{)}^2\end{array}$$

$$\begin{array}{r}f([C\text{sc}ii]{}_{158}^{\mathrm{PDR}})=[C\text{sc}ii]{}_{158}^{\mathrm{PDR}}/[C\text{sc}ii]{}_{158}=0.82(\pm 0.01)\\ +0.41(\pm 0.04)log(S_{63}/S_{158})\end{array}$$

$$\begin{array}{rl}log([O\text{sc}i]{}_{63}/[C\text{sc}ii]{}_{158}^{\mathrm{PDR}})=& -0.16+0.71log{S}_{63}/S_{158}\\ & +0.09(log{S}_{63}/S_{158}{)}^2\end{array}$$

$$\begin{array}{r}log(S_{63}/S_{158})=0.22+1.34logx-0.25(logx{)}^2\end{array}$$

$$L_{\mathrm{FIR}}\approx 2G\pi (D/2{)}^2{\varphi }_{\mathrm{A}}$$

Full text

1. Introduction

One of the most fundamental processes studied in virtually any field of physics is the dissipation of energy (Thomson, 1874). In particular, in extra-galactic astrophysics, investigating how interstellar gas cools down is crucial to our understanding of galaxy formation and evolution, since gravity is only able to collapse structures when they are sufficiently cold.

Thirty years ago, data obtained with the Kuiper Airborne Observatory (KAO) revealed that the far-infrared (FIR) spectra of nearby galaxies were populated with some of the most intense emission lines observed across the electromagnetic spectrum, indicating that they are very efficient cooling channels of the interstellar medium (ISM) (Watson et al., 1984; Stacey et al., 1991; Lord et al., 1996). A decade later, the Infrared Space Observatory (ISO; Kessler et al., 1996) increased the number of galaxies with FIR fine-structure line detections to the dozens (Malhotra et al., 1997, 2001; Luhman et al., 1998, 2003). But twenty more years needed to pass until the Herschel Space Observatory (Herschel hereafter; Pilbratt et al., 2010) re-opened a new window into the FIR universe providing spectroscopic data for significantly larger samples of nearby, intermediate and high-redshift galaxies (e.g., Graciá-Carpio et al., 2011; Díaz-Santos et al., 2013, 2014; Rigopoulou et al., 2014; Magdis et al., 2014; Brisbin et al., 2015; Cormier et al., 2015; Rosenberg et al., 2015; Malhotra et al., 2017).

The most important fine-structure lines emitted by atomic species in the $\sim$ 50–200 $\mu$m wavelength range are: [N iii]57.3$\,\mu$m  ([N iii]57), [O i]63.2$\,\mu$m  ([O i]63), [O iii]88.4$\,\mu$m  ([O iii]88), [N ii]121.9$\,\mu$m  ([N ii]122), [O i]145.5$\,\mu$m  ([O i]145), [C ii]157.7$\,\mu$m  ([C ii]158) and [N ii]205.2$\,\mu$m  ([N ii]205). Each of them originates from a different phase of the ISM, with C+ and O emission mostly arising from the neutral and molecular medium within dense photo-dissociation regions (PDRs) surrounding newly formed massive stars (Tielens & Hollenbach, 1985; Hollenbach & Tielens, 1997, and references therein), and N*+,++* and O*++* emission being produced by warm ionized gas, both in the diffuse medium and dense (H ii) regions.

Of particular importance among star-forming galaxies are the so-called luminous IR galaxies (LIRGs; $L_{\rm IR}$ = 1011-12 $L_{\odot}$). LIRGs cover the entire evolutionary merger sequence, ranging from isolated galaxies, to early interacting systems, to advanced mergers. They exhibit enhanced star formation rates (SFRs) and specific SFRs (SSFRs = SFR/$M_{\star}$), consequence of the funneling of large amounts of gas and dust towards their nuclei due to the loss of angular momentum during the dynamical interaction. And while the presence of active galactic nuclei (AGN) in LIRGs is frequent (Petric et al., 2011; Alonso-Herrero et al., 2012), their contribution to the bolometric luminosity of the hosts is still very limited in comparison to ultra-luminous IR galaxies (ULIRGs; $L_{\rm IR}$ $\geq$ 1012 $L_{\odot}$ Veilleux et al., 2009). Therefore, nearby LIRGs are a key galaxy population bridging the gap between normal, Milky-Way (MW) type star-forming galaxies and the most extreme, AGN-dominated (quasar-like) systems (Sanders & Mirabel, 1996). This diversity is also reflected on the fact that they cover the entire transition between main-sequence (MS) galaxies and starbursts (Díaz-Santos et al., 2013).

At high redshfit, LIRGs dominate the obscured star formation activity between z $\sim$ 1–3 (e.g., Berta et al., 2011; Murphy et al., 2011; Magnelli et al., 2011), and a number of works have already shown that local LIRGs (and not ULIRGs) are probably the closest local analogs of this high-z, IR-bright galaxy population (Desai et al., 2007; Pope et al., 2008; Menéndez-Delmestre et al., 2009; Díaz-Santos et al., 2010b; Stacey et al., 2010). Thus, a comprehensive study of the physical properties of low-redshift LIRGs, and specifically of their ISM, is critical for our understanding of the evolution of galaxies and AGN across cosmic time.

To this end, we have performed a systematic study of the most important FIR cooling lines of the ISM in a complete, flux-limited sample of nearby LIRGs, the Great Observatories All-sky LIRG Survey (GOALS; Armus et al., 2009), using Herschel and its Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al., 2010) as well as the Fourier-transform spectrometer (FTS) of the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al., 2010). We combine these observations with mid-IR (MIR) spectroscopy previously obtained with the Infrared Spectrograph (IRS; Houck et al., 2004) on board the Sptizer Space Telescope (Spitzer hereafter; Werner et al., 2004) to provide a panchromatic view of the heating and cooling of the ISM in LIRGs across a wide range of integrated properties such as IR luminosity, compactness, dust temperature, AGN activity, and merger stage.

The paper is organized as follows: In §2 we present the LIRG sample and the new Herschel spectroscopy, as well as the ancillary observations used in this work. In §3 we describe the processing and analysis of the data. The basic results are presented in §4. The ISM properties derived for the LIRG sample are discussed throughout §5 in relation to specific emission line ratios. In particular, in §5.1.1 and §5.1.3 we describe how the fractional contribution of PDRs to the [C ii]158  emission shapes the observed trend between the [C ii]158  deficit and the FIR color of LIRGs. Section §5.1.2 presents the electron densities found for the ionized gas phase of the ISM derived from the nitrogen lines and discusses the implications. In §5.2.1 we present a link between the [O i]63/[C ii]${}^{\rm PDR}_{158}$  and the kinetic temperature of the PDR gas. We explore PDR covering factors and metallicity variations in §5.2.2, and confront FIR emission line ratios involving oxygen and nitrogen emission lines to photo-ionization models of H ii  regions in §5.2.3. We derive the average PDR properties for each galaxy in the sample in §5.3 and show the existence of a critical break in the PDR conditions as a function of luminosity surface density in §5.4. We discuss the physical implications of this results in §5.5. The summary of the results and conclusions are given in §6.

2. Sample and Observations

2.1. The GOALS Sample

The Great Observatories All-sky LIRG Survey (GOALS; Armus et al. 2009) encompasses the complete sample of 202 LIRG and ULIRG systems contained in the IRAS Revised Bright Galaxy Sample (RBGS; Sanders et al., 2003), which in turn is also a complete, flux-limited sample of 629 galaxies with IRAS $S_{60\,\mu m}\,>\,5.24$ Jy and Galactic latitudes $|b|\,>\,5\,^{\circ}$. There are 180 LIRGs and 22 ULIRGs in GOALS and their median redshift is z = 0.0215 (or $\sim\,95.2$ Mpc), with the closest galaxy being at z = 0.0030 (15.9 Mpc; NGC 2146) and the farthest at z = 0.0918 (400 Mpc; IRAS 07251-0248). To date, there are many published works that have already exploited the science content of multi-wavelength data obtained mostly from space-born facilities including GALEX UV (Howell et al., 2010), HST optical and near-IR imaging (Haan et al., 2011; Kim et al., 2013), along with Chandra X-ray (Iwasawa et al., 2011), as well as Spitzer/IRS (Díaz-Santos et al., 2010b; Díaz-Santos et al., 2011; Petric et al., 2011; Stierwalt et al., 2013, 2014; Inami et al., 2013) and Herschel/PACS and SPIRE spectroscopy (Díaz-Santos et al., 2013, 2014; Zhao et al., 2013; Lu et al., 2014, 2015; Zhao et al., 2016b). Moreover, a number of ground-based observatories such as VLA, CARMA and ALMA have also been used to observe the GOALS sample (e.g., Murphy et al., 2013a; Murphy, 2013b; Xu et al., 2014, 2015; Zhao et al., 2016a, among others).

The RBGS, and therefore the GOALS sample, were defined based on IRAS observations. However, the higher angular resolution achieved by Spitzer allowed us to spatially disentangle galaxies within the same LIRG system into separate components. From the more than 290 individual galaxies in GOALS, not all were observed by Herschel. In systems with two or more galactic nuclei, minor companions with MIPS 24$\,\mu$m flux density ratios smaller than 1:5 with respect to the brightest galaxy were not targeted due to their small contribution to the total IR luminosity of the system. The method used to calculate the IR luminosities of individual galaxies ($L_{\rm IR\,[8-1000\,\mu m]}$, as defined in Sanders & Mirabel 1996)111$F_{\rm IR\,[8-1000\,\mu m]}\,=\,1.8\,\times\,10^{-14}\,(13.48\,S_{\rm 12\mu m}\,+5.16\,S_{\rm 25\mu m}\,+2.58\,S_{\rm 60\mu m}\,+\,S_{\rm 100\mu m})$ [W m*-2*], with $S_{\nu}$ in [Jy]. $L_{\rm IR}$ = $4\pi\,D_{\rm L}^{2}\,F_{\rm IR}$. The luminosity distances, $D_{\rm L}$, are taken from Armus et al. (2009). is described in Díaz-Santos et al. (2013) (see end of section 3.2 for further details on the scaling of the IR luminosities used in this work).

2.2. Herschel/PACS Observations

We have obtained FIR spectroscopic observations for 200 LIRG systems from GOALS using the Integral Field Spectrometer (IFS) of the PACS instrument on board Herschel (IRASF08339+6517 and IRASF09111–1007 were not observed). Since some targets contain multiple components, there are 241 individual galaxies with available spectra in at least one emission line. PACS/IFS range spectroscopy of the [O i]63, [O iii]88, [N ii]122  and [C ii]158  fine-structure emission lines was obtained for 239, 161, 75 and 239 individual sources, respectively. Most of the data were collected as part of our OT1 and OT2 programs (OT1_larmus_1, OT2_larmus_1; P.I.: L. Armus) accounting more than 200 hours of observing time in total. Additional observations that are publicly available in the Herschel archive were included from various projects. The main programs from where these complementary data were gathered are: KPGT_esturm_1 (P.I.: E. Sturm), KPOT_pvanderw_1 (P.I.: P. van der Werf) and OT1_dweedman_1 (P.I.: D. Weedman).

The IFS on PACS is able to perform simultaneous spectroscopy in the $51-73$ or $70-105\,\mu$m (3rd and 2nd orders, respectively; “blue” camera) and the $102-210\,\mu$m (1st order; “red” camera) ranges. The IFU is composed of a 5 $\times$ 5 array of individual detectors (spaxels) each with a field of view (FoV) of $\sim\,$9.4″  on a side, for a total of 47″ $\times$ 47″. The physical size of the PACS FoV at the median distance of our LIRG sample is $\sim$ 20 kpc on a side. The number of spectral elements in each pixel is 16, which are rearranged together via an image slicer over two 16 $\times$ 25 Ge:Ga detector arrays (blue and red cameras). The spectral range selected for the observations was scanned several times, increasing the spectral resolution up to at least Nyquist sampling.

While we requested line maps for some LIRGs of the sample (from two to a few raster positions depending on the target), pointed (one single raster) chop-nod observations were taken for the majority of galaxies. For those galaxies with maps, only one raster position was used to obtain the line fluxes presented in this work. For a more detailed discussion on how the observations were set up, we refer the reader to Díaz-Santos et al. (2013).

2.3. Herschel/SPIRE Observations

In addition to the PACS/IFS spectra, we obtained observations of the [N ii]205  emission line using the SPIRE Fourier Transform Spectrometer (FTS) for 121 galaxies in the GOALS sample (Lu et al., 2017, OT1_nlu_1; P.I.: N. Lu;). These observations were part of a broader project whose primarily aim is to study the dense and warm molecular gas in LIRGs (see Lu et al., 2014, 2015). Details about the data processing as well as the results concerning the [N ii]205  line emission are presented in Zhao et al. (2013, 2016b). See section 3.2 for further details about how these data are used.

2.4. Spitzer/IRS Observations

As part of the Spitzer GOALS legacy program, all galaxies observed with Herschel/PACS have available Spitzer/IRS low resolution, R $\sim\,60-120$ (SL module: $5.2-14.5\,\mu$m, and LL module: $14-38\,\mu$m), as well as medium resolution, R $\sim\,600$ (SH module: $9.9-19.6\,\mu$m, and LH module: $18.7-37.2\,\mu$m), slit spectroscopy. The IRS spectra were extracted with the Spitzer IRS Custom Extraction (SPICE) software222http://irsa.ipac.caltech.edu/data/SPITZER/docs/dataanalysistools/tools/spice/, using the standard extraction aperture and a point source calibration mode. The projected angular sizes of the apertures on the sky are 3.7″ $\times$ 12″  at 9.8$\,\mu$m in SL, 10.6″ $\times$ 35″  at 26$\,\mu$m in LL, 4.7″ $\times$ 15.5″  at 14.8$\,\mu$m in SH, and 11.1″ $\times$ 36.6″  at 28$\,\mu$m in LH. Thus, the area covered by the SL and LL apertures is approximately equivalent (within a factor of $\sim\,2$) to that of an individual spaxel and a 3 $\times$ 3 spaxel box of the IFS in PACS, respectively. Likewise, the SH and LH apertures are approximately equivalent to one spaxel and a 3 $\times$ 3 aperture box in PACS, respectively. Aside from the line fluxes, which were presented in Stierwalt et al. (2014), other observables derived from the IRS data that we use in this work are: the strength of the 9.7$\,\mu$m silicate feature, $S_{\rm 9.7\mu m}$333The silicate strength is defined as $S_{\rm 9.7\mu m}$ = ln ($f_{\rm peak}/f_{\rm cont}$), where $f_{\rm peak}$ is the flux density at the peak absorption (or emission) close to 9.7 $\mu$m, and $f_{\rm cont}$ is the flux density of the continuum emission measured outside of the feature, interpolated at the wavelength of the peak. Thus, negative values indicate the feature appears in absorption while positive values indicate it is emission., and the equivalent width (EW) of the 6.2$\,\mu$m PAH, both of which were presented in Stierwalt et al. (2013). We refer the reader to these works for further details about the reduction, extraction and calibration of the IRS spectra as well as for the main results derived from the analysis.

3. Herschel/PACS Data Reduction and Analysis

3.1. Data Processing

The Herschel Interactive Processing Environment (HIPE; v13.0) application was used to retrieve the raw data from the Herschel Science Archive (HSA444http://herschel.esac.esa.int/Science_Archive.shtml) as well as to process them. We used the script for “LineScan” observations (also valid for “range” mode) to reduce our spectra. We processed the data from level 0 up to level 2 using the following steps: Flag and reject saturated data, perform initial calibrations, flag and reject “glitches”, compute the differential signal of each on-off pair of data-points for each chopper cycle, divide by the relative spectral response function, divide by the response, convert frames to PACS cubes, and correct for flat-fielding. Next, for each camera (red or blue), HIPE builds the wavelength grid, for which we chose a final rebinning with an oversample = 2, and an upsample = 1 that corresponds to a Nyquist sampling. The spectral resolution achieved for each line was derived directly from the data and is $\sim$ 82 km s*-1*  for [O i]63, $\sim$ 120 km s*-1*  for [O iii]88, $\sim$ 287 km s*-1*  for [N ii]122, and $\sim$ 235 km s*-1*  for [C ii]158. The final processing steps were: flag and reject remaining outliers, rebin all selected cubes on consistent wavelength grids and, finally, average the nod-A and nod-B rebinned cubes (all cubes at the same raster position are averaged). This is the final science-grade product currently possible for single raster observations. From this point on, the analysis of the spectra was performed using in-house developed idl routines.

3.2. Data Analysis

To obtain the line flux of a source we use an iterative procedure described in detail in Díaz-Santos et al. (2013). The only difference with respect to the method explained in that work is that, due to the highly variable profile of the lines and underlying continuum, for some sources we had to slightly modify the spectral range over which the final continuum-subtracted spectrum is integrated. This was done on a case by case basis over the entire galaxy sample and for each line and spaxel, to ensure that the selected range was correct. The uncertainty associated with the line flux is calculated as the standard deviation of the final fitted underlying continuum integrated over the same wavelength range as the line. That is, uncertainties for all quantities used across this analysis are based on the individual spectrum of each line, therefore reflecting the uncertainties associated with –and measured directly from– the data. Absolute photometric uncertainties, which can be as high as $\sim\,$11–12%555http://herschel.esac.esa.int/Docs/PACS/html/ch04s10.html, are not taken into account except for the analysis performed in section 5.3 (the version of the calibration files used in this work was PACS_CAL_69_0)666http://herschel.esac.esa.int/twiki/bin/view/Public/PacsCalTreeHistory.

We extracted the line fluxes for our LIRGs from the spectra using three different apertures: (a) the spaxel at which the 63$\,\mu$m continuum emission of the galaxy peaks (hereafter referred as “central” spaxel); (b) in a 3 $\times$ 3 spaxel box centered on the central spaxel as defined in (a), limited by the PACS FoV, and (c) the total FoV (5 $\times$ 5 spaxel box). In order to recover the total flux of the source from the spectrum extracted from the central spaxel (method (a)), we apply a point-source aperture correction (which varies as a function of wavelength). We note that this correction will only recover the total flux in sources that are unresolved by PACS. For extended sources the fluxes obtained in this manner will be lower limits to the integrated flux of the galaxy.

Table 1 presents the measurements made available in electronic form for each line ([O i]63, [O iii]88, [N ii]122, and [C ii]158). The line and continuum fluxes provided are those obtained using the (a) method as well as the “best” spatially-integrated value for each individual galaxy. This is defined as the highest flux value that maximizes the S/N ratio of the measurement (line or continuum). In other words, we select, among the three methods described above, the one that has the lowest noise in the measurement while still accounting for the entire flux of the galaxy enclosed by the PACS FoV. As an additional constraint we require that no other source is contained within the aperture used to represent the integrated galaxy flux. Note, however, that some galaxies have companions at distances $\lesssim$ 9.4″  (a PACS spaxel). These objects are marked in the tables and figures. As mentioned in section 2.4, the angular size of a PACS spaxel is roughly similar to that of the aperture used to extract the Spitzer/IRS spectra. We note that the Spitzer and Herschel pointings usually coincide within $\lesssim$ 2″. For further details regarding the pointing accuracy and centering of a source within a given spaxel we refer the reader to Díaz-Santos et al. (2013).

In order to estimate the fractional contribution of the PDR component to the total [C ii]158  emission based on the [C ii]158/[N ii]205  ratio in our LIRGs in section 5.1.1, we extracted an additional set of [C ii]158  spectra using a circular aperture with a diameter equal to the angular size of the SPIRE beam at 205 $\mu$m ($\approx$ 17″), to which we applied an aperture correction based on the PACS spectrometer beam efficiency maps (v6) provided in the PACS instrument and calibration webpages777http://herschel.esac.esa.int/twiki/bin/view/Public/PacsCalibrationWeb, after rebinning them to the regular 5 $\times$ 5 spaxel FoV. We note that we do not apply the correction factors based on the FIR color of galaxies provided in Lu et al. (2017) to the [N ii]205  emission. We use instead the original, point-source calibrated fluxes since we do not use those corrections in our Herschel/PACS data.

The individual FIR luminosities of galaxies belonging to a LIRG system formed by of two or more components were derived in a similar manner to the method used to calculate their $L_{\rm IR}$  (section 2.1). In order to obtain the $L_{\rm IR}$  and $L_{\rm FIR}$  of a given galaxy for the different extraction appertures described above, we scaled the integrated IRAS FIR luminosity of the system ($L_{\rm FIR\,[42.5-122.5\,\mu m]}$, as defined in Helou et al. 1985)888$F_{\rm FIR\,[42.5-122.5\,\mu m]}\,=\,1.26\,\times\,10^{-14}\,(2.58\,S_{\rm 60\mu m}\,+\,S_{\rm 100\mu m})$ [W m*-2*], with $S_{\nu}$ in [Jy]. $L_{\rm FIR}$ = $4\pi\,D_{\rm L}^{2}\,F_{\rm FIR}$. The luminosity distances, $D_{\rm L}$, are taken from Armus et al. (2009). and $L_{\rm IR}$  with the ratio between the continuum flux density of each individual galaxy evaluated at 63$\,\mu$m in the PACS spectrum (measured in the same aperture as the line) and the total IRAS 60$\,\mu$m flux density.

4. Results

4.1. The Contribution of AGN

Following the formulation described in Veilleux et al. (2009), we calculate the potential contribution of an AGN to the MIR and bolometric luminosities of each galaxy in GOALS (see also Petric et al., 2011) employing up to five Spitzer/IRS diagnostics, depending on the availability of data: the [Ne v]14.3/[Ne ii]12.8  and [O iv]25.9/[Ne ii]12.8  emission line ratios, the 6.2$\,\mu$m PAH  EW (see also Armus et al., 2007), the $S_{\rm 30}$/$S_{\rm 15}$  dust continuum slope999We modified the reference value for the $S_{\rm 30}$/$S_{\rm 15}$  ratio of a pure starburst/H ii  source from 22.4 (see Table 9 in Veilleux et al., 2009) to 10 in order to reflect more accurately the actual distribution of $S_{\rm 30}$/$S_{\rm 15}$  ratios seen in the GOALS sample. We also adapt the PAH EW diagnostic, which in (Veilleux et al., 2009) is developed for the 7.7 $\mu$m feature, to be used with the 6.2$\,\mu$m PAH, and modify the reference value for pure starbursts and its bolometric correction factor accordingly. We further assume there is no PAH destruction due to the AGN when calculating its fractional contribution using this method (Díaz-Santos et al., 2008; Díaz-Santos et al., 2010a; Esquej et al., 2014)., and the Laurent diagram, which is based on a decomposition of the MIR spectrum of galaxies in three individual components: AGN, PAH and H ii  emissions (see Laurent et al., 2000, for a detailed explanation of this diagnostic). These indicators are based on how the central AGN modifies the MIR line and continuum spectrum of a normal star-forming galaxy –either through the ionization of the surrounding gas to higher states and/or via the heating of dust at higher temperatures than star-formation–, and provide an estimate of its fractional contribution to the MIR emission, $\alpha_{\rm AGN}^{\rm MIR}$. Once these values are known, corrections based on SED templates of pure starbursting and AGN-powered sources can be applied to derive the fractional contribution of AGN to the bolometric luminosity of galaxies, $\alpha_{\rm AGN}^{\rm bol}$  (Veilleux et al., 2009, their Table 10, and discussion therein).

While individually each of these diagnostics has its own particular limitations, the combination of all of them allows for a reasonable quantification of the AGN’s average fractional luminosity contribution, <$$\alpha_{\rm AGN}^{\rm MIR}$$> and <$$\alpha_{\rm AGN}^{\rm bol}$$>. Because there is a significant number of non-detections in one or more of the relevant MIR features employed to calculate <$$\alpha_{\rm AGN}^{\rm MIR}$$> or <$$\alpha_{\rm AGN}^{\rm bol}$$>, we used the Astronomy SURVival analysis package, asurv (Lavalley et al., 1992, v1.3), which adopts the maximum-likelihood Kaplan-Meier (KM) estimator to compute the mean value of univariate distributions containing censored data (Feigelson & Nelson, 1985). Table 2 presents the $\alpha_{\rm AGN}^{\rm MIR}$  values for each MIR diagnostic, as well as <$$\alpha_{\rm AGN}^{\rm MIR}$$> and <$$\alpha_{\rm AGN}^{\rm bol}$$>.

4.2. FIR Line Deficits

Using the best galaxy-integrated measurements described in section 3.2, we present in Figure 1 the [O i]63, [O iii]88, [N ii]122, [C ii]158  and [N ii]205  emission line deficits101010The word “deficit” refers to the deficiency of a given line flux when compared to the dust continuum emission in a galaxy. This term was historically coined to express the decrease in the line (gas) cooling efficiency with respect to that of the dust in IR luminous galaxies. –expressed as the line-to-FIR continuum flux ratio– for the entire GOALS LIRG sample as a function of the $S_{\nu 63\,\mu m}/S_{\nu 158\,\mu m}$ ($S_{\rm 63}$/$S_{\rm 158}$) continuum flux density ratio, which is a first order tracer of the average dust temperature in galaxies, $T_{\rm dust}$. The dynamic range in both x- and y-axes is the same in all panels to facilitate the comparison. We have fitted the data using a functional form of the type:

[TABLE]

where $\epsilon_{0}$ denotes a limiting line/FIR ratio for sources with cold FIR colors and no deficits, and $\delta$ is the $S_{\rm 63}$/$S_{\rm 158}$  at which the line/FIR ratio has been reduced by a factor of e with respect to $\epsilon_{0}$. The $\epsilon_{0}$ parameter can be understood as the nominal cooling efficiency of each line with respect to that of big dust grains, representative of normal star-forming galaxies. The parameters obtained from the fits to the entire GOALS sample can be found in Table 3 and the fits are displayed in Figure 1 as solid black lines, with the associated 1 $\sigma$  dispersion around the fit shown as dotted lines. While the choice of this particular functional form is arbitrary, it provides a better description of the trends than a power-law fit, as it is known that the line to FIR ratios do not increase indefinitely at the lower end of the $T_{\rm dust}$  distribution and IR luminosities (e.g., Malhotra et al., 1997, 2001; Brauher et al., 2008; De Looze et al., 2014; Cormier et al., 2015; Herrera-Camus et al., 2015). That is, the line to FIR ratios level off at low $T_{\rm dust}$  reflecting a cooling efficiency “ceiling” (which depends on each line) of the gas in PDRs, H ii  regions and the diffuse ISM, with respect to the energy dissipated by dust in thermal equilibrium in normal galaxies.

We convert the $S_{\rm 63}$/$S_{\rm 158}$  ratios into dust temperatures (see upper x-axis in Figure 1) by assuming that the observed FIR continuum emission is produced by a single-temperature modified black body (mBB) with a fixed emissivity index $\beta$ = 1.8 and whose emission is optically thin. This is a reasonable approximation for spectral energy distribution (SED) fits that do not include data at $\lambda$ $\lesssim$ 60 $\mu$m. We also provide a practical equation that relates both quantities, $S_{\rm 63}$/$S_{\rm 158}$  and $T_{\rm dust}$, using the following approximation:

[TABLE]

We note that this equation is only valid for the dust temperatures and $S_{\rm 63}$/$S_{\rm 158}$  ratios spanned by the galaxies in the GOALS sample. That is 21 $\leq$ $T_{\rm dust}$ $\leq$ 48 K or 0.1 $\leq$ $S_{\rm 63}$/$S_{\rm 158}$ $\leq$ 4. The error in $T_{\rm dust}$  obtained from this expression is $\lesssim$ 0.5 K for the dynamic ranges mentioned.

Figure 1 shows that there is a common trend for most lines to show stronger deficits as the average $T_{\rm dust}$  becomes warmer –including the two [N ii]  lines, which arise from the ionized medium. The [O i]63  deficit shows a decline of approximately an order of magnitude and a large scatter. The [N ii]122, [C ii]158  and [N ii]205, exhibit stronger deficits, of up to two orders of magnitude, and tighter trends. For the [O iii]88  line, although there may be a deficit at the highest $S_{\rm 63}$/$S_{\rm 158}$  values ($\gtrsim$ 2), the exponential fit yields a result that is statistically indistinguishable from a flat trend. Binning the data and obtaining the median values for each bin provides the same result.

A possible interpretation is that the dispersion in the [O iii]88/FIR ratio as a function of $T_{\rm dust}$  may be related to the location where the line and dust emission originate within the star-forming regions. As we argue in sections 5.1.3 and 5.5 (see also Díaz-Santos et al., 2013), most of the energy reprocessed by dust may be arising from grains in front of the PDRs, mixed with the ionized gas, and at low optical depths into the molecular cloud. Low ionization lines, that originate from the PDR itself or the outer edge of the H ii  region ([N ii]122, [N ii]205  and [C ii]158) would show less scatter around the mean. On the other hand, lines arising from highly ionized gas within the H ii  regions, like [O iii]88, and lines for which part of the emission may be arising from deeper into the molecular cloud, like [O i]63, would show larger dispersion around the trend average (see $\sigma$/$\epsilon_{0}$ in Table 3). The potential contribution of an AGN could also introduce additional dispersion in the line deficits, specially at high $S_{\rm 63}$/$S_{\rm 158}$  (Fischer et al., 2014; González-Alfonso et al., 2015) or when the AGN overwhelmingly dominates the MIR emission (Díaz-Santos et al., 2014), but we do not find galaxies with large fractional MIR or bolometric AGN contributions to be distributed with a significantly larger scatter than star-formation dominated sources, based on the Spitzer/IRS spectral diagnostics described in section 4.1.

Figure 2 presents the same line deficits but as a function of the luminosity surface density, $\Sigma_{\rm IR}$, defined as the ratio of effective luminosity, $L_{\rm IR,eff}$ = $L_{\rm IR}$/2, divided by the effective area of the source that contains half of its luminosity measured at 70 $\mu$m, $\pi\rm R^{2}_{\rm 70\mu m,eff}$  for galaxies with available FIR size measurements. These correlations are overall tighter than with $T_{\rm dust}$  (including [O iii]88), suggesting a closer physical connection between the cause(s) that give rise to the line deficits, and the concentration of dust-reprocessed energy –or IR “compactness”– of LIRGs (see discussion in section 5.4). The scatter in the trends of those lines that have a PDR origin, [O i]63  and [C ii]158, is especially small and remarkably constant in relative terms at any $\Sigma_{\rm IR}$. Moreover, the trend is followed by nearly all LIRGs regardless of their FIR color or $T_{\rm dust}$. We have fitted these correlations with a second-order polynomial function:

[TABLE]

The best fit parameters are presented in Table 4. Note that in Figure 2 the fits have been set to the maximum value of the respective quadratic equation below the wavelength at which the maximum is reached, such that the ratio remains constant.

5. Discussion

5.1. [C ii]158  Emission from Ionized Gas and PDRs

Because of the low ionization potential of the Carbon atom (11.26 eV), the [C ii]158  emission line can be produced not only in regions of ionized gas ([C ii]${}^{\rm ion}_{158}$) but also in the dense, neutral ISM ([C ii]${}^{\rm PDR}_{158}$). However, due to the low critical density of the transition when it is collisionally excited by free electrons or protons ($n^{\rm cr,[CII]}_{\rm e}$ $\simeq$ 44 cm*-3*, at $T$ = 8000 K; Goldsmith et al., 2012), the [C ii]158  is rapidly thermalized in mildly dense, ionized environments. Thus, unless the volume filling factor of the diffuse medium is very high (see section 5.1.2), most of the [C ii]158  emission –especially in actively star-forming galaxies– is expected to arise from dense PDRs surrounding young, massive stars (Hollenbach & Tielens, 1997), where the [C ii]158  is collisionally excited by neutral and molecular hydrogen, with $n_{\rm H}^{cr}$,[CII] $\simeq$ 3.0 $\times$ 103 cm*-3*  and $n_{\rm H_{2}}^{cr}$,[CII] $\simeq$ 6.1 $\times$ 103 cm*-3*, at $T$ = 100 K (Goldsmith et al., 2012).

5.1.1 Dense PDRs

We can use the [N ii]205  line to estimate the amount of [C ii]158  emission produced in the ionized phase of the ISM, [C ii]${}^{\rm ion}_{158}$  (e.g., Oberst et al., 2006; Beirão et al., 2012; Croxall et al., 2012; Kapala et al., 2015). The layer within the Strömgren sphere where nitrogen is singly ionized ranges between 29.60 and 14.53 eV, close to that where [C ii]${}^{\rm ion}_{158}$  also originates, 24.38 to 13.6 eV. In addition, given the similar $n_{\rm e}^{\rm cr}$  and $E_{\rm ul}$/$k_{\rm B}$  of both transitions, the [C ii]${}^{\rm ion}_{158}$/[N ii]205  ratio is roughly constant and depends weakly on $n_{\rm e}$, the intensity of ionizing field, $q$, and the kinetic temperature of the gas, $T^{\rm kin}_{\rm gas}$.

The photo-ionization models presented in Oberst et al. (2006) predict a [C ii]${}^{\rm ion}_{158}$/[N ii]205 $\simeq$ 3 $\pm$ 0.5 for a range of $n_{\rm e}$ up to $\sim$ 103 cm*-3*  (see section 5.1.2). For convenience, we use this constant ratio to estimate [C ii]${}^{\rm ion}_{158}$  and subtract it from the total [C ii]158  flux, which yields [C ii]${}^{\rm PDR}_{158}$  (= [C ii]158–[C ii]${}^{\rm ion}_{158}$). Figure 3 (upper panel) shows the fraction of [C ii]  arising from PDRs, f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}{158}}) = [C ii]${}^{\rm PDR}_{158}$/[C ii]158, as a function of $S_{\rm 63}$/$S_{\rm 158}$  and $T_{\rm dust}$  for the entire GOALS sample. In Díaz-Santos et al. (2013) we showed that the presence of an AGN in LIRGs does not play a role in the decreasing of the [C ii]158/FIR  ratio as a function of $L_{\rm IR}$  or $T_{\rm dust}$. To show that it does not have an affect on f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}{158}}) either, we color-code Figure 3 a function of the fractional contribution of the AGN to the MIR emission of each galaxy, based on the 6.2 $\mu$m PAH EW diagnostic (see section 4.1). Neither sources with large MIR AGN fractions nor those with the largest bolometric contributions (<$$\alpha_{\rm AGN}^{\rm MIR}$$> $\geq$ 0.5) are found systematically in a different region of the parameter space than star formation dominated galaxies – or contribute to increase the scatter of the correlation.

A fit to the data provides:

[TABLE]

with a dispersion of 0.07. For reference, the location of the Milky Way is also shown in the figure, assuming a luminosity-weighted average $T_{\rm dust}$ = 25 ($\pm$ 5) K and considering that PDR emission accounts for the remaining amount of [C ii]  that is not associated with the ionized medium, which is estimated to be between $\sim$ 1/3 and 1/2 (Goldsmith et al., 2015). While the scatter is large, we can identify a broad trend for galaxies with warmer $T_{\rm dust}$  to show larger f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}_{158}}) (black solid line), increasing from $\sim$ 60 %, close to the MW value, to nearly 95 %. That is, there is a larger PDR contribution to the total [C ii]158  emission in warmer systems, indicating that, even though the [C ii]158  line shows a large deficit with respect to the FIR emission in progressively more luminous, warmer galaxies, most of the extra [C ii]  produced in them originates in dense PDRs.

This scenario is supported by the trend presented in Figure 3 (bottom panel), which shows also a positive correlation of f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}{158}}) with the average hardness of the radiation field seen by the ionized gas, as traced by the [O iii]88/[N ii]122  line ratio (see section 5.2.3). Therefore, f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}{158}}) increases in environments associated with recent episodes of massive star formation. This is also consistent with H ii  regions in the LIRG nuclei to be more enshrouded (optically and geometrically thicker; see discussion in section 5.4) than those of evolved star-forming complexes where the stellar winds from massive stars and supernovae have already cleared out most of the dust from the star formation sites (Blitz & Shu, 1980; Larson, 1981). That is, galaxies with more evolved H ii  regions (i.e., posterior to experience a starburst event; or simply having more modest star-formation, like the MW) have more of the [C ii]158  emission arising from the ionized gas (likely from the low density ISM, see below). We also note that there is no clear trend between f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}_{158}}) and the electron density of the ionized gas, $n_{\rm e}$, as traced by the [N ii]122/[N ii]205  ratio (see color code in the figure). We discuss further implications of the trend seen in Figure 3 (top) in section 5.1.3.

5.1.2 Ionized Gas

As shown in Figure 3, even though the [C ii]${}^{\rm PDR}_{158}$  dominates the total [C ii]158  emission, the [C ii]${}^{\rm ion}_{158}$  contribution is not negligible. This ionized component can subsequently originate from both the diffuse ([C ii]${}^{\rm diff,{\sc HII}}_{158}$) and dense medium ([C ii]${}^{\rm ion,{\sc HII}}_{158}$). Inami et al. (2013) used Spitzer/IRS high resolution spectroscopy to probe modestly ionized gas within the H ii  regions of the LIRG nuclei via the [S iii]18.7$\,\mu$m/[S iii]33.5$\,\mu$m  line ratio111111Both lines have $n_{\rm e}^{\rm cr}$ $\gtrsim$ 2000 cm*-3*., and calculated the average $n_{\rm e}$  for most galaxies to be typically $\sim\,$ 100 to a few hundred cm*-3*, with a median of $\sim$ 300 cm*-3*. Nearly 30 % of the galaxies for which both lines are detected, though, show ratios consistent with the ionized medium being in the low density limit ($n_{\rm e}$ $\lesssim$ 100 cm*-3*). We note that the layer within the Strömgren sphere where Sulfur is doubly ionized is located between 34.79 and 23.34 eV, while C2+ transitions to C*+* at $h\nu$ $<$ 24.38 eV. Considering that the electron density should at least remain constant, if not increase towards the denser PDR region ($n_{\rm e}$  of a few 103 cm*-3*  have been found in H ii  regions using optical emission lines of singly ionized sulfur; Osterbrock, 1989), this means that the $n_{\rm e}$  derived using the [S iii]18.7/[S iii]33.5  ratio likely represents a lower limit to the density of the volume from where [C ii]${}^{\rm ion}_{158}$  and singly ionized nitrogen emission arises within H ii  regions. And because the densities derived from the Sulfur lines are significantly larger than the [C ii]${}^{\rm ion}_{158}$  and [N ii]205  critical densities, this implies that the [C ii]${}^{\rm ion,{\sc HII}}_{158}$  emission has likely been thermalized, therefore suggesting that most [C ii]${}^{\rm ion}_{158}$  and [N ii]205  emission is produced in the diffuse ISM, with a modest contribution from the dense ionized phase – unless the diffuse medium is extremely thin (low volume density) or its average $T^{\rm kin}_{\rm gas}$  very low.

We can also use the [N ii]122/[N ii]205  ratio in combination with the models from Oberst et al. (2006) to derive the average electron densities in our LIRG sample. As mentioned above, the region where emission from singly ionized nitrogen atoms originates in an H ii  region largely overlaps with that of the [C ii]${}^{\rm ion,{\sc HII}}_{158}$  emission. Figure 4 shows the distribution of [N ii]122/[N ii]205  ratios and $n_{\rm e}$  for those galaxies with available measurements of both lines. We find densities between $\sim$ 20 and 100 cm*-3*, with a median value of 41 cm*-3*  and mean of 45 cm*-3*. These values are very similar to what has been found by other studies of normal and starbursting galaxies. For instance, Zhao et al. (2016b) find $n_{\rm e}$ = 22 cm*-3*  for a subsample of GOALS LIRGs using ISO data, and Herrera-Camus et al. (2016) find $n_{\rm e}$ = 30 cm*-3*  for spatially resolved regions of 21 nearby, normal star-forming galaxies selected from the Herschel KINGFISH and Beyond the Peak Herschel surveys. In the Milky Way, the average value measured by Goldsmith et al. (2015) with Herschel/PACS is 29 cm*-3*.

As we noted above, the $n_{\rm e}$  values we find based on the nitrogen line ratio are smaller than those inferred from the sulfur lines, suggesting that both [N ii]205  and [C ii]${}^{\rm ion}_{158}$  transitions have been thermalized in the dense H ii  regions and thus mostly originate from the diffuse ISM. In this case, if we assume a $n_{\rm e}$ $\simeq$ 500 cm*-3*  and a $n_{\rm e}$ $\lesssim$ 1 cm*-3*  for the dense and diffuse ionized medium respectively, the mean value of the [N ii]122/[N ii]205  ratio implies an average volume filling factor, $\phi_{\rm V}$, of $\lesssim$ 5 % for the H ii  regions with respect to the overall volume of ionized emitting gas. Of course, this is only a very rough estimate, and $\phi_{\rm V}$  will vary significantly depending on the compactness of the galaxy.

A number of reasons have been proposed to explain the observed deficit of FIR lines of purely ionized species – mostly singly and double ionized nitrogen and oxygen ([N ii]122, [N ii]205, [N iii]57, [O iii]52  and [O iii]88). A boost in the average intensity of the radiation field has been suggested by several authors, who have been able to reproduce the low [N ii]122/FIR ratios observed in warm LIRGs using the Cloudy spectral synthesis code (Ferland et al., 1998) by increasing the ionization parameter from log(U) $\sim$ –2.5 to –1.5 (Voit, 1992; Abel et al., 2009; Graciá-Carpio et al., 2011; Fischer et al., 2014). However, the deficits of the higher ionization species, [N iii]57  [O iii]52, and [O iii]88  could not be explained by this effect. In view of the discussion above, our results indicate that a progressive thermalization of those lines with the lowest critical densities could also contribute to the observed deficits. Indeed, we find that the [N ii]205  line, which has the lowest $n_{\rm e}^{\rm cr}$  among our set of ionized lines, presents the largest decline in the [N ii]205/FIR ratio, of almost two orders of magnitude; followed by the [N ii]122, and [O iii]88  lines, which have increasingly higher $n_{\rm e}^{\rm cr}$.

5.1.3 Dust Temperature and PDR fraction as Drivers of the [C ii]158  Deficit

The total luminosity of an optically-thin modified black-body (mBB) increases as a function of its temperature as $T^{4+\beta}$, where $\beta$ is the emissivity index of the emitting material. Thus, an increase of the average $T_{\rm dust}$  in a LIRG from $\sim$ 25 to $\sim$ 50 K can boost its $L_{\rm FIR}$  by at least a factor of $\sim$ 50 (for $\beta$ = 1.8; Planck Collaboration et al. 2011)121212Note that the peak emission of a modified black-body around these temperatures is always within the wavelength range used to calculate $L_{\rm FIR}$, and therefore it can be approximately scaled using the Stephan-Boltzan law.. As we argued in Díaz-Santos et al. (2013), the [C ii]158  deficit is likely caused by an increase of the amount of dust heated to high temperatures in the ionized region –close to the PDR front– of dust-bounded star-forming regions (see e.g., Abel et al., 2009; Paladini et al., 2012), where the density of material is expected to increase significantly (Draine, 2011). Furthermore, the deficit seems to be restricted to the nuclear region of galaxies, where the starburst is ongoing (Díaz-Santos et al., 2014; Smith et al., 2016). Thus, most of the FIR dust continuum emitted at the BB peak would be associated with active, still dust-bounded star-forming regions. On the other hand, and in light of our findings in the previous sections, we expect the [N ii]205  and [C ii]${}^{\rm ion}_{158}$  emission to originate mostly from low density, “fossil” H ii  regions and diffuse ionized gas not associated with OB stars. Then, if the [C ii]${}^{\rm ion}_{158}$  and FIR emissions truly arise from different ISM phases, the [C ii]${}^{\rm ion}_{158}$/FIR  ratio should be mostly controlled by the $L_{\rm FIR}$  boosting, due to the increase of $T_{\rm dust}$  in progressively more active starbursts.

The top left panel of Figure 5 shows [C ii]${}^{\rm ion}_{158}$/FIR  as a function of $T_{\rm dust}$  (derived from the $S_{\rm 63}$/$S_{\rm 158}$  ratio) as well as a fit in log-log space to the data (dashed line; see the figure caption for the precise equation). Under the assumptions discussed above, this trend indeed suggests that the rise in the average radiation field intensity produced by young, massive stars born during the ongoing starburst is mostly transferred into an increase of $T_{\rm dust}$, which enhances the [C ii]158  deficit through the boosting of the FIR continuum. If we assume that the FIR emission scales with \hbox{T_{\rm dust}}^{4+\beta} (where $\beta$ = 1.8), the suppression seen in the [C ii]${}^{\rm ion}_{158}$/FIR  ratio of the warmest LIRGs can be explained in terms of a factor of almost $\sim$ 30 excess in the IR continuum emission rather than a deficit in the observed line flux. This scaling is also able to recover reasonably well the shape of the trend seen in Figure 5. The fit to the data is even better if we consider the BB emission to be optically thick below 100 $\mu$m – an assumption frequently made when fitting the SED of heavily dust-obscured sources (e.g., Wilson et al., 2014).

However, the increase in the dust cooling efficiency131313The term “dust cooling efficiency” is used here to refer to the emitted luminosity per dust particle (or per dust mass, assuming that there are no significant variations in the dust grain composition and/or size distribution among our galaxies) due to the increase in temperature. alone cannot be the only physical parameter involved in the variation of the line-to-FIR ratios because: (1) that would imply that none of the energy produced by massive stars at the core of star-forming regions would be transformed into gas heating, and (2) there is still a significant dispersion of a factor of $\sim$ 3 or larger in the [C ii]158/FIR  ratio at a given $T_{\rm dust}$  (see Figure 1). The top right panel of Figure 5 shows the [C ii]158/FIR  ratio after correcting the FIR emission for the increase due to hotter $T_{\rm dust}$  using the fit to the [C ii]${}^{\rm ion}_{158}$/FIR  ratio described above (top left panel). After this “IR excess” has been taken into account, we can see a trend for [C ii]158  emission to rise as a function of $T_{\rm dust}$, suggesting that part of the heating from the stellar radiation field is indeed transferred to the gas, increasing the cooling carried out through the [C ii]158  (and other) line(s). Note that the dispersion also increases with $T_{\rm dust}$  as a result of the PDR contribution to the total [C ii]158  budget (see color coding in the figure).

Figure 5 (bottom) conveys the same information but in the form of the classical [C ii]158/FIR  ratio, which is displayed as a function of $S_{\rm 63}$/$S_{\rm 158}$  and $T_{\rm dust}$. Note that the $S_{\rm 63}$/$S_{\rm 158}$  ratio shown here is approximately equal to the $S_{\rm 70}$/$S_{\rm 160}$ ratio that can be measured with Herschel/PACS or SOFIA/FIFI-LS in spatially resolved nearby galaxies, or at high redshifts with a combination of ALMA bands. Regions of similar f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}_{158}}) contribution are marked with the same line styles as in the top right panel. This figure shows again that the dispersion observed in the [C ii]${}^{\rm ion}_{158}$/FIR  ratio as a function of $T_{\rm dust}$  is caused by the different contributions of the PDR component to the total [C ii]158  emission. This is consistent with the scenario put forward in section 5.1.1 regarding a larger fraction of a given dust mass ($M_{\rm dust}$) being heated up to higher temperatures by the absorption of ionizing photons in the most compact, luminous starbursts. Because the cooling efficiency of dust increases as $T_{\rm dust}$4+β, a significant amount of dust present within the H ii  regions, close to the PDRs (see Draine, 2011), could be able to reprocess most of the energy emitted by the stars before it reaches the surrounding neutral/molecular medium. Indeed, the temperatures measured from our data are typical of PDR surfaces (i.e., at $A_{V}$ $\sim$ 0) exposed to moderate to high radiation fields ($G$ $\gtrsim$ 102$G_{\rm 0}$; Tielens 2005; see section 5.3). In addition, a fraction of the energy is transferred via the photo-electric effect into gas heating, as shown in the top right panel, which cools down, increasing the [C ii]${}^{\rm PDR}_{158}$  emission.

Recently, using Herschel/HIFI, Goicoechea et al. (2015) have studied the variation of the [C ii]158/FIR  ratio on $\sim$ 1 pc scales across the Orion molecular cloud 1 (OMC1) and the region surrounding the Trapezium cluster, which is being ionized by the intense UV field of massive O stars. They find that there is a broad correlation between [C ii]158/FIR  and the column density of dust through the molecular cloud –tighter than with $L_{\rm FIR}$/$M_{\rm gas}$  (see Graciá-Carpio et al., 2011)–, and conclude that the [C ii]158 -emitting column relative to the total dust column along each line of sight is responsible for the variations of [C ii]158/FIR  observed through the cloud. This trend is similar to that found for the nuclear emission of LIRGs, which show a correlation between [C ii]158/FIR  and (1) the strength of the 9.7$\mu$m silicate absorption feature, $S_{\rm 9.7\mu m}$, probing the opacity, and thus the column density and total mass of the MIR-absorbing dust (Díaz-Santos et al., 2013); and (2) the opacity of the free-free emission in H ii  regions as measured through the spectral index between 1.5 and 6 GHz continuum emission (Barcos-Muñoz et al., 2017). However, although these results are in agreement at face value, the argument made for the case of the OMC1 by Goicoechea et al. (2015) –in which the [C ii]158  deficit would be entirely caused by an increase in the dust column density– is probably not the complete answer for the case of the kpc-scale, integrated emission of LIRGs. In terms of total $M_{\rm dust}$, the decrease of more than an order of magnitude in [C ii]158/FIR  would require a similar increase in the amount of dust mass in the galaxy, something that is not seen even in the most extreme ULIRGs, which have similar $M_{\rm dust}$  as dusty, normal star-forming (MS) galaxies (da Cunha et al., 2010). Also, we note that the column density fit used in Goicoechea et al. (2015) assume that there are no dust temperature variations. On the other hand, Lombardi et al. (2014) presented column density and effective dust temperature maps based on NIR, Herschel and Planck data for the entire Orion complex at a lower angular resolution than Goicoechea et al. (2015). Considering both works, temperature variations are in the range $\sim$ 15–35 K, where regions with relatively high temperatures emit significantly higher fluxes, even if the optical depth is substantially lower.

In summary, while it is possible that the column density in LIRGs may indeed increase due to a compactification of the star-forming regions (see end of next section), this interpretation is still compatible with the picture put forward above, where a given amount of dust mass would be accumulated closer to the heating source, thus also increasing $T_{\rm dust}$  and the corresponding radiative cooling. Disentangling between the two effects is, however, extremely challenging in galaxies where it is not possible to spatially resolve the different temperature components of the dust emission.

5.2. FIR Oxygen Line Emission

5.2.1 [O i]63/[C ii]${}^{\rm PDR}_{158}$  as a Probe of the Gas Kinetic Temperature

The [O i]63  line is an important coolant of the warm, neutral ISM. A substantial fraction of the [O i]63  emission is believed to originate in PDRs, mostly at optical depths similar to those of the locations responsible for the C+ emission, $A_{V}$ $\lesssim$ 4 – although neutral oxygen may exist deeper in the molecular cloud before it combines into CO at $A_{V}$ $\sim$ 10 (Hollenbach & Tielens, 1997; Kaufman et al., 1999). While [O i]63  can be as weak as $\sim$ 20 % of the [C ii]158  luminosity emitted in a normal star-forming galaxy, its relative contribution to the ISM cooling can become increasingly important in the dense and warm environments of extreme starbursting galaxies ($n^{\rm cr,[OI]_{63}}_{\rm H}$ $\simeq$ 106 cm*-3*), even surprassing the [C ii]158  line in the warmest LIRGs (Rosenberg et al., 2015). Using ISO/LWS data, Malhotra et al. (2001) found a correlation between the [O i]63/[C ii]158  ratio and the FIR color of star-forming galaxies. Our Herschel/PACS data and Spitzer/IRS observations allow us to explore this correlation in depth for the entire GOALS sample; a 4-fold increase in the number of galaxies available relative to previous studies.

Because of the neutral gas phase origin of the oxygen fine-structure emission, we use the [C ii]${}^{\rm PDR}_{158}$  component alone to compare it with the [O i]63  line. For those galaxies without [N ii]205  observations we use equation 4 to estimate their f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}_{158}}) (capped at 90 %) from the $S_{\rm 63}$/$S_{\rm 158}$  ratio. Figure 6 shows that [O i]63/[C ii]${}^{\rm PDR}_{158}$  is tightly correlated with the FIR $S_{\rm 63}$/$S_{\rm 158}$  continuum ratio, with a scatter of only 0.13 dex. We have used the statistical equilibrium radiative transfer code RADEX (van der Tak et al., 2007), to model the variation of the [O i]63/[C ii]${}^{\rm PDR}_{158}$  ratio as a function of the kinetic temperature of the gas ($T^{\rm kin}_{\rm gas}$), assuming $T^{\rm kin}_{\rm gas}$ = $\eta$ $T_{\rm dust}$. We have further assumed that the excitation is dominated by collisions with neutral hydrogen atoms, a gas density of $n_{\rm H}$ = 103 cm*-3*  (see section 5.3), and a column density of $N_{\rm H}$ = 1016 cm*-2*. In this regime, and for lower volume densities and columns up to $N_{\rm H}$ $\simeq$ 1018 cm*-2*, both lines are optically thin and also $n_{\rm H}$ $<$ $n_{\rm H}^{cr}$. Despite the simplicity of the approach, the RADEX model reproduces the observed trend between [O i]63/[C ii]${}^{\rm PDR}_{158}$  and $T_{\rm dust}$  remarkably well (black solid line in Figure 6; the uncertainty around the trend is 1$\sigma$ = 0.30, dotted lines) with a best fit for $\eta$ = 2.06 $\pm$ 0.03. The value of $\eta$ decreases to 2.02 and 1.60 for $N_{\rm H}$ = 1017 and 1018 cm*-2*  respectively, since at higher columns the gas does not need to be as warm to produce the same [O i]63/[C ii]${}^{\rm PDR}_{158}$  ratios. With any of these scalings we infer gas kinematic temperatures similar or lower than the excitation energy of the [C ii]158  fine-structure transition, $E_{\rm ul,[CII]_{158}}$/$k_{\rm B}$ = 92 K. This is in contrast to suggestions that [C ii]158  deficit is mainly due to thermal saturation of the [C ii]158  emission (Muñoz & Oh, 2016).

The trend can be fitted as a function of $S_{\rm 63}$/$S_{\rm 158}$  using a second order polynomial function:

[TABLE]

Inversely, one can obtain $S_{\rm 63}$/$S_{\rm 158}$  as a function of the line ratio:

[TABLE]

with $x$ = [O i]63/[C ii]${}^{\rm PDR}_{158}$, and using equation 2 obtain $T_{\rm dust}$  from $S_{\rm 63}$/$S_{\rm 158}$, which is interchangeable with $T^{\rm kin}_{\rm gas}$  assuming the best scaling factor obtained from the model fit to our data, $T^{\rm kin}_{\rm gas}$ = 2.06 $T_{\rm dust}$. We caution that these equations assume that the [O i]63  line is optically thin, and are only valid approximately over the range spanned by our galaxy sample, i.e., 0.25 $\lesssim$ [O i]63/[C ii]${}^{\rm PDR}_{158}$ $\lesssim$ 3.5 and 28 $\lesssim$ $T_{\rm dust}$ $\lesssim$ 76 K, or equivalently over the ranges shown in Figure 6.

Some of the warmest, most luminous galaxies seem to deviate from the correlation, showing systematically lower [O i]63/[C ii]${}^{\rm PDR}_{158}$  ratios than those predicted by the trend. The most likely explanation is that in these LIRGs the [O i]63  emission becomes optically thick due to an increase of the in situ gas column density and/or foreground absorption by cold gas across the line of sight (e.g., Poglitsch et al., 1996). This is in agreement with the fact that most galaxies falling below the fitted correlation show large 9.7 $\mu$m silicate strengths as well, $S_{\rm 9.7\mu m}$ $\lesssim$ –1.5 (see color-coding in Figure 6), indicating a similar increase of the dust opacity towards the hot, MIR emitting background source. Optically thick [O i]63  has been already reported by studies with available information of the [O i]145  emission line (Malhotra et al., 2001; Farrah et al., 2013; Rosenberg et al., 2015). Since [O i]145  is optically thin, a value of the [O i]145/[O i]63  ratio in excess of $\sim$ 0.1 implies some degree of optical thickness in [O i]63.

There are only 12 galaxies in GOALS ($\sim$ 5% of the sample) displaying [O i]63/[C ii]${}^{\rm PDR}_{158}$  ratios lower than half the value of the fit shown in Figure 6 (this threshold is denoted by the dotted-dashed line), and five of them show clear absorption in the [O i]63  line profile (IRASF08572+3915, NGC 4418, Arp 220, IRASF17207–0014 and IRAS17578–0400), in some cases showing P-Cygni or inverse P-Cygni profiles indicating outflowing or inflowing gas (González-Alfonso et al., 2012; Falstad et al., 2015, 2017). Figure 12 in the Appendix shows the Herschel/PACS spectra of these galaxies, where signs of [O i]63  absorption are visually noticeable. Additionally, a small fraction of galaxies in GOALS with spectra covering the wavelength range of the OH doublet at 119 $\mu$m, OH119  (IRAS07251–0248, IRASF15250+3609 and IRASF17207–0014) show signs of strong molecular absorption, which is characteristic of FIR-thick, warm, compact objects (González-Alfonso et al., 2015). Also, in the ULIRGs sample from Farrah et al. (2013), 4 out of the 5 galaxies that are below the [O i]63/[C ii]${}^{\rm PDR}_{158}$  threshold also present signs of OH119  absorption. In these highly obscured sources, it is possible that an AGN contributes to the dust heating (Spoon et al., 2013), but it is impossible to determine the exact nature of the buried source due to its extreme optical thickness.

5.2.2 [O iii]88/[O i]63  and the Ionized/Neutral Gas Filling Factors

In contrast to the [O i]63  line, [O iii]88  emission originates from ionized gas with lower densities ($n_{\rm e}^{\rm cr}$ $\sim$ 500 cm*-3*). Thus, the two lines clearly probe distinct phases of the ISM, providing a first order approximation to the ratio of volume filling factors, $\phi_{\rm V}$, between the warm ionized and neutral gas, assuming the lines are optically thin and that the ratio of temperatures between both phases is relatively constant across galaxies. As pointed out by Cormier et al. (2015), the [O iii]88/[O i]63  ratio is between 2–10 times higher in low metallicity galaxies than in normal star-forming galaxies, which suggest a correspondingly higher $\phi_{\rm V}$  of ionized gas with respect to PDR emission in the former.

In Figure 7 we show [O iii]88/[O i]63  as a function of the FIR $S_{\rm 63}$/$S_{\rm 158}$  ratio for the entire GOALS sample as open circles. We also include the low metallicity galaxies from the Dwarf Galaxy Sample (DGS) presented in Cormier et al. (2015). To obtain the $S_{\rm 63}$/$S_{\rm 158}$  ratio for the low metallicity sample, which only have published data of the $S_{\rm 70}$/$S_{\rm 100}$ ratio, we use the same method described in section 4.2 to obtain $T_{\rm dust}$  and then interpolate the fitted mBB emission at 63 and 158 $\mu$m to obtain the $S_{\rm 63}$/$S_{\rm 158}$  ratio for each galaxy. Figure 7 suggests that there is no clear correlation between [O iii]88/[O i]63  and the FIR color in LIRGs. The median [O iii]88/[O i]63  ratio of our sample is 0.46 $\pm$ 0.30. Low metallicity galaxies populate a region of the parameter space with higher $T_{\rm dust}$  ($\gtrsim$ 30 K) and larger $\phi_{\rm V}$  of ionized gas ($\gtrsim$ 5 $\times$) than dusty, metal-rich systems. Figure 7 also indicates that galaxies with [O iii]88/[O i]63  ratios in excess of $\simeq$ 1 have on average larger [Ne iii]15.6/[Ne ii]12.8  emission line ratios as well ($\gtrsim$ 0.5), regardless of the nature of the sample. This is in agreement with a picture in which a larger volume of the ISM is ionized, which is likely caused by a combination of less effective dust cooling (due to a lack of metals) and by the presence of harder ionization fields, rising the gas temperature and the ionization state. In addition, in dusty galaxies the PDR covering factor is likely higher than in normal or low metallicity sources, and thus more energy is captured by neutral and molecular material and cooled via the [O i]63  line emission.

We note that there are a few LIRGs lying in the same region of the diagnostic diagram as low metallicity galaxies, with [O iii]88/[O i]63 $\geq$ 1. However, most of them have $S_{\rm 9.7\mu m}$ $\lesssim$ –1.5 and are outliers in the [O i]63/[C ii]${}^{\rm PDR}_{158}$  vs. $S_{\rm 63}$/$S_{\rm 158}$  correlation (sources below the threshold marked by the dotted-dashed line in Figure 6), indicating a significant optical thickness of the [O i]63  line, which would artificially boost the oxygen line ratio to higher values (see pluses in Figure 7).

5.2.3 [O iii]88/[N ii]122: A Tracer of the Ionization Field Hardness?

The similar critical densities of the [O iii]88  and [N ii]122  emission lines ($n_{\rm e}^{\rm cr}$ $\sim$ 500 cm*-3*  and $\sim$ 300 cm*-3*, respectively) but different ionization potentials ($\sim$ 35.1 eV and 14.5 eV, respectively) make their ratio a good tracer of the average hardness of the radiation field in a galaxy, or equivalently, the mean, luminosity-weighted effective temperature of the massive stars born in the starburst, $T_{\rm eff}$. Indeed, within the context of star formation, [O iii]88  emission is mostly powered by O stars ($T_{\rm eff}$ $\gtrsim$ 30000 K (e.g., Ferkinhoff et al., 2011). Since O stars are short lived, with lifespans of around a few Myr ($\lesssim$ 10 Myr for an O6 type), $T_{\rm eff}$  rapidly drops with age as the stellar population evolves after a burst of star formation. Alternatively, if a powerful AGN exists, the line flux of highly ionized species can be significantly boosted by its narrow line region. However, we have explored possible correlations between the [O iii]88/[N ii]122  ratio and the AGN fractional contribution to the MIR and/or bolometric luminosity of LIRGs and found no significant trend with any of the diagnostics described in section 4.1, suggesting that these lines originate predominantly from star formation (except maybe in a few particular cases, where the AGN contribution could be significant).

We have used itera, the idl Tool for Emission-line Ratio Analysis (Groves et al., 2008), which is based on the Starburst99 stellar evolution synthesis models (Leitherer et al., 1999) and the shock and photo-ionization mappings iii code (Allen et al., 2008), to explore the properties of the dense ionized gas surrounding newly formed stars as a function of the main physical parameters that describe the emission line nebula, such as $n_{\rm e}$, q141414q $\equiv$ Q_{0}/4\pi R^{2}$$n_{\rm e} = $Uc$ where q and U are different definitions of the ionization parameter, $Q_{0}$ is the number of ionizing Lyman-continuum photons, R is the radius of the Strömgren sphere, and $n_{\rm e}$  is the electron density in the ionized nebula. Both q and U are proportional to $G$/$n_{\rm H}$, if a fixed stellar population and size for the star-forming region are assumed., metallicity (in units of $Z_{\odot}$) and the age of the starburst. Figure 8 shows the [O iii]88/[N ii]122  as a function of [Ne iii]15.6/[Ne ii]12.8  for the LIRGs with available Herschel and Spitzer measurements. For the ITERA models, we have assumed that the starburst has occurred instantaneously (single stellar population models), which is a likely scenario for local LIRGs, and selected the stellar atmospheres from Levesque et al. (2010) with high mass loss. We have adjusted the [O iii]88/[N ii]122  ratio predicted by the models down by 0.18 dex to roughly account for the difference in oxygen abundances between these models and the values that we use in section 5.3 to model the PDR emission of the same galaxies with the models from Kaufman et al. (1999) (O/H = 3 $\times$ 10*-4*). We note that this only affects the y-axis in Figure 8. We explore the ITERA models for a range of parameters: q = 2–40 $\times$ 107 cm s*-1*, starburst ages = 1–10 Myr, and metallicities of $Z_{\odot}$, and 2 $Z_{\odot}$.

As we can see, the models reproduce well the range of line ratios observed in our sample, although some of the LIRGs with the lowest [O iii]88/[N ii]122  seem to fall below the model grid. Note that two of these galaxies (IRASF 05189-2524 and NGC 1068) harbor an AGN with <$$\alpha_{\rm AGN}^{\rm bol}$$> $\geq$ 0.5 (flagged with crosses in Figure 8). Sources with [O iii]88/[N ii]122 $\lesssim$ 1 can only be reproduced by models with 2 $Z_{\odot}$, starburst ages of $\simeq$ 1 Myr and the lowest ionization parameter q = 2 $\times$ 107 cm s*-1*  (solid green line). We note however, that this region of the parameter space could be reached at slightly higher metallicities than 2 $Z_{\odot}$, but such models are not currently available in itera. The remaining bulk of the sample is well described by models with solar or super-solar metallicity and starburst ages ranging from $\sim$ 2 to 5 Myr, in agreement with Inami et al. (2013). The starburst age is only slightly degenerate with q, and variations in the [O iii]88/[N ii]122  ratio can be mostly attributed to a change in the ionization parameter (regardless of Z, except for ages $\lesssim$ 2 Myr). Particularly, values of [O iii]88/[N ii]122 $\gtrsim$ 5 (or lower at [Ne iii]15.6/[Ne ii]12.8 $\lesssim$ 0.5), cannot be reproduced by any model of any age unless q $>$ 4 $\times$ 107 cm s*-1*, and $\gtrsim$ 8 $\times$ 107 cm s*-1*  for [O iii]88/[N ii]122  close to $\simeq$ 10.

We note that even though higher [O iii]88/[N ii]122  can be reached at sub-solar metallicities and could explain the location of the sources in the DGS studied by Cormier et al. (2015) (filled stars in Figure 8), this is not likely the reason for the increase in the line ratio observed in dusty systems like LIRGs, but rather an increase in q. Thus, in dusty star-forming galaxies the [O iii]88/[N ii]122  ratio probes mostly the ionizing parameter q for starbursts of a few Myr (i.e., when these emission lines are actually detected), with [O iii]88/[N ii]122  scaling roughly linearly from $\simeq$ 1 to 10 for values of q ranging from 2 $\times$ 107 to 2 $\times$ 108 cm s*-1*. This relation can be refined if a measurement of the [Ne iii]15.6/[Ne ii]12.8  ratio is also available.

5.3. PDR Modeling

The [C ii]158  and [O i]63  fine-structure lines are two of the main coolants of the neutral ISM (Malhotra et al., 1997; Rosenberg et al., 2015) and therefore they are the most useful to constrain its physical conditions. We have used the PDR Toolkit (PDRT) wrapper (Pound & Wolfire, 2008) to derive the main parameters that govern the PDR emission: the intensity of the UV inter-stellar radiation field (ISRF) that impinges the PDR surface, $G$, measured in units of the local Galactic value ($G_{\rm 0}$ = 1.6 $\times$ 10*-3* erg s*-1* cm*-2*; Habing 1968), and the volume density of the neutral gas, $n_{\rm H}$. The PDRT is based on the models developed by Kaufman et al. (1999, 2006), which calculate the temperature structure and line emission from a PDR, based on the formalism for a one-dimensional semi-infinite slab (Tielens & Hollenbach, 1985). We assume the case in which the PDRs are illuminated only on one side. This choice is the most reasonable for starburst galaxies, especially LIRGs, where most of the newly-formed massive stars are likely still embedded in their parent molecular clouds (e.g., Díaz-Santos et al., 2007). That is, while a generic molecular cloud in the disk of a LIRG may be seeing an isotropic ISRF from older stellar populations, most of the line flux will be overwhelmingly dominated by the energy cooled in warm PDRs associated with the most recent burst. Strictly, these models can only be compared to individual star-forming regions. However, they can also provide useful information regarding the typical conditions of the PDRs in a given galaxy when only spatially integrated measurements are available.

We use the integrated [C ii]${}^{\rm PDR}_{158}$, [O i]63  and FIR fluxes as input to the PDRT models, excluding sources with limits in any of these quantities. Note that the [C ii]  fluxes used here are those obtained after subtracting the ionized gas contribution from the line emission151515To obtain the [C ii]${}^{\rm PDR}_{158}$  flux for galaxies without a measurement of the [N ii]205  line we use equation 4 to estimate f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}_{158}}) (capped at 90 %) from the FIR $S_{\rm 63}$/$S_{\rm 158}$  continuum ratio. The errors are fully propagated based on the dispersion of the correlation. (see section 5.1.1); a critical step needed to obtain consistent results from the PDR analysis. The data of each galaxy are compared, via $\chi^{2}$  minimization, to the ([C ii]${}^{\rm PDR}_{158}$+[O i]63)/FIR and [C ii]${}^{\rm PDR}_{158}$/[O i]63  ratios predicted by the PDR model for each combination of parameters, –0.5 $\leq$ log ($G$/$G_{\rm 0}$) $\leq$ 6.5 and 1 $\leq$ log ($n_{\rm H}$/cm*-3*) $\leq$ 7, in 0.25 dex steps. The flux uncertainties employed to calculate the $\chi^{2}$  include a 11 % error due to the absolute photometric uncertainty, added in quadrature to the measured error in the actual spectra. To obtain the map of the probability distribution for each $G$  and $n_{\rm H}$  combination averaged over entire GOALS sample we combine the $\chi_{\nu}^{2}$  maps of each galaxy using a weighting of $e^{-\chi_{\nu[G_{\rm 0},n_{\rm H}]}^{2}/2}$.

We present the results in Figure 9 separately for sub-LIRGs ($L_{\rm IR}$ $<$ 1011 $L_{\odot}$) (left panel), LIRGs (central panel) and ULIRGs (right panel). We exclude from this analysis those galaxies whose bolometric luminosity is dominated by an AGN (<$$\alpha_{\rm AGN}^{\rm bol}$$> $\geq$ 0.5), sources where there is a mismatch between the spatial location of the lines and/or continuum emission peaks, galaxies with companions within the aperture used to measure their fluxes, and sources where the [O i]63  line shows signs of self-absorption (see section 5.2.1). We do not attempt to correct for the latter effect since the models already account for the possible optical thickness of the line within the PDR, and self-absorption due to intervening cold material along the line(s) of sight is completely unconstrained. Observations at significantly higher spectral resolution than our PACS data, such as those obtained for the Milky-Way with Herschel/HIFI, would be needed in order to correct for this effect (Pineda et al., 2013; Gerin et al., 2015; Langer et al., 2016).

As we can see, the location of the main cloud of best fitting parameters (purple color, $\chi_{\nu}^{2}$ $\sim$ 1) show that most star-forming LIRGs have average ISRF intensities in the range of $\sim$ 101-3.5 $G_{\rm 0}$  and gas densities $n_{\rm H}$ $\sim$ 1–103 cm*-3*. In this regime of the parameter space of relatively low $G$  and $n_{\rm H}$  where the [C ii]  line dominates the cooling ([O i]63/[C ii]${}^{\rm PDR}_{158}$ $\lesssim$ 1), both parameters are degenerate161616This is more clear when looking at the parameter space of the $\chi_{\nu}^{2}$  map of an individual galaxy. In Figure 9 we show the averaged map, <$$\chi_{\nu}^{2}$$>, which is biased towards the bottom (best value) of the $\chi_{\nu}^{2}$  distribution for each individual source due to the weighting. Thus the overall <$$\chi_{\nu}^{2}$$> map is effectively the accumulation of the best $\chi_{\nu}^{2}$  values for the entire galaxy sample., and the relevant quantity setting the observed line ratio is $G$/$n_{\rm H}$  (see also Malhotra et al., 2001), which is approximately proportional to the temperature of the gas. This is because the numerator and the denominator of the ([C ii]${}^{\rm PDR}_{158}$+[O i]63)/FIR ratio mainly depend on the $n_{\rm H}$  and $G$, respectively. In turn, the [O i]63/[C ii]${}^{\rm PDR}_{158}$  ratio is sensitive to $G$, but also particularly to the gas temperature as we showed in section 5.2.1 ($T^{\rm kin}_{\rm gas}$ $\propto$ $T_{\rm dust}$).

Because in a one-side illuminated cloud the (optically thin) FIR intensity radiated back is $I_{\rm FIR}$ $\simeq$ 2 $\times$ $G$ /4$\pi$ [erg s*-1* cm*-2* sr*-1*] (Kaufman et al., 1999), the integrated FIR luminosity produced by a galaxy containing such typical PDR regions is:

[TABLE]

where $D$ is the angular diameter of the overall starburst, and $\phi_{\rm A}$  is the beam (area) filling factor of the PDR-emitting regions within it. We note that there is a secondary region of the parameter space, around $G$ $\simeq$ 1 $G_{\rm 0}$  and $n_{\rm H}$ $\simeq$ 104-5 cm*-3*, where low $\chi_{\nu}^{2}$  values are also found. While this range of parameters is equally probable based on the input data, such a low $G$  is rather unlikely. Even assuming $\phi_{\rm A}$ = 1 (but see section 5.4), a $G$ = 1 $G_{\rm 0}$  implies that in order to generate the total $L_{\rm FIR}$  emitted by a LIRG, a few 1011 $L_{\odot}$, the size of the starbursting region would have to be $D$ $\gtrsim$ 10 kpc (increasing as $\phi_{\rm A}^{-1/2}$). While this might be possible for a few objects, the typical sizes of LIRGs measured in the MIR and FIR are on average significantly smaller than that, with diameters of just a few kpc, or even below 1 kpc in the case of ULIRGs (Díaz-Santos et al., 2010b; Rujopakarn et al., 2011; Lutz et al., 2016; Chu et al., 2017). Thus, we only consider solutions with $G$ $>$ 100.5 $G_{\rm 0}$.

The best fitted $G$  values obtained through the PDR modeling are several orders of magnitude lower than 106 $G_{\rm 0}$, which implies that grain charging should not be a major cause of the reduced photoelectric heating efficiency of the gas (see Kaufman et al. 1999 for details). In other words, PAH molecules should remain mostly neutral on average, even in the warmest LIRGs. This is in agreement with results based on ISO data showing that the [C ii]158  to PAH6.2  (or overall 5–10 $\mu$m wavelength range) emission ratio remains roughly constant in normal star-forming galaxies and ULIRGs (Luhman et al., 2003; Helou et al., 2001, respectively). We note that in those works the ionized gas component was not subtracted from the total [C ii]158  emission when compared to the PAH6.2  (see section 5.1.1) and while AGN identification was performed, its contribution to the IR or bolometric power output of the galaxies was not taken into account. Both issues probably contribute to the large dispersion seen in the [C ii]158/PAH6.2  ratio across each sample ($\sim$ 0.2 dex). The values of our best fitted $G$  and the implication of having mostly a neutral PAH population, are also in agreement with the fact that the PAH7.7/PAH11.3  ratio –a tracer of the ionization state of the PAH molecules (Draine & Li, 2001)– does not vary significantly in LIRGs ($\lesssim$ 30%; Stierwalt et al., 2014).

Interestingly, within the context of the PDR models used here, the $G$/$n_{\rm H}$  ratios of around one third of our LIRGs (see section 5.4) are larger than the threshold ($G$/$n_{\rm H}$ $\simeq$ 2 cm3) at which the radiation pressure on the dust starts driving grains through the gas at velocities greater than the average turbulent velocity of the gas assumed by the model ($v_{\rm drift}$ $\simeq$ $v_{\rm turb}$ = 1.5 km s*-1*). This suggests that galaxies with $G$/$n_{\rm H}$  above this critical value (dotted line in Figure 9) could have PDRs that may be dynamically unstable (Kaufman et al., 1999, see also Draine 2011 for a discussion regarding dust drift velocities within H ii  regions).

The purple regions in the left panel of Figure 9 show that sub-LIRGs have $G$  values spanning the lower end of the parameter space of LIRGs (see center panel), but have similar gas densities. On the other hand, the most luminous galaxies in GOALS (ULIRGs, $L_{\rm IR}$ $\geq$ 1012 $L_{\odot}$; right panel) have the highest $G$  values, $\simeq$ 10 $\times$ the average of the sample, but also the same range of gas densities, $n_{\rm H}$ $\sim$ 1–103 cm*-3*. Therefore, there is a progression for $G$  to increase with $L_{\rm IR}$, but not an equivalent increase in the PDR gas density, as it would be expected from more compact environments associated with merger-driven star-formation. The GOALS ULIRG sample agrees reasonably well with the results obtained for the local ULIRGs in the HERUS sample studied by Farrah et al. (2013), who obtained similar $G$  values, between $\simeq$ 3 $\times$ 102-3 $G_{\rm 0}$, albeit with an upper limit for the derived gas densities, $n_{\rm H}$ $\lesssim$ 2 $\times$ 102 cm*-3*, a factor of a few lower than our estimates.

In section 5.2.1 we estimated $T^{\rm kin}_{\rm gas}$ $\simeq$ 50-100 K, significantly lower than the PDR surface temperatures favored by the models, $T_{\rm PDR}$ $\simeq$ 150–1000 K (see contours in Figure 9). A possible explanation for this discrepancy is that most of the line emission arises from a region deep in the cloud – far removed from the warm PDR surface. In this case, thermal saturation of the [C ii]158  line would occur at higher $G$/$n_{\rm H}$  than expected in the models, closer to the upper limit that we measure in the most luminous sources, around $G$/$n_{\rm H}$ $\simeq$ 10 cm3  and $T^{\rm kin}_{\rm gas}$ $\simeq$ 100 K $\approx$ $\Delta$E/$k_{\rm B}$ = 92 K. We note, however, that despite the discrepancy in absolute values there is a broad correlation between $T_{\rm PDR}$  and $T^{\rm kin}_{\rm gas}$.

5.4. Probing the ISM Structure Using $G$/$n_{\rm H}$:

A Characteristic Break in $\Sigma_{\rm IR}$

If most of the IR luminosity of a galaxy is produced in PDRs, the local ISRF, traced by $G$, and the galaxy-integrated $L_{\rm FIR}$  are tied, via equation 7, through the overall size of the starburst, D, and the area filling factor of the PDR-emitting star-forming regions within the source, $\phi_{\rm A}$171717We note that area filling factors larger than unity are possible if multiple PDRs are stacked along the same line of sight.. This means that if we know the actual size of the starburst in a galaxy, we can estimate $\phi_{\rm A}$. Figure 10 shows the ratio of $L_{\rm FIR,eff}$  to $G$  as a function of the effective physical area of our galaxies measured at 70 $\mu$m, where the starburst emission peaks. Within the uncertainties, we find that all LIRGs fall below the dashed line marking $\phi_{\rm A}$ = 1. However, $\phi_{\rm A}$  can be as small as $\simeq$ 10*-3* for some of the largest galaxies.

In Figure 11 we show the $G$/$n_{\rm H}$  ratio as a function of the luminosity surface density for our LIRG sample. We can see that there is a clear correlation between both quantities, but only above $\Sigma_{\rm IR}^{*}$ $\simeq$ 5 $\times$ 1010 L${{}_{\odot}}\,$kpc*-2*  (see also Elbaz et al., 2011). Below this threshold, an increase in $\Sigma_{\rm IR}$  is not followed by an increase of $G$/$n_{\rm H}$. This implies that while the number density of star-forming regions within LIRGs in this dynamic range is progressively larger, the average PDR properties do not vary significantly. This is clear from the color-coding of the data, which indicates the $\phi_{\rm A}$  of each galaxy. Increasing $\Sigma_{\rm IR}$  from $\sim$ 109 to $\sim$ 5 $\times$ 1010 L${{}_{\odot}}\,$kpc*-2*  can be accomplished by an equal increase in $\phi_{\rm A}$, from 10*-3* to a few 10*-2*, with $G$/$n_{\rm H}$  remaining at a relatively constant value of $\simeq$ 0.32 $G_{\rm 0}$ cm3  (see dotted line). Indeed, $G$/$n_{\rm H}$  ratios of the order of unity or lower are also found in the nuclei and extended disks of normal, nearby star-forming galaxies, like M 82 or NGC 891 (e.g., Contursi et al., 2013; Hughes et al., 2015). Above $\Sigma_{\rm IR}^{*}$, even though $\phi_{\rm A}$  still keeps increasing, the upturn in $\Sigma_{\rm IR}$  is mostly driven by a rise in $G$/$n_{\rm H}$, pointing to a change in the physical conditions of the PDRs towards more intense ISRFs for a given gas density. In fact, most of the increase in $G$/$n_{\rm H}$  at high $\Sigma_{\rm IR}$  is due to an increase in $G$, while variations in $n_{\rm H}$  are the source of scatter in the trend at a any $\Sigma_{\rm IR}$. For reference, the solid line in Figure 11 shows a linear dependence (not a fit) of the form: $G$/$n_{\rm H}$ [$G_{\rm 0}$ cm3] = 1.5 $\times$ 10*-11* $\times$ $\Sigma_{\rm IR}$ [L${{}_{\odot}}\,$kpc*-2*].

5.5. Physical Interpretation:

Young, Compact, Dusty Star-Forming Regions

It has been proposed that the deficit of PDR ([C ii], [O i]) line emission with respect to that of the IR continuum could be understood in terms of the geometry of the star-forming regions. This essentially represents the change of the ratio between the PDR the surface area and volume of the dust-emitting region – a ratio that would decrease linearly as the star-forming region grows in size. This would be further accentuated as the starburst evolves and H ii  regions begin to overlap, sharing a common PDR envelope, creating effectively a giant, single star-forming region – a “giant Orion”. This picture does not seem consistent with our results for two reasons: (1) the areal filling factor of PDR emission we measure in the majority of LIRGs is much less than unity as might be expected from a single, large H ii  region filling the central starburst, and (2) we estimate a volume filling factor of dense ionized gas of $\lesssim$ 5%, implying that it mostly originates from compact H ii  regions.

The connection between $\Sigma_{\rm IR}$  and $G$/$n_{\rm H}$  shown in Figure 11 relates a kpc-scale property of LIRGs to the physics of the individual PDRs. In other words, the fact that both quantities scale with each other above $\Sigma_{\rm IR}^{*}$  indicates that in this regime, the luminosity surface density of the entire starburst increases with the average energy per particle density at the surface of PDRs inside the few pc-scale star-forming regions. A joint, progressive increase of $\Sigma_{\rm IR}$  and $G$/$n_{\rm H}$  can be accomplished by increasing the fraction of younger, more massive stars per star-forming region181818Here we do not imply a change in the initial mass function, but rather refer to younger star-forming regions in which massive stars, with orders of magnitude higher light-to-mass ratios than solar mass stars, have still not disappeared (ages less than a few Myr; Inami et al., 2013).. Indeed, the presence of a larger number of younger stars with harder ionization fields would provide the necessary boost in $G$  and $L_{\rm IR}$  to produce the observed correlation seen in Figure 11, above $\Sigma_{\rm IR}^{*}$. Besides, the spatial distribution of gas and dust within these regions must be also different. Very young, less-developed star-forming regions have not had time to expel their dusty envelopes and thus have a larger fraction of the surrounding material closer to the newly born stars. At the same time, their radial density structure peaks at the position of the ionization front, close to the PDR, in contrast to dustless H ii  regions where the density profile is rather constant (Draine, 2011). Note that in this case the volume and column density of gas and dust may remain constant, as only a redistribution of the intervening material along the radial direction is needed. However, $T_{\rm dust}$, and thus $G$, must increase since dust grains are overall closer to the heating source.

Our results indicate that this dusty phase, in which we are able to see a significant number of star-forming regions that are likely matter-bounded and still embedded in their molecular cocoon (e.g., Kawamura et al., 2009; Miura et al., 2012; Whitmore et al., 2014), seems to be associated primarily with the period when the overall starburst region is also most compact (above $\Sigma_{\rm IR}^{*}$). A mechanism often proposed to produce this compaction of the gas and dust in galaxies are major mergers (Sanders et al., 1988; Hopkins et al., 2008). While there is not a clear sequence where galaxies with increasingly larger $\Sigma_{\rm IR}$  and $G$/$n_{\rm H}$  are found to be in progressively more advanced interacting stages, we find that most ($\sim$ 80 %) late-stage mergers (Stierwalt et al., 2013) have $\Sigma_{\rm IR}$ $>$ 5 $\times$ 1010 L${{}_{\odot}}\,$kpc*-2*, and lie along the slope in Figure 11. This lack of a simple correlation with merger stage is consistent with the idea that star-forming clumps and stellar clusters can form and be quickly destroyed by powerful stellar feedback in gas-rich starbursts and mergers on relatively short (10–20 Myr) timescales (Whitmore et al., 2014; Oklopčić et al., 2017; Linden et al., 2017), significantly shorter than the timescale of the entire merging process (over several hundred Myr). Thus, rather than a single or small number of giant, long-lived, and relatively normal H ii  regions, it may be more accurate to think of the nuclear starburst in a LIRG as a collection of dense, dusty and young star-forming regions that are individually short-lived but constantly replenished, like “fireworks” triggered by a galactic merger.

6. Conclusions

We obtained Herschel/PACS spectroscopy of the main FIR fine-structure emission lines for a sample of $\sim$ 240 galaxies in GOALS, a complete flux-limited and luminosity-limited sample of all 60$\,\mu$m-selected LIRGs systems detected in the nearby universe (z $<$ 0.09). We combined these observations with Herschel/SPIRE and MIR Spitzer/IRS spectroscopic data to derive the main physical characteristics of the neutral and ionized gas in dense PDRs and H ii  regions and provide a comprehensive view of the heating and cooling of the ISM as a function of galaxy-integrated properties. We have found the following results:

• •

The FIR lines with the most pronounced “deficits” (greatest decline of the line flux with respect to FIR continuum emission) as a function of $S_{\rm 63}$/$S_{\rm 158}$, or equivalently $T_{\rm dust}$, are the [C ii]158, [N ii]122  and [N ii]205  emission lines, which show a decrease of a factor $\sim$ 50 in the warmest LIRGs. The [O i]63  line displays a smaller deficit (a factor of $\sim$ 10). The [O iii]88  line shows no deficit with $T_{\rm dust}$  but rather an increasing scatter at high $S_{\rm 63}$/$S_{\rm 158}$. However, all lines show an increasing deficit with rising infrared luminosity surface density, $\Sigma_{\rm IR}$.

• •

We use the [C ii]158/[N ii]205  ratio to derive the contribution of photo-dissociation regions (PDRs) to the total [C ii]158  emission, [C ii]${}^{\rm PDR}_{158}$, as well as its fractional contribution, f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}{158}}) = [C ii]${}^{\rm PDR}_{158}$/[C ii]158. We find that f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}{158}}) broadly increases as a function of $S_{\rm 63}$/$S_{\rm 158}$  from $\sim$ 60 %, a value typical of Milky Way-like, to up to 95 % in the warmest LIRGs. f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}_{158}}) also increases as a function of the hardness of the ionizing radiation field as traced by the [O iii]88/[N ii]122  ratio.

• •

We find that the scatter in the [C ii]158-to-FIR  ratio as a function of $S_{\rm 63}$/$S_{\rm 158}$  is correlated with f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}{158}}). LIRGs with higher f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}{158}}) show progressively smaller [C ii]158  deficits at a given $T_{\rm dust}$, indicating that while warmer LIRGs do in fact have larger [C ii]158  deficits than normal star-forming galaxies, those with higher f(\hbox{[C\,{\sc ii}]{}^{\rm PDR}_{158}}) show less pronounced ones.

• •

Most LIRGs show [N ii]122/[N ii]205  ratios compatible with electron densities of the ionized gas in the range $n_{\rm e}$ $\simeq$ 101-2cm*-3*, with a median of 41 cm*-3*. Given the higher $n_{\rm e}$  values obtained from the MIR sulfur line ratios, this suggests that most of the [N ii]205  emission (and [C ii]${}^{\rm ion}_{158}$) is thermalized in the H ii  regions and thus mostly originate from the diffuse ionized gas phase. Assuming a simple, two-component ISM model with dense ($n_{\rm e}$ $\simeq$ 500 cm*-3*) and diffuse ($n_{\rm e}$ $\lesssim$ 1 cm*-3*) ionized gas regions, the measured [N ii]122/[N ii]205  ratio implies an average volume filling factor $\phi_{\rm V}$ $\leq$ 5 % for dense H ii  regions with respect to the overall volume of ionized emitting gas in the majority of LIRGs.

• •

We model the [O i]63/[C ii]${}^{\rm PDR}_{158}$  ratio and find that it is well correlated with the kinetic temperature of the gas in PDRs when the [O i]63  emission is optically thin. Galaxies with [O i]63/[C ii]${}^{\rm PDR}_{158}$  ratios lower than expected by the model tend to display strong absorption of the 9.7 $\mu$m silicate feature, indicating large dust opacities and suggesting that [O i]63  may be optically thick or self-absorbed in about 5 % of the sample.

• •

We find that the [O iii]88/[N ii]122  ratio is a good tracer of the ionization parameter, q, which ranges from 2 $\times$ 107 to 2 $\times$ 108 cm s*-1*  proportionally with [O iii]88/[N ii]122 $\sim$ 1–10. Moreover, most LIRGs are compatible with solar or super-solar metallicities and starburst ages $\simeq$ 2–5 Myr old, in agreement with previous results.

• •

We derive the intensity of the UV radiation field heating the PDRs, $G$, and the volume density of their neutral gas, $n_{\rm H}$. We find values in the range $G$ $\sim$ 101-3.5 $G_{\rm 0}$, increasing with $L_{\rm IR}$, and $n_{\rm H}$ $\sim$ 1–103cm*-3*  regardless of the luminosity of galaxies. The $G$/$n_{\rm H}$, ratio ranges between $\simeq$ 0.1 and 50 $G_{\rm 0}$ cm3  for the entire sample, with ULIRGs populating the upper end of the distribution at $G$/$n_{\rm H}$ $\gtrsim$ 2 $G_{\rm 0}$ cm3. We use the scaling relation between $G$  and $L_{\rm IR}$  to estimate the area filling factor of the PDRs within the overall area of the starburst region, $\phi_{\rm A}$, and find $\phi_{\rm A}$ $\simeq$ 1–10*-3*.

• •

The $G$/$n_{\rm H}$  ratio shows a distinct transition at $\Sigma_{\rm IR}^{*}$ $\simeq$ 5 $\times$ 1010 L${{}_{\odot}}\,$kpc*-2*, with the ratio being constant below this value, and increasing linearly with $\Sigma_{\rm IR}$  above it. We suggest that below this threshold, the increase in $\Sigma_{\rm IR}$  is purely driven by an increase of the number density of star-forming regions in the galaxy ($\propto$ $\phi_{\rm A}$), without changing the average PDR physical conditions. Above $\Sigma_{\rm IR}^{*}$, $G$/$n_{\rm H}$  also starts to increase, signaling a departure from the typical PDR properties of normal star-forming galaxies towards more intense/harder radiation fields and compact geometries typical of starbursting sources.

We note that the correlation between $\Sigma_{\rm IR}$  and $G$/$n_{\rm H}$  links a macroscopic property of galaxies with the local, small-scale physics of their PDR regions, and indicates that there is a threshold above which the average radiation field intensity per particle density in individual star-forming regions scales linearly with the overall energy density of the entire starburst region. We discuss the possibility that the observed correlation above $\Sigma_{\rm IR}^{*}$  arises not only from a boost in the radiation field intensity provided by young massive stars but also due to finding these star-forming regions in an earlier, dustier phase of their evolution. While we do not find a clear trend between the position of galaxies along the $G$/$n_{\rm H}$  vs. $\Sigma_{\rm IR}$  correlation and the interaction stage of the systems, we find that most late-stage mergers ($\sim$ 80%) are above $\Sigma_{\rm IR}^{*}$  and show enhanced $G$/$n_{\rm H}$  ratios, probably reflecting the fast, “fireworks”-like process of star-cluster formation and disruption during the merger.

Acknowledgments

We thank the referee for her/his useful comments. T.D-S. would like to thank M. Wolfire, J. Pineda, C. Ferkinhoff, D. Brisbin and N. Scoville for stimulating discussions about PDR physics and models. T.D-S. acknowledges support from ALMA-CONICYT project 31130005 and FONDECYT regular project 1151239. G.C.P. was supported by a FONDECYT Postdoctoral Fellowship (No. 3150361). N.L. acknowledges support from the NSFC grant No. 11673028. K.I. acknowledges support by the Spanish MINECO under grant AYA2016-76012-C3-1-P and MDM-2014-0369 of ICCUB (Unidad de Excelencia ’María de Maeztu’). This work was carried out in part at the Jet Propulsion Laboratory, which is operated for NASA by the California Institute of Technology. This work is based on observations made with the Herschel Space Observatory, an European Space Agency Cornerstone Mission with science instruments provided by European-led Principal Investigator consortia and significant participation from NASA. The Spitzer Space Telescope is operated by the Jet Propulsion Laboratory, California Institute of Technology, under NASA contract 1407. 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, and of NASA’s Astrophysics Data System (ADS) abstract service. T.D-S. wants to thank the NASA Herschel Science Center (NHSC), and in particular D. Shupe, for providing access to their computer cluster with which most of the Herschel/PACS data-sets were processed. Part of this work was carried out at the Aspen Center for Physics, which is supported by the National Science Foundation grant PHY-1066293.