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

We present an analysis of , [O iii]88, [N ii]122, and 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, , and as a function of FIR color and infrared luminosity surface density, . The median electron density of the ionized gas in LIRGs, based on the [N ii]122/[N ii]205 ratio, is = 41 cm−3. We find that the dispersion in the deficit of LIRGs is attributed to a varying fractional contribution of photodissociation regions (PDRs) to the observed emission, f( ) = / , which increases from ∼60% to ∼95% in the warmest LIRGs. The / ratio is tightly correlated with the PDR gas kinetic temperature in sources where is not optically thick or self-absorbed. For each galaxy, we derive the average PDR hydrogen density, , and intensity of the interstellar radiation field, G, in units of and find G/ ratios of ∼0.1–50 cm3, with ULIRGs populating the upper end of the distribution. There is a relation between G/ and , showing a critical break at ≃ 5 × 1010 L⊙ kpc−2. Below , G/ remains constant, ≃0.32 cm3, and variations in are driven by the number density of star-forming regions within a galaxy, with no change in their PDR properties. Above , G/ increases rapidly with , 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.


INTRODUCTION
One of the most fundamental processes studied in virtually any field of physics is the dissipation of energy (Thomson 1874). In particular, in extra-galactic astrophysics, investigating how interstellar gas cools down is crucial to our understanding of galaxy formation and evo- † Contact email: tanio.diaz@mail.udp.cl lution, 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 farinfrared (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. 1997Luhman et al. 1998Luhman et al. , 2003. But twenty more years needed to pass until the Herschel Space Observatory (Herschel hereafter; Pilbratt et al. 2010) re-opened a new window into the FIR universe providing spectroscopic data for significantly larger samples of nearby, intermediate and high-redshift galaxies (e.g., Graciá-Carpio et al. 2011;Díaz-Santos et al. 2013Rigopoulou et al. 2014;Magdis et al. 2014;Brisbin et al. 2015;Cormier et al. 2015;Rosenberg et al. 2015;Malhotra et al. 2017).
The most important fine-structure lines emitted by atomic species in the ∼ 50-200 µm wavelength range are: [N iii] Each of them originates from a different phase of the ISM, with C + and O emission mostly arising from the neutral and molecular medium within dense photo-dissociation regions (PDRs) surrounding newly formed massive stars (Tielens & Hollenbach 1985;Hollenbach & Tielens 1997, and references therein), and N +,++ and O ++ emission being produced by warm ionized gas, both in the diffuse medium and dense (H ii) regions.
Of particular importance among star-forming galaxies are the so-called luminous IR galaxies (LIRGs; L IR = 10 11−12 L ⊙ ). 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 ⋆ ), consequence of the funneling of large amounts of gas and dust towards their nuclei due to the loss of angular momentum during the dynamical interaction. And while the presence of active galactic nuclei (AGN) in LIRGs is frequent (Petric et al. 2011;Alonso-Herrero et al. 2012), their contribution to the bolometric luminosity of the hosts is still very limited in comparison to ultra-luminous IR galaxies (ULIRGs; L IR ≥ 10 12 L ⊙ 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 (quasarlike) systems (Sanders & Mirabel 1996). This diversity is also reflected on the fact that they cover the entire transition between main-sequence (MS) galaxies and starbursts (Díaz-Santos et al. 2013).
At high redshfit, LIRGs dominate the obscured star formation activity between z ∼ 1-3 (e.g., Berta et al. 2011;Murphy et al. 2011;Magnelli et al. 2011), and a number of works have already shown that local LIRGs (and not ULIRGs) are probably the closest local analogs of this high-z, IR-bright galaxy population (Desai et al. 2007;Pope et al. 2008;Menéndez-Delmestre et al. 2009;Díaz-Santos et al. 2010b;Stacey et al. 2010). Thus, a comprehensive study of the physical properties of lowredshift LIRGs, and specifically of their ISM, is critical for our understanding of the evolution of galaxies and AGN across cosmic time.
To this end, we have performed a systematic study of the most important FIR cooling lines of the ISM in a complete, flux-limited sample of nearby LIRGs, the Great Observatories All-sky LIRG Survey (GOALS; Armus et al. 2009), using Herschel and its Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) as well as the Fourier-transform spectrometer (FTS) of the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010). We combine these observations with mid-IR (MIR) spectroscopy previously obtained with the Infrared Spectrograph (IRS; Houck et al. 2004) on board the Sptizer Space Telescope (Spitzer hereafter; Werner et al. 2004) to provide a panchromatic view of the heating and cooling of the ISM in LIRGs across a wide range of integrated properties such as IR luminosity, compactness, dust temperature, AGN activity, and merger stage. The paper is organized as follows: In §2 we present the LIRG sample and the new Herschel spectroscopy, as well as the ancillary observations used in this work. In §3 we describe the processing and analysis of the data. The basic results are presented in §4. The ISM properties derived for the LIRG sample are discussed throughout §5 in relation to specific emission line ratios. In particular, in §5.1.1 and §5.1.3 we describe how the fractional contribution of PDRs to the [C ii] 158 emission shapes the observed trend between the [C ii] 158 deficit and the FIR color of LIRGs. Section §5.1.2 presents the electron densities found for the ionized gas phase of the ISM derived from the nitrogen lines and discusses the implications. In §5.2.1 we present a link between the [O i] 63 /[C ii] PDR 158 and the kinetic temperature of the PDR gas. We explore PDR covering factors and metallicity variations in §5.2.2, and confront FIR emission line ratios involving oxygen and nitrogen emission lines to photo-ionization models of H ii regions in §5.2.3. We derive the average PDR properties for each galaxy in the sample in §5.3 and show the existence of a critical break in the PDR conditions as a function of luminosity surface density in §5.4. We discuss the physical implications of this results in §5.5. The summary of the results and conclusions are given in §6.

