<figure><img src="image_1.png" /></figure>

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

T. Díaz-Santos¹, L. Armus², V. Charmandaris³,⁴, N. Lu⁵,⁶, S. Stierwalt⁷,⁸, G. Stacey⁹, S. Malhotra¹⁰, P. P. van der Werf¹¹, J. H. Howell², G. C. Privon¹²,¹³, J. M. Mazzarella¹⁴, P. F. Goldsmith¹⁵, E. J. Murphy⁸, L. Barcos-Muñoz⁸,¹⁶, S. T. Linden⁷,⁸, H. Inami¹⁷, K. L. Larson², A. S. Evans⁷,⁸, P. Appleton¹⁴,¹⁸, K. Iwasawa¹⁹,²⁰, S. Lord²¹, D. B. Sanders²², and J. A. Surace²

¹ Núcleo de Astronomía de la Facultad de Ingeniería, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile; tanio.diaz@mail.udp.cl

² Infrared Processing and Analysis Center, MC 314-6, Caltech, 1200 E. California Blvd., Pasadena, CA 91125, USA

³ Institute for Astronomy, Astrophysics, Space Applications & Remote Sensing, National Observatory of Athens, GR-15236, Penteli, Greece

⁴ University of Crete, Department of Physics, GR-71003, Heraklion, Greece

⁵ China-Chile Joint Center for Astronomy (CCJCA), Camino El Observatorio 1515, Las Condes, Santiago, Chile

⁶ National Astronomical Observatories, Chinese Academy of Sciences (CAS), Beijing 100012, China

⁷ Department of Astronomy, University of Virginia, P.O. Box 400325, Charlottesville, VA 22904, USA

⁸ National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

⁹ Department of Astronomy, Cornell University, Ithaca, NY 14853, USA

¹⁰ Astrophysics Science Division, Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA

¹¹ Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands

¹² Departamento de Astronomía, Universidad de Concepción, Casilla 160-C, Concepción, Chile

¹³ Instituto de Astrofísica, Facultad de Física, Pontificia Universidad Católica de Chile, Casilla 306, Santiago 22, Chile

¹⁴ Infrared Processing and Analysis Center, MC 100-22, Caltech, 1200 E. California Blvd., Pasadena, CA 91125, USA

¹⁵ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

¹⁶ Joint ALMA Observatory, Alonso de Córdova 3107, Vitacura, Santiago, Chile

¹⁷ Centre de Recherche Astrophysique de Lyon, Universite de Lyon, Universite Lyon 1, CNRS, Observatoire de Lyon, 9 avenue Charles Andre, Saint-Genis Laval Cedex F-69561, France

¹⁸ NASA Herschel Science Center, IPAC, California Institute of Technology, MS 100-22, Cech, Pasadena, CA 91125, USA

¹⁹ Institut de Cincies del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Marti i Franques 1, E-08028 Barcelona, Spain

²⁰ ICREA, Pg. Lluís Companys 23, E-08010 Barcelona, Spain

²¹ The SETI Institute, 189 Bernardo Avenue, Suite 100, Mountain View, CA 94043, USA

²² Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

Received 2017 February 9; revised 2017 July 16; accepted 2017 July 17; published 2017 August 29

## Abstract

