# Gas and dust cooling along the major axis of M 33
(HerM33es)

Carsten Kramer, Thomas Nikola, Sibylle Anderl, Frank Bertoldi, Médéric Boquien, Jonathan Braine, Christof Buchbender, Françoise Combes, Christian Henkel, Israel Hermelo, et al.

► To cite this version:

Carsten Kramer, Thomas Nikola, Sibylle Anderl, Frank Bertoldi, Médéric Boquien, et al. Gas and dust cooling along the major axis of M 33 (HerM33es): Herschel/PACS [C II] and [O I] observations. Astronomy and Astrophysics - A&A, 2020, 639, pp.A61. 10.1051/0004-6361/201936754. hal-02893678

HAL Id: hal-02893678

https://hal.science/hal-02893678

Submitted on 8 Jul 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Gas and dust cooling along the major axis of M 33 (HerM33es)

## Herschel/PACS [C II] and [O I] observations*,**

Carsten Kramer¹,², Thomas Nikola³, Sibylle Anderl⁴,⁵, Frank Bertoldi⁶, Médéric Boquien⁷, Jonathan Braine⁸,
Christof Buchbender⁹, Françoise Combes¹⁰, Christian Henkel¹¹,¹², Israel Hermelo¹³, Frank Israel¹⁴,
Monica Relaño¹⁵,²⁰, Markus Röllig⁹, Karl Schuster¹, Fatemeh Tabatabaei¹⁶,¹⁹, Floris van der Tak¹⁷,¹⁸,
Simon Verley¹⁵,²⁰, Paul van der Werf¹⁴, Martina Wiedner¹⁰, and Emmanuel M. Xilouris²¹

*(Affiliations can be found after the references)*

Received 22 September 2019 / Accepted 5 May 2020

## ABSTRACT

**Context.** M 33 is a gas rich spiral galaxy of the Local Group. Its vicinity allows us to study its interstellar medium (ISM) on linear scales corresponding to the sizes of individual giant molecular clouds.

**Aims.** We investigate the relationship between the two major gas cooling lines and the total infrared (TIR) dust continuum.

**Methods.** We mapped the emission of gas and dust in M 33 using the far-infrared lines of [C II] and [O I](63 µm) and the total infrared continuum. The line maps were observed with the PACS spectrometer on board the Herschel Space Observatory. These maps have 50 pc resolution and form a ~370 pc wide stripe along its major axis covering the sites of bright H II regions, but also more quiescent arm and inter-arm regions from the southern arm at 2 kpc galacto-centric distance to the south out to 5.7 kpc distance to the north. Full-galaxy maps of the continuum emission at 24 µm from Spitzer/MIPS, and at 70 µm, 100 µm, and 160 µm from Herschel/PACS were combined to obtain a map of the TIR.

**Results.** TIR and [C II] intensities are correlated over more than two orders of magnitude. The range of TIR translates to a range of far ultraviolet (FUV) emission of $G_{0,obs} \sim 2$ to 200 in units of the average Galactic radiation field. The binned [C II]/TIR ratio drops with rising TIR, with large, but decreasing scatter. The contribution of the cold neutral medium to the [C II] emission, as estimated from VLA H I data, is on average only 10%. Fits of modified black bodies to the continuum emission were used to estimate dust mass surface densities and total gas column densities. A correction for possible foreground absorption by cold gas was applied to the [O I] data before comparing it with models of photon dominated regions. Most of the ratios of [C II]/[O I] and ([C II]+[O I])/TIR are consistent with two model solutions. The median ratios are consistent with one solution at $n \sim 2 \times 10^2 \text{ cm}^{-3}$, $G_0 \sim 60$, and a second low-FUV solution at $n \sim 10^4 \text{ cm}^{-3}$, $G_0 \sim 1.5$.

**Conclusions.** The bulk of the gas along the lines-of-sight is represented by a low-density, high-FUV phase with low beam filling factors ~1. A fraction of the gas may, however, be represented by the second solution.

**Key words.** galaxies: ISM – galaxies: individual: M 33 – infrared: galaxies – infrared: ISM

## 1. Introduction

The strongest cooling line of the interstellar medium (ISM) in galaxies is usually the [C II](158 µm) line (e.g., Brauher et al. 2008). It is an important extinction-free tracer of star formation. Only in regions of high density and high temperature can the [O I] 63 µm line become the dominant coolant (Kaufman et al. 1999). Mapping the emission of [C II] and [O I] in nearby galaxies at high angular resolutions allows us to study their variation with the local environment, covering not only giant molecular clouds (GMCs) and sites of massive star formation at different galacto-centric distances, but also the diffuse inter-arm gas.

Galactic [C II] emission originates from various ISM phases. On global scales in the Milky Way, 30% of the [C II] emission stems from molecular gas which is bright in CO, that is from photon dominated regions (PDRs), while 25% arises from CO-dark molecular gas, 25% from atomic gas, and another 20%

from diffuse, ionized gas (Pineda et al. 2014). These contributions vary depending on the environment and the location in the Milky Way. In a recent study of ten active star-forming regions in the nearby galaxies NGC 3184 and NGC 628, Abdullah et al. (2017) showed that dense PDRs are the dominant [C II] emitters contributing ~70%, with other important contributions from the warm ionized medium (WIM), while the atomic gas only contributes with less than 5% on average. They also conclude that the relative strengths of all components vary significantly, depending on the physical properties of the gas.

In general, [C II] emission is well correlated with the star formation rate (SFR) in galaxies. The correlation between SFR tracers derived from Hα and 24 µm emission or from the total infrared continuum (TIR), which is integrated between 1 µm and 1 mm wavelength, and [C II] holds well on kiloparsec scales (Herrera-Camus et al. 2015). However, on scales of ~50 pc, one might expect to find a larger scatter as star-forming regions can be spatially distinguished from the diffuse emission. The TIR traces the emission of different grain size populations, from PAHs to large grains, and in a variety of environments and phases. It measures the dust heating caused by FUV photons, but also by soft optical photons of insufficient energy to overcome the grain work function ($E < 6$ eV, Tielens & Hollenbach 1985) and any Coulomb potential to eject electrons, which heat the

<footnote>* Maps of TIR, [C II], [O I] shown in Figs. 2 and 3 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/639/A61</footnote>

<footnote>** Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.</footnote>
gas. Kapala et al. (2015) find that the [CII]/TIR ratio increases with galacto-centric distance in M 31, and Kapala et al. (2017) show that this is caused by changes in the relative hardness of the absorbed stellar radiation field, which are caused by varying stellar populations, dust opacity, and galaxy metallicity, while the photo-electric heating efficiency, the balance of the photo-electric heating rate, and the grain FUV absorption rate (Tielens 2008) may stay constant.

The gas metallicity has a profound impact on the thermal balance of the ISM and hence on the star formation. Low metallicity environments imply lower dust abundance and a drop of the dust-to-gas ratio (Draine et al. 2014), allowing FUV photons of newborn OB stars to penetrate more deeply into molecular clouds.

M 33, also known as the Triangulum galaxy, is a gas-rich Sc flocculent galaxy. With a total baryonic mass of ~5 × 10¹⁰ M☉ (Corbelli 2003) it is the third most massive galaxy in the Local Group, after the Milky Way and M 31 which are about a factor 20 more massive. Its proximity of 840 kpc (Galleti et al. 2004) (10'' ≈ 40.7 pc) makes it one of the nearest disk galaxies. It is only moderately inclined by 56° (Zaritsky et al. 1989). Its metallicity is about half solar (Lin et al. 2017; Toribio San Cipriano et al. 2016; Magrini et al. 2010), similar to that of the Large Magellanic Cloud (Hunter et al. 2007). These properties make M 33 an object of choice to study the interplay of low-metallicity ISM and star formation on local and global scales.

M 33 harbours several of the brightest giant H II regions of the Local Group, including NGC 604. Higdon et al. (2003) observed the six brightest H II regions of M 33 using ISO/LWS and its 70'' beam (280 pc). They concluded that more than half of the observed [C II] arises from PDRs, and that the ionized gas lines can be modeled as arising from a single H II component within their beams. A cut along the major axis of M 33 was studied by Kramer et al. (2013) using ISO/LWS data taken at galacto-centric distances from –8 kpc to +8 kpc, comparing these data with maps of CO and H I, integrated intensities. They suggested that the fraction of [C II] emission stemming from the cold neutral medium traced by H I rises from only 15% in the inner ±4 kpc of M 33 to 40% in the outer parts. They also found the [C II]/TIR ratio to rise from ~0.5% in the inner galaxy to about 4% in the outer parts, at distances beyond 4 kpc.

Here, we present maps of [CII], [OI] (63 µm), and the TIR taken along the major axis of M 33 using Herschel/PACS within the framework of the HerM33es key project (Kramer et al. 2010). Due to its proximity, the spatial resolution provided by Herschel, 12'' for the [CII] line corresponding to ~50 pc, allows us to probe the gas and dust at the scale of individual GMCs (Gratier et al. 2017; Tabatabaei et al. 2014; Xilouris et al. 2012; Boquien et al. 2011). At these scales we expect the global galactic conditions affecting GMCs to be constant and the conditions of these clouds to be exclusively altered by the local star formation. It is then possible to study the star formation and state of the ISM locally, but within the broader galactic context. These maps include among others the southern arm, the nuclear region, and several H II regions, among them BCLMP 691 and BCLMP 302. The PACS data of the latter region had been presented by Mookerjea et al. (2011) and they are included in the present work. The two H II regions BCLMP 302 and BCLMP 691 were also studied within the HerM33es project by Mookerjea et al. (2016) and Braine et al. (2012), respectively, using velocity-resolved spectra of [CII] obtained with HIFI in combination with spectra of H I, and CO, discussing the presence of CO-dark molecular gas, and more generally, trying to disentangle the contribution of the different ISM phases along the lines-of-sight.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** M 33 map of the total far-infrared (TIR) continuum in color together with [CII] contours, both at 12'' resolution. TIR units are erg s⁻¹ cm⁻² sr⁻¹. [CII] contours are 3.4, 7.4, 16.0, 33.6, 73.0 in units of 10⁻⁶ erg s⁻¹ cm⁻² sr⁻¹. A polygon marks the outer edges of the regions mapped in [CII]. Coordinates are in RA and Dec (eq. J2000). The nucleus at 1:33:50.9, 30:39:35.8 (J2000) (Skrutskie et al. 2006) is marked, together with a few prominent giant H II regions. Ellipses delin-uate galacto-centric distances of 2, 4, and 6 kpc.</figcaption></figure>