The GOALS Sample
The Great Observatories All-sky LIRG Survey (GOALS; Armus et al. 2009) encompasses the complete sample of 202 LIRG and ULIRG systems contained in the IRAS Revised Bright Galaxy Sample (RBGS; Sanders et al. 2003), which in turn is also a complete, flux-limited sample of 629 galaxies with IRAS S 60 µ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 multiwavelength data obtained mostly from space-born facilities including GALEX UV (Howell et al. 2010), HST optical and near-IR imaging (Haan et al. 2011;Kim et al. 2013), along with Chandra X-ray (Iwasawa et al. 2011), as well as Spitzer /IRS (Díaz-Santos et al. 2010b, 2011Petric et al. 2011;Stierwalt et al. 2013Stierwalt et al. , 2014Inami et al. 2013) and Herschel /PACS and SPIRE spec-troscopy (Díaz-Santos et al. 2013Zhao et al. 2013;Lu et al. 2014Lu et al. , 2015Zhao et al. 2016b). Moreover, a number of ground-based observatories such as VLA, CARMA and ALMA have also been used to observe the GOALS sample (e.g., Murphy et al. 2013a;Murphy 2013b;Xu et al. 2014Xu et al. , 2015Zhao et al. 2016a, among others).
The RBGS, and therefore the GOALS sample, were defined based on IRAS observations. However, the higher angular resolution achieved by Spitzer allowed us to spatially disentangle galaxies within the same LIRG system into separate components. From the more than 290 individual galaxies in GOALS, not all were observed by Herschel. In systems with two or more galactic nuclei, minor companions with MIPS 24 µm flux density ratios smaller than 1:5 with respect to the brightest galaxy were not targeted due to their small contribution to the total IR luminosity of the system. The method used to calculate the IR luminosities of individual galaxies (L IR [8−1000 µm] , as defined in Sanders & Mirabel 1996)  Most of the data were collected as part of our OT1 and OT2 programs (OT1 larmus 1, OT2 larmus 1; P.I.: L. Armus) accounting more than 200 hours of observing time in total. Additional observations that are publicly available in the Herschel archive were included from various projects. The main programs from where these complementary data were gathered are: KPGT esturm 1 (P.I.: E. Sturm), KPOT pvanderw 1 (P.I.: P. van der Werf) and OT1 dweedman 1 (P.I.: D. Weedman).
The IFS on PACS is able to perform simultaneous spectroscopy in the 51−73 or 70−105 µm (3rd and 2nd orders, respectively; "blue" camera) and the 102 − 210 µm (1st 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).

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

Spitzer/IRS Observations
As part of the Spitzer GOALS legacy program, all galaxies observed with Herschel /PACS have available Spitzer /IRS low resolution, R ∼ 60 − 120 (SL module: 5.2 − 14.5 µm, and LL module: 14 − 38 µm), as well as medium resolution, R ∼ 600 (SH module: 9.9 − 19.6 µm, and LH module: 18.7 − 37.2 µm), slit spectroscopy. The IRS spectra were extracted with the Spitzer IRS Custom Extraction (SPICE) software 3 , 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 LL apertures is approximately equivalent (within a factor of ∼ 2) to that of an individual spaxel and a 3 × 3 spaxel box of the IFS in PACS, respectively. Likewise, the SH and LH apertures are approximately equivalent to one spaxel and a 3 × 3 aperture box in PACS, respectively. Aside from the line fluxes, which were presented in Stierwalt et al. (2014), other observables derived from the IRS data that we use in this work are: the strength of the 9.7 µm silicate feature, S 9.7µm 4 , and the equivalent width (EW) of the 6.2 µm PAH, both of which were presented in Stierwalt et al. (2013). We refer the reader to these works for further details about the reduction, extraction and calibration of the IRS spectra as well as for the main results derived from the analysis.

Data Processing
The Herschel Interactive Processing Environment (HIPE; v13.0) application was used to retrieve the raw data from the Herschel Science Archive (HSA 5 ) as well as to process them. We used the script for "LineScan" observations (also valid for "range" mode) to reduce our spectra. We processed the data from level 0 up to level 2 using the following steps: Flag and reject saturated data, perform initial calibrations, flag and reject "glitches", compute the differential signal of each on-off pair of data-points for each chopper cycle, divide by the relative spectral response function, divide by the response, convert frames to PACS cubes, and correct for flat-fielding. Next, for each camera (red or blue), HIPE builds the wavelength grid, for which we chose a final rebinning with an oversample = 2, and an upsample = 1 that corresponds to a Nyquist sampling. The spectral resolution achieved for each line was derived directly from the data and is ∼ 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: 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.

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% 6 , 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) 7 .
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 as "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.  158 ). The line and continuum fluxes provided are those obtained using the (a) method as well as the "best" spatially-integrated value for each individual galaxy. This is defined as the highest flux value that maximizes the S/N ratio of the measurement (line or continuum). In other words, we select, among the three methods described above, the one that has the lowest noise in the measurement while still accounting for the entire flux of the galaxy enclosed by the PACS FoV. As an additional constraint we require that no other source is contained within the aperture used to represent the integrated galaxy flux. Note, however, that some galaxies have companions at distances 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 PACS spectrometer beam efficiency maps (v6) provided in the PACS instrument and calibration webpages 8 , after rebinning them to the regular 5 × 5 spaxel FoV. We note that we do not apply the correction factors based on the FIR color of galaxies provided in Lu et al. (2017) to the [N ii] 205 emission. We use instead the original, point-source calibrated fluxes since we do not use those corrections in our Herschel /PACS data.
The individual FIR luminosities of galaxies belonging to a LIRG system formed by of two or more components were derived in a similar manner to the method used to calculate their L IR (section 2.1). In order to obtain the L IR and L FIR of a given galaxy for the different extraction appertures described above, we scaled the integrated IRAS FIR luminosity of the system (L FIR [42.5−122.5 µm] , as defined in Helou et al. 1985) 9 and L IR with the ratio 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.