We present an analysis of [O I]$_{63}$, [O III]$_{88}$, [N II]$_{122}$, and [C II]$_{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. We find pronounced declines ("deficits") of line-to-FIR continuum emission for [N II]$_{122}$, [O I]$_{63}$, and [C II]$_{158}$ as a function of FIR color and infrared luminosity surface density, $\Sigma_{IR}$. The median electron density of the ionized gas in LIRGs, based on the [N II]$_{122}$/[N II]$_{205}$ ratio, is $n_e = 41 \text{ cm}^{-3}$. We find that the dispersion in the [C II]$_{158}$ deficit of LIRGs is attributed to a varying fractional contribution of photodissociation regions (PDRs) to the observed [C II]$_{158}$ emission, $f([C II]_{158}^{\text{PDR}}) = [C II]_{158}^{\text{PDR}}/[C II]_{158}$, which increases from ~60% to ~95% in the warmest LIRGs. The [O I]$_{63}$/$[C II]_{158}^{\text{PDR}}$ ratio is tightly correlated with the PDR gas kinetic temperature in sources where [O I]$_{63}$ is not optically thick or self-absorbed. For each galaxy, we derive the average PDR hydrogen density, $n_H$, and intensity of the interstellar radiation field, $G$, in units of $G_0$ and find $G/n_H$ ratios of ~0.1–50 $G_0 \text{ cm}^{-3}$, with ULIRGs populating the upper end of the distribution. There is a relation between $G/n_H$ and $\Sigma_{IR}$, showing a critical break at $\Sigma_{IR}^* \approx 5 \times 10^{10} L_{\odot} \text{ kpc}^{-2}$. Below $\Sigma_{IR}^*$, $G/n_H$ remains constant, $\approx 0.32 G_0 \text{ cm}^{-3}$, and variations in $\Sigma_{IR}$ are driven by the number density of star-forming regions within a galaxy, with no change in their PDR properties. Above $\Sigma_{IR}^*$, $G/n_H$ increases rapidly with $\Sigma_{IR}$, signaling a departure from the typical PDR conditions found in normal star-forming galaxies toward more intense/harder radiation fields and compact geometries typical of starbursting sources.

**Key words:** galaxies: evolution – galaxies: ISM – galaxies: nuclei – galaxies: starburst – infrared: galaxies

**Supporting material:** machine-readable tables

## 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 extragalactic 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) reopened a new window into the FIR universe, providing spectroscopic data for significantly larger samples of nearby, intermediate-redshift, and high-redshift galaxies (e.g., Graciá-Carpio et al. 2011; Díaz-Santos et al. 2013, 2014; Magdis et al. 2014; Rigopoulou 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 ~50–200 µm wavelength range are [N III] 57.3 µm ([N III]$_{57}$), [O I] 63.2 µm ([O I]$_{63}$), [O III] 88.4 µm ([O III]$_{88}$), [N II] 121.9 µm ([N II]$_{122}$), [O I] 145.5 µm ([O I]$_{145}$), [C II] 157.7 µm ([C II]$_{158}$), and [N II] 205.2 µ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 photodissociation 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, in both the diffuse medium and dense (H II) regions.

Of particular importance among star-forming galaxies are the so-called luminous IR galaxies (LIRGs; $L_{IR} = 10^{11-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<sub>☉</sub>), as a consequence of the funneling of large amounts of gas and dust toward their nuclei due to the loss of angular momentum during the dynamical interaction. And while the presence of active galactic nuclei (AGNs) 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 ultraluminous IR galaxies (ULIRGs; $L_{IR} \ge 10^{12} 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 in the fact that they cover the entire transition between main-sequence (MS) galaxies and starbursts (Díaz-Santos et al. 2013).

At high redshift, LIRGs dominate the obscured star formation activity between $z \sim 1$ and 3 (e.g., Berta et al. 2011; Magnelli et al. 2011; Murphy 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 AGNs 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 Spitzer 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 Section 2 we present the LIRG sample and the new Herschel spectroscopy, as well as the ancillary observations used in this work. In Section 3 we describe the processing and analysis of the data. The basic results are presented in Section 4. The ISM properties derived for the LIRG sample are discussed throughout Section 5 in relation to specific emission-line ratios. In particular, in Sections 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 Section 5.2.1 we present a link between the [O I]$_{63}$/[C II]$_{158}^{PDR}$ and the kinetic temperature of the PDR gas. We explore PDR covering factors and metallicity variations in Section 5.2.2 and confront FIR emission line ratios involving oxygen and nitrogen emission lines to photoionization models of H II regions in Section 5.2.3. We derive the average PDR properties for each galaxy in the sample in Section 5.3 and show the existence of a critical break in the PDR conditions as a function of luminosity surface density in Section 5.4. We discuss the physical implications of these results in Section 5.5. The summary of the results and conclusions are given in Section 6.

## 2. Sample and Observations

### 2.1. The GOALS Sample

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°$. There are 180 LIRGs and 22 ULIRGs in GOALS, and their median redshift is $z = 0.0215$ (or ~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), and Chandra X-ray (Iwasawa et al. 2011), as well as Spitzer/IRS (Díaz-Santos et al. 2010b, 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, 2016b; Lu et al. 2014, 2015). 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 2013; Murphy et al. 2013; Xu et al. 2014, 2015; Zhao et al. 2016a, among others).

The RBGS, and therefore the GOALS sample, was 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 µ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 and FIR luminosities of individual galaxies ($L_{IR[8–1000micron]}$, as defined in Sanders & Mirabel 1996$^{23}$, and $L_{FIR[42.5–122.5micron]}$, as defined in Helou et al. 1985$^{23}$) is described at the end of Section 3.2.

## 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 (IRASF 08339+6517 and IRASF 09111–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 for more than 200 hr 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 µm (third and second orders, respectively; “blue” camera) and the 102–210 µm (first order; “red” camera) ranges. The IFU is composed of a 5 × 5 array of individual detectors (spaxels) each with a field of view (FOV) of ~9′′4 on a side, for a total of 47′′ × 47′′. The physical size of the PACS FOV at the median distance of our LIRG sample is ~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 × 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 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 primary 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 µm; LL module: 14–38 µm), and medium-resolution, $R \sim 600$ (SH module: 9.9–19.6 µm; LH module: 18.7–37.2 µm), slit spectroscopy. The IRS spectra were extracted with the Spitzer IRS Custom Extraction (SPICE) software,$^{24}$ using the standard extraction aperture and a point-source calibration mode. The projected angular sizes of the apertures on the sky are 3′′7 × 12′′ at 9.8 µm in SL, 10′′6 × 35′′ at 26 µm in LL, 4′′7 × 15′′5 at 14.8 µm in SH, and 11′′1 × 36′′6 at 28 µm in LH. Thus, the area covered by the SL and SH, and the LL and LH apertures is approximately equivalent (within a factor of ~2) to that of an individual spaxel and a 3 × 3 spaxel box of the PACS/IFS, respectively. Aside from the line fluxes, which are available in Stierwalt et al. (2014), other observables measured from the IRS data used in this work are the strength of the 9.7 µm silicate feature,$^{25}$ $S_{9.7µm}$, and the equivalent width (EW) of the 6.2 µm PAH, both of which were provided 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 (HSA$^{26}$), 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, calculate 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 ~82 km s$^{-1}$ for [O I]$_{63}$, ~120 km s$^{-1}$ for [O III]$_{88}$, ~287 km s$^{-1}$ for [N II]$_{122}$, and ~235 km s$^{-1}$ for [C II]$_{158}$. The final processing steps were as follows: 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.

<footnote>$^{23}$ $F_{IR[8–1000µm]} = 1.8 \times 10^{-14}$ (13.48 $S_{12µm} + 5.16$ $S_{25µm} + 2.58$ $S_{60µm} + S_{100µm}$) [W m$^{-2}$], with $S_ν$ in Jy. $L_{IR} = 4π D_L^2 F_{IR}$. The luminosity distances, $D_L$, are taken from Armus et al. (2009).</footnote>

<footnote>$^{24}$ http://irsa.ipac.caltech.edu/data/SPITZER/docs/dataanalystools/tools/spice/</footnote>

<footnote>$^{25}$ The silicate strength is defined as $S_{9.7µm} = \ln(f_{peak}/f_{cont})$, where $f_{peak}$ is the flux density at the peak absorption (or emission) close to 9.7 µm, and $f_{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 that the feature appears in absorption, while positive values indicate that it is emission.</footnote>

<footnote>$^{26}$ http://herschel.esac.esa.int/Science_Archive.shtml</footnote>
## 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 ~11%–12%,<fnref n="27" /> 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).<fnref n="28" />

We extracted the line fluxes for our LIRGs from the spectra using three different apertures: (a) the spaxel at which the 63 µm continuum emission of the galaxy peaks (hereafter referred to as the “central” spaxel); (b) in a 3 × 3 spaxel box centered on the central spaxel as defined in (a), limited by the PACS FOV; and (c) the total FOV (5 × 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 method (a), as well as the “best” spatially integrated value for each individual galaxy. The latter is defined as the highest flux value that maximizes the signal-to-noise ratio of the measurements (line and continuum). In other words, we select, among the three methods described above, the one that provides the lowest noise in the measurements while still accounting for the entire line and continuum fluxes 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 ≲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 ≲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 µm (≈17′′), to which we applied an aperture correction based on the

<table><caption>**Table 1**
Table Content of Emission-line and Continuum Data</caption><thead><tr><th>Column No.</th><th>Quantity</th><th>Units</th></tr></thead><tbody><tr><td>1</td><td>ID</td><td>...</td></tr><tr><td>2</td><td>Galaxy name</td><td>...</td></tr><tr><td>3, 4</td><td>Raster (x, y)</td><td>pixels</td></tr><tr><td>5, 6</td><td>Central spaxel (x, y)</td><td>pixels</td></tr><tr><td>7, 8</td><td>R.A. and decl. of (4, 5)</td><td>...</td></tr><tr><td>9</td><td>Dist. to [O I]<sub>63</sub> cont.</td><td>arcsec</td></tr><tr><td>10</td><td>Central line flux</td><td>×10<sup>-17</sup> W m<sup>-2</sup></td></tr><tr><td>11</td><td>Central cont. flux dens.</td><td>Jy</td></tr><tr><td>12</td><td>Companion galaxy?</td><td>...</td></tr><tr><td>13, 14, 15</td><td>R-C, R-L, C-L</td><td>...</td></tr><tr><td>16</td><td>Best line flux</td><td>×10<sup>-17</sup> W m<sup>-2</sup></td></tr><tr><td>17</td><td>Best cont. flux dens.</td><td>Jy</td></tr><tr><td>18</td><td>Best measurement</td><td>...</td></tr><tr><td>19</td><td>AOR</td><td>...</td></tr><tr><td>20</td><td>Program</td><td>...</td></tr></tbody></table>

**Note.** The data table for each FIR emission line is available in the electronic edition of this paper. The columns include the line and continuum fluxes observed with Herschel/PACS, as well as a number of measurements: (1) Identification number. (2) Name of the galaxy. (3, 4) The raster used to obtain the galaxy measurements. (5, 6) Reference central spaxel, defined as the closest spaxel to the 63 µm continuum peak of the galaxy within the 5 × 5 PACS FOV in the raster (cols. (3) and (4)) of the AOR (col. (19)). (7, 8) Right ascension and declination of the reference spaxel. (9) Distance from the reference spaxel to the 63 µm continuum peak. (10) Flux and uncertainty of the line measured from the reference spaxel (method (a)). Uncertainties in cols. (10), (11), (16), and (17) do not include the absolute photometric calibration uncertainty; negative values in cols. (10)–(11) indicate upper limits. (11) Flux density and uncertainty of the continuum under the line measured in the reference spaxel (method (a)). (12) Does a companion galaxy exist within the reference spaxel? 1 = yes; 0 = no. (13, 14, 15) Do the spaxels of the reference (R), continuum (C), and line (L) emission peaks coincide among each other? 1 = yes; 0 = no. (16) Galaxy-integrated flux and uncertainty of the line measured from the best aperture (see col. (18)). Negative values in cols. (16)–(17) indicate upper limits. (17) Galaxy-integrated flux density and uncertainty of the continuum under the line measured in the best aperture (see col. (18)). (18) Best measurement type: methods (a), (b), or (c) (see text for details). (19) AOR ID of the data set used. (20) Program ID of the AOR. The IR luminosities and luminosity surface densities of the galaxies, $L_{IR}$ and $\Sigma_{IR}$ (the latter defined as $(L_{IR}/2)/\pi R_{70 \mu m, eff}^2$), can be found at http://goals.ipac.caltech.edu.

(This table is available in its entirety in machine-readable form.)

PACS spectrometer beam efficiency maps (v6) provided on the PACS instrument and calibration Web pages,<fnref n="29" /> after rebinning them to the regular 5 × 5 spaxel FOV. We note that we do not apply the correction factors to the [N II]205 µm emission based on the FIR color of galaxies provided in Lu et al. (2017). We use instead the original, point-source calibrated fluxes since we do not use those corrections in our Herschel/PACS data.

In order to obtain the $L_{IR}$ and $L_{FIR}$ of individual galaxies belonging to a LIRG system formed by two or more components, for the different extraction apertures described above, we scaled the integrated IRAS IR and FIR luminosities of the system with the ratio<fnref n="30" /> between the continuum flux density of each individual galaxy evaluated at 63 µm in the PACS spectrum (measured in the same aperture as the line) and the total IRAS 60 µm flux density.

<footnote n="27">²⁷ http://herschel.esac.esa.int/Docs/PACS/html/ch04s10.html</footnote>

<footnote n="28">²⁸ http://herschel.esac.esa.int/twiki/bin/view/Public/PacsCalTreeHistory</footnote>

<footnote n="29">²⁹ http://herschel.esac.esa.int/twiki/bin/view/Public/PacsCalibrationWeb</footnote>

<footnote n="30">³⁰ $F_{FIR}$ [42.5–122.5 µm] = 1.26 × 10⁻¹⁴ (2.58 $S_{60µm}$ + $S_{100µm}$) [W m⁻²], with $S_ν$ in [Jy]. $L_{FIR} = 4π D_L^2 F_{FIR}$. The luminosity distances, $D_L$, are taken from Armus et al. (2009).</footnote>
# 4. Results

## 4.1. MIR and Bolometric AGN Contributions

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 [NeV]14.3/[NeI]12.8 and [OIV]25.9/[NeII]12.8 emission line ratios, the 6.2 µm PAH EW (see also Armus et al. 2007), the S30/S15 dust continuum slope<fnref n="31" />, 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—through the ionization of the surrounding gas to higher states and/or via the heating of dust to higher temperatures than star formation—and provide an estimate of its fractional contribution to the MIR emission, αMIRAGN. Once these values are known, corrections based on spectral energy distribution (SED) templates of pure starbursts and AGN-powered sources can be applied to derive the fractional contribution of AGNs to the bolometric luminosity of galaxies, αbolAGN (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, $\langle \alpha_{AGN}^{MIR} \rangle$ and $\langle \alpha_{AGN}^{bol} \rangle$. Because there are a significant number of nondetections in one or more of the relevant MIR features employed to calculate $\langle \alpha_{AGN}^{MIR} \rangle$ or $\langle \alpha_{AGN}^{bol} \rangle$, 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 αMIRAGN values for each MIR diagnostic, as well as $\langle \alpha_{AGN}^{MIR} \rangle$ and $\langle \alpha_{AGN}^{bol} \rangle$.

## 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 deficits<fnref n="32" />— expressed as the line-to-FIR continuum flux ratio—for the entire GOALS LIRG sample as a function of the Sν63 µm/Sν158 µm (S63/S158) continuum flux density ratio, which is a first-order tracer of the average dust temperature in galaxies, Tdust. 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:

$$L_{\text{line}}/L_{\text{FIR}} = \epsilon_0 e^{-(S_{63}/S_{158})/\delta}, \quad (1)$$

where $\epsilon_0$ denotes a limiting line/FIR ratio for sources with cold FIR colors and no deficits, and $\delta$ is the S63/S158 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 fits to the data are displayed in Figure 1 as solid black lines, with the associated 1$\sigma$ dispersion around the fit shown as dotted lines. The derived parameters can be found in Table 3. 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 Tdust 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 Tdust, 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 S63/S158 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 blackbody (mBB) with a fixed emissivity index $\beta = 1.8$ and whose emission is optically thin. This is a reasonable approximation for SED fits that do not include data at $\lambda \lesssim 60$ µm. We also provide a practical equation that relates both quantities, S63/S158 and Tdust, using the following approximation:

$$T_{\text{dust}} = 20.24 + 14.54(S_{63}/S_{158}) - 3.75(S_{63}/S_{158})^2 \\ + 0.46(S_{63}/S_{158})^3. \quad (2)$$

We note that this equation is only valid for the dust temperatures and S63/S158 ratios spanned by the galaxies in the GOALS sample, that is, $21 \text{ K} \le T_{\text{dust}} \le 48 \text{ K}$ or $0.1 \le S_{63}/S_{158} \le 4$. The error in Tdust obtained from this expression is $\lesssim 0.5 \text{ 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 Tdust 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 S63/S158 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 I]63/FIR and [O III]88/FIR ratios as a function of Tdust may be related to the location where the line emission and dust emission originate, respectively, 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 also originate from the PDR itself or the outer edge of the H II region ([N II]122,

<footnote n="31">³¹ We modified the reference value for the S30/S15 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 S30/S15 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 µm feature, to be used with the 6.2 µm PAH, and we modify the reference value for pure starbursts and its bolometric correction factor accordingly. We further assume that there is no PAH destruction due to the AGN when calculating its fractional contribution using this method (Díaz-Santos et al. 2008, 2010a; Esquej et al. 2014).</footnote>

<footnote n="32">³² The 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.</footnote>
<table><caption>**Table 2**
Table Content of AGN Fractions</caption><thead><tr><th rowspan="3">ID<br>(1)</th><th rowspan="2">[Ne V]<sub>14.3</sub>/<br>[Ne II]<sub>12.8</sub><br>MIR<br>(2)</th><th rowspan="2">[O IV]<sub>25.9</sub>/<br>[Ne II]<sub>12.8</sub><br>MIR<br>(3)</th><th rowspan="2">6.2 µm PAH<br>EW<br>MIR<br>(4)</th><th rowspan="2">S<sub>30</sub>/S<sub>15</sub><br>Ratio<br>MIR<br>(5)</th><th rowspan="2">Laurent<br>Diagram<br>MIR<br>(6)</th><th colspan="2">⟨α<sub>AGN</sub>⟩</th></tr><tr><th colspan="2">KM Estimator</th></tr><tr><th></th><th></th><th></th><th></th><th></th><th></th><th>MIR<br>(7)</th><th>Bol.<br>(8)</th></tr></thead><tbody><tr><td>0</td><td>-0.01</td><td>0.01 ± 0.01</td><td>0.23 ± 0.01</td><td>0.38 ± 0.01</td><td>0.60 ± 0.05</td><td>0.25 ± 0.10</td><td>0.09 ± 0.04</td></tr><tr><td>1</td><td>-0.02</td><td>-0.03</td><td>0.40 ± 0.03</td><td>0.01 ± 0.01</td><td>0.58 ± 0.05</td><td>0.20 ± 0.11</td><td>0.04 ± 0.02</td></tr><tr><td>2</td><td>-0.01</td><td>0.01 ± 0.01</td><td>0.04 ± 0.01</td><td>0.39 ± 0.02</td><td>0.35 ± 0.08</td><td>0.16 ± 0.08</td><td>0.07 ± 0.05</td></tr><tr><td>3</td><td>-0.03</td><td>0.02 ± 0.01</td><td>0.81 ± 0.01</td><td>0.63 ± 0.01</td><td>0.72 ± 0.04</td><td>0.44 ± 0.16</td><td>0.23 ± 0.09</td></tr><tr><td>4</td><td>0.01 ± 0.01</td><td>0.03 ± 0.01</td><td>0.27 ± 0.01</td><td>0.32 ± 0.01</td><td>0.57 ± 0.05</td><td>0.24 ± 0.09</td><td>0.09 ± 0.03</td></tr></tbody></table>

Note. The columns include the fractional AGN contributions to the emission of LIRGs and the associated uncertainties, based on different *Spitzer* / IRS diagnostics. (1) Identification number to match with Table 1. (2–6) MIR AGN fractions based on each individual diagnostic (see text for details). Negative values indicate upper limits in the case of the [Ne V]14.3/[Ne II]12.8 and [O IV]25.9/[Ne II]12.8 ratios (cols. (2) and (3)), and lower limits in the case of the 6.2 µm PAH EW, S30/S15 ratio, and Laurent diagram diagnostics (cols. (4)–(6)). (7, 8) Average MIR and bolometric AGN fractions based on all diagnostics. The bolometric fractions based on individual diagnostics can be found at http://goals.ipac.caltech.edu.

(This table is available in its entirety in machine-readable form.)

[N II]205, and [C II]158), would show less scatter around the deficit's trend. 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, especially at high S63/S158 (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). However, 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. An additional source of scatter could arise from the fact that the most embedded regions in ULIRGs may be optically thick in both gas and dust (Sakamoto et al. 2008; Fischer et al. 2014; González-Alfonso et al. 2015). In some of these sources, the integrated line deficits could instead reflect the conditions of the surrounding, less obscured medium, thus disguising them as normal galaxies with large line-to-FIR ratios.

Figure 2 presents the same line deficits but as a function of the luminosity surface density, $\Sigma_{IR}$, defined as the ratio of effective luminosity, $L_{IR,eff} = L_{IR}/2$, divided by the effective area (containing half of the source's luminosity) measured at 70 µm, $\pi R_{70/\mu m,eff}^2$, for galaxies with available FIR size measurements. These correlations are overall tighter than with $T_{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 Sections 5.4 and 5.5). 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_{IR}$. Moreover, the trend is followed by nearly all LIRGs regardless of their FIR color, $T_{dust}$, or AGN contribution. We have fitted these correlations with a second-order polynomial function:

$$ \log(L_{\text{line}}/L_{\text{IR}}) = \alpha_0 + \alpha_1 \log \Sigma_{\text{IR}} + \alpha_2 (\log \Sigma_{\text{IR}})^2. \quad (3) $$

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. [CII]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]158<sup>ion</sup>) but also in the dense, neutral ISM ([C II]158<sup>PDR</sup>). However, due to the low critical density of the transition when it is collisionally excited by free electrons or protons ($n_e^{cr,[CII]} \sim 44 \text{ 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 ($n_H^{cr,[CII]} \sim 3.0 \times 10^3 \text{ cm}^{-3}$ and $n_{H_2}^{cr,[CII]} \sim 6.1 \times 10^3 \text{ 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]158<sup>ion</sup> (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]158<sup>ion</sup> also originates, 24.38–13.6 eV. In addition, given the similar $n_e^{cr}$ and $E_{ul}/k_B$ of both transitions, the [C II]158<sup>ion</sup>/[N II]205 ratio is roughly constant and depends weakly on $n_e$, the intensity of the ionizing field, $q$, and the kinetic temperature of the gas, $T_{gas}^{kin}$.

The photoionization models presented in Oberst et al. (2006) predict a [C II]158<sup>ion</sup>/[N II]205 $\sim 3 \pm 0.5$ for a range of $n_e$ up to $\sim 10^3 \text{ cm}^{-3}$ (see Section 5.1.2). For convenience, we use this constant ratio to estimate [C II]158<sup>ion</sup> and subtract it from the total [C II]158 flux, which yields [C II]158<sup>PDR</sup> (=[C II]158–[C II]158<sup>ion</sup>). Figure 3 (top panel) shows the fraction of [C II] arising from PDRs, $f([C II]158^{PDR}) = [C II]158^{PDR}/[C II]158$, as a function of S63/S158 and $T_{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_{IR}$ or $T_{dust}$. To show that it does not have an
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Line deficits of the [O I] 63.2 µm, [O III] 88.4 µm, [N II] 121.9 µm, [C II] 157.7 µm, and [N II] 205.2 µm emission lines, defined as the line-to FIR[42.5–122.5µm] flux ratio as a function of the FIR S<sub>63</sub>/S<sub>158</sub> continuum ratio for the entire GOALS sample (open circles; see legend in the middle right panel). The dynamic range of the x- and y-axes is the same for all panels. The data shown here represent the best integrated emission values available for each galaxy and line, as described in Section 3.2. The top x-axis is the T<sub>dust</sub> of a modified blackbody (with β = 1.8) that has an S<sub>63</sub>/S<sub>158</sub> ratio equal to the values shown in the bottom x-axis. The data are color-coded as a function of the total IR luminosity, L<sub>IR[8–1000µm]</sub>, measured within the same aperture used to obtain the line and continuum emissions. Sources marked with black dots are galaxies where there is a mismatch between the location where the line and/or continuum emissions peak (see Table 1). Galaxies with companions within the aperture used to measure their integrated flux are also marked with black dots. In addition to the GOALS sample, we also show the ULIRG sample from Farrah et al. (2013) (open squares), which populates better the high-T<sub>dust</sub> regime. Galaxies with an AGN contributing more than 50% to their bolometric luminosity, 〈α<sup>bol</sup><sub>AGN</sub>〉 ≥ 0.5 (see Section 4.1), are marked with black crosses. The solid lines show fits to the data using Equation (1), with the dotted lines representing the 1σ dispersion of the data with respect to the best fit. The dispersion is calculated in the same way for all the fits performed throughout the paper. Neither galaxies marked as black dots nor AGN-dominated sources are used for the fits (here or in any other fit presented throughout the paper).</figcaption></figure>
<table><caption>**Table 3**
Best Parameters from the Line Deficit versus S<sub>63</sub>/S<sub>158</sub> Fits</caption><thead><tr><th>Line</th><th>ε<sub>0</sub><br>(× 10<sup>-3</sup>)</th><th>δ<br>(e-fold)</th><th>1σ Disp.<br>(× 10<sup>-3</sup>)</th></tr></thead><tbody><tr><td>[O I]<sub>63</sub></td><td>4.37 ± 0.93</td><td>1.15 ± 0.30</td><td>1.32</td></tr><tr><td>[O III]<sub>88</sub></td><td>0.94 ± 0.17</td><td>N/A</td><td>0.94</td></tr><tr><td>[N II]<sub>122</sub></td><td>1.27 ± 0.16</td><td>0.67 ± 0.08</td><td>0.11</td></tr><tr><td>[C II]<sub>158</sub></td><td>14.0 ± 0.9</td><td>0.68 ± 0.04</td><td>1.69</td></tr><tr><td>[N II]<sub>205</sub></td><td>1.26 ± 0.16</td><td>0.50 ± 0.04</td><td>0.07</td></tr></tbody></table>

Note. The parameters $ε_0$ and $δ$ correspond to the fits of the line deficits (line-to-FIR ratios) as a function of the S<sub>63</sub>/S<sub>158</sub> ratio for the entire GOALS sample using Equation (1). The fits are presented in Figure 1 (solid lines).

affect on $f([C \text{II}]_{158}^{\text{PDR}})$ either, we color-code Figure 3 as a function of the fractional contribution of the AGN to the MIR emission of each galaxy, based on the 6.2 µm PAH EW diagnostic (see Section 4.1). Neither sources with large MIR AGN fractions nor those with the largest bolometric contributions ($\langle\alpha_{\text{AGN}}^{\text{MIR}}\rangle \ge 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

$$f([C \text{II}]_{158}^{\text{PDR}}) = [C \text{II}]_{158}^{\text{PDR}} / [C \text{II}]_{158} = 0.82(\pm 0.01) + 0.41(\pm 0.04)\log(S_{63}/S_{158}), \quad (4)$$

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_{\text{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 ~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_{\text{dust}}$ to show larger $f([C \text{II}]_{158}^{\text{PDR}})$ (black solid line), increasing from ~60%, close to the MW value, to nearly 95%. That is, there is a larger PDR contribution to the total [C II] emission in warmer systems, indicating that, even though the [C II] line shows a larger 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([C \text{II}]_{158}^{\text{PDR}})$ with the average hardness of the radiation field seen by the dense ionized gas, as traced by the [O III]$_{88}$/[N II]$_{122}$ line ratio (see Section 5.2.3). Therefore, $f([C \text{II}]_{158}^{\text{PDR}})$ increases in environments associated with recent episodes of massive star formation. This is also consistent with H II regions in the LIRG nuclei being 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 experiencing a starburst event, or simply having more modest star formation, like the MW) have more of the [C II] 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([C \text{II}]_{158}^{\text{PDR}})$ and the electron density of the ionized gas, $n_e$, as traced by the [N II]$_{122}$/[N II]$_{205}$ ratio (see color coding in the bottom panel of Figure 3). 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 \text{II}]_{158}^{\text{PDR}}$ dominates the total [C II] emission, contribution from the ionized gas, $[C \text{II}]_{158}^{\text{ion}}$, is not negligible. This ionized component can subsequently originate from both the diffuse ($[C \text{II}]_{158}^{\text{ion,diff}}$) and dense medium ($[C \text{II}]_{158}^{\text{ion,HII}}$). 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 µm/$[S \text{III}]_{33.5}$ µm line ratio$^{33}$ and calculated the average $n_e$ for most galaxies to be typically ~100 to a few hundred cm$^{-3}$, with a median of ~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_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 $C^{2+}$ transitions to $C^+$ at $\text{h}\nu < 24.38$ eV. Considering that the electron density should at least remain constant, if not increase toward the denser PDR region ($n_e$ of a few times $10^3$ cm$^{-3}$ have been found in H II regions using optical emission lines of singly ionized sulfur; Osterbrock 1989), this means that the $n_e$ derived using the $[S \text{III}]_{18.7}/[S \text{III}]_{33.5}$ ratio likely represents a lower limit to the density of the volume from where $[C \text{II}]_{158}^{\text{ion,HII}}$ 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 \text{II}]_{158}^{\text{ion}}$ and $[N \text{II}]_{205}$ critical densities, this implies that the $[C \text{II}]_{158}^{\text{ion,HII}}$ emission has likely been thermalized, therefore suggesting that most of the $[C \text{II}]_{158}^{\text{ion}}$ and $[N \text{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_{\text{gas}}^{\text{kin}}$ is very low.

We can also use the $[N \text{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 \text{II}]_{158}^{\text{ion,HII}}$ emission. Figure 4 shows the distribution of $[N \text{II}]_{122}$/[N II]$_{205}$ ratios and $n_e$ for those galaxies with available measurements of both lines. We find densities between ~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 starburst galaxies. For instance, Zhao et al. (2016b) find $n_e = 22$ cm$^{-3}$ for a subsample of GOALS LIRGs using ISO data, and Herrera-Camus et al. (2016) find $n_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_e$ values we find based on the nitrogen line ratio are smaller than those inferred from the sulfur lines, suggesting that both $[N \text{II}]_{205}$ and $[C \text{II}]_{158}^{\text{ion}}$ have been thermalized in the dense H II regions and thus mostly

<footnote>$^{33}$ Both lines have $n_e^{\text{cr}} \gtrsim 2000$ cm$^{-3}$.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Line deficits of the [O I]$_{63}$, [O III]$_{88}$, [N II]$_{122}$, [C II]$_{158}$, and [N II]$_{205}$ emission lines, defined as the ratio of line luminosity to $L_{IR[8-1000\mu m]}$ as a function of the IR luminosity surface density, defined as $\Sigma_{IR} = (L_{IR}/2)/\pi R_{70\mu m,eff}^2$, for galaxies in GOALS with available measurements of their FIR sizes (taken from Lutz et al. 2016; open circles). Symbols are as in Figure 1. We note that while the $L_{IR}$ used in the y-axis represents the total IR luminosity of the galaxy (calculated in the same aperture as the line luminosity), the value used in the x-axis is the effective luminosity, $L_{IR,eff}$, where $L_{IR,eff} = L_{IR}/2$, since the measured sizes refer to the half-light radii of the sources. The solid lines represent a fit to the data using a second-order polynomial (see Equation (3)). The best-fit parameters are tabulated in Table 4. 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. For reference, in the panel showing the [C II]$_{158}$ deficit we also plot the best fit found by Lutz et al. (2016) using the same parameterization (dashed line; within the errors of our fit) and the best fit originally found by Díaz-Santos et al. (2013) using galaxy sizes measured in the MIR with *Spitzer* and assuming a linear log-log relation (dot-dashed line).</figcaption></figure>

originate from the diffuse ISM. In this case, if we assume an $n_e$ $\simeq 500 \text{ cm}^{-3}$ and an $n_e \lesssim 1 \text{ 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_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_V$ will vary significantly depending on the compactness of the galaxy.
<table><caption>**Table 4**
Best Parameters from the Line Deficit versus $\Sigma_{IR}$ Fits</caption><thead><tr><th>Line</th><th>α<sub>0</sub></th><th>α<sub>1</sub></th><th>α<sub>2</sub></th><th>1σ Disp. (dex)</th><th>Max. Ratio (×10<sup>-3</sup>)</th></tr></thead><tbody><tr><td>[O I]<sub>63</sub></td><td>-15.83</td><td>2.711</td><td>-0.1412</td><td>0.24</td><td>1.52</td></tr><tr><td>[O III]<sub>88</sub></td><td>-8.78</td><td>1.272</td><td>-0.0721</td><td>0.32</td><td>0.69</td></tr><tr><td>[N II]<sub>122</sub></td><td>-7.66</td><td>1.023</td><td>-0.0613</td><td>0.22</td><td>0.40</td></tr><tr><td>[C II]<sub>158</sub></td><td>-11.47</td><td>2.043</td><td>-0.1145</td><td>0.17</td><td>4.39</td></tr><tr><td>[N II]<sub>205</sub></td><td>-8.24</td><td>1.165</td><td>-0.0719</td><td>0.20</td><td>0.30</td></tr></tbody></table>

**Note.** The parameters correspond to the fits of the line deficits (line-to-$L_{IR}$ ratios) as a function of $\Sigma_{IR}$ for galaxies in GOALS with available FIR sizes, using Equation (3). The fits are presented in Figure 2 (solid lines). The last column presents the maximum line-to-$L_{IR}$ ratio of the distribution.

A number of reasons have been proposed to explain the observed deficit of FIR lines of purely ionized species—mostly singly and doubly 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 the 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_e^{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_e^{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 blackbody (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_{dust}$ in a LIRG from $\sim 25$ to $\sim 50$ K can boost its $L_{FIR}$ by at least a factor of $\sim 50$ (for $\beta = 1.8$; Planck Collaboration et al. 2011).<fnref n="34" /> 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. 2017). 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]158 emission to originate mostly from low-density, “fossil” H II

<figure><img src="image_4.png" /></figure>

**Figure 3.** Fraction of PDR contribution to the total [C II]158 emission, $f([C II]_{158}^{PDR})$, as a function of (1) the FIR $S_{63}/S_{158}$ continuum ratio (top) and (2) the [O III]88/[N II]122 emission line ratio (bottom), for the LIRGs in GOALS with available measurements in the relevant lines. Top panel: position of the Milky Way (MW) is shown for reference as a black square. Bottom panel: LIRGs with measurements of the [N II]122 and [N II]205 lines. Galaxies without available data in any of the $N^+$ lines are shown as small black open circles. The remaining symbols are as in Figure 1. In the top panel, galaxies are color-coded as a function of the AGN fractional contribution to their MIR emission, as estimated using the 6.2 µm PAH EW diagnostic (see Section 4.1). The solid black line represents a fit to the data (see Equation (4)). In the bottom panel, galaxies are color-coded as a function of their [N II]122/[N II]205 ratio, showing that $f([C II]_{158}^{PDR})$ does not depend on $n_e$ (see Accurso et al. 2017). The correlation (black solid line) is described by the following fit: $f([C II]_{158}^{PDR}) = -0.757(\pm 0.012) + 0.176(\pm 0.020)\log([O III]_{88}/[N II]_{122})$, with a dispersion of 0.04.

regions and diffuse ionized gas not associated with OB stars. Then, if the [C II]158$^{ion}$ and FIR emissions truly arise from different ISM phases, the [C II]158$^{ion}$/FIR ratio should be mostly controlled by the $L_{FIR}$ boosting, due to the increase of $T_{dust}$ in progressively more active starbursts.

The top left panel of Figure 5 shows [C II]158$^{ion}$/FIR as a function of $T_{dust}$ (derived from the $S_{63}/S_{158}$ ratio), as well as a fit in log-log space to the data (dashed line; see the figure caption for the fitted equation). Under the assumptions discussed above, this trend indeed suggests that the rise in the radiation field intensity caused by the ongoing starburst is mostly transferred into an increase of $T_{dust}$, which enhances the [C II]158

<footnote n="34">³⁴ Note that the peak emission of a modified blackbody around these temperatures is always within the wavelength range used to calculate $L_{FIR}$, and therefore it can be approximately scaled using the Stephan-Boltzmann law.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** [C II]158/[N II]205 (circles) and [N II]122/[N II]205 (diamonds) line ratios as a function of the electron density, $n_e$, for the LIRGs in GOALS with available data on both nitrogen lines, based on the emission-line models from Oberst et al. (2006) (solid line). The dotted line represents the typical [C II]158/[N II]205 ratio for ionized gas that we use to calculate [C II]158 (see text). Galaxies are color-coded as a function of the [S III]18.7/[S III]33.5 line ratio measured with the long-high (LH) resolution module of Spitzer/IRS (Inami et al. 2013), which traces denser ($n_e \simeq 10^{3-4.5} \text{ cm}^{-3}$) and higher ionized gas (34.79 > $h\nu$/eV > 23.34) than that probed by the nitrogen lines. Galaxies with lower/upper limits in [S III]18.7/[S III]33.5 are marked as upward/downward-pointing semicircles. Higher [S III]18.7/[S III]33.5 ratios imply higher densities. Most LIRGs show [N II]122/[N II]205 ratios compatible with $n_e \simeq 10^{1-2} \text{ cm}^{-3}$, with a median value of 41 cm$^{-3}$ and mean of 45 cm$^{-3}$. The remaining symbols are as in Figure 1.</figcaption></figure>

deficit through the boosting of the FIR continuum. If we assume that the FIR emission scales with $T_{\text{dust}}^{4+\beta}$ (where $\beta = 1.8$), the suppression seen in the [C II]158/FIR ratio of the warmest LIRGs can be explained in terms of a factor of ~30 excess in the IR luminosity 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 µ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 efficiency<fnref n="35" /> 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 the starburst would be transformed into gas heating, and (2) there is still a significant dispersion, of a factor of ~3 or larger, in the [C II]158/FIR ratio at a given $T_{\text{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_{\text{dust}}$ using the fit to the [C II]158/FIR ratio described above (top left panel). After this “IR excess” has been taken into account, we can see a trend for the [C II]158 emission to rise as a function of $T_{\text{dust}}$, suggesting that part of the heating 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_{\text{dust}}$ as a result of the varying

PDR contribution to the total [C II]158 budget (see color coding in the figure).

The bottom panel of Figure 5 conveys the same information but in the form of the classical [C II]158/FIR ratio, which is displayed as a function of $S_{63}/S_{158}$ and $T_{\text{dust}}$. Note that the $S_{63}/S_{158}$ ratio shown here is approximately equal to the $S_{70}/S_{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$ ([C II]158)PDR 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]158/FIR ratio as a function of $T_{\text{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, in which a larger fraction of the dust mass ($M_{\text{dust}}$) in the most compact, luminous starbursts, would be heated up to higher temperatures by the absorption of ionizing photons. In other words, because the cooling efficiency of dust increases as $T_{\text{dust}}^{4+\beta}$, 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 10^2 G_0$; Tielens 2005; see Section 5.3). In addition, a fraction of the energy is transferred via the photoelectric effect into gas heating, as shown in the top right panel, which cools down, increasing the [C II]158PDR emission.

Recently, using Herschel/HIFI, Goicoechea et al. (2015) have studied the variation of the [C II]158/FIR ratio on ~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_{\text{FIR}}/M_{\text{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 µm silicate absorption feature, $S_{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 kiloparsec-scale, integrated emission of LIRGs. In terms of total $M_{\text{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_{\text{dust}}$ to 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) assumes that there are no dust temperature variations. On the other hand, Lombardi et al. (2014) presented column density and effective

<footnote n="35">³⁵ The 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.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The quantities shown in the y-axes of all panels, $[C \text{ \textbar{} II}]_{158}^{\text{ion}}$/FIR (top left), $[C \text{ \textbar{} II}]_{158}$/FIR/$[C \text{ \textbar{} II}]_{158}^{\text{ion}}$-fit (top right), and $[C \text{ \textbar{} II}]_{158}$/FIR (bottom), are presented as a function of the FIR $S_{63}/S_{158}$ continuum ratio for the entire GOALS sample. Top left: $[C \text{ \textbar{} II}]_{158}^{\text{ion}}$ deficit. The dashed black line is a log-log fit to the data where log($[C \text{ \textbar{} II}]_{158}^{\text{ion}}$/FIR) = -3.343(±0.021) -1.82(±0.11)log($S_{63}/S_{158}$), and with a dispersion of 0.20 dex. Top right: $[C \text{ \textbar{} II}]_{158}$ deficit divided by the fit to the $[C \text{ \textbar{} II}]_{158}^{\text{ion}}$/FIR ratio, with a diagnostic grid indicating different contributions of the PDR component to the total $[C \text{ \textbar{} II}]_{158}$ emission, $f([C \text{ \textbar{} II}]_{158}^{\text{PDR}})$. The dashed line indicates a ~0% PDR fraction, set by the fit to the data in the top left panel. Additional increasing contributions of $[C \text{ \textbar{} II}]_{158}^{\text{PDR}}$ = [1, 3, 6, 9, 19] × $[C \text{ \textbar{} II}]_{158}^{\text{ion}}$ are shown with dotted lines, from ~50% to 95%. Bottom: same data and grid as in the top right panel, but presented in the more typical form of the $[C \text{ \textbar{} II}]_{158}$ deficit plot. All panels are color-coded as a function of $f([C \text{ \textbar{} II}]_{158}^{\text{PDR}})$. Galaxies with lower/upper limits in $f([C \text{ \textbar{} II}]_{158}^{\text{PDR}})$ are marked as upward/downward-pointing semicircles. The remaining symbols are as in Figure 1. Small open black circles are sources with unavailable $f([C \text{ \textbar{} II}]_{158}^{\text{PDR}})$ due to lack of measurement of the $[N \text{ \textbar{} II}]_{205}$ line.</figcaption></figure>

dust temperature maps based on near-IR, 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 of ~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 owing to a compactification of the star-forming regions (see end of next section), this interpretation is also 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_{\text{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. $[\text{O I}]_{63}/[\text{C II}]_{158}^{\text{PDR}}$ as a Probe of the Gas Kinetic Temperature

The $[\text{O I}]_{63}$ line is an important coolant of the warm, neutral ISM. A substantial fraction of the $[\text{O I}]_{63}$ emission is believed to originate in PDRs, mostly at optical depths similar to those of the locations responsible for the $\text{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;
<figure><img src="image_7.png" /><figcaption>**Figure 6.** [O I]$_{63}$/[C II]$_{158}^{PDR}$ ratio as a function of the FIR $S_{63}/S_{158}$ continuum ratio for the entire GOALS sample (open circles) and the LIRGs from Farrah et al. (2013) (open squares), using the central spaxel aperture to focus on the galaxy nuclei and enhance the $T_{dust}$ dynamic range. Galaxies are color-coded based on the strength of the 9.7 µm silicate absorption feature, $S_{9.7\mu m}$. The remaining symbols are as in Figure 1. The black line is a fit to the data using the radiative transfer code RADEX. The code calculates the population levels of the fine-structure transitions of the O I and C II species as a function of $T_{gas}^{kin}$ (see upper y-axis), which we have parameterized to depend linearly on $T_{dust}$ (see the text and Equations (5) and (6)). The dotted lines are the 1 σ dispersion around the trend. The remarkably tight correlation between [O I]$_{63}$/[C II]$_{158}^{PDR}$ and $T_{dust}$ (or $T_{gas}^{kin}$) only breaks down at low $S_{9.7\mu m}$, albeit only in the most extremely luminous, warm LIRGs and ULIRGs (the dot-dashed line represents 0.5 times the fitted trend, a threshold below which the [O I]$_{63}$ emission may be becoming optically thick or be self-absorbed). As can be seen, the majority of galaxies show $T_{gas}^{kin}$ lower than the temperatures of the energy levels of the emission lines considered ($E_{ul,[CII]158}/k_B = 92$ K, $E_{ul,[OII]63}/k_B = 228$ K).</figcaption></figure>

Kaufman et al. 1999). While [O I]$_{63}$ can be as weak as ~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_H^{cr,[OII]63} \simeq 10^6$ cm$^{-3}$), even surpassing 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 fourfold 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]$_{158}^{PDR}$ 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$ ([C II]$_{158}^{PDR}$) (capped at 90%) from the $S_{63}/S_{158}$ ratio. Figure 6 shows that [O I]$_{63}$/[C II]$_{158}^{PDR}$ is tightly correlated with the FIR $S_{63}/S_{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]$_{158}^{PDR}$ ratio as a function of the kinetic temperature of the gas ($T_{gas}^{kin}$), assuming $T_{gas}^{kin} = \eta T_{dust}$. We have further assumed that the excitation is dominated by collisions with neutral hydrogen atoms, a gas density of $n_H = 10^3$ cm$^{-3}$ (see Section 5.3), and a column density of $N_H = 10^{16}$ cm$^{-2}$. In this regime, and for lower volume densities and columns up to $N_H \simeq 10^{18}$ cm$^{-2}$, both lines are optically thin and also $n_H < n_H^{cr}$.

Despite the simplicity of the approach, the RADEX model reproduces the observed trend between [O I]$_{63}$/[C II]$_{158}^{PDR}$ and $T_{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_H = 10^{17}$ and $10^{18}$ 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]$_{158}^{PDR}$ ratios. With any of these scalings we infer gas kinematic temperatures similar to or lower than the excitation energy of the [C II]$_{158}$ fine-structure transition, $E_{ul,[CII]158}/k_B = 92$ K. This is in contrast to suggestions that the [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_{63}/S_{158}$ using a second-order polynomial function:

$$ \log([\text{O I}]_{63}/[\text{C II}]_{158}^{\text{PDR}}) = -0.16 + 0.71 \log S_{63}/S_{158} + 0.09(\log S_{63}/S_{158})^2. \quad (5) $$

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

$$ \log(S_{63}/S_{158}) = 0.22 + 1.34 \log x - 0.25(\log x)^2, \quad (6) $$

with $x = [\text{O I}]_{63}/[\text{C II}]_{158}^{\text{PDR}}$, and use Equation (2) to obtain $T_{dust}$ from $S_{63}/S_{158}$. Then, assuming the best scaling factor obtained from the model fit to our data, $T_{gas}^{kin} = 2.06 T_{dust}$, $T_{gas}^{kin}$ can be estimated. 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 [\text{O I}]_{63}/[\text{C II}]_{158}^{\text{PDR}} \lesssim 3.5$ and $25 \lesssim T_{dust} \lesssim 50$ K, or equivalently over the $S_{63}/S_{158}$ and $T_{gas}^{kin}$ 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]$_{158}^{PDR}$ 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 µm silicate strengths as well, $S_{9.7\mu m} \lesssim -1.5$ (see color coding in Figure 6), indicating a similar increase of the dust opacity toward 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 ~0.1 implies some degree of optical thickness in [O I]$_{63}$.

There are only 12 galaxies in GOALS (~5% of the sample) displaying [O I]$_{63}$/[C II]$_{158}^{PDR}$ ratios lower than half the value of the fit shown in Figure 6 (this threshold is denoted by the dotted line), and five of them show clear absorption in the [O I]$_{63}$ line profile (IRASF 08572+3915, NGC 4418, Arp 220, IRASF 17207–0014, and IRAS 17578–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 7 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
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Herschel/PACS spectra (black line) of the five galaxies with [O I]$_{63}$/[C II]$_{158}$ PDR ratios lower than half the value of the fit to the correlation in Figure 6 (dot-dashed line) that show clear signs of absorption in their [O I]$_{63}$ line profiles: IRASF 08572+3915, NGC 4418, Arp 220, IRASF 17207–0014, and IRAS 17578–0400. The continuum emission was estimated using the wavelength ranges delimited by the two pairs of vertical dotted lines at the sides of the emission line (see Section 3 for details). The Gaussian fit to the line uses the wavelength range delimited by the two vertical dashed lines, which are also employed for integrating the line flux.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** [O III]88/[O I]63 ratio as a function of the FIR S<sub>63</sub>/S<sub>158</sub> ratio for the GOALS sample with available measurements of the lines (open circles), and the DGS from Cormier et al. (2015) (open stars). All galaxies are color-coded as a function of the *Spitzer/IRS* short-high (SH) [Ne III]15.6/[Ne II]12.8 ratio, or as black open symbols when there is no information in any of the neon lines. The solid and dotted lines are the median value for the GOALS sample and the median absolute deviation, respectively. The remaining symbols are as in Figure 1. There is no correlation between [O III]88/[O I]63 and T<sub>dust</sub> for most of the GOALS galaxies. The DGS shows higher [O III]88/[O I]63 and [Ne III]15.6/[Ne II]12.8 ratios, indicative of larger φ<sub>V</sub> of ionized gas (with respect to PDRs), and younger starbursts and/or more intense radiation fields. A few LIRGs show also high [O III]88/[O I]63 ratios. However, in most sources this is due to a significant optical thickness of the [O I]63 line (see galaxies marked with plus signs in this figure, based on the threshold [dot-dashed line] shown in Figure 6).</figcaption></figure>

doublet at 119 µm, OH<sub>119</sub> (IRAS 07251–0248, IRASF 15250+3609, and IRASF 17207–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 ULIRG sample from Farrah et al. (2013), four out of the five galaxies that are below the [O I]63/[C II]<sup>PDR</sup> threshold also present signs of OH<sub>119</sub> 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, which arises from the dense
neutral medium, the [O III]88 emission originates from ionized
gas with lower densities (n_e^cr ~ 500 cm⁻³). Thus, the two lines
clearly probe distinct phases of the ISM, providing a first-order
approximation to the ratio of volume filling factors, φ_V,
between the warm ionized gas and the neutral gas, assuming
that 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 and 10 times higher in low-metallicity galaxies than
in normal star-forming galaxies, which suggests a correspond-
ingly higher φ_V of ionized gas with respect to PDR emission in
the former.

In Figure 8 we show [O III]88/[O I]63 as a function of the FIR S<sub>63</sub>/S<sub>158</sub> 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<sub>63</sub>/S<sub>158</sub> ratio for the low-metallicity sample, which
only have published data of the S<sub>70</sub>/S<sub>100</sub> ratio, we use the same
method described in Section 4.2 to obtain T<sub>dust</sub> and then
interpolate the fitted mBB emission at 63 and 158 µm to obtain
the S<sub>63</sub>/S<sub>158</sub> ratio for each galaxy. Figure 8 suggests that there is
no clear correlation between [O III]<sub>88</sub>/[O I]<sub>63</sub> and the FIR color
in LIRGs. The median [O III]<sub>88</sub>/[O I]<sub>63</sub> ratio of our sample is
0.46 ± 0.30. Low-metallicity galaxies populate a region of the
parameter space with higher T<sub>dust</sub> (≥30 K) and larger φ<sub>V</sub> of
ionized gas (≥5×) than dusty, metal-rich systems. Figure 8
also indicates that galaxies with [O III]<sub>88</sub>/[O I]<sub>63</sub> ratios in excess
of ≃1 have on average larger [Ne III]15.6/[Ne II]12.8 emission
line ratios as well (≥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, likely due to a combination of
less effective dust cooling (caused by a lack of metals) and to
the presence of harder ionization fields, raising the gas
temperature and its ionization state. Conversely, 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 ≥ 1. However, most of them have S<sub>9.7</sub> µm ≲
−1.5 and are outliers in the [O I]63/[C II]<sup>PDR</sup> versus S<sub>63</sub>/S<sub>158</sub>
correlation described in Section 5.2.1 (sources below the
threshold marked by the dot-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 plus
signs in Figure 8).

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_e^cr ~ 500 and ~300 cm⁻³, respectively) but
different ionization potentials (~35.1 eV and 14.5 eV, respec-
tively) 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<sub>eff</sub>. Indeed, within the context of
star formation, [O III]88 emission is mostly powered by O stars
(T<sub>eff</sub> ≳ 30,000 K; e.g., Ferkinhoff et al. 2011). Since O stars
are short-lived, with life spans of around a few Myr (≈10 Myr
for an O6 type), T<sub>eff</sub> rapidly drops with age as the stellar
population evolves after the burst of star formation. Alter-
natively, 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 dominant).

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 photoionization 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
<figure><img src="image_10.png" /><figcaption>**Figure 9.** [O III]88/[N II]122 ratio as a function of the [Ne III]15.6/[Ne II]12.8 ratio for the GOALS sample (red filled circles). The solid lines represent MAPPINGS III models for varying metallicity and age of the starburst (shown in different colors) and ionization parameter *q* (different line styles). Models with similar ages but different *q* are linked with dotted lines. Most LIRGs are compatible with solar- or supersolar-metallicity stellar populations ≈2–5 Myr old, except for the ones with the lowest [O III]88/[N II]122 ratios, and expand the entire range in *q*, from 2 × 10⁷ to 2 × 10⁸ cm s⁻¹. Galaxies from the DGS from Cormier et al. (2015) are shown as small, black filled stars.</figcaption></figure>

nebula, such as $n_e$, $q$, metallicity<fnref n="36" /> (in units of $Z_\odot$), and the age of the starburst. Figure 9 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 9. We explore the ITERA models for a range of parameters: $q = (2-40) \times 10^7 \text{ cm s}^{-1}$, starburst ages = 1–10 Myr, and metallicities of 1 and $2Z_\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 $\langle \alpha_{\text{AGN}}^{\text{bol}} \rangle \ge 0.5$ (flagged with crosses in Figure 9). Sources with [O III]88/[N II]122 $\lesssim$ 1 can only be reproduced by models with $2Z_\odot$, starburst ages of $\simeq$1 Myr, and the lowest ionization parameter $q = 2 \times 10^7 \text{ cm s}^{-1}$ (solid green line). We note, however, that this region of the parameter space could be reached at slightly higher metallicities than $2Z_\odot$, but such models are not currently available in ITERA. The remaining bulk of the sample is well described by models with solar- or supersolar-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 starting at a slightly lower limit 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 10^7 \text{ cm s}^{-1}$, or $\gtrsim 8 \times 10^7 \text{ cm s}^{-1}$ when [O III]88/[N II]122 is close to $\simeq$10. We note that even though higher [O III]88/[N II]122 values can be reached at subsolar metallicities and could explain the location of the sources in the DGS studied by Cormier et al. (2015) (filled stars in Figure 9), 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 10^7$ to $2 \times 10^8 \text{ 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 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 interstellar radiation field (ISRF) that impinges the PDR surface. $G$, measured in units of the local Galactic value ($G_0 = 1.6 \times 10^{-3} \text{ erg s}^{-1} \text{ cm}^{-2}$; Habing 1968), and the volume density of the neutral gas, $n_\text{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 starburst. 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]158$^{PDR}$, [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 emission<fnref n="37" /> (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]158$^{PDR}$+[O I]63)/FIR

<footnote n="36">³⁶ $q \equiv Q_0/4\pi R^2 n_e = U_c$, 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_e$ is the electron density in the ionized nebula. Both $q$ and $U$ are proportional to $G/n_\text{H}$, if a fixed stellar population and size for the star-forming region are assumed.</footnote>

<footnote n="37">³⁷ To obtain the [C II]158$^{PDR}$ flux for galaxies without a measurement of the [N II]205 line, we use Equation (4) to estimate $f([C II]158^{PDR})$ (capped at 90%) from the FIR $S_{63}/S_{158}$ continuum ratio. The errors are fully propagated based on the dispersion of the correlation.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Results from the PDRT model fitting for three subsamples of galaxies from GOALS: sub-LIRGs (left), LIRGs (middle), and ULIRGs (right). The figures show the average probability distribution for each subsample, $\langle\chi^2\rangle$, as a function of G and $n_H$. See text for details about what galaxies were used in the modeling, as well as how the $\chi^2$ maps of the galaxies were averaged. The combined best-fit values are the regions of the parameter space colored purple. The input data for the models are the measured $[C \text{ \textbar} \text{II}]_{158}^{\text{PDR}}$, $[\text{O I}]_{63}$, and FIR fluxes of each galaxy. The dashed line denotes a $G/n_H = 2 \text{ cm}^3$, the threshold at which $v_{\text{drift}} \simeq v_{\text{turb}}$ (Kaufman et al. 1999). The solid lines indicate the temperature of the PDR surface, $T_{\text{PDR}}$, for a given combination of G and $n_H$ as derived from the model.</figcaption></figure>

and $[C \text{ \bar{I}}]_{158}^{\text{PDR}}/[\text{O I}]_{63}$ ratios predicted by the PDR model for each combination of parameters, $-0.5 \le \log(G/G_0) \le 6.5$ and $1 \le \log(n_H/cm^{-3}) \le 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_H$ combination averaged over the entire GOALS sample, we stack the $\chi^2$ maps of each galaxy using a weighting of $e^{-\chi^2/[G_0 n_H]}$.

We present the results in Figure 10 separately for sub-LIRGs ($L_{\text{IR}} < 10^{11} L_{\odot}$) (left panel), LIRGs (middle panel), and ULIRGs (right panel). We exclude from this analysis those galaxies whose bolometric luminosity is dominated by an AGN ($\langle\alpha_{\text{AGN}}^{\text{bol}}\rangle \ge 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 $[\text{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 MW 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^2 \sim 1$) shows that most star-forming LIRGs have average ISRF intensities in the range of $\sim 10^{1-3.5} G_0$ and gas densities $n_H \sim 10-10^3 cm^{-3}$. In this regime of the parameter space of relatively low G and $n_H$, where the $[C \text{ \bar{I}}]_{158}^{\text{PDR}} \lesssim 1$, both parameters are degenerate<fnref n="38" /> and the relevant quantity setting the observed line ratio is $G/n_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 \text{ \bar{I}}]_{158}^{\text{PDR}} + [\text{O I}]_{63})/\text{FIR ratio}$ mainly depend on $n_H$ and G, respectively. In turn, the $[\text{O I}]_{63}/[C \text{ \bar{I}}]_{158}^{\text{PDR}}$ ratio is sensitive to G, but also particularly to the gas temperature, as we showed in Section 5.2.1 ($T_{\text{gas}}^{\text{kin}} \propto T_{\text{dust}}$).

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

$$L_{\text{FIR}} \approx 2 G \pi (D/2)^2 \phi_{\text{A}}, \quad (7)$$

where D is the angular diameter of the overall starburst and $\phi_A$ is the beam (areal) 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_0$ and $n_H \simeq 10^{4-5} cm^{-3}$, where low $\chi^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_A = 1$ (but see Section 5.4), a $G = 1 G_0$ implies that in order to generate the total $L_{\text{FIR}}$ emitted by a LIRG, a few times $10^{11} L_{\odot}$, the size of the starbursting region would have to be $D \gtrsim 10$ kpc (increasing as $\phi_A^{-1/2}$). While this might be possible in 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 kiloparsecs, 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 > 10^{0.5} G_0$.

The best-fitted G values obtained through the PDR modeling are several orders of magnitude lower than $10^6 G_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 \text{ \bar{I}}]_{158}/\text{PAH}_{6.2}$ (or overall 5–10 µm wavelength range) emission ratio remains roughly constant in normal star-forming galaxies and ULIRGs (respectively, Luhman et al. 2003; Helou et al. 2001). We note that in those works the ionized gas component was not subtracted from the total $[C \text{ \bar{I}}]_{158}$ emission when compared to the $\text{PAH}_{6.2}$ (see Section 5.1.1) and that while AGN identification was

<footnote n="38">³⁸ This is more clear when looking at the parameter space of the $\chi^2$ map of an individual galaxy. In Figure 10 we show the averaged map, $\langle\chi^2\rangle$, which is biased toward the bottom (best value) of the $\chi^2$ distribution for each individual source, due to the weighting. Thus, the overall $\langle\chi^2\rangle$ map is effectively the accumulation of the best $\chi^2$ values for the entire galaxy sample.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** $L_{\text{FIR,eff}}/G$ ratio (where $L_{\text{FIR,eff}} = L_{\text{FIR}}/2$) as a function of the effective (half-light) area of the LIRGs in GOALS. The sizes have been taken from Lutz et al. (2016). Galaxies are color-coded as a function of the $G/n_H$ ratio. Dashed and dotted lines indicate different values of the area filling factor of the PDR-emitting regions within the galaxy, $\phi_A$, as derived from Equation (7).</figcaption></figure>

performed, its contribution to the IR or bolometric power output of the galaxies was not taken into account. Both issues probably contributed to the large dispersion seen in the [C II]$_{158}$/PAH$_{6.2}$ ratio across each sample (~0.2 dex). The values of our best-fitted G and the indication that the PAH population is mostly neutral on average are also in agreement with the fact that the PAH$_{7.7}$/PAH$_{11.3}$ ratio—a tracer of the ionization state of the PAH molecules (Draine & Li 2001)—does not vary significantly in LIRGs (≈30%; Stierwalt et al. 2014).

Interestingly, within the context of the PDR models used here, the $G/n_H$ ratios derived for around one-third of our LIRGs (see Section 5.4) are larger than the threshold ($G/n_H \simeq 2 \text{ cm}^{-3}$) at which the radiation pressure on dust starts driving grains through the gas at velocities greater than the average turbulent velocity of the gas assumed by the model ($v_{\text{drift}} \simeq v_{\text{turb}} = 1.5 \text{ km s}^{-1}$). This suggests that galaxies with $G/n_H$ above this critical value (dotted line in Figure 10) 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 10 show that sub-LIRGs have G values spanning the lower end of the parameter space of LIRGs (see middle panel), but similar gas densities. On the other hand, the most luminous galaxies in GOALS (ULIRGs, $L_{\text{IR}} \ge 10^{12} 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_H \sim 10-10^3 \text{ cm}^{-3}$. Therefore, there is a progression for G to increase with $L_{\text{IR}}$, but not an equivalent increase in the PDR gas density, as 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 10^{2-3} G_0$, albeit with an upper limit for the derived gas densities, $n_H \lesssim 2 \times 10^2 \text{ cm}^{-3}$, a factor of a few lower than our estimates.

In Section 5.2.1 we estimated $T_{\text{gas}}^{\text{kin}} \simeq 50-100 \text{ K}$, significantly lower than the PDR surface temperatures favored by the models, $T_{\text{PDR}} \simeq 150-1000 \text{ K}$ (see contours in Figure 10). 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_H$ than expected in the models, closer to the upper limit that we measure in the most luminous sources, around $G/n_H \simeq 10 \text{ cm}^3$ and $T_{\text{gas}}^{\text{kin}} \simeq 100 \text{ K} \approx \Delta E/k_B = 92 \text{ K}$. We note, however, that despite the discrepancy in absolute values there is a broad correlation between $T_{\text{PDR}}$ and $T_{\text{gas}}^{\text{kin}}$.

### 5.4. Probing the ISM Structure Using $G/n_H$: A Characteristic Break in $\Sigma_{\text{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_{\text{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_A$.<fnref n="39" /> This means that if we know the actual size of the starburst in a galaxy, we can estimate $\phi_A$. Figure 11 shows the ratio of $L_{\text{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_A = 1$. However, $\phi_A$ can be as small as $\simeq 10^{-3}$ for some of the largest galaxies.

In Figure 12 we show the $G/n_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_{\text{IR}}^* \simeq 5 \times 10^{10} L_{\odot} \text{ kpc}^{-2}$ (see also Elbaz et al. 2011). Below this threshold, an increase in $\Sigma_{\text{IR}}$ is not followed by an increase of $G/n_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_A$ of each galaxy. Increasing $\Sigma_{\text{IR}}$ from $\sim 10^9$ to $\sim 5 \times 10^{10} L_{\odot} \text{ kpc}^{-2}$ can be accomplished by an equal increase in $\phi_A$, from $10^{-3}$ to a few times $10^{-2}$, with $G/n_H$ remaining at a relatively constant value of $\simeq 0.32 G_0 \text{ cm}^3$ (see dotted line). Similarly, $G/n_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 M82 or NGC 891 (e.g., Contursi et al. 2013; Hughes et al. 2015). Above $\Sigma_{\text{IR}}^*$, even though $\phi_A$ still keeps increasing, the upturn in $\Sigma_{\text{IR}}$ is mostly driven by a rise in $G/n_H$, pointing to a change in the physical conditions of the PDRs toward more intense ISRFs for a given gas density. In fact, most of the increase in $G/n_H$ at high $\Sigma_{\text{IR}}$ is due to an increase in G, while variations in $n_H$ are the source of scatter in the trend at any $\Sigma_{\text{IR}}$. For reference, the dashed line in Figure 12 shows a linear dependence (not a fit) of the form $G/n_H[G_0 \text{ cm}^3] = 1.5 \times 10^{-11} \times \Sigma_{\text{IR}}[L_{\odot} \text{ kpc}^{-2}]$.

### 5.5. Physical Interpretation: Young, Compact, Dusty Star-forming Regions

It has been proposed that the deficit of PDR line emission ([C II], [O I]) with respect to that of the IR dust continuum could be understood in terms of the geometry of the star-forming regions within a starburst. This essentially represents the change of the ratio between the PDR surface area and the volume of the dust-emitting regions—a ratio that, in this context, would decrease linearly as the star-forming regions grow in size. This would be further accentuated as the starburst evolves and H II regions begin to overlap, sharing a common PDR envelope and creating effectively a giant, single star-forming region—a “giant Orion.”

<footnote n="39">³⁹ We note that areal filling factors larger than unity are possible if multiple PDRs are stacked along the same line of sight.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** G/n<sub>H</sub> ratio as a function of the effective Σ<sub>IR</sub> for the GOALS sample (where Σ<sub>IR</sub> = L<sub>IR,eff</sub>/πR<sup>2</sup><sub>eff,70µm</sub> and L<sub>IR,eff</sub> = L<sub>IR</sub>/2). The galaxies are color-coded as a function of the PDR area filling factor, φ<sub>A</sub>. Below Σ<sup>*</sup><sub>IR</sub> ≃ 5 × 10<sup>10</sup> L<sub>☉</sub> kpc<sup>-2</sup>, G/n<sub>H</sub> remains constant, ranging between ≃0.2 and 0.6 cm<sup>3</sup> (the dotted line shows the median value, 0.32 ± 0.21 G<sub>0</sub> cm<sup>3</sup>), suggesting that the boosting in Σ<sub>IR</sub> is driven by a higher φ<sub>A</sub>, i.e., a larger number of star-forming regions per unit area. Above that threshold, G/n<sub>H</sub> increases linearly with Σ<sub>IR</sub> (dashed line; where G/n<sub>H</sub> [G<sub>0</sub> cm<sup>3</sup>] = 1.5 × 10<sup>-11</sup> × Σ<sub>IR</sub> [L<sub>☉</sub> kpc<sup>-2</sup>]), indicating a change in the physical properties of the PDRs. At a fixed Σ<sub>IR</sub> the scatter in G/n<sub>H</sub> is mostly due to variations in n<sub>H</sub>, with lower-n<sub>H</sub> galaxies displaying higher G/n<sub>H</sub> ratios.</figcaption></figure>

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, which is what might be expected from a single, large H II region filling the central starburst, and (2) we estimate the volume filling factor of the dense ionized gas to be <5% (see Section 5.1.2), implying that it should mostly originate from compact H II regions.

The connection between Σ<sub>IR</sub> and G/n<sub>H</sub> shown in Figure 12 relates a kiloparsec-scale property of LIRGs to the physics of individual PDRs. In other words, the fact that both quantities scale with each other above Σ<sup>*</sup><sub>IR</sub> 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 parsec-scale star-forming regions. A joint, progressive increase of Σ<sub>IR</sub> and G/n<sub>H</sub> can be accomplished by increasing the fraction of younger, more massive stars per star-forming region.<fnref n="40" /> Indeed, the presence of a larger number of younger stars with intense/harder ionization fields would provide the necessary boost in G and L<sub>IR</sub> to produce the correlation observed in Figure 12, above Σ<sup>*</sup><sub>IR</sub>. Besides, the spatial distribution of gas and dust within these regions must be also different. Very young, less developed star-forming regions probably 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 may peak 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<sub>dust</sub>, 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 Σ<sup>*</sup><sub>IR</sub>). 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 Σ<sub>IR</sub> and G/n<sub>H</sub> are found to be in progressively more advanced interacting stages, we find that most (~80%) late-stage mergers (Stierwalt et al. 2013) have Σ<sub>IR</sub> > 5 × 10<sup>10</sup> L<sub>☉</sub> kpc<sup>-2</sup> and lie along the slope in Figure 12. 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; Linden et al. 2017; Oklopčić 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 ~240 galaxies in

<footnote n="40">⁴⁰ Here 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).</footnote>
GOALS, a complete flux-limited and luminosity-limited sample of all 60 µm-selected LIRG 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:

1. 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_{63}/S_{158}$, or equivalently $T_{\text{dust}}$, are the $[\text{C II}]_{158}$, $[\text{N II}]_{122}$, and $[\text{N II}]_{205}$ emission lines, which show a decrease of a factor of ~50 in the warmest LIRGs. The $[\text{O I}]_{63}$ line displays a smaller deficit (a factor of ~10). The $[\text{O III}]_{88}$ line shows no deficit with $T_{\text{dust}}$, but rather an increasing scatter at high $S_{63}/S_{158}$. However, all lines show an increasing deficit with rising infrared luminosity surface density, $\Sigma_{\text{IR}}$.

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

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

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

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

6. We find that the $[\text{O III}]_{88}/[\text{N II}]_{122}$ ratio is a good tracer of the ionization parameter, $q$, which ranges from $2 \times 10^7$ to $2 \times 10^8 \text{ cm s}^{-1}$ proportionally with $[\text{O III}]_{88}/[\text{N II}]_{122} \sim 1-10$. Moreover, most LIRGs are compatible with solar or supersolar metallicities and

starburst ages $\simeq 2-5$ Myr, in agreement with previous results.

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

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

We note that the correlation between $\Sigma_{\text{IR}}$ and $G/n_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_{\text{IR}}^*$ arises not only from a boost in the radiation field intensity provided by young massive stars but also as a result of 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_H$ versus $\Sigma_{\text{IR}}$ correlation and the interaction stage of the systems, we find that most late-stage mergers ($\sim 80\%$) are above $\Sigma_{\text{IR}}^*$ and show enhanced $G/n_H$ ratios, probably reflecting the fast, “fireworks”-like process of star cluster formation and disruption during the merger.

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 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, a 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.

## ORCID iDs

T. Díaz-Santos https://orcid.org/0000-0003-0699-6083

V. Charmandaris https://orcid.org/0000-0002-2688-1956

N. Lu https://orcid.org/0000-0002-8948-1044

S. Stierwalt https://orcid.org/0000-0002-2596-8531

S. Malhotra https://orcid.org/0000-0002-9226-5350

P. P. van der Werf https://orcid.org/0000-0001-5434-5942

G. C. Privon https://orcid.org/0000-0003-3474-1125

P. F. Goldsmith https://orcid.org/0000-0002-6622-8396

E. J. Murphy https://orcid.org/0000-0001-7089-7325

L. Barcos-Muñoz https://orcid.org/0000-0003-0057-8892

S. T. Linden https://orcid.org/0000-0002-1000-6081

P. Appleton https://orcid.org/0000-0002-7607-8766

K. Iwasawa https://orcid.org/0000-0002-4923-3281

D. B. Sanders https://orcid.org/0000-0002-1233-9998

J. A. Surace https://orcid.org/0000-0001-7291-0087

## References

Abel, N. P., Dudley, C., Fischer, J., Satyapal, S., & van Hoof, P. A. M. 2009, ApJ, 701, 1147

Accurso, G., Saintonge, A., Bisbas, T. G., & Viti, S. 2017, MNRAS, 464, 3315

Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20

Alonso-Herrero, A., Pereira-Santaella, M., Rieke, G. H., & Rigopoulou, D. 2012, ApJ, 744, 2

Armus, L., Charmandaris, V., Bernard-Salas, J., et al. 2007, ApJ, 656, 148

Armus, L., Mazzarella, J. M., Evans, A. S., et al. 2009, PASP, 121, 559

Barcos-Muñoz, L., Leroy, A. K., et al. 2017, ApJ, 843, 117

Beirão, P., Armus, L., Helou, G., et al. 2012, ApJ, 751, 144

Berta, S., Magnelli, B., Nordon, R., et al. 2011, A&A, 532, A49

Blitz, L., & Shu, F. H. 1980, ApJ, 238, 148

Brauher, J. R., Dale, D. A., & Helou, G. 2008, ApJS, 178, 280

Brisbin, D., Ferkinhoff, C., Nikola, T., et al. 2015, ApJ, 799, 13

Chu, J. K., Sanders, D. B., Larson, K. L., et al. 2017, ApJS, 229, 25

Contursi, A., Poglitsch, A., Grácia Carpio, J., et al. 2013, A&A, 549, A118

Cormier, D., Madden, S. C., Lebouteiller, V., et al. 2015, A&A, 578, A53

Croxall, K. V., Smith, J. D., Wolfire, M. G., et al. 2012, ApJ, 747, 81
da Cunha, E., Charmandaris, V., Díaz-Santos, T., et al. 2010, A&A, 523, A78
De Looze, I., Cormier, D., Lebouteiller, V., et al. 2014, A&A, 568, A62
Desai, V., Armus, L., Spoon, H. W. W., et al. 2007, ApJ, 669, 810
Díaz-Santos, T., Alonso-Herrero, A., Colina, L., et al. 2008, ApJ, 685, 211

Díaz-Santos, T., Alonso-Herrero, A., Colina, L., et al. 2010a, ApJ, 711, 328

Díaz-Santos, T., Alonso-Herrero, A., Colina, L., Ryder, S. D., & Knapen, J. H. 2007, ApJ, 661, 149

Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2014, ApJL, 788, L17

Díaz-Santos, T., Armus, L., Charmandaris, V., et al. 2013, ApJ, 774, 68

Díaz-Santos, T., Charmandaris, V., Armus, L., et al. 2011, ApJ, 741, 32

Díaz-Santos, T., Charmandaris, V., Armus, L., et al. 2010b, ApJ, 723, 993

Draine, B. T. 2011, ApJ, 732, 100

Draine, B. T., & Li, A. 2001, ApJ, 551, 807

Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119

Esquej, P., Alonso-Herrero, A., González-Martín, O., et al. 2014, ApJ, 780, 86

Falstad, N., González-Alfonso, E., Aalto, S., et al. 2015, A&A, 580, A52

Falstad, N., González-Alfonso, E., Aalto, S., & Fischer, J. 2017, A&A, 597, A105

Farrah, D., Lebouteiller, V., Spoon, H. W. W., et al. 2013, ApJ, 776, 38

Feigelson, E. D., & Nelson, P. I. 1985, ApJ, 293, 192

Ferkinhoff, C., Brisbin, D., Nikola, T., et al. 2011, ApJL, 740, L29

Ferland, G. J., Korista, K. T., Verner, D. A., et al. 1998, PASP, 110, 761

Fischer, J., Abel, N. P., González-Alfonso, E., et al. 2014, ApJ, 795, 117

Gerin, M., Ruaud, M., Goicoechea, J. R., et al. 2015, A&A, 573, A30
Goicoechea, J. R., Teyssier, D., Etxaluze, M., et al. 2015, ApJ, 812, 75

Goldsmith, P. F., Langer, W. D., Pineda, J. L., & Velusamy, T. 2012, ApJS, 203, 13

Goldsmith, P. F., Yıldız, U. A., Langer, W. D., & Pineda, J. L. 2015, ApJ, 814, 133

González-Alfonso, E., Fischer, J., Graciá-Carpio, J., et al. 2012, A&A, 541, A4

González-Alfonso, E., Fischer, J., Sturm, E., et al. 2015, ApJ, 800, 69

Graciá-Carpio, J., Sturm, E., Hailey-Dunsheath, S., et al. 2011, ApJL, 728, L7

Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3

Groves, B., Dopita, M. A., Sutherland, R. S., et al. 2008, ApJS, 176, 438

Haan, S., Surace, J. A., Armus, L., et al. 2011, AJ, 141, 100

Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421

Helou, G., Malhotra, S., Hollenbach, D. J., Dale, D. A., & Contursi, A. 2001, ApJL, 548, L73

Helou, G., Soifer, B. T., & Rowan-Robinson, M. 1985, ApJL, 298, L7

Herrera-Camus, R., Bolatto, A., Smith, J. D., et al. 2016, ApJ, 826, 175

Herrera-Camus, R., Bolatto, A. D., Wolfire, M. G., et al. 2015, ApJ, 800, 1

Hollenbach, D. J., & Tielens, A. G. G. M. 1997, ARA&A, 35, 179

Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356

Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 18

Howell, J. H., Armus, L., Mazzarella, J. M., et al. 2010, ApJ, 715, 572

Hughes, T. M., Foyle, K., Schirm, M. R. P., et al. 2015, A&A, 575, A17

Inami, H., Armus, L., Charmandaris, V., et al. 2013, ApJ, 777, 156

Iwasawa, K., Sanders, D. B., Teng, S. H., et al. 2011, A&A, 529, A106

Kapala, M. J., Sandstrom, K., Groves, B., et al. 2015, ApJ, 798, 24

Kaufman, M. J., Wolfire, M. G., & Hollenbach, D. J. 2006, ApJ, 644, 283

Kaufman, M. J., Wolfire, M. G., Hollenbach, D. J., & Luhman, M. L. 1999, ApJ, 527, 795

Kawamura, A., Mizuno, Y., Minamidani, T., et al. 2009, ApJS, 184, 1

Kessler, M. F., Steinz, J. A., Anderegg, M. E., et al. 1996, A&A, 315, L27

Kim, D.-C., Evans, A. S., Vavilkin, T., et al. 2013, ApJ, 768, 102

Langer, W. D., Goldsmith, P. F., & Pineda, J. L. 2016, A&A, 590, A43

Larson, R. B. 1981, MNRAS, 194, 809

Laurent, O., Mirabel, I. F., Charmandaris, V., et al. 2000, A&A, 359, 887

Lavalley, M., Isobe, T., & Feigelson, E. 1992, in ASP Conf. Ser. 25, Astronomical Data Analysis Software and Systems I, ed. D. M. Worrall, C. Biemesderfer, & J. Barnes (San Francisco, CA: ASP), 245

Leitherer, C., Schaefer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Levesque, E. M., Kewley, L. J., & Larson, K. L. 2010, AJ, 139, 712

Linden, S. T., Evans, A. S., Rich, J., et al. 2017, ApJ, 843, 91

Lombardi, M., Bouy, H., Alves, J., & Lada, C. J. 2014, A&A, 566, A45

Lord, S. D., Hollenbach, D. J., Haas, M. R., et al. 1996, ApJ, 465, 703

Lu, N., Zhao, Y., Xu, C. K., et al. 2015, ApJL, 802, L11

Lu, N., Zhao, Y., Díaz-Santos, T., et al. 2017, ApJS, 230, 1

Lu, N., Zhao, Y., Xu, C. K., et al. 2014, ApJL, 787, L23

Luhman, M. L., Satyapal, S., Fischer, J., et al. 1998, ApJL, 504, L11

Luhman, M. L., Satyapal, S., Fischer, J., et al. 2003, ApJ, 594, 758

Lutz, D., Berta, S., Contursi, A., et al. 2016, A&A, 591, A136

Magdis, G. E., Rigopoulou, D., Hopwood, R., et al. 2014, ApJ, 796, 63

Magnelli, B., Elbaz, D., Chary, R. R., et al. 2011, A&A, 528, A35

Malhotra, S., Helou, G., Stacey, G., et al. 1997, ApJL, 491, L27

Malhotra, S., Kaufman, M. J., Hollenbach, D., et al. 2001, ApJ, 561, 766

Malhotra, S., Rouchs, J. E., Finkelstein, K., et al. 2017, ApJ, 835, 110

Menéndez-Delmestre, K., Blain, A. W., Smail, I., et al. 2009, ApJ, 699, 667

Miura, R. E., Kohno, K., Tosaki, T., et al. 2012, ApJ, 761, 37

Muñoz, J. A., & Oh, S. P. 2016, MNRAS, 463, 2085

Murphy, E. J. 2013, ApJ, 777, 58

Murphy, E. J., Chary, R.-R., Dickinson, M., et al. 2011, ApJ, 732, 126

Murphy, E. J., Stierwalt, S., Armus, L., Condon, J. J., & Evans, A. S. 2013, ApJ, 768, 2

Oberst, T. E., Parshley, S. C., Stacey, G. J., et al. 2006, ApJL, 652, L125

Oklopčić, A., Hopkins, P. F., Feldmann, R., et al. 2017, MNRAS, 465, 952

Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Mill Valley, CA: University Science Books)

Paladini, R., Umana, G., Veneziani, M., et al. 2012, ApJ, 760, 149

Petric, A. O., Armus, L., Howell, J., et al. 2011, ApJ, 730, 28

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Pineda, J. L., Langer, W. D., Velusamy, T., & Goldsmith, P. F. 2013, A&A, 554, A103

Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2011, A&A, 536, A25

Poglitsch, A., Herrmann, F., Genzel, R., et al. 1996, ApJL, 462, L43
Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Pope, A., Chary, R., Alexander, D. M., et al. 2008, ApJ, 675, 1171

Pound, M. W., & Wolfire, M. G. 2008, in ASP Conf. Ser. 394, Astronomical Data Analysis Software and Systems XVII, ed. R. W. Argyle, P. S. Bunclark, & J. R. Lewis (San Francisco, CA: ASP), 654

Rigopoulou, D., Hopwood, R., Magdis, G. E., et al. 2014, ApJL, 781, L15

Rosenberg, M. J. F., van der Werf, P. P., Aalto, S., et al. 2015, ApJ, 801, 72

Rujopakarn, W., Rieke, G. H., Eisenstein, D. J., & Juneau, S. 2011, ApJ, 726, 93

Sakamoto, K., Wang, J., Wiedner, M. C., et al. 2008, ApJ, 684, 957

Sanders, D. B., Mazzarella, J. M., Kim, D.-C., Surace, J. A., & Soifer, B. T. 2003, AJ, 126, 1607

Sanders, D. B., & Mirabel, I. F. 1996, ARA&A, 34, 749

Sanders, D. B., Soifer, B. T., Elias, J. H., et al. 1988, ApJ, 325, 74

Smith, J. D. T., Croxall, K., Draine, B., et al. 2017, ApJ, 834, 5

Spoon, H. W. W., Farrah, D., Lebouteiller, V., et al. 2013, ApJ, 775, 127

Stacey, G. J., Geis, N., Genzel, R., et al. 1991, ApJ, 373, 423

Stacey, G. J., Hailey-Dunsheath, S., Ferkinhoff, C., et al. 2010, ApJ, 724, 957

Stierwalt, S., Armus, L., Charmandaris, V., et al. 2014, ApJ, 790, 124

Stierwalt, S., Armus, L., Surace, J. A., et al. 2013, ApJS, 206, 1

Thomson, W. 1874, Natur, 9, 441

Tielens, A. G. G. M. 2005, The Physics and Chemistry of the Interstellar Medium (Cambridge: Cambridge Univ. Press)

Tielens, A. G. G. M., & Hollenbach, D. 1985, ApJ, 291, 722

van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

Veilleux, S., Rupke, D. S. N., Kim, D., et al. 2009, ApJS, 182, 628

Voit, G. M. 1992, ApJ, 399, 495

Watson, D. M., Genzel, R., Townes, C. H., Werner, M. W., & Storey, J. W. V. 1984, ApJL, 279, L1

Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1

Whitmore, B. C., Brogan, C., Chandar, R., et al. 2014, ApJ, 795, 156

Wilson, C. D., Rangwala, N., Glenn, J., et al. 2014, ApJL, 789, L36

Xu, C. K., Cao, C., Lu, N., et al. 2014, ApJ, 787, 48

Xu, C. K., Cao, C., Lu, N., et al. 2015, ApJ, 799, 11

Zhao, Y., Lu, N., Xu, C. K., et al. 2016a, ApJ, 820, 118

Zhao, Y., Lu, N., Xu, C. K., et al. 2016b, ApJ, 819, 69

Zhao, Y., Lu, N., Xu, C. K., et al. 2013, ApJL, 765, L13