## 2. Observations

### 2.1. PACS spectroscopy

The region mapped with Herschel in [CII] and [OI] covers a radial strip along the major axis of M 33, which is, with a few gaps about 35' long, and 1.5' wide (Figs. 1–3). The total area covered is 38.5 arcmin². Data were taken with PACS (Poglitsch et al. 2010) on Herschel (Pilbratt et al. 2010). [NII](122 µm) data were also recorded but suffer from baseline instabilities and are not discussed here. Most data were taken in unchopped mode with an observation of an off-source position at the beginning and end of each observation. The off position for the two northern most observations was at 1h35m28s87, +31°46′1.92″. For the other observations, an off position at 1h30m20s90, +30°38′25.44″ was used. The only data taken in wavelength
<figure><img src="image_3.png" /></figure>

**Fig. 2.** Central part of the PACS maps of M 33. *From left to right:* intensities of total infrared continuum (TIR), [C II], [O I](63 µm), [C II]/TIR, and [C II]/[O I] intensity ratios. Intensities are in units of erg s⁻¹ cm⁻² sr⁻¹. Offsets are in arcminutes relative to the nucleus, after de-rotation of the position angle by 22.5°. The TIR contour is at 2 × 10⁻³ erg s⁻¹ cm⁻² sr⁻¹, the threshold intensity used in the correlation plots to distinguish TIR-bright and TIR-weak regions (Figs. 4–6, 8, 5). The color wedge is the same as in Fig. 1. [C II] contours are the same as in Fig. 1. [O I] contours are 1.4, 3.0, 6.3, 13.8 in units of 10⁻⁶ erg s⁻¹ cm⁻² sr⁻¹. The map of [C II]/TIR is overlayed with contours at the 0.6% and 1% level of the ratio. The map of [C II]/[O I] is overlayed with the contour of 3.5. Long tic marks in the color wedges indicate contour levels.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Three northern most regions of M 33 covered with PACS. **Top row:** northernmost region. **Middle row:** northern region. **Bottom row:** BCLMP 691. **From left to right:** intensities of total infrared continuum (TIR), [C II], [O I](63 µm), [C II]/TIR, and [C II]/[O I] intensity ratios. Intensities are in units of erg s⁻¹ cm⁻² sr⁻¹. Offsets are in arcminutes relative to the nucleus, the same as in Fig. 2. The TIR color wedge is the same as in Fig. 1. Contours of TIR, [C II], [O I], [C II]/TIR, and [C II]/[O I] are the same as in Fig. 2. The dashed [C II] contour for BCLMP 691 corresponds to a constant intensity of 22.4 × 10⁻⁶ erg s⁻¹ cm⁻² sr⁻¹ derived from the best-fitting y-intercept in the correlation plot of log[C II]/TIR vs. logTIR, for a slope of -1 (Fig. 9).</figcaption></figure>

switching mode early-on in the observation campaign are those of the H II region BCLMP 302 at 2 kpc galacto-centric distance in the north of M 33 (Mookerjea et al. 2011). These data are included here for completeness.