The Contribution of AGN
Following the formulation described in Veilleux et al. (2009), we calculate the potential contribution of an AGN to the MIR and bolometric luminosities of each galaxy in GOALS (see also Petric et al. 2011) employing up to five Spitzer /IRS diagnostics, depending on (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-4) of the AOR (col. 19); (7,8) Right ascension and declination of the reference spaxel; (9) Distance from the 63 µm continuum peak to the reference spaxel; (10) Flux and uncertainty of the line measured from the reference spaxel ((a) method. Uncertainty does not account for absolute photometric calibration uncertainty in cols. (10,11,16,17)). Negative values in cols. The IR luminosities and luminosity surface densities of the galaxies, LIR and ΣIR -the latter defined as (LIR/2)/πR 2 70µm,eff -, can be found at: http://goals.ipac.caltech.edu the availability of data: the [Ne v] 14.3 /[Ne ii] 12.8 and [O iv] 25.9 /[Ne ii] 12.8 emission line ratios, the 6.2 µm PAH EW (see also Armus et al. 2007), the S 30 /S 15 dust continuum slope 10 , and the Laurent diagram, which is based on a decomposition of the MIR spectrum of galaxies in three individual components: AGN, PAH and H ii emissions (see Laurent et al. 2000, for a detailed explanation of this diagnostic). These indicators are based on how the central AGN modifies the MIR line and continuum spectrum of a normal star-forming galaxy -either through the ionization of the surrounding gas to higher states and/or via the heating of dust at higher temperatures than starformation-, and provide an estimate of its fractional contribution to the MIR emission, α MIR AGN . Once these val- 10 We modified the reference value for the S 30 /S 15 ratio of a pure starburst/H ii source from 22.4 (see Table 9 in Veilleux et al. 2009) to 10 in order to reflect more accurately the actual distribution of S 30 /S 15 ratios seen in the GOALS sample. We also adapt the PAH EW diagnostic, which in (Veilleux et al. 2009) is developed for the 7.7 µm feature, to be used with the 6.2 µm PAH, and modify the reference value for pure starbursts and its bolometric correction factor accordingly. We further assume there is no PAH destruction due to the AGN when calculating its fractional contribution using this method ( ues are known, corrections based on SED templates of pure starbursting and AGN-powered sources can be applied to derive the fractional contribution of AGN to the bolometric luminosity of galaxies, α bol AGN (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, <α MIR AGN > and <α bol AGN >. Because there is a significant number of non-detections in one or more of the relevant MIR features employed to calculate <α MIR AGN > or <α bol AGN >, we used the Astronomy SURVival analysis package, asurv (Lavalley et al. 1992, v1.3), which adopts the maximumlikelihood Kaplan-Meier (KM) estimator to compute the mean value of univariate distributions containing censored data (Feigelson & Nelson 1985). Table 2  -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 (S 63 /S 158 ) continuum flux density ratio, which is a first order tracer of the average dust temperature in galaxies, T dust . The dynamic range in both x-and y-axes is the same in all panels to facilitate the comparison. We have fitted the data using a functional form of the type: where ǫ 0 denotes a limiting line/FIR ratio for sources with cold FIR colors and no deficits, and δ is the S 63 /S 158 at which the line/FIR ratio has been reduced by a factor of e with respect to ǫ 0 . The ǫ 0 parameter can be understood as the nominal cooling efficiency of each line with respect to that of big dust grains, representative of normal star-forming galaxies. The parameters obtained from the fits to the entire GOALS sample can be found in Table 3 and the fits are displayed in Figure 1 as solid black lines, with the associated 1 σ dispersion around the fit shown as dotted lines. While the choice of this particular functional form is arbitrary, it provides a better description of the trends than a power-law fit, as it is known that the line to FIR ratios do not increase indefinitely at the lower end of the T dust distribution and IR luminosities (e.g., Malhotra et al. 1997Malhotra et al. , 2001Brauher et al. 2008;De Looze et al. 2014;Cormier et al. 2015;Herrera-Camus et al. 2015). That is, the line to FIR ratios level off at low T dust reflecting a cooling efficiency "ceiling" (which depends on each line) of the gas in PDRs, H ii regions and the diffuse ISM, with respect to the energy dissipated by dust in thermal equilibrium in normal galaxies. We convert the S 63 /S 158 ratios into dust temperatures (see upper x-axis in Figure 1) by assuming that the observed FIR continuum emission is produced by a singletemperature modified black body (mBB) with a fixed emissivity index β = 1.8 and whose emission is optically thin. This is a reasonable approximation for spectral energy distribution (SED) fits that do not include data at λ 60 µm. We also provide a practical equation that relates both quantities, S 63 /S 158 and T dust , using the following approximation: We note that this equation is only valid for the dust temperatures and S 63 /S 158 ratios spanned by the galaxies in the GOALS sample. That is 21 ≤ T dust ≤ 48 K or 0.1 ≤ S 63 /S 158 ≤ 4. The error in T dust obtained from this expression is 0.5 K for the dynamic ranges mentioned. Figure 1 shows that there is a common trend for most lines to show stronger deficits as the average T dust becomes warmer -including the two [N ii] lines, which arise from the ionized medium.  ( 2), the exponential fit yields a result that is statistically indistinguishable from a flat trend. Binning the data and obtaining the median values for each bin provides the same result.
A possible interpretation is that the dispersion in the [O iii] 88 /FIR ratio as a function of T dust may be related to the location where the line and dust emission originate within the star-forming regions. As we argue in sections 5.1.3 and 5.5 (see also Díaz-Santos et al. 2013), most of the energy reprocessed by dust may be arising from grains in front of the PDRs, mixed with the ionized gas, and at low optical depths into the molecular cloud. Low ionization lines, that originate from the PDR itself or the outer edge of the H ii region ( Table 3). The potential contribution of an AGN could also introduce additional dispersion in the line deficits, specially at high S 63 /S 158 (Fischer et al. 2014;González-Alfonso et al. 2015) or when the AGN overwhelmingly dominates the MIR emission (Díaz-Santos et al. 2014), but we do not find galaxies with large fractional MIR or bolometric AGN contributions to be distributed with a significantly larger scatter than star-formation dominated sources, based on the Spitzer /IRS spectral diagnostics described in section 4.1. Figure 2 presents the same line deficits but as a function of the luminosity surface density, Σ IR , defined as the ratio of effective luminosity, L IR,eff = L IR /2, divided by the effective area of the source that contains half of its luminosity measured at 70 µm, πR 2 70µm,eff 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 section 5.4). The scatter in the trends of those lines that have a PDR origin, [O i] 63 and [C ii] 158 , is especially small and remarkably constant in relative terms at any Σ IR . Moreover, the trend is followed by nearly all LIRGs regardless of their FIR color or T dust . We have fitted these correlations with a secondorder polynomial function: 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.  . 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 dust of a modified black body (with β = 1.8) that has a S 63 /S 158 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 IR [8−1000 µm] , 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 populate better the high T dust regime. Galaxies with an AGN contributing more than 50 % to their bolometric luminosity, <α bol AGN > ≥ 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 in 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). Table 3 Best Parameters From The Line Deficit vs. S63/S158 Fits Note. -The parameters ǫ0 and δ correspond to the fits of the line deficits as a function of the S63/S158 ratio for the entire GOALS sample using equation 1. The fits are presented in Figure 1 (solid lines). Note. -The parameters correspond to the fits of the line deficits as a function of ΣIR for galaxies in GOALS with available FIR sizes, using equation 3. The fits are presented in Figure 2 (solid lines).
ionized environments. Thus, unless the volume filling factor of the diffuse medium is very high (see section 5.1.2), most of the [C ii] 158 emission -especially in actively star-forming galaxies-is expected to arise from dense PDRs surrounding young, massive stars (Hollenbach & Tielens 1997), where the [C ii] 158 is collisionally excited by neutral and molecular hydrogen, with n cr

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] ion 158 (e.g., Oberst et al. 2006;Beirão et al. 2012;Croxall et al. 2012;Kapala et al. 2015). The layer within the Strömgren sphere where nitrogen is singly ionized ranges between 29.60 and 14.53 eV, close to that where [C ii] ion 158 also originates, 24.38 to 13.6 eV. In addition, given the similar n cr e and E ul /k B of both transitions, the [C ii] ion 158 /[N ii] 205 ratio is roughly constant and depends weakly on n e , the intensity of ionizing field, q, and the kinetic temperature of the gas, T kin gas . The photo-ionization models presented in Oberst et al. (2006) Figure 3 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 (<α MIR AGN > ≥ 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: 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 dust = 25 (± 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 dust to show larger f ([C ii] PDR 158 ) (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] 158 emission in warmer systems, indicating that, even though the [C ii] 158 line shows a large deficit with respect to the FIR emission in progressively more luminous, warmer galaxies, most of the extra [C ii] produced in them originates in dense PDRs. This scenario is supported by the trend presented in Figure 3 (bottom panel), which shows also a positive correlation of f ([C ii] PDR 158 ) with the average hardness of the radiation field seen by the ionized gas, as traced by the 158 ) increases in environments associated with recent episodes of massive star formation. This is also consistent with H ii regions in the LIRG nuclei to be more enshrouded (optically and geometrically thicker; see discussion in section 5.4) than those of evolved starforming complexes where the stellar winds from massive stars and supernovae have already cleared out most of the dust from the star formation sites (Blitz & Shu 1980;Larson 1981). That is, galaxies with more evolved H ii regions (i.e., posterior to experience a starburst event; or simply having more modest star-formation, like the MW) have more of the [C ii] 158 emission arising from the ionized gas (likely from the low density ISM, see below). We also note that there is no clear trend between f ([C ii] PDR 158 ) and the electron density of the ionized gas, n e , as traced by the [N ii] 122 /[N ii] 205 ratio (see color code in the figure). We discuss further implications of the trend seen in Figure 3 (top) in section 5.1.3.

Ionized Gas
As shown in Figure 3,  ratio as a function of the IR luminosity surface density, defined as Σ IR = (L IR /2)/πR 2 70µm,eff , 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 parametrization (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 (dotted-dashed line). 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 iii]33.5 µm line ratio 12 , 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 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 hν < 24.38 eV. Considering that the electron density should at least remain constant, if not increase towards the denser PDR region (n e of a few 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 iii] 18.7 /[S iii] 33.5 ratio likely represents a lower limit to the density of the volume from where [C ii] ion 158 and singly ionized nitrogen emission arises within H ii regions. And because the densities derived from the Sulfur lines are significantly larger than the [C ii] ion 158 and [N ii] 205 critical densities, this implies that the [C ii] ion,HII 158 emission has likely been thermalized, therefore suggesting that most [C ii] ion 158 and [N ii] 205 emission is produced in the diffuse ISM, with a modest contribution from the dense ionized phaseunless the diffuse medium is extremely thin (low volume density) or its average T kin gas very low. We can also use the [N ii] 122 /[N ii] 205 ratio in combination with the models from Oberst et al. (2006) to derive the average electron densities in our LIRG sample. As mentioned above, the region where emission from singly ionized nitrogen atoms originates in an H ii region largely overlaps with that of the [C ii] ion,HII 158 emission. Figure 4 shows the distribution of [N ii] 122 /[N ii] 205 ratios and n e for those galaxies with available measure-ments 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 starbursting 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 ii] 205 and [C ii] ion 158 transitions have been thermalized in the dense H ii regions and thus mostly originate from the diffuse ISM. In this case, if we assume a n e ≃ 500 cm −3 and a n e 1 cm −3 for the dense and diffuse ionized medium respectively, the mean value of the [N ii] 122 /[N ii] 205 ratio implies an average volume filling factor, φ V , of 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 φ V will vary significantly depending on the compactness of the galaxy.
A number of reasons have been proposed to explain the observed deficit of FIR lines of purely ionized species -mostly singly and double ionized nitrogen and oxygen ([N ii] 122 . 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 ) ∼ -2.5 to -1.5 (Voit 1992;Abel et al. 2009;Graciá-Carpio et al. 2011;Fischer et al. 2014 The total luminosity of an optically-thin modified black-body (mBB) increases as a function of its temperature as T 4+β , where β is the emissivity index of the emitting material. Thus, an increase of the average T dust in a LIRG from ∼ 25 to ∼ 50 K can boost its L FIR by at least a factor of ∼ 50 (for β = 1.8; Planck Collaboration et al. 2011) 13 . 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 13 Note that the peak emission of a modified black-body around these temperatures is always within the wavelength range used to calculate L FIR , and therefore it can be approximately scaled using the Stephan-Boltzan law.
the ionized region -close to the PDR front-of dustbounded star-forming regions (see e.g., Abel et al. 2009;Paladini et al. 2012), where the density of material is expected to increase significantly (Draine 2011). Furthermore, the deficit seems to be restricted to the nuclear region of galaxies, where the starburst is ongoing (Díaz-Santos et al. 2014;Smith et al. 2016). Thus, most of the FIR dust continuum emitted at the BB peak would be associated with active, still dust-bounded star-forming regions. On the other hand, and in light of our findings in the previous sections, we expect the [N ii] 205 and [C ii] ion 158 emission to originate mostly from low density, "fossil" H ii regions and diffuse ionized gas not associated with OB stars. Then, if the [C ii] ion 158 and FIR emissions truly arise from different ISM phases, the [C ii] ion 158 /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] ion 158 /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 precise equation). Under the assumptions discussed above, this trend indeed suggests that the rise in the average radiation field intensity produced by young, massive stars born during the ongoing starburst is mostly transferred into an increase of T dust , which enhances the [C ii] 158 deficit through the boosting of the FIR continuum. If we assume that the FIR emission scales with T dust 4+β (where β = 1.8), the suppression seen in the [C ii] ion 158 /FIR ratio of the warmest LIRGs can be explained in terms of a factor of almost ∼ 30 excess in the IR continuum emission rather than a deficit in the observed line flux. This scaling is also able to recover reasonably well the shape of the trend seen in Figure 5. The fit to the data is even better if we consider the BB emission to be optically thick below 100 µman 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 14 alone cannot be the only physical parameter involved in the variation of the line-to-FIR ratios because: (1) that would imply that none of the energy produced by massive stars at the core of star-forming regions would be transformed into gas heating, and (2) there is still a significant dispersion of a factor of ∼ 3 or larger in the [C ii] 158 /FIR ratio at a given T 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 dust using the fit to the [C ii] ion 158 /FIR ratio described above (top left panel). After this "IR excess" has been taken into account, we can see a trend for [C ii] 158 emission to rise as a function of T dust , suggesting that part of the heating from the stellar radiation field is indeed transferred to the gas, increasing the cooling carried out through the [C ii] 158 (and other) line(s). Note that the dispersion also increases with T dust as a result of the PDR contribution to the total [C ii] 158 budget (see color coding in the figure).  158 /FIR ratio as a function of T 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 regard-ing a larger fraction of a given dust mass (M dust ) being heated up to higher temperatures by the absorption of ionizing photons in the most compact, luminous starbursts. Because the cooling efficiency of dust increases as T dust 4+β , a significant amount of dust present within the H ii regions, close to the PDRs (see Draine 2011), could be able to reprocess most of the energy emitted by the stars before it reaches the surrounding neutral/molecular medium. Indeed, the temperatures measured from our data are typical of PDR surfaces (i.e., at A V ∼ 0) exposed to moderate to high radiation fields (G 10 2 G 0 ; Tielens 2005; see section 5.3). In addition, a fraction of the energy is transferred via the photo-electric effect into gas heating, as shown in the top right panel, which cools down, increasing the [C ii] PDR 158 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 FIR /M 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µ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 kpcscale, integrated emission of LIRGs. In terms of total M 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 dust as dusty, normal star-forming (MS) galaxies (da Cunha et al. 2010). Also, we note that the column density fit used in Goicoechea et al. (2015) assume that there are no dust temperature variations. On the other hand, Lombardi et al. (2014) presented column density and effective dust temperature maps based on NIR, Herschel and Planck data for the entire Orion complex at a lower angular resolution than Goicoechea et al. (2015). Considering both works, temperature variations are in the range ∼ 15-35 K, where regions with relatively high temperatures emit significantly higher fluxes, even if the optical depth is substantially lower.
In summary, while it is possible that the column density in LIRGs may indeed increase due to a compactification of the star-forming regions (see end of next section), this interpretation is still compatible with the picture put forward above, where a given amount of dust mass would be accumulated closer to the heating source, thus also increasing T 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.  (Rosenberg et al. 2015). Using ISO /LWS data, Malhotra et al. (2001) found a correlation between the [O i] 63 /[C ii] 158 ratio and the FIR color of star-forming galaxies. Our Herschel /PACS data and Spitzer /IRS observations allow us to explore this correlation in depth for the entire GOALS sample; a 4-fold increase in the number of galaxies available relative to previous studies.
Because of the neutral gas phase origin of the oxygen fine-structure emission, we use the [ ratio as a function of the kinetic temperature of the gas (T kin gas ), assuming T kin gas = η 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 ≃ 10 18 cm −2 , both lines are optically thin and also n H < n cr H . Despite the simplicity of the approach, the RADEX model reproduces the observed trend between and T dust remarkably well (black solid line in Figure 6; the uncertainty around the trend is 1σ = 0.30, dotted lines) with a best fit for η = 2.06 ± 0.03. The value of η 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] PDR 158 ratios. With any of these scalings we infer gas kinematic temperatures similar or lower than the excitation energy of the [C ii] 158 fine-structure transition, E ul,[CII]158 /k B = 92 K. This is in contrast to suggestions that [C ii] 158 deficit is mainly due to thermal saturation of the [C ii] 158 emission and T dust (or T kin gas ) only breaks down at low S 9.7µm , albeit only in the most extremely luminous, warm LIRGs and ULIRGs (the dotted-dashed line marks 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 kin gas lower than the temperatures of the energy levels of the emission lines considered (E ul,[CII] 158 /k B = 92 K, E ul,[OI] 63 /k B = 228 K). (Muñoz & Oh 2016).
The trend can be fitted as a function of S 63 /S 158 using a second order polynomial function: 3.5 and 28 T dust 76 K, or equivalently over the ranges shown in Figure 6.
Some of the warmest, most luminous galaxies seem to deviate from the correlation, showing systematically lower [O i] 63 /[C ii] PDR 158 ratios than those predicted by the trend. The most likely explanation is that in these LIRGs the [O i] 63 emission becomes optically thick due to an increase of the in situ gas column density and/or foreground absorption by cold gas across the line of sight (e.g., Poglitsch et al. 1996). This is in agreement with the fact that most galaxies falling below the fitted correlation show large 9.7 µm silicate strengths as well, S 9.7µm -1.5 (see color-coding in Figure 6), indicating a similar increase of the dust opacity towards the hot, MIR emitting background source. Optically thick [O i] 63 has been already reported by studies with available information of the [O i] 145 emission line Farrah et al. 2013;Rosenberg et al. 2015 (González-Alfonso et al. 2012;Falstad et al. 2015Falstad et al. , 2017. Figure 12 in the Appendix shows the Herschel /PACS spectra of these galaxies, where signs of [O i] 63 absorption are visually noticeable. Additionally, a small fraction of galaxies in GOALS with spectra covering the wavelength range of the OH doublet at 119 µm, OH 119 (IRAS07251-0248, IRASF15250+3609 and IRASF17207-0014) show signs of strong molecular absorption, which is characteristic of FIR-thick, warm, compact objects . Also, in the ULIRGs sample from Farrah et al. (2013), 4 out of the 5 galaxies that are below the [O i] 63 /[C ii] PDR 158 threshold also present signs of OH 119 absorption. In these highly obscured sources, it is possible that an AGN contributes to the dust heating ), but it is impossible to determine the exact nature of the buried source due to its extreme optical thickness.

[O iii]88/[O i]63 and the Ionized/Neutral Gas Filling Factors
In contrast to the [O i] 63 line, [O iii] 88 emission originates from ionized gas with lower densities (n cr e ∼ 500 cm −3 ). Thus, the two lines clearly probe distinct phases of the ISM, providing a first order approximation to the ratio of volume filling factors, φ V , between the warm ionized and neutral gas, assuming the lines are optically thin and that the ratio of temperatures between both phases is relatively constant across galaxies. As pointed out by Cormier et al. (2015), the [O iii] 88 /[O i] 63 ratio is between 2-10 times higher in low metallicity galaxies than in normal star-forming galaxies, which suggest a correspondingly higher φ V of ionized gas with respect to PDR emission in the former.
In Figure 7   the low metallicity sample, which only have published data of the S 70 /S 100 ratio, we use the same method described in section 4.2 to obtain T dust and then interpolate the fitted mBB emission at 63 and 158 µm to obtain the S 63 /S 158 ratio for each galaxy. Figure 7  ratio of our sample is 0.46 ± 0.30. Low metallicity galaxies populate a region of the parameter space with higher T dust ( 30 K) and larger φ V of ionized gas ( 5 ×) than dusty, metal-rich systems. Figure 7 also indicates that galaxies with [O iii] 88 /[O i] 63 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, which is likely caused by a combination of less effective dust cooling (due to a lack of metals) and by the presence of harder ionization fields, rising the gas temperature and the ionization state. In addition, in dusty galaxies the PDR covering factor is likely higher than in normal or low metallicity sources, and thus more energy is captured by neutral and molecular material and cooled via the [O i] 63 line emission.
We note that there are a few LIRGs lying in the same region of the diagnostic diagram as low metallicity galaxies, with [O iii] 88 /[O i] 63 ≥ 1. However, most of them have S 9.7µm -1.5 and are outliers in the vs. S 63 /S 158 correlation (sources be-  The similar critical densities of the [O iii] 88 and [N ii] 122 emission lines (n cr e ∼ 500 cm −3 and ∼ 300 cm −3 , respectively) but different ionization potentials (∼ 35.1 eV and 14.5 eV, respectively) make their ratio a good tracer of the average hardness of the radiation field in a galaxy, or equivalently, the mean, luminosityweighted effective temperature of the massive stars born in the starburst, T eff . Indeed, within the context of star formation, [O iii] 88 emission is mostly powered by O stars (T eff 30000 K (e.g., Ferkinhoff et al. 2011). Since O stars are short lived, with lifespans of around a few Myr ( 10 Myr for an O6 type), T eff rapidly drops with age as the stellar population evolves after a burst of star formation. Alternatively, if a powerful AGN exists, the line flux of highly ionized species can be significantly boosted by its narrow line region. However, we have explored possible correlations between the [O iii] 88 /[N ii] 122 ratio and the AGN fractional contribution to the MIR and/or bolometric luminosity of LIRGs and found no significant trend with any of the diagnostics described in section 4.1, suggesting that these lines originate predominantly from star formation (except maybe in a few particular cases, where the AGN contribution could be significant).
We have used itera, the idl Tool for Emissionline Ratio Analysis , which is based on the Starburst99 stellar evolution synthesis mod-els (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 nebula, such as n e , q 15 , metallicity (in units of Z ⊙ ) and the age of the starburst. Figure 8 shows  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 × 10 −4 ). We note that this only affects the y-axis in Figure 8. We explore the ITERA models for a range of parameters: q = 2-40 × 10 7 cm s −1 , starburst ages = 1-10 Myr, and metallicities of Z ⊙ , and 2 Z ⊙ .
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 <α bol AGN > ≥ 0.5 (flagged with crosses in Figure 8). Sources with [O iii] 88 /[N ii] 122 1 can only be reproduced by models with 2 Z ⊙ , starburst ages of ≃ 1 Myr and the lowest ionization parameter q = 2 × 10 7 cm s −1 (solid green line). We note however, that this region of the parameter space could be reached at slightly higher metallicities than 2 Z ⊙ , but such models are not currently available in itera. The remaining bulk of the sample is well described by models with solar or super-solar metallicity and starburst ages ranging from ∼ 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 2 Myr can be reached at sub-solar metallicities and could explain the location of the sources in the DGS studied by Cormier et al. (2015) (filled stars in Figure 8), this is not likely the reason for the increase in the line ratio observed in dusty systems like LIRGs, but rather an increase in q. Thus, in dusty star-forming galaxies the [O iii] 88 /[N ii] 122 ratio probes mostly the ionizing parameter q for starbursts of a few Myr (i.e., when these emission lines are actually detected), with [O iii] 88 /[N ii] 122 scaling roughly linearly from ≃ 1 to 10 for values of q ranging from 2 × 10 7 to 2 × 10 8 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.

PDR Modeling
The [C ii] 158 and [O i] 63 fine-structure lines are two of the main coolants of the neutral ISM (Malhotra et al. 1997;Rosenberg et al. 2015) and therefore they are the most useful to constrain its physical conditions. We have used the PDR Toolkit (PDRT) wrapper (Pound & Wolfire 2008) to derive the main parameters that govern the PDR emission: the intensity of the UV inter-stellar radiation field (ISRF) that impinges the PDR surface, G, measured in units of the local Galactic value (G 0 = 1.6 × 10 −3 erg s −1 cm −2 ; Habing 1968), and the volume density of the neutral gas, n H . The PDRT is based on the models developed by Kaufman et al. (1999Kaufman et al. ( , 2006, which calculate the temperature structure and line emission from a PDR, based on the formalism for a one-dimensional semi-infinite slab (Tielens & Hollenbach 1985). We assume the case in which the PDRs are illuminated only on one side.
This choice is the most reasonable for starburst galaxies, especially LIRGs, where most of the newly-formed massive stars are likely still embedded in their parent molecular clouds (e.g., Díaz-Santos et al. 2007). That is, while a generic molecular cloud in the disk of a LIRG may be seeing an isotropic ISRF from older stellar populations, most of the line flux will be overwhelmingly dominated by the energy cooled in warm PDRs associated with the most recent burst. Strictly, these models can only be compared to individual star-forming regions. However, they can also provide useful information regarding the typical conditions of the PDRs in a given galaxy when only spatially integrated measurements are available.
We use the integrated [C ii] PDR 158 , [O i] 63 and FIR fluxes as input to the PDRT models, excluding sources with limits in any of these quantities. Note that the [C ii] fluxes used here are those obtained after subtracting the ionized gas contribution from the line emission 16 (see section 5.1.1); a critical step needed to obtain consistent results from the PDR analysis. ratios predicted by the PDR model for each combination of parameters, -0.5 ≤ log (G/G 0 ) ≤ 6.5 and 1 ≤ log (n H /cm −3 ) ≤ 7, in 0.25 dex steps. The flux uncertainties employed to calculate the χ 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 entire GOALS sample we combine the χ 2 ν maps of each galaxy using a weighting of e −χ 2 ν[G 0 ,n H ] /2 . We present the results in Figure 9 separately for sub-LIRGs (L IR < 10 11 L ⊙ ) (left panel), LIRGs (central panel) and ULIRGs (right panel). We exclude from   (Kaufman et al. 1999). The solid lines indicate the temperature of the PDR surface, T PDR , for a given combination of G and n H as derived from the model. this analysis those galaxies whose bolometric luminosity is dominated by an AGN (<α bol AGN > ≥ 0.5), sources where there is a mismatch between the spatial location of the lines and/or continuum emission peaks, galaxies with companions within the aperture used to measure their fluxes, and sources where the [O i] 63 line shows signs of self-absorption (see section 5.2.1). We do not attempt to correct for the latter effect since the models already account for the possible optical thickness of the line within the PDR, and self-absorption due to intervening cold material along the line(s) of sight is completely unconstrained. Observations at significantly higher spectral resolution than our PACS data, such as those obtained for the Milky-Way with Herschel /HIFI, would be needed in order to correct for this effect (Pineda et al. 2013;Gerin et al. 2015;Langer et al. 2016).
As we can see, the location of the main cloud of best fitting parameters (purple color, χ 2 ν ∼ 1) show that most star-forming LIRGs have average ISRF intensities in the range of ∼ 10 1−3.5 G 0 and gas densities n H ∼ 1-10 3 cm −3 . In this regime of the parameter space of relatively low G and n H where the [C ii] line dominates the cooling ( 1), both parameters are degenerate 17 , 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 ( ratio is sensitive to G, but also particularly to the gas temperature as we showed in section 5.2.1 (T kin gas ∝ T dust ). Because in a one-side illuminated cloud the (optically 17 This is more clear when looking at the parameter space of the χ 2 ν map of an individual galaxy. In Figure 9 we show the averaged map, <χ 2 ν >, which is biased towards the bottom (best value) of the χ 2 ν distribution for each individual source due to the weighting. Thus the overall <χ 2 ν > map is effectively the accumulation of the best χ 2 ν values for the entire galaxy sample. (Kaufman et al. 1999), the integrated FIR luminosity produced by a galaxy containing such typical PDR regions is:

thin) FIR intensity radiated back is
where D is the angular diameter of the overall starburst, and φ A is the beam (area) filling factor of the PDR-emitting regions within it. We note that there is a secondary region of the parameter space, around G ≃ 1 G 0 and n H ≃ 10 4−5 cm −3 , where low χ 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 φ A = 1 (but see section 5.4), a G = 1 G 0 implies that in order to generate the total L FIR emitted by a LIRG, a few 10 11 L ⊙ , the size of the starbursting region would have to be D 10 kpc (increasing as φ −1/2 A ). While this might be possible for a few objects, the typical sizes of LIRGs measured in the MIR and FIR are on average significantly smaller than that, with diameters of just a few kpc, or even below 1 kpc in the case of ULIRGs (Díaz-Santos et al. 2010b;Rujopakarn et al. 2011;Lutz et al. 2016;Chu et al. 2017). Thus, we only consider solutions with G > 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 ii] 158 to PAH 6.2 (or overall 5-10 µm wavelength range) emission ratio remains roughly constant in normal star-forming galaxies and ULIRGs (Luhman et al. 2003;Helou et al. 2001, respectively). We note that in those works the ionized gas component was not subtracted from the total [C ii] 158 emission when compared to the PAH 6.2 (see section 5.1.1) and while AGN identification was performed, its contribution to the IR or bolometric power output of the galaxies was not taken into account. Both issues probably contribute to the large dispersion seen in the [C ii] 158 /PAH 6.2 ratio across each sample (∼ 0.2 dex). The values of our best fitted G and the implication of having mostly a neutral PAH population, are also in agreement with the fact that the 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 of around one third of our LIRGs (see section 5.4) are larger than the threshold (G/n H ≃ 2 cm 3 ) at which the radiation pressure on the dust starts driving grains through the gas at velocities greater than the average turbulent velocity of the gas assumed by the model (v drift ≃ v turb = 1.5 km s −1 ). This suggests that galaxies with G/n H above this critical value (dotted line in Figure 9) could have PDRs that may be dynamically unstable (Kaufman et al. 1999, see also Draine 2011 for a discussion regarding dust drift velocities within H ii regions).
The purple regions in the left panel of Figure 9 show that sub-LIRGs have G values spanning the lower end of the parameter space of LIRGs (see center panel), but have similar gas densities. On the other hand, the most luminous galaxies in GOALS (ULIRGs, L IR ≥ 10 12 L ⊙ ; right panel) have the highest G values, ≃ 10 × the average of the sample, but also the same range of gas densities, n H ∼ 1-10 3 cm −3 . Therefore, there is a progression for G to increase with L IR , but not an equivalent increase in the PDR gas density, as it would be expected from more compact environments associated with merger-driven star-formation. The GOALS ULIRG sample agrees reasonably well with the results obtained for the local ULIRGs in the HERUS sample studied by Farrah et al. (2013), who obtained similar G values, between ≃ 3 × 10 2−3 G 0 , albeit with an upper limit for the derived gas densities, n H 2 × 10 2 cm −3 , a factor of a few lower than our estimates.
In section 5.2.1 we estimated T kin gas ≃ 50-100 K, significantly lower than the PDR surface temperatures favored by the models, T PDR ≃ 150-1000 K (see contours in Figure 9). A possible explanation for this discrepancy is that most of the line emission arises from a region deep in the cloud -far removed from the warm PDR surface. In this case, thermal saturation of the [C ii] 158 line would occur at higher G/n H than expected in the models, closer to the upper limit that we measure in the most luminous sources, around G/n H ≃ 10 cm 3 and T kin gas ≃ 100 K ≈ ∆E/k B = 92 K. We note, however, that despite the discrepancy in absolute values there is a broad correlation between T PDR and T kin gas .
5.4. Probing the ISM Structure Using G/n H : A Characteristic Break in Σ IR If most of the IR luminosity of a galaxy is produced in PDRs, the local ISRF, traced by G, and the galaxyintegrated L 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, φ A 18 . This means that if we know the actual size of the starburst in a galaxy, we can estimate φ A . Figure 10 shows the ratio of L FIR,eff to G as a function of the effective physical area of our galaxies measured at 70 µm, where the starburst emission peaks. Within the uncertainties, we find that all LIRGs fall below the dashed line marking φ A = 1. However, φ A can be as small as ≃ 10 −3 for some of the largest galaxies.
In Figure 11 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 Σ * IR ≃ 5 × 10 10 L ⊙ kpc −2 (see also Elbaz et al. 2011). Below this threshold, an increase in Σ 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 φ A of each galaxy. Increasing Σ IR from ∼ 10 9 to ∼ 5 × 10 10 L ⊙ kpc −2 can be accomplished by an equal increase in φ A , from 10 −3 to a few 10 −2 , with G/n H remaining at a relatively constant value of ≃ 0.32 G 0 cm 3 (see dotted line). Indeed, 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 M 82 or NGC 891 (e.g., Contursi et al. 2013;Hughes et al. 2015). Above Σ * IR , even though φ A still keeps increasing, the upturn in Σ IR is mostly driven by a rise in G/n H , pointing to a change in the physical conditions of the PDRs towards more intense IS-RFs for a given gas density. In fact, most of the increase in G/n H at high Σ IR is due to an increase in G, while variations in n H are the source of scatter in the trend at a any Σ IR . For reference, the solid line in Figure 11 shows a linear dependence (not a fit) of the form: Figure 11. The G/n H ratio as a function of the effective Σ IR for the GOALS sample (where Σ IR = L IR,eff /πR 2 eff,70µm and L IR,eff = L IR /2). The galaxies are color-coded as a function of the PDR area filling factor, φ A . Below Σ * IR ≃ 5 × 10 10 L ⊙ kpc −2 , G/n H remains constant, ranging between ≃ 0.2-0.6 cm 3 (the dotted line shows the median value, 0.32 ± 0.21 G 0 cm 3 ), suggesting that the boosting in Σ IR is driven by a higher φ A , i.e., a larger number of star-forming regions per unit area. Above that threshold, G/n H increases linearly with Σ IR (dashed line; where G/n H [G 0 cm 3 ] = 1.5 × 10 −11 × Σ IR [L ⊙ kpc −2 ]) indicating a change in the physical properties of the PDRs. At a fixed Σ IR the scatter in G/n H is mostly due to variations in n H , with lower n H galaxies displaying higher G/n H ratios.
G/n H [G 0 cm 3 ] = 1.5 × 10 −11 × Σ IR [L ⊙ kpc −2 ]. 5.5. Physical Interpretation: Young, Compact, Dusty Star-Forming Regions It has been proposed that the deficit of PDR ([C ii], [O i]) line emission with respect to that of the IR continuum could be understood in terms of the geometry of the star-forming regions. This essentially represents the change of the ratio between the PDR the surface area and volume of the dust-emitting region -a ratio that would decrease linearly as the star-forming region grows in size. This would be further accentuated as the starburst evolves and H ii regions begin to overlap, sharing a common PDR envelope, creating effectively a giant, single star-forming region -a "giant Orion". This picture does not seem consistent with our results for two reasons: (1) the areal filling factor of PDR emission we measure in the majority of LIRGs is much less than unity as might be expected from a single, large H ii region filling the central starburst, and (2) we estimate a volume filling factor of dense ionized gas of 5%, implying that it mostly originates from compact H ii regions.
The connection between Σ IR and G/n H shown in Figure 11 relates a kpc-scale property of LIRGs to the physics of the individual PDRs. In other words, the fact that both quantities scale with each other above Σ * IR indicates that in this regime, the luminosity surface density of the entire starburst increases with the average energy per particle density at the surface of PDRs inside the few pc-scale star-forming regions. A joint, progressive increase of Σ IR and G/n H can be accomplished by increasing the fraction of younger, more massive stars per star-forming region 19 . Indeed, the presence of a larger number of younger stars with harder ionization fields would provide the necessary boost in G and L IR to produce the observed correlation seen in Figure 11, above Σ * IR . Besides, the spatial distribution of gas and dust within these regions must be also different. Very young, less-developed star-forming regions have not had time to expel their dusty envelopes and thus have a larger fraction of the surrounding material closer to the newly born stars. At the same time, their radial density structure peaks at the position of the ionization front, close to the PDR, in contrast to dustless H ii regions where the density profile is rather constant (Draine 2011). Note that in this case the volume and column density of gas and dust may remain constant, as only a redistribution of the intervening material along the radial direction is needed. However, T dust , and thus G, must increase since dust grains are overall closer to the heating source.
Our results indicate that this dusty phase, in which we are able to see a significant number of star-forming regions that are likely matter-bounded and still embedded in their molecular cocoon (e.g., Kawamura et al. 2009;Miura et al. 2012;Whitmore et al. 2014), seems to be associated primarily with the period when the overall starburst region is also most compact (above Σ * IR ). A mechanism often proposed to produce this compaction of the gas and dust in galaxies are major mergers (Sanders et al. 1988;Hopkins et al. 2008). While there is not a clear sequence where galaxies with increasingly larger Σ IR and G/n H are found to be in progressively more advanced interacting stages, we find that most (∼ 80 %) late-stage mergers (Stierwalt et al. 2013) have Σ IR > 5 × 10 10 L ⊙ kpc −2 , and lie along the slope in Figure 11. This lack of a simple correlation with merger stage is consistent with the idea that starforming clumps and stellar clusters can form and be quickly destroyed by powerful stellar feedback in gasrich starbursts and mergers on relatively short (10-20 Myr) timescales (Whitmore et al. 2014;Oklopčić et al. 2017;Linden et al. 2017), significantly shorter than the timescale of the entire merging process (over several hundred Myr). Thus, rather than a single or small number of giant, long-lived, and relatively normal H ii regions, it may be more accurate to think of the nuclear starburst in a LIRG as a collection of dense, dusty and young starforming regions that are individually short-lived but constantly replenished, like "fireworks" triggered by a galactic merger.

CONCLUSIONS
We obtained Herschel /PACS spectroscopy of the main FIR fine-structure emission lines for a sample of ∼ 240 galaxies in GOALS, a complete flux-limited and luminosity-limited sample of all 60 µm-selected LIRGs systems detected in the nearby universe (z < 0.09). We combined these observations with Herschel /SPIRE and MIR Spitzer /IRS spectroscopic data to derive the main physical characteristics of the neutral and ionized gas in dense PDRs and H ii regions and provide a comprehensive view of the heating and cooling of the ISM as a function of galaxy-integrated properties. We have found the following results: • The FIR lines with the most pronounced "deficits" (greatest decline of the line flux with respect to FIR continuum emission) as a function of S 63 /S 158 , or equivalently • 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 ∼ 10 1−3.5 G 0 , increasing with L IR , and n H ∼ 1-10 3 cm −3 regardless of the luminosity of galaxies. The G/n H , ratio ranges between ≃ 0.1 and 50 G 0 cm 3 for the entire sample, with ULIRGs populating the upper end of the distribution at G/n H 2 G 0 cm 3 . We use the scaling relation between G and L IR to estimate the area filling factor of the PDRs within the overall area of the starburst region, φ A , and find φ A ≃ 1-10 −3 .
• The G/n H ratio shows a distinct transition at Σ * IR ≃ 5 × 10 10 L ⊙ kpc −2 , with the ratio being constant below this value, and increasing linearly with Σ IR above it. We suggest that below this threshold, the increase in Σ IR is purely driven by an increase of the number density of star-forming regions in the galaxy (∝ φ A ), without changing the average PDR physical conditions. Above Σ * IR , G/n H also starts to increase, signaling a departure from the typical PDR properties of normal star-forming galaxies towards more intense/harder radiation fields and compact geometries typical of starbursting sources.
We note that the correlation between Σ 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 Σ * IR arises not only from a boost in the radiation field intensity provided by young massive stars but also due to finding these star-forming regions in an earlier, dustier phase of their evolution. While we do not find a clear trend between the position of galaxies along the G/n H vs. Σ IR correlation and the interaction stage of the systems, we find that most late-stage mergers (∼ 80%) are above Σ * IR and show enhanced G/n H ratios, probably reflecting the fast, "fireworks"-like process of star-cluster formation and disruption during the merger.