The strip consists of 21 individual observations (ObsIds, Table A.1). The footprint of the PACS array is 47" × 47", resolved into 5×5 pixels. Each spatial pixel (spaxel) has a size of 9.4". For each center position, a 3 × 3 raster was observed with a step-size of 24", resulting in a total coverage of 95" × 95" per center position. The half power beamwidths (HPBWs) for the [O I] and [C II] line observations are 9.5" and 11.5", respectively. Line widths in M 33 are not resolved: the instrumental spectral resolution is about 90 km s⁻¹ and 240 km s⁻¹, respectively (PACS Observer's Manual V2.3), far larger than line widths observed in M 33 in [C II] (Mookerjea et al. 2016; Braine et al. 2012), CO or H I (Druard et al. 2014; Gratier et al. 2017).

PACS data reduction of the unchopped scans was done using HIPE version 15.0.1 (Ott 2010). The data were calibrated using the PACS calibration tree version 78. The data pipeline was used, including a correction for transients (glitches) and including spectral flat fielding, to obtain level 2 data products for each of the individual 21 ObsIds. Stepping through the pipeline and displaying intermediate results indicated that the “transient correction” is necessary to obtain good final data products. At this

point, the “Off” data products were kept separate from the “On” data products. After running this initial pipeline, a script “Spectoscopy: Mosaic multiple observations” from the menu “PACS Useful scripts” was used to combine the individual ObsIds. This script was modified slightly to average all Off-positions together, subtract them from each On-position, and to create the final data cube. The output is a single FITS file that contains the spectra at each spaxel of the combined map on a 3" grid. Figure B.1 shows a few exemplary spectra.

Python scripts were developed to derive line integrated intensities and to estimate the baseline noise. Polynomials of up to third order were fit to the spectra, masking the edges, which suffer from increased noise and artifacts, and also the region around the expected line position from H I data (Warner et al. 1973). The script determined the best fitting polynomial, and subtracted it. Integrated line intensities were derived by summing over a narrow wavelength range, 2.8 times the resolution, centered on the H I velocities at each PACS spaxel.

The baseline noise was estimated by calculating the standard deviation (root mean square, rms) within the wavelength range used to fit the baseline. The spectral noise varies with spaxel. The median baseline noise, that is, the statistical error, are (1.6±1.1)×10⁻⁷ erg s⁻¹ cm⁻² sr⁻¹ for [C II], (3 ± 7) × 10⁻⁷ erg s⁻¹ cm⁻² sr⁻¹ for [O I](63 µm). For each position the local 1σ uncertainty
was estimated by adding in quadrature the local noise levels and the relative calibration errors of 15% for [C II] and [O I], respectively.

## 2.2. Total infrared continuum and atomic hydrogen

The integrated intensity of the continuum emission between wavelengths of 1 µm and 1 mm, the TIR, has been estimated for each pixel from a weighted sum of the MIPS/Spitzer 24 µm and PACS 70 µm, 100 µm, and 160 µm maps. Before combining these maps, all images were convolved to a common resolution of 12" using the dedicated kernels provided by Aniano et al. (2011) and regridded to a common reference frame with a pixel size of 3" projected onto the [C II] PACS positions. The PACS 70 µm map (Boquien et al. 2015) was rescaled to MIPS fluxes by multiplying with the MIPS scaling factor (mc = 1.03) and dividing by the color correction factor (cc = 0.98) which were interpolated from Table 3 of the PACS report by Müller et al. (2011) assuming a modified blackbody with an emissivity index of β = 1.5 and an average dust temperature of 35 K (Tabatabaei et al. 2014): $S_{70}^{\text{MIPS}} = mc/cc \cdot S_{70}^{\text{PACS}}$. Next, the TIR was calculated using weighting factors from Eq. (1) and Table 1 of Boquien et al. (2011), and using all data larger than zero:

$$ \log S_{\text{TIR}} = 0.220 \cdot \log S_{24} + 0.202 \cdot \log S_{70}^{\text{MIPS}} + 0.243 \cdot \log S_{100} + 0.311 \cdot \log S_{160} + 1.361. \quad (1) $$

In a similar approach, Galametz et al. (2013) used spatially resolved Herschel maps of the KINGFISH sample of nearby galaxies to derive weighting factors to estimate the TIR. Galametz et al. (2013, cf. their Fig. 6) find that the determined calibration coefficients are very similar to those measured in M 33 by Boquien et al. (2011) using Spitzer and Herschel maps.

The median statistical error of the TIR in M 33 as measured in emission free areas, is $(2.0 \pm 1.0) \times 10^{-5} \text{ erg s}^{-1} \text{ cm}^{-2} \text{ sr}^{-1}$. For each position, the local $1\sigma$ uncertainty was estimated by adding in quadrature the statistical error and the relative calibration error of 25%, assuming 20% error on the 160 µm data, and ~10% error on the 100, 70, and 24 µm data, respectively (Boquien et al. 2011).

The TIR map (Fig. 1) resulting from Eq. (1), shows the flocculent spiral structure of M 33, a number of prominent giant H II regions, the arm and interarm regions, together with the drop of intensities with increasing galacto-centric radius.

To estimate the contribution of the cold neutral medium to the [C II] emission, we used a VLA map of M 33 in the atomic hydrogen 21 cm line at 12" resolution (Gratier et al. 2010), matching the resolutions of the other data.

# 3. Analysis

## 3.1. [C II] emission and the total infrared continuum

The maps of TIR and [C II] (Figs. 1–3) show a striking similarity. Both tracers vary by more than two orders of magnitude in intensity. They are bright in the inner star forming arms and gradually fainter in the outskirts at several kiloparsec radial distance. A correlation plot of TIR against [C II] shows their tight correlation (Fig. 4). The results of linear least-squares fits to log I(TIR) against log I([C II]) are shown in Fig. 4 and Table 1. The scatter is larger at fainter intensities and larger than the statistical errors of the individual positions. To gain more insight, we have subdivided the observed positions between those of TIR intensities above and below a given value, tracing the dense, warm,

<figure><img src="image_5.png" /><figcaption>Fig. 4. Total infrared continuum intensity (TIR) versus [C II] intensities in M 33. Intensities are given in erg s$^{-1}$ cm$^{-2}$ sr$^{-1}$. Each point corresponds to one position on a 12" grid at 12" resolution, with an intensity 3 times above the local rms, for [C II], [O I], and TIR. Typical $1\sigma$ errorbars are shown. Colors correspond to the four different regions of the inner and outer galaxy, and above and below a TIR threshold, as described in Sect. 3.1. Straight lines delineate the results of unweighted linear least-squares fits to all data (drawn), the inner galaxy (dashed), and the outer galaxy (dotted) (Table 1).</figcaption></figure>

star-forming spiral arms on the one hand, and the diffuse interarm dust on the other hand. A second subdivision distinguishes between the inner and outer galaxy. The TIR threshold and the radial boundary have no specific physical meaning and were selected to emphasize differences between the four regions. We chose a TIR threshold of $2 \times 10^{-3} \text{ erg s}^{-1} \text{ cm}^{-2} \text{ sr}^{-1}$ (shown as contour in Figs. 2 and 3) and a radial boundary of 2.93 kpc (12') to separate the inner and outer galaxy (cf. Verley et al. 2009). The fit slopes and intercepts of the inner and outer regions of the galaxy do not differ within the errors (Table 1). We find indications for a steepening of the slope for the TIR-bright outer regions (Fig. 4).

As the dust grains emitting the observed TIR emission are heated by the stellar radiation field, the TIR can be used to estimate the FUV energy flux $G_{0,\text{obs}}$, which is given in units of the average radiation field in the solar neighborhood, the Habing field of $1.6 \times 10^{-3} \text{ erg cm}^{-2} \text{ s}^{-1}$ (Habing 1968):

$$ G_{0,\text{obs}} = C_1 C_2 \frac{\text{TIR}}{1.6 \times 10^{-3}/(4\pi)} \quad (2) $$

with TIR in units of erg s$^{-1}$ cm$^{-2}$ sr$^{-1}$. The factor of $C_1 = 0.5$ approximately takes into account the absorption of visible
<table><caption>**Table 1.** Results of unweighted linear least-squares fits to logTIR = b + m × log[C II] (Fig. 4) with TIR and [C II] in units of erg s⁻<fnref n="1" /> cm⁻² sr⁻<fnref n="1" />.</caption><thead><tr><th>Region</th><th>m ± σ</th><th>b ± σ</th><th>r</th></tr></thead><tbody><tr><td>All</td><td>0.99 ± 0.11</td><td>2.15 ± 0.52</td><td>0.91</td></tr><tr><td>Inner</td><td>0.95 ± 0.13</td><td>2.01 ± 0.60</td><td>0.90</td></tr><tr><td>Outer</td><td>0.96 ± 0.29</td><td>1.96 ± 1.50</td><td>0.87</td></tr></tbody></table>

**Notes.** The Pearson correlation coefficient *r* serves as a measure of the linear correlation between TIR and [C II].

photons by grains (Kaufman et al. 1999; Tielens & Hollenbach 1985).

The factor $C_2$ corrects for the fraction of FUV photons leaking out of the galaxy. To estimate the FUV attenuation in M 33, Boquien et al. (2015) combined GALEX FUV maps with 24 µm maps, following Kennicutt et al. (2009). The typical attenuation in M 33 is around 0.6 mag in the FUV band, consistent within the scatter with 0.53 mag previously found by Verley et al. (2009). In the low metallicity environment of M 33, about half of the FUV photons are hence on average absorbed by dust and reradiated in the TIR, while the other half escapes, that is, $C_2 = 2$.

In M 33 at 50 pc resolution, the estimated FUV field $G_{0,obs}$ ranges between ~2 and 200 (Fig. 5). The selected TIR threshold of $2 \times 10^{-3}$ erg s⁻<fnref n="1" /> cm⁻² sr⁻<fnref n="1" /> corresponds to a FUV field of $G_{0,obs} = 15$. Further below, we used the observed ([C II]+[O I])/TIR and [C II]/[O I] ratios together with PDR models of given density and FUV field to improve on this estimate.

The average [C II]/TIR ratio in M 33 is (0.64 ± 0.21)% (Table 2). Figure 5 (left) shows the data at the individual positions, together with unweighted binned averages. For low TIR values the [C II]/TIR ratio shows a large scatter with a high averaged value of (1.1 ± 0.4)% at log(TIR) = -3.75, with the TIR intensities in units of erg s⁻<fnref n="1" /> cm⁻² sr⁻<fnref n="1" />. The binned ratios drop smoothly with TIR, reaching (0.5 ± 0.1)% at high TIR values of log(TIR) = -1.75. The binned averages are consistent with the results of a linear least-squares fit (Fig. 5).

To explore the [C II]/TIR variation further, we study its variation with galacto-centric distance (Fig. 5, right). The [C II]/TIR ratio stays almost constant within a galacto-centric radius of ~4 kpc (cf. Verley et al. 2009). The two northern regions show rising [C II]/TIR ratios, as had already been seen in the ISO/LWS cut (Kramer et al. 2013, hereafter K2013). The northernmost region shows a mean ratio that is a factor 2 higher than in the inner galaxy, 1.3%, with an increased scatter of 0.4% (Table 3). The corresponding map (Fig. 3) shows a small region of low [C II]/TIR ratios, where TIR is just above the threshold of $2 \times 10^{-3}$ erg s⁻<fnref n="1" /> cm⁻² sr⁻<fnref n="1" />, but surrounded by an extended region of high ratios >1% where TIR is below the threshold and [C II] is also weak.

The ISO/LWS cut along the major axis of M 33 (K2013) shows a flat [C II]/FIR<fnref n="1" /> distribution in the inner galaxy, and an abrupt rise of the average ratio in the northern and southern outskirts, at galacto-centric distances beyond ~4.5 kpc out to ~7.5 kpc. Using the FIR/TIR correction factors listed for each position in Table A.1 of K2013, the resulting [C II]/TIR values are consistent within the $1\sigma$ scatter with the ratios determined here. The Herschel/PACS data presented here improve on this work due to their high sensitivities at much better angular resolution, mapping the GMCs also perpendicular to the major axis of M 33.

## 3.2. Emission of [C II], [O I] 63 µm, and the TIR

The critical densities of the [C II](158 µm) and [O I](63 µm) lines for excitation by collisions with H are $3 \times 10^3$ cm⁻³ and $5 \times 10^5$ cm⁻³, respectively, with upper energy levels $E/k_B$ of 92 K and 228 K, respectively (Kaufman et al. 1999). In contrast to the [C II] line, the [O I](63 µm) line is expected to be excited almost exclusively in the dense, warm interface regions of UV illuminated molecular clouds. The distribution of these two emission lines reflects these strong differences in excitation requirements. The intensities of [C II] and [O I] emission are well correlated near the peaks and ridges of the spiral arms of M 33, as is seen in the maps of the central region and of BCLMP 691 (Figs. 2 and 3) and in the correlation plot (Fig. 6, cf. Table 4). At lower intensity levels of [C II] and [O I], in the more diffuse inter-arm regions, the correlation between the two gas tracers is weak. This is seen most prominently in the maps of the two northern most regions (Fig. 3).

As expected, [C II] is stronger than [O I] at most positions. The observed [C II]/[O I](63 µm) ratio ranges between ~0.3 and 20 with a median of 4.5 ± 2.6. The observed ([C II]+[O I](63 µm)/TIR ratios vary between 0.3 and 2.9% with a median of 0.8 ± 0.3% (Fig. 7). Table 5 lists the average values and the scatter for the four regions defined in Sect. 3.1. The TIR bright positions all show ratios of [C II]/[O I] > 1 and ([C II]+[O I])/TIR < 2%, while the positions where [O I] intensities exceed those of [C II], and ([C II]+[O I])/TIR > 2% all lie in the TIR weak regions.

PDR models may allow one to estimate the local excitation conditions, the local densities and the local FUV fields. In order to compare the observations with the predictions of PDR models, we need to consider other components of the ISM that may contribute to the emission. A part of the [C II] emission may stem from the cold neutral medium (CNM) or from a diffuse ionized medium. The [O I](63 µm) line may become optically thick and be affected by foreground absorption. Below, we discuss these possibilities.

### 3.2.1. [C II] from the cold neutral medium

To estimate the contribution of the cold neutral medium (CNM) to the observed [C II] emission, we used H I VLA data, at all positions at which [C II] has been detected, and at almost the same angular resolution. Following the approach of K2013, we first corrected the H I line intensities for the contribution from the surfaces of PDRs assuming a typical $G_0/n$ ratio of $10^{-3}$. Next, the [C II] emission from the remaining neutral gas, the CNM, was estimated assuming a given fractional abundance of C⁺/H in the low-metallicity environment of M 33, optically thin H I and [C II] emission, and a density and temperature which are typical for diffuse atomic clouds (Table 6). While some regions at larger galacto-centric distances show enhanced CNM fractions, as already seen by K2013, the average CNM contribution is only ~10% (Fig. B.2), less than the estimated [C II] calibration error. We do not subtract this small contribution from the CNM from the observed [C II] intensities to estimate the [C II] emission from PDRs.

### 3.2.2. [C II] from the ionized medium

Another part of the [C II] emission may not arise from the neutral gas of PDRs modeled by Kaufman et al. (2006), but from the diffuse, ionized gas. Observations of the [N II](122 µm, 205 µm) lines would allow us to determine electron densities and the

<footnote n="1">¹ K2013 discussed the far-infrared emission (FIR), integrated between 42.5 and 122.5 µm.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Left: [C II]/TIR intensity ratio versus TIR and corresponding FUV field $G_{0,obs}$ in M 33 (see Eq. (2) in Sect. 3.1). The large pink points and errorbars show binned ratios (Table 2). Errorbars represent the standard deviation of the data. The dashed line shows the result of an unweighted least-squares fit to $log([C II]/TIR) = b + m \times logTIR$ and all ratios, which gives $m = -0.16 \pm 0.10$ and $b = -2.64 \pm 0.26$, with correlation coefficient $r = -0.39$. For a few of the individual points $1\sigma$ errorbars are shown. Right: [C II]/TIR versus projected galacto-centric radius in M 33. The large pink points and errorbars show ratios binned over selected ranges of radii (Table 3).</figcaption></figure>

<table><caption>**Table 2.** Binned [C II]/TIR ratios shown in Fig. 5.</caption><thead><tr><th>log TIR</th><th>[C II]/TIR<br>%</th></tr></thead><tbody><tr><td>All</td><td>0.64 ± 0.23</td></tr><tr><td>-3.75</td><td>1.13 ± 0.36</td></tr><tr><td>-3.25</td><td>0.85 ± 0.39</td></tr><tr><td>-2.75</td><td>0.66 ± 0.20</td></tr><tr><td>-2.25</td><td>0.56 ± 0.13</td></tr><tr><td>-1.75</td><td>0.49 ± 0.10</td></tr></tbody></table>

**Notes.** The width of the TIR bins is 0.5 dex. Errors of the individual points are ignored.

importance of the ionized gas (Oberst et al. 2006; Croxall et al. 2012, 2017). Parkin et al. (2013, 2014) studied maps of [C II], [O I], and TIR emission in M 51 and Centaurus A. PDR models better fit the observed [C II]/[O I] and [C II]+[O I]/TIR ratios after correcting for the fraction of [C II] emission originating from the ionized medium using [N II](122 µm, 205 µm). They estimate that 70–80% of [C II] arise in ionized gas of the nucleus and center of M 51, 50% in its arm and interarm regions. In Centaurus A the ionized fraction is 10–20%. Hughes et al. (2015) use the [N II](205 µm) line to estimate fractions of ~16–64% in NGC 891. They also explore the use of an empirical relation

<table><caption>**Table 3.** [C II]/TIR ratios averaged over different galacto-centric distances (cf. Fig. 5).</caption><thead><tr><th colspan="2">Radial range</th><th>[C II]/TIR</th></tr><tr><th>(arcsec)</th><th>(kpc)</th><th>(%)</th></tr></thead><tbody><tr><td>-700–700</td><td>-2.8–2.8</td><td>0.62 ± 0.19</td></tr><tr><td>750–870</td><td>3.1–3.5</td><td>0.63 ± 0.23</td></tr><tr><td>1090–1210</td><td>4.4–4.9</td><td>0.77 ± 0.31</td></tr><tr><td>1320–1410</td><td>5.4–5.7</td><td>1.30 ± 0.41</td></tr></tbody></table>

**Notes.** Errors of the individual points are ignored.

between [N II] and 24 µm emission, to construct a more complete map of the [C II] fraction from the ionized gas. For M 33, we do not attempt to correct the [C II] intensities for the contribution from the ionized gas as observations of the [N II] lines are missing.

### 3.2.3. TIR

A fraction of the TIR emission may stem from non-PDR phases of the ISM like the CNM. We did not try to correct the TIR emission for these contributions before using the PDR models. However, the correlation between TIR and H I intensities is poor in M 33 (Fig. B.3). An unweighted linear least-squares fit gives
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Observed intensities of [O I](63 µm) and [C II] in M 33 with representative 1σ errorbars. The straight black line corresponds to a ratio of 1. The dashed black line corresponds to the average [C II]/[O I] ratio (Table 5). The red and blue dashed lines are the results of unweighted least squares fits to the TIR-bright inner and outer regions, respectively (Table 4).</figcaption></figure>

a correlation coefficient of $r = 0.45$, which is much lower than for the TIR-[C II] relation. This indicates that most of the TIR emission stems from PDR regions, as the CNM contribution to the [C II] emission is also low.

### 3.2.4. Self-absorbed [O I](63 µm) emission

The [O I](63 µm) line is expected to have a higher opacity than the [C II] line and hence may be affected by self-absorption caused by cold foreground clouds of atomic oxygen, as has been seen in velocity resolved spectra and with narrow beams toward bright background sources in the Milky Way (Schneider et al. 2018; Gerin et al. 2015; Leurini et al. 2015; Karska et al. 2014; Lis et al. 2001; Timmermann et al. 1996). Spectra of the [O I](63 µm) line taken toward ULIRGs are often heavily self-absorbed (Rosenberg et al. 2015). Israel et al. (2017) also used PACS/Herschel to observe the center of Centaurus A. Using PDR models allowed them to compare the observed [O I](63 µm) intensity with the intensity predicted from modeling the emission of other FIR lines. In the circumnuclear disk Israel et al. (2017) find high optical depths of the [O I] line of 1.0–1.5. While these PDR models take into account optical depth effects of the spectral lines emerging from the simulated slab of gas, they do not consider possible absorption by foreground gas of different excitation conditions.

The intrinsic line widths observed in M 33 in [C II], CO and H I are far smaller than the velocity resolution of PACS/Herschel (~90 km s⁻¹). Using HIFI/Herschel HerM33es observations, Mookerjea et al. (2016) and Braine et al. (2012) find full-widths at half-maximum (FWHM) of [C II] of 7 to 17 km s⁻¹ in three regions along the major axis of M 33: BCLMP 302, BCLMP 691, and the nucleus. Druard et al. (2014) and Gratier et al. (2017) studied the variation of FWHMs derived from maps

<table><caption>**Table 4.** Results of unweighted linear least squares fits to log[O I] = b + m × log[C II] (Fig. 6).</caption><thead><tr><th>Region</th><th>m ± σ</th><th>b ± σ</th><th>r</th></tr></thead><tbody><tr><td>Bright, inner</td><td>0.84 ± 0.01</td><td>-1.43 ± 0.86</td><td>0.74</td></tr><tr><td>Bright, outer</td><td>1.49 ± 0.25</td><td>1.65 ± 5.65</td><td>0.75</td></tr></tbody></table>

of the complete galaxy in CO 2–1 taken with the IRAM 30 m telescope and H I taken with the VLA. The line widths averaged in radial bins of 1 kpc drop with distance, out to 7.5 kpc: from about 8 to 5 km s⁻¹ for CO, and from about 17 to 12 km s⁻¹ for H I. As the [O I](63 µm) line is expected to arise from the dense, warm cloud interfaces, its linewidths should be similar to those of CO, and smaller than those of [C II], which may also trace other, more diffuse phases.

Assuming low densities in the cold, foreground layers of gas, the majority of the oxygen atoms are in their ground state, allowing to derive an upper limit of the [O I] opacities. The column density of atomic oxygen can then be approximated in terms of the opacity $τ_0$ at the center velocity of the J = 1 – 2 [O I] 63 µm line by

$$N(\text{OI}) = 2 \times 10^{17} \tau_0 \Delta v_{\text{FWHM}} \quad (3)$$

in cm⁻², with the FWHM line width in km s⁻¹ (Vastel et al. 2002; Liseau et al. 2006). For a typical line width of 7 km s⁻¹ (see above) and an average oxygen abundance of 4 × 10⁻⁴ in M 33 (see below), the line center opacity exceeds 1 for $N([\text{OI}]) = 1.4 \times 10^{18} \text{cm}^{-2}$ and $N(\text{H}) > 3 \times 10^{21} \text{cm}^{-2}$. Following Crawford et al. (1985), the opacity is more generally written as function of local density $n$ and temperature $T$ of the [O I] two level system:

$$\tau_0 = \frac{\lambda^3 A_{\text{ul}}}{8\pi\Delta v_{\text{FWHM}}} \left[ \left(1 + \frac{n_{\text{cr}}}{n}\right) \exp(228/T) - 1 \right] \times \left[ \frac{\frac{3}{5} \exp(-228/T)}{1 + \frac{3}{5} \exp(-228/T) + \frac{n_{\text{cr}}}{n}} \right] N(\text{OI}) \quad (4)$$

with the Einstein A-coefficient $A_{ul} = 8.46 \times 10^{-5} \text{ s}^{-1}$, the critical density $n_{cr} = 4.7 \times 10^5 \text{ cm}^{-3}$ for collisions with H-atoms, $h\nu/k_B = 228 \text{ K}$, and the ratio of statistical weights $g_u/g_l = 3/5$. The resulting hydrogen column density for a line center opacity of the [O I] 63 µm line of 1 agrees within 10% with the result from Eq. (3) for $10^3 \text{ cm}^{-3} \le n \le 10^5 \text{ cm}^{-3}$ and $20 \text{ K} \le 400 \text{ K}$. This is also consistent with the results of RADEX radiative transfer modeling (van der Tak et al. 2007).

We used the dust emission to estimate total hydrogen column densities at the positions observed in M 33 by fitting a single modified black body (MBB) to the 70, 100, 160 µm fluxes, while keeping $\beta = 1.5$ constant, deriving the dust temperature and dust mass surface density. To find the best-fitting SED for each pixel on the map, a $\chi^2$ function was minimized using the Levenberg-Marquardt algorithm (Xilouris et al. 2012). For a constant gas-to-dust ratio of 150 (Kramer et al. 2010), this implies typical hydrogen column densities of $2.4 \times 10^{21} \text{ cm}^{-2}$ per beam, implying moderate [O I] opacities of 0.7. Hydrogen column densities peak at $\sim 10^{22} \text{ cm}^{-2}$ ($A_V \sim 10 \text{ mag}$) (Fig. B.4), with corresponding [O I] optical depths of $\sim 3$.

The standard PDR models of Kaufman et al. (2006), Pound & Wolfie (2008), used below, assume homogeneous slabs of an optical extinction of $A_V = 10 \text{ mag}$. The models compute a simultaneous solution for the chemistry, the thermal balance,
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** [C II]/[O I](63 µm) ratio plotted against the ([C II]+[O I](63 µm)/TIR ratio for M 33. [C II], [O I], and TIR intensities are from positions with intensities above the local 3σ value. Colors correspond to the four different regions of the inner and outer galaxy, and above and below a TIR threshold, as described in Sect. 3.1. *Left*: observed ratios with typical 1σ errors. *Right*: ratios with corrected [O I] intensities, as described in Sect. 3.2.4. In addition, median and standard deviations are shown for the four regions (errorbars in color) (cf. Table 5).</figcaption></figure>

<table><caption>**Table 5.** Median [C II]/[O I] and ([C II]+[O I])/TIR ratios and rms values as observed and after correction for [O I] self absorption.</caption><thead><tr><th rowspan="3">Region</th><th rowspan="3">[C II]/[O I]</th><th rowspan="3">([C II]+[O I])/TIR</th><th colspan="4">Solutions of PDR models</th></tr><tr><th colspan="2">Moderate solution</th><th colspan="2">Low-FUV solution</th></tr><tr><th>%</th><th>n<br>(cm<sup>-3</sup>)</th><th>G<sub>0</sub></th><th>n<br>(cm<sup>-3</sup>)</th><th>G<sub>0</sub></th></tr></thead><tbody><tr><td colspan="7" style="text-align: left; font-weight: bold;">Observed ratios</td></tr><tr><td>All</td><td>4.51 ± 2.58</td><td>0.75 ± 0.31</td><td>60</td><td>20</td><td>10<sup>4</sup></td><td>1</td></tr><tr><td>Bright, inner</td><td>4.79 ± 2.43</td><td>0.68 ± 0.19</td><td>60</td><td>20</td><td>10<sup>4</sup></td><td>1</td></tr><tr><td>Bright, outer</td><td>3.71 ± 3.68</td><td>0.72 ± 0.19</td><td>30</td><td>20</td><td>10<sup>4</sup></td><td>1</td></tr><tr><td>Weak, inner</td><td>3.83 ± 2.54</td><td>0.91 ± 0.27</td><td>100</td><td>30</td><td>10<sup>4</sup></td><td>1</td></tr><tr><td>Weak, outer</td><td>3.01 ± 2.36</td><td>1.07 ± 0.54</td><td>300</td><td>60</td><td>10<sup>4</sup></td><td>1</td></tr><tr><td colspan="7" style="text-align: left; font-weight: bold;">Ratios with corrected [O I]</td></tr><tr><td>All</td><td>3.0 ± 2.1</td><td>0.82 ± 0.53</td><td>200</td><td>60</td><td>10<sup>4</sup></td><td>1.5</td></tr><tr><td>Bright, inner</td><td>3.06 ± 1.92</td><td>0.76 ± 0.2</td><td>200</td><td>60</td><td>10<sup>4</sup></td><td>1.5</td></tr><tr><td>Bright, outer</td><td>2.78 ± 2.69</td><td>0.77 ± 0.18</td><td>200</td><td>70</td><td>10<sup>4</sup></td><td>1.5</td></tr><tr><td>Weak, inner</td><td>3.27 ± 2.28</td><td>0.94 ± 0.89</td><td>300</td><td>60</td><td>10<sup>4</sup></td><td>1.5</td></tr><tr><td>Weak, outer</td><td>2.46 ± 2.02</td><td>1.12 ± 0.55</td><td>600</td><td>80</td><td>10<sup>4</sup></td><td>2.0</td></tr><tr><td>#1</td><td>4.59 ± 0.97</td><td>0.76 ± 0.25</td><td>900</td><td>20</td><td>5 × 10<sup>3</sup></td><td>3.0</td></tr></tbody></table>

**Notes.** Solutions of the PDR models of Pound & Wolfire (2008), Kaufman et al. (2006) are shown for the median ratio of all positions and for position #1 with one of the highest observed [C II] intensities $1.2 \times 10^{-4}$ erg s$^{-1}$ cm$^{-2}$ sr$^{-1}$.

and also the radiative transfer. As already said above, optical depths effects of the emergent line emission are taken into account. From these models, it is known that the [O I](63 µm) line emission becomes optically thick, with opacities of several, over the entire parameter space of n, G₀ sampled by the models. These models do, however, not consider absorption of the emission of warm background gas by colder foreground gas, as would be expected for GMCs which are internally heated by star-formation.

To estimate the possible effect of foreground absorption, we modeled the reduction in line center brightness temperatures assuming for each of the two source components beam filling factors of 1 and

$$T_R = J_v(T_{ex}) [1 - \exp[-\tau_0]] \quad (5)$$

with $J_v(T_{ex}) = 228 (\exp(\frac{228}{T_{ex}}) - 1)^{-1}$. The two source components are added, allowing for foreground absorption:

$$T_R = T_{R,fg} + \exp(-\tau_{0,fg}) T_{R,bg}. \quad (6)$$

Using in addition Eq. (4) for the opacity of the [O I] line, we furthermore assumed for both components FWHM linewidths of 7 km s$^{-1}$ and a density of $n = 10^4$ cm$^{-3}$. We assume that the background GMCs lie on average in the mid-plane of the galaxy and considered only half of the total column density in the following. In the absence of high-resolution spectra, which would allow for fitting free parameters to the observed line profiles (e.g., Graf et al. 1993), we assumed here ad-hoc that 80% of this half total column density is in the foreground components at 15 K, while 20% is in a background component at 200 K. This results in a reduction of the emerging line intensity by a factor of 3.3 for the highest total column densities of $10^{22}$ cm$^{-2}$ observed. The median corrections are $1.32 \pm 1.8$. The median corresponds to a total column density of $\sim 1.17 \times 10^{21}$ cm$^{-2}$. A smaller fraction of foreground gas would lead to lower reduction factors and vice versa. We take these estimates as rough, first estimates
<table><caption>**Table 6.** Assumptions to estimate the fraction of [C II] emission from the cold neutral medium (CNM).</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>N(H I, PDR)</td><td>3.25 × 10<sup>20</sup> cm<sup>-2</sup></td></tr><tr><td>X(C<sup>+</sup>)</td><td>5.9 × 10<sup>-5</sup></td></tr><tr><td>n<sub>CNM</sub></td><td>100 cm<sup>-3</sup></td></tr><tr><td>T<sub>CNM</sub></td><td>80 K</td></tr></tbody></table>

**Notes.** N(H I, PDR) is the estimated column density of atomic hydrogen at the surfaces of the PDRs. X(C$^+$) is the fractional abundance of C$^+$, n$_{CNM}$ and T$_{CNM}$ are density and temperature of the CNM.

and apply them to the data before comparing them with PDR models.

Figure 7 shows the observed intensity ratios and the ratios after correction of the [O I] intensities. After the correction, the [C II]/[O I](63 µm) ratio of all positions ranges between 0.13 and 13.7 with a median of 3.0 ± 2.1, and the ([C II]+[O I](63 µm))/TIR ratio ranges between 0.3 and 11.6% with a median of 0.8 ± 0.5%. The median ratios do not differ much between the different regions in M 33 (Table 5).

### 3.2.5. PDR models

To better understand the observed ratios of [C II]/[O I] and ([C II]+[O I])/TIR, corrected for [O I] self-absorption, we compare them with the predictions of PDR models (Kaufman et al. 2006; Pound & Wolfire 2008). These assume FUV illuminated slabs of an optical extinction of $A_V = 10$ mag and solar metallicities, for a range of local densities and FUV fields. The FUV field $G_0$ used in the PDR models is the local FUV field heating the slabs of dust and gas, leading to emission of the TIR, [C II], and [O I]. The ratio of the local $G_0$ from the models and the $G_{0,obs}$ from the observed TIR emission gives the beam filling factor of the emitting clouds. Beam filling factors in the observed ratios cancel out under the reasonable first order assumption that the three tracers ([C II], [O I], TIR) have the same filling factor. We ignore here that the excitation conditions for the [O I] line indicate a somewhat smaller beam filling factor than for the other two tracers.

For many of the observed ratios, the PDR models do not provide a unique $G_0$, n solution. For these positions, two solutions exist, one at high density and low FUV field, and another one at moderate values of $G_0$ and n. An example is shown in Fig. B.6, which exhibits contours of the average ratios of all observed positions (Table 5) as a function of density and FUV field of the Kaufman PDR models. There exist two best-fitting solutions, which are both consistent with the ratios of [C II], [O I], and TIR.

Solutions with moderate n, $G_0$ exist for the bulk of the ratios in M 33, after correcting for [O I] self-absorption (Fig. 8 upper panel). The ratios span the range of $log(G_0) < 3$ and $1 \le \log(n/cm^{-3}) \le 4.25$. The average ratios of all four regions can be explained by FUV-fields of $n \sim 2 \times 10^2$ cm$^{-3}$ and $G_0 \sim 60$ (Fig. B.6, Table 5). The observed FUV fields $G_{0,obs}$ range between 2 and 200, which indicates beam filling factors of about 1. The average ratios of the four regions do not differ much and lead to similar best fitting n, $G_0$ solutions. Even after correcting the [O I] intensities, some of the ratios still lie outside of the parameter space spanned by the moderate solutions: for example the ratios with ([C II]+[O I])/TIR ~ 0.8% and [C II]/[O I] > 4. Corrections of the [C II] intensities for contributions from the diffuse, ionized gas may, however, move all points

to lower [C II]/[O I] ratios and lower ([C II]+[O I])/TIR ratios, and into this parameter space.

The low-FUV solutions of the PDR models do, on the other hand, cover the entire range of ratios (Fig. 8 lower panel). The ratios lie in the regime $0 \le \log G_0 \le 0.75$ and $3.5 \le \log(n/cm^{-3}) \le 4.5$. The average of all ratios is best fit by $n \sim 10^4$ cm$^{-3}$, $G_0 \sim 1.5$ (Fig. B.6, Table 5). Given the range of observed $G_{0,obs}$, this would indicate beam filling factors $\Phi_{G0}$ of between ~1 for the regions with the lowest TIR and ~100 for the active arm regions with the highest observed TIR. While a fraction of the PDRs along the lines-of-sight may be represented by the low-FUV solution, this cannot be the sole solution. Such a high number of PDRs along the lines-of-sight seems not consistent with the observed peak optical extinctions of $A_V \sim 10$ mag, which resemble those of a single PDR model slab. However, only the surface columns of these models of about 1 mag, or of a hydrogen column density of ~$10^{21}$ cm$^{-2}$, emit [C II] (cf. Fig. 2 in Kaufman et al. 1999) and much of the remaining, deeper layers may not exist. High filling factors also seem not consistent with the relatively short lines-of-sight in M 33 with its moderate inclination of only 56°. Comparing the observed [C II] intensities with those predicted by the PDR models also gives high beam filling factors for the low-FUV solutions. Table 5 lists the ratios for a position with particularly strong [C II] intensity. The best fitting low-FUV solution for this position is $n = 5 \times 10^3$ cm$^{-3}$, $G_0 = 3$ (Fig. B.6). For this solution, the ratio of observed to modeled [C II] intensity is $\Phi_{CII} \sim 10$, which again seems difficult to reconcile with the models. For the moderate solution at the [C II] peak position, on the other hand, $\Phi_{CII}$ is ~1.

The degeneracy of the PDR model solutions has been discussed, for example by Hughes et al. (2015), Parkin et al. (2013), Kramer et al. (2005), Higdon et al. (2003), and Malhotra et al. (2001) for a variety of normal galaxies, who all prefer the moderate solution, often argueing that the low-FUV solution would require far too many PDR slabs along the lines-of-sight to be consistent with the observed intensities.

However, it is also clear that GMCs and their OB associations illuminating gas and dust show structure on a wide range of scales, and that the PDR model slab of constant volume and column density, illuminated by a constant FUV field is only a first order approximation. The inner parts of the GMCs and their OB associations and H II regions are illuminated by high FUV-fields, while the less dense, colder, outer, and more remote, diffuse parts of the GMCs are illuminated by much lower FUV-fields. We therefore cannot exclude that both solutions of the PDR models are at play along the lines-of-sight in M 33.

Detailed modeling of the spectral energy distributions (SEDs) may shed light on the relative fractions of gas and dust which exhibit the two solutions of the Kaufman PDR models. Aniano et al. (2012) define PDRs as those regions which are heated by starlight intensities which are a factor 100 or more higher than the ISRF in the solar neighborhood (Mathis et al. 1983)$^2$, $U > 100$, tracing massive star forming regions heated by OB stars. They use Draine & Li (2007) dust models, which assume a distribution of radiation fields between a minimum and a maximum value, fitting them to the observed SEDs between 3.6 µm and 500 µm at each position of the galaxies NGC 628 and NGC 6946, creating maps of the PDR fraction. They find that only a fraction of 12% to 14% of the TIR emission stems from such PDRs. More than 85% of the TIR emission in these two galaxies stems from the diffuse ISM heated by low

<footnote>$^2$ The Mathis et al. (1983) estimate for the ISRF has $G_0 = 1.14$ (e.g., Draine 2011).</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Diagnostic diagram of the [C II]/[O I] (63 µm) ratio plotted against the ([C II]+[O I] (63 µm)/TIR ratio for M 33, using the corrected [O I] intensities (see Sect. 3.2.4). Median and standard deviation of the four regions are shown in color, and of all data (star and large errorbars). Colors correspond to the four different regions as described in Sect. 3.1: the TIR bright, inner region (red), the bright, outer region (blue), the weak, inner region (green), and the weak, outer region (black). Superimposed is a grid of constant hydrogen nuclei density n (dashed contours) and FUV field strength G₀ (solid contours) from PDR models (Kaufman et al. 1999, 2006). We note that for many ratios there are two solutions for n, G₀. Upper panel: solutions with moderate n, G₀. Lower panel: low-FUV solutions.</figcaption></figure>

starlight intensities. The latter regions may resemble the regions in M 33 which are characterized by the low-FUV solution of the Kaufman PDR models.

### 3.3. Metallicities

To first order, one may expect that low metallicities lead to a rise of [C II]/TIR ratios, as a lowered dust abundance leads to deeper penetration lengths of UV photons, increasing the [C II] surface layers (Bolatto et al. 1999). This may increase [C II] intensities, while the thermal dust flux stays about constant (Israel et al. 1996). Indeed, low metallicities have been proposed to explain the observed high [C II]/TIR ratios found for example in dwarf galaxies (Cormier et al. 2015) and in the outskirts of spirals (e.g., Kapala et al. 2015, 2017).

M 33 exhibits about half solar metallicities together with shallow gradients of dropping metallicities with increasing galacto-centric radius (Magrini et al. 2010). Toribio San

Cipriano et al. (2016) measured radial abundance gradients of O/H and C/H from observations of H II regions in M 33. From optical recombination lines they find:

$$12 + \log(O/H) = 8.76 - (0.33 \pm 0.13) \times R/R_{25} \quad (7)$$

$$12 + \log(C/H) = 8.64 - (0.61 \pm 0.11) \times R/R_{25}. \quad (8)$$

The value of the 25th mag B-band isophotal radius ($R_{25}$) is 28 arcmin (6.84 kpc). This is the first work to discuss the C/H abundance gradient of M 33 based on more than just two H II regions (cf. the discussion in Toribio San Cipriano et al. 2016). The observed C/H gradient is about twice as steep as the O/H gradient, similar to what is found in other small galaxies with subsolar metallicities. The observed metallicity gradient may contribute to the observed indications of a rise of the [C II]/TIR ratio in the outer regions of M 33 (Fig. 5). However, the scatter of the observed [C II]/TIR ratio at a given radial distance in M 33 (Fig. 5) and hence about constant metallicity, indicates that other mechanisms can become important.
<figure><img src="image_10.png" /><figcaption>Fig. 9. [C II]/TIR intensity ratio versus TIR for the northern H II region BCLMP 691 in M 33. The blue drawn line shows the result of an unweighted linear least-squares fit to the TIR-bright data of BCLMP 691 varying only the y-intercept, while keeping the slope fixed at -1. The black dashed line shows the result of an unweighted linear least-squares fit to all M 33 data (cf. Fig. 5). A few typical 1σ errorbars are shown.</figcaption></figure>

Here, we take as one example the H II region BCLMP 691, which is located at a galacto-centric distance of 3.3 kpc in the far north of M 33 (Braine et al. 2012) at about constant metallicity (Eq. (8)). For a small region of ~0.9 arcmin² above the TIR threshold, BCLMP 691 exhibits a marked drop of the observed [C II]/TIR ratio from the outskirts of the H II region to its center and the peak of TIR emission by more than a factor of 3, from ~1% to 0.3% (Figs. 3 and 9).

For optically thin emission, the TIR equals

$$ \text{TIR} = \int B_v(T_d)\tau_d dv = \sum_d \kappa_0 \int B_v(T_d)(v/v_0)^\beta dv \quad (9) $$

with the total dust mass surface density $\Sigma_d$, the line-of-sight (l.o.s.) weighted averaged dust temperature $T_d$, the Planck function $B_v$, the opacity $\tau_d$, the average grain cross-section per gram $\kappa_0$ at frequency $v_0$, and the l.o.s. averaged dust emissivity index $\beta$, assuming optically thin emission.

Figure 9 shows the observed log [C II]/TIR vs. log TIR in BCLMP 691. A linear least squares-fit to the positions above the TIR threshold (Fig. 3), keeping the slope fixed to -1 results in a correlation coefficient of $r = -0.84$, indicating a good correlation. This slope is consistent with a constant slab of [C II] emission of $2.2 \times 10^{-5}$ erg s⁻¹ cm⁻² sr⁻¹, the dashed contour in the corresponding [C II] map (Fig. 3). The [C II] map indeed exhibits a roughly constant emission over the TIR-bright region,

indicating a constant column density and excitation temperature. MBB fits to the SEDs within the small area around the peak of TIR emission in BCLMP 691, where TIR intensities are above the threshold, show that dust temperatures stay fairly constant (23.2±1.9 K) with a median value of the estimated errors of $T_{\text{dust}}$ of 1.3 K. On the other hand, the dust mass surface densities vary by a factor of ~4 between 200 and $800 M_{\odot}$ beam⁻¹ (Fig. B.5), a variation which is a factor 3.5 larger than the median of the estimated errors, which is $170 M_{\odot}$ beam⁻¹. The observed steep drop of [C II]/TIR with TIR in a region of about constant metallicity is hence naturally explained.

## 4. Summary and conclusions

The emission lines of [C II] and [O I] (63 µm) were mapped along the major axis of the Local Group galaxy M 33. These maps have a width of ~370 pc and cover a region of 38 arcmin² at 50 pc resolution, allowing to resolve arm and interarm regions. The southern most region lies at 2 kpc galacto-centric distance and, located at the other side of the galaxy, the northern most region lies at 5.7 kpc distance. These maps much improve on the 1-dimensional [C II] cut at 280 pc resolution, which had been observed along the major axis of M 33 using ISO/LWS (Kramer et al. 2013).

We combined full-galaxy maps at 24 µm, 70 µm, 100 µm, and 160 µm to construct a map of the total infrared continuum emission (TIR), integrated between 1 µm and 1000 µm wavelength using the kernels provided by Aniano et al. (2011) and the weighting factors derived by Boquien et al. (2011). The observed range of TIR intensities translates to a range of FUV fluxes of $G_{0,\text{obs}} \sim 2$ to 200 in units of the average Galactic radiation field.

We find that the TIR and [C II] intensities are tightly correlated over two orders of magnitude. The average [C II]/TIR ratio of 0.64 ± 0.23% is not significantly higher than the average [C II]/TIR ratio of 0.48 ± 0.21% found in a sample of 54 nearby galaxies by Smith et al. (2017). [C II]/TIR ratios observed in the two northern most regions at 4.5 and 5.5 galacto-centric distances show increasing average ratios of 0.8 and 1.3, respectively. The large-scale variation of the [C II]/TIR ratios is consistent with the ISO/LWS observations. The resolution and extent of the PACS map allows to distinguish between the diffuse interarm regions of the inner and outer galaxy, and the arm regions. The [C II]/TIR ratio averaged over bins of 0.5 dex decreases with increasing TIR from 1.1 ± 0.4% for regions with weak TIR emission to 0.5 ± 0.1% for the arm-regions with highest TIR, while the scatter of binned averages of the [C II]/TIR ratios decreases as well.

The drop of [C II]/TIR ratios toward sites of massive star forming regions, where TIR peaks, is most clearly visible in the [C II]/TIR map of one of the northern H II regions, BCLMP 691 at 3.3 kpc galacto-centric distance, for positions where TIR is bright. In this case, the drop of [C II]/TIR ratios is consistent with a [C II] surface layer of constant intensity, which is independent of TIR. The rise of TIR is caused by a rise of dust mass surface densities, for about constant dust temperatures and emissivities, as is confirmed by the results of modified black body (MBB) fits. For this H II region, the observed steep drop of [C II]/TIR with TIR in a region of about constant metallicity is hence naturally explained. However, this does not rule out that on larger scales the drop of metallicities with galacto-centric distance observed in M 33 (Magrini et al. 2010) is decisive to explain the observed drop of [C II]/TIR with TIR on these scales.

In M 33, [C II] intensities are stronger than those of the [O I] (63 µm) line at almost all positions. The [C II]/[O I] ratio
varies between ~0.2 and 20, with a mean of 4.5 ± 2.6. At first glance, the [O I] and [C II] maps of the inner galaxy resemble each other. However, closer inspection shows that [O I] and [C II] are only correlated in the TIR-brightest regions while there is only little correlation of both tracers in the TIR weak regions. The observed ([C II]+[O I])/TIR ratio varies between ~0.3% and 3%, with an average of 0.75 ± 0.3%.

Data of atomic hydrogen, taken with the VLA at the same resolution as the [C II] data, were used to estimate the contribution of the cold neutral medium (CNM) to the observed [C II] emission, following the approach of K2013. As anticipated from this work, the averaged contribution is only ~10%. We did not correct the [C II] emission for this minor contribution. Modified black bodies were fit to the continuum emission at 70, 100, 160 µm to derive dust temperatures and surface densities, which were used to estimate total gas and O column densities, and the optical depths of the [O I](63 µm) line. A simple model of cold foreground and warm background gas then allowed us to estimate that for the highest H I column densities of 10<sup>22</sup> cm<sup>-2</sup>, [O I](63 µm) line intensities are reduced by a factor of 3.3, while the median correction factor is 1.3 ± 1.8. The [O I] data were corrected for this effect. An additional correction of the [C II] data, for diffuse, ionized gas, was not attempted here.

The observed [C II]/[O I] and ([C II]+[O I])/TIR ratios were corrected for [O I] foreground absorption, and then compared with standard PDR models (Kaufman et al. 2006). Averages of all observed ratios are similar to the averages of the four individual regions. They all show two solutions of the PDR models, a moderate solution with *n* ~ 2 × 10<sup>2</sup> cm<sup>-3</sup>, G<sub>0</sub> ~ 60, and a low-FUV solution with *n* ~ 10<sup>4</sup> cm<sup>-3</sup>, G<sub>0</sub> ~ 1.5. The bulk of the observed positions can be modeled by a moderate solution. This solution implies low beam filling factors of ~1. The low-FUV solution, on the other hand, cannot be the sole solution for all gas along the lines of sight, as it would imply very high beam filling factors ≫1, which are inconsistent with the observed FUV fields, the [C II] intensities, and the total column densities.

**Acknowledgements.** We thank the anonymous referee for insightful comments which helped to improve the paper, Mark Wolfire and Alessandra Contursi for helpful discussion, and the NHSC team at IPAC for their support in the data reduction process. M.R. and S.V. acknowledge support by the research projects AYA2014-53506-P and AYA2017-84897-P from the Spanish Ministerio de Economía y Competitividad, from the European Regional Development Funds (FEDER) and the Junta de Andalucía (Spain) grants FQM108. This study has been partially financed by the Consejería de Conocimiento, Investigación y Universidad, Junta de Andalucía and European Regional Development Fund (ERDF), ref. SOMM17/6105/UGR. FST thanks the Spanish Ministry of Economy and Competitiveness (MINECO) for support under grant number AYA2016-76219-P.

## References

Abdullah, A., Brandl, B. R., Groves, B., et al. 2017, ApJ, 842, 4

Aniano, G., Draine, B. T., Gordon, K. D., & Sandstrom, K. 2011, PASP, 123, 1218

Aniano, G., Draine, B. T., Calzetti, D., et al. 2012, ApJ, 756, 138

Bolatto, A. D., Jackson, J. M., & Ingalls, J. G. 1999, ApJ, 513, 275

Boquien, M., Calzetti, D., Combes, F., et al. 2011, AJ, 142, 111

Boquien, M., Calzetti, D., Aalto, S., et al. 2015, A&A, 578, A8

Braine, J., Gratier, P., Kramer, C., et al. 2012, A&A, 544, A55

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

Corbelli, E. 2003, MNRAS, 342, 199

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

Crawford, M. K., Genzel, R., Townes, C. H., & Watson, D. M. 1985, ApJ, 291, 755

Croxall, K. V., Smith, J. D., Wolfire, M. G., et al. 2012, ApJ, 747, 81

Croxall, K. V., Smith, J. D., Pellegrini, E., et al. 2017, ApJ, 845, 96

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton: Princeton University Press)

Draine, B. T., & Li, A. 2007, ApJ, 657, 810

Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172

Druard, C., Braine, J., Schuster, K. F., et al. 2014, A&A, 567, A118

Galametz, M., Kennicutt, R. C., Calzetti, D., et al. 2013, MNRAS, 431, 1956

Galleti, S., Bellazzini, M., & Ferraro, F. R. 2004, A&A, 423, 925

Gerin, M., Ruaud, M., Goicoechea, J. R., et al. 2015, A&A, 573, A30

Graf, U. U., Eckart, A., Genzel, R., et al. 1993, ApJ, 405, 249

Gratier, P., Braine, J., Rodriguez-Fernandez, N. J., et al. 2010, A&A, 522, A3

Gratier, P., Braine, J., Schuster, K., et al. 2017, A&A, 600, A27

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

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

Higdon, S. J. U., Higdon, J. L., van der Hulst, J. M., & Stacey, G. J. 2003, ApJ, 592, 161

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

Hunter, I., Dufton, P. L., Smartt, S. J., et al. 2007, A&A, 466, 277

Israel, F. P., Maloney, P. R., Geis, N., et al. 1996, ApJ, 465, 738

Israel, F. P., Güsten, R., Meijerink, R., Requena-Torres, M. A., & Stutzki, J. 2017, A&A, 599, A53

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

Kapala, M. J., Groves, B., Sandstrom, K., et al. 2017, ApJ, 842, 128

Karska, A., Herpin, F., Bruderer, S., et al. 2014, A&A, 562, A45

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

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

Kennicutt, R. C., Jr., Hao, C. N., Calzetti, D., et al. 2009, ApJ, 703, 1672

Kramer, C., Mookerjea, B., Bayet, E., et al. 2005, A&A, 441, 961

Kramer, C., Buchbender, C., Xilouris, E. M., et al. 2010, A&A, 518, L67

Kramer, C., Abreu-Vicente, J., García-Burillo, S., et al. 2013, A&A, 553, A114 (K2013)

Leurini, S., Wyrowski, F., Wiesemeyer, H., et al. 2015, A&A, 584, A70

Lin, Z., Hu, N., Kong, X., et al. 2017, ApJ, 842, 97

Lis, D. C., Keene, J., Phillips, T. G., et al. 2001, ApJ, 561, 823

Liseau, R., Justtanont, K., & Tielens, A. G. G. M. 2006, A&A, 446, 561

Magrini, L., Stanghellini, L., Corbelli, E., Galli, D., & Villaver, E. 2010, A&A, 512, A63

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

Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 500, 259

Mookerjea, B., Kramer, C., Buchbender, C., et al. 2011, A&A, 532, A152

Mookerjea, B., Israel, F., Kramer, C., et al. 2016, A&A, 586, A37

Müller, T., Okumura, K., & Klaas, U. 2011, PACS Photometer Passbands and Colour Correction Factors for Various Source SEDs, PICC-ME-TN-038, Tech. rep.

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

Ott, S. 2010, in Astronomical Data Analysis Software and Systems XIX, eds. Y. Mizumoto, K. I. Morita, & M. Ohishi, ASP Conf. Ser., 434, 139

Parkin, T. J., Wilson, C. D., Schirm, M. R. P., et al. 2013, ApJ, 776, 65

Parkin, T. J., Wilson, C. D., Schirm, M. R. P., et al. 2014, ApJ, 787, 16

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

Pineda, J. L., Langer, W. D., & Goldsmith, P. F. 2014, A&A, 570, A121

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Pound, M. W., & Wolfire, M. G. 2008, in Astronomical Data Analysis Software and Systems XVII, eds. R. W. Argyle, P. S. Bunclark, & J. R. Lewis, ASP Conf. Ser., 394, 654

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

Schneider, N., Röllig, M., Simon, R., et al. 2018, A&A, 617, A45

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

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

Tabatabaei, F. S., Braine, J., Xilouris, E. M., et al. 2014, A&A, 561, A95

Tielens, A. G. G. M. 2008, ARA&A, 46, 289

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

Timmermann, R., Bertoldi, F., Wright, C. M., et al. 1996, A&A, 315, L281

Toribio San Cipriano, L., García-Rojas, J., Esteban, C., Bresolin, F., & Peimbert, M. 2016, MNRAS, 458, 1866

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

Vastel, C., Polehampton, E. T., Baluteau, J. P., et al. 2002, ApJ, 581, 315

Verley, S., Corbelli, E., Giovanardi, C., & Hunt, L. K. 2009, A&A, 493, 453

Warner, P. J., Wright, M. C. H., & Baldwin, J. E. 1973, MNRAS, 163, 163

Xilouris, E. M., Tabatabaei, F. S., Boquien, M., et al. 2012, A&A, 543, A74

Zaritsky, D., Elston, R., & Hill, J. M. 1989, AJ, 97, 97

<footnote n="1">¹ Institut de Radioastronomie Millimétrique (IRAM), 300 rue de la Piscine, 38406 Saint Martin d'Hères, France
e-mail: kramer@iram.fr</footnote>

<footnote n="2">² IRAM, Av. Divina Pastora 7, 18012 Granada, Spain</footnote>

<footnote n="3">³ Cornell University, Ithaca, NY 14852, USA</footnote>
⁴ Frankfurter Allgemeine Zeitung, Hellerhofstraße 2-4, 60327 Frankfurt am Main, Germany

⁵ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

⁶ Argelander Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany

⁷ Unidad de Astronomía, Universidad de Antofagasta, Av. Angamos 601, Antofagasta 1270300, Chile

⁸ Laboratoire d'Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, Allée Georoy Saint-Hilaire, 33615 Pessac, France

⁹ KOSMA, I. Physikalisches Institut, Universität zu Köln, Zülpicher Straße 77, 50937 Köln, Germany

¹⁰ Observatoire de Paris, LERMA, College de France, CNRS, PSL Univ., Sorbonne University, UPMC, Paris, France

¹¹ Max Planck Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

¹² Department of Astronomy, King Abdulaziz University, PO Box 80203, 21589 Jeddah, Saudi Arabia

¹³ Instituto de Astrofísica de Andalucía (IAA-CSIC), CAHA, Glorieta de la Astronomía s/n, 18008 Granada, Spain

¹⁴ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands

¹⁵ Dept. Física Teórica y del Cosmos, Universidad de Granada, 18012 Granada, Spain

¹⁶ Institute for Research in Fundamental Sciences-IPM, Larak Garden, 19395-5531 Tehran, Iran

¹⁷ SRON Netherlands Institute for Space Research, Landleven 12, 9747 AD Groningen, The Netherlands

¹⁸ Kapteyn Astronomical Institute, University of Groningen, Groningen, The Netherlands

¹⁹ Instituto de Astrofísica de Canarias, Vía L'actea S/N, 38205 La Laguna, Spain

²⁰ Instituto Universitario Carlos I de Física Teórica y Computacional, Facultad de Ciencias, Universidad de Granada, 18071 Granada, Spain

²¹ Institute for Astronomy, Astrophysics, Space Applications & Remote Sensing, National Observatory of Athens, P. Penteli, 15236 Athens, Greece
## Appendix A: Observing parameters of PACS spectroscopy observations

<table><caption>**Table A.1.** Observing parameters of all 22 fields mapped in M 33 with Herschel/PACS in line spectroscopy mode.</caption><thead><tr><th>Target-name</th><th>RA<br>(eq. 2000)</th><th>Dec</th><th>Integr. time<br>(s)</th><th>ObsID</th><th>Obs. day</th></tr></thead><tbody><tr><td>M 33-1350N</td><td>1<sup>h</sup>34<sup>m</sup>31<sup>s</sup>15</td><td>+31<sup>d</sup>00<sup>m</sup>23<sup>s</sup>90</td><td>6777</td><td>1342213044</td><td>616</td></tr><tr><td>M 33-1150N</td><td>1<sup>h</sup>34<sup>m</sup>25<sup>s</sup>19</td><td>+30<sup>d</sup>57<sup>m</sup>19<sup>s</sup>20</td><td>6779</td><td>1342213043</td><td>616</td></tr><tr><td>M 33-800N</td><td>1<sup>h</sup>34<sup>m</sup>14<sup>s</sup>80</td><td>+30<sup>d</sup>51<sup>m</sup>55<sup>s</sup>80</td><td>6781</td><td>1342213045</td><td>616</td></tr><tr><td>M 33-600N</td><td>1<sup>h</sup>34<sup>m</sup>08<sup>s</sup>80</td><td>+30<sup>d</sup>48<sup>m</sup>51<sup>s</sup>00</td><td>6779</td><td>1342213046</td><td>616</td></tr><tr><td>M 33-500N</td><td>1<sup>h</sup>34<sup>m</sup>06<sup>s</sup>12</td><td>+30<sup>d</sup>47<sup>m</sup>19<sup>s</sup>26</td><td>3876</td><td>1342189072</td><td>238</td></tr><tr><td>M 33-400N</td><td>1<sup>h</sup>34<sup>m</sup>02<sup>s</sup>90</td><td>+30<sup>d</sup>45<sup>m</sup>46<sup>s</sup>20</td><td>6777</td><td>1342212583</td><td>609</td></tr><tr><td>M 33-350N</td><td>1<sup>h</sup>34<sup>m</sup>01<sup>s</sup>42</td><td>+30<sup>d</sup>45<sup>m</sup>00<sup>s</sup>10</td><td>6777</td><td>1342212582</td><td>609</td></tr><tr><td>M 33-250N</td><td>1<sup>h</sup>33<sup>m</sup>58<sup>s</sup>45</td><td>+30<sup>d</sup>43<sup>m</sup>27<sup>s</sup>70</td><td>6777</td><td>1342212579</td><td>609</td></tr><tr><td>M 33-200N</td><td>1<sup>h</sup>33<sup>m</sup>57<sup>s</sup>00</td><td>+30<sup>d</sup>42<sup>m</sup>41<sup>s</sup>50</td><td>6777</td><td>1342212540</td><td>608</td></tr><tr><td>M 33-150N</td><td>1<sup>h</sup>33<sup>m</sup>55<sup>s</sup>48</td><td>+30<sup>d</sup>41<sup>m</sup>55<sup>s</sup>30</td><td>6775</td><td>1342212538</td><td>608</td></tr><tr><td>M 33-100N</td><td>1<sup>h</sup>33<sup>m</sup>54<sup>s</sup>00</td><td>+30<sup>d</sup>41<sup>m</sup>09<sup>s</sup>10</td><td>6775</td><td>1342212536</td><td>608</td></tr><tr><td>M 33-50N</td><td>1<sup>h</sup>33<sup>m</sup>52<sup>s</sup>51</td><td>+30<sup>d</sup>40<sup>m</sup>22<sup>s</sup>90</td><td>6775</td><td>1342212534</td><td>608</td></tr><tr><td>M 33</td><td>1<sup>h</sup>33<sup>m</sup>51<sup>s</sup>02</td><td>+30<sup>d</sup>39<sup>m</sup>36<sup>s</sup>70</td><td>6775</td><td>1342213047</td><td>616</td></tr><tr><td>M 33-50S</td><td>1<sup>h</sup>33<sup>m</sup>49<sup>s</sup>54</td><td>+30<sup>d</sup>38<sup>m</sup>50<sup>s</sup>60</td><td>6775</td><td>1342212535</td><td>608</td></tr><tr><td>M 33-100S</td><td>1<sup>h</sup>33<sup>m</sup>48<sup>s</sup>10</td><td>+30<sup>d</sup>38<sup>m</sup>04<sup>s</sup>30</td><td>6775</td><td>1342212537</td><td>608</td></tr><tr><td>M 33-150S</td><td>1<sup>h</sup>33<sup>m</sup>46<sup>s</sup>58</td><td>+30<sup>d</sup>37<sup>m</sup>18<sup>s</sup>20</td><td>6773</td><td>1342212539</td><td>608</td></tr><tr><td>M 33-200S</td><td>1<sup>h</sup>33<sup>m</sup>45<sup>s</sup>10</td><td>+30<sup>d</sup>36<sup>m</sup>31<sup>s</sup>90</td><td>6773</td><td>1342212578</td><td>609</td></tr><tr><td>M 33-250S</td><td>1<sup>h</sup>33<sup>m</sup>43<sup>s</sup>61</td><td>+30<sup>d</sup>35<sup>m</sup>45<sup>s</sup>80</td><td>6773</td><td>1342212580</td><td>609</td></tr><tr><td>M 33-300S</td><td>1<sup>h</sup>33<sup>m</sup>42<sup>s</sup>10</td><td>+30<sup>d</sup>34<sup>m</sup>59<sup>s</sup>50</td><td>6773</td><td>1342212581</td><td>609</td></tr><tr><td>M 33-400S</td><td>1<sup>h</sup>33<sup>m</sup>39<sup>s</sup>20</td><td>+30<sup>d</sup>33<sup>m</sup>27<sup>s</sup>20</td><td>6773</td><td>1342212584</td><td>609</td></tr><tr><td>M 33-450S</td><td>1<sup>h</sup>33<sup>m</sup>37<sup>s</sup>68</td><td>+30<sup>d</sup>32<sup>m</sup>41<sup>s</sup>00</td><td>6773</td><td>1342212585</td><td>609</td></tr><tr><td>M 33-500S</td><td>1<sup>h</sup>33<sup>m</sup>36<sup>s</sup>20</td><td>+30<sup>d</sup>31<sup>m</sup>54<sup>s</sup>80</td><td>6773</td><td>1342212586</td><td>609</td></tr></tbody></table>
**Appendix B: Additional figures**

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

**Fig. B.1.** PACS spectra of [C II] (blue) and [O I](63 µm) (red) at selected positions in M 33. All spectra are shown at their native angular resolution, before convolution to 12" resolution. The velocities have been corrected for the systemic velocity of M 33 of –180 km s⁻¹. Dashed, vertical lines mark the residual velocity of the H I line (Warner et al. 1973). Spectra are shown after subtraction of polynomial baselines of up to 3rd order. The two red (blue) fields mark the velocity range, which was used to calculate the baseline rms for the [O I] ([C II]) line. Baselines at more extreme velocities are not shown. [C II] integrated intensities were calculated over the inner range between the blue fields (±538 km s⁻¹), and the integrated intensities for the [O I] line were calculated over the spectral range (in white) between the red fields. The double-peak structure of the weak [O I] line seen in the right spectrum of the *lower row* is attributed to baseline noise.
<figure><img src="image_12.png" /><figcaption>**Fig. B.2.** Fraction of [C II] emission from the cold neutral medium (CNM), estimated from H I emission.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. B.3.** TIR versus H I intensities in M 33. Intensities are given in erg s⁻¹ cm⁻² sr⁻¹. Typical 1σ errorbars are shown. For H I we assume an calibration error of 15% (cf. K2013). The drawn line shows the result of an unweighted linear least-squares fit to log TIR = b + m × log HI, which gives a slope m = 0.48 ± 0.11, y-intercept b = 3.07 ± 1.32, and correlation coefficient r = 0.45.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. B.4.** Total hydrogen column densities derived from fits of modified black bodies to the 70 µm, 100 µm, and 160 µm data at each position, assuming a constant gas-to-dust ratio (see Sect. 3.2.4).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. B.5.** Dust mass surface density vs. dust temperature, derived from fits of modified black bodies, for the TIR bright positions in the H II region BCLMP 691. The thick dashed line shows the median dust temperature and thin dashed lines correspond to the rms scatter.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. B.6.** Contours of the ratios of [O I]/[C II], and ([C II]+[O I])/TIR, with [O I] corrected (Table 5), as a function of cloud density *n* and incident FUV flux *G*₀ using the PDR model of Kaufman et al. (2006). **Left:** median ratios of all positions. **Right:** ratios at a [C II] peak position.</figcaption></figure>