Gas and dust cooling along the major axis of M33 (HerM33es) -- Herschel/PACS [CII] and [OI] observations

M33 is a gas rich spiral galaxy of the Local Group. We investigate the relationship between the two major gas cooling lines and the total infrared (TIR) dust continuum. We mapped the emission of gas and dust in M33 using the far-infrared lines of [CII] and [OI](63um) and the TIR. The line maps were observed with Herschel/PACS. These maps have 50pc resolution and form a ~370pc wide stripe along its major axis covering the sites of bright HII regions, but also more quiescent arm and inter-arm regions from the southern arm at 2kpc galacto-centric distance to the south out to 5.7kpc distance to the north. Full-galaxy maps of the continuum emission at 24um from Spitzer/MIPS, and at 70um, 100um, and 160um from PACS were combined to obtain a map of the TIR. TIR and [CII] intensities are correlated over more than two orders of magnitude. The range of TIR translates to a range of far ultraviolet (FUV) emission of G0,obs~2 to 200 in units of the average Galactic radiation field. The binned [CII]/TIR ratio drops with rising TIR, with large, but decreasing scatter. Fits of modified black bodies (MBBs) 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 [OI] data before comparing it with models of photon dominated regions (PDRs). Most of the ratios of [CII]/[OI] and ([CII]+[OI])/TIR are consistent with two model solutions. The median ratios are consistent with one solution at n~2x10^2 cm-3, G0~60, and and a second low-FUV solution at n~10^4 cm-3, G0~1.5. 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.


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 COdark molecular gas, 25% from atomic gas, and another 20% from diffuse, ionized gas (Pineda et al. 2014). These contribu-Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.
Maps of TIR, [C ii], [O i] shown in Figures 2,3 are available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/. tions 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 gas. Kapala et al. (2015) find that the [C ii]/TIR ratio increases Article number, page 1 of 19 arXiv:2005.03683v2 [astro-ph.GA] 12 May 2020 A&A proofs: manuscript no. herm33es-pacsspec 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 photoelectric 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 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 [C ii], [O i] (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 [C ii] 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 [C ii] 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.  4, 7.4, 16.0, 33.6, 73.0 in units of 10 −6 erg s −1 cm −2 sr −1 . A polygon marks the outer edges of the regions mapped in [C ii]. Coordinates are in R.A. 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 delineate galacto-centric distances of 2, 4, and 6 kpc.

PACS spectroscopy
The region mapped with Herschel in [C ii] and [O i] 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,2,3). The total area covered is 38.5 arcmin 2 . Data were taken with PACS (Poglitsch et al. 2010) on Herschel (Pilbratt et al. 2010). [N ii](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 1 h 35 m 28.87 s , +31 • 46 1.92 . For the other observations, an off position at 1 h 30 m 20.90 s , +30 • 38 25.44 was used. The only data taken in wavelength 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 stepsize 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 −1 and 240 km s −1 , 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. 2014Gratier 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 −7 erg s −1 cm −2 sr −1 for [C ii], (3 ± 7) · 10 −7 erg s −1 cm −2 sr −1 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.

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 MIPS 70 = mc/cc · S PACS 70 . Next, the TIR was calculated using weighting factors from Equation 1 and Table 1 of Boquien et al. (2011), and using all data larger than zero: log S TIR = 0.220 · log S 24 + 0.202 · log S MIPS 70 + 0.243 · log S 100 + 0.311 · log S 160 + 1.361.
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 ± 1.0) · 10 −5 erg s −1 cm −2 sr −1 . For each position, the local 1 σ 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.

[C ii] emission and the total infrared continuum
The maps of TIR and [C ii] (Figs. 1, 2, 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, 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 · 10 −3 erg s −1 cm −2 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,obs , which is given in units of  the average radiation field in the solar neighborhood, the Habing field of 1.6 · 10 −3 erg cm −2 s −1 (Habing 1968): 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 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 ( 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 (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 · 10 −3 erg s −1 cm −2 sr −1 , but surrounded by an  Table 1). 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 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σ 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.
log  (Fig. 7). Table 5 lists the average values and the scatter for the four regions defined in Section 3. PDR models may allow one to estimate the local excitation conditions, the local densities and the local FUV fields. In order  (Table 3).
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.

[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.

[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 importance of the ionized gas (Oberst et al. 2006;Croxall et al. 2012Croxall et al. , 2017. Parkin et al. (2013aParkin et al. ( , 2014 (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 · 10 −4 erg s −1 cm −2 sr −1 .

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 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 N(H i,PDR) 3.25 · 10 20 cm −2 X(C + ) 5.9 · 10 −5 n CNM 100 cm −3 T CNM 80 K emission stems from PDR regions, as the CNM contribution to the [C ii] emission is also low.

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  in cm −2 , with the FWHM line width in km s −1 (Vastel et al. 2002;Liseau et al. 2006). For a typical line width of 7 km s −1 (see above) and an average oxygen abundance of 4 · 10 −4 in M 33 (see below), the line center opacity exceeds 1 for N([O i]) = 1.4 · 10 18 cm −2 and N(H)> 3 · 10 21 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: with the Einstein A-coefficient A ul = 8.46 · 10 −5 s −1 , the critical density n cr = 4.7 · 10 5 cm −3 for collisions with H-atoms, hν/k B = 228 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 cm −3 ≤ n ≤ 10 5 cm −3 and 20 K ≤ 400 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 β = 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 χ 2 function was minimized using the Levenberg-Marquardt algorithm (Xilouris et al. 2012). For a constant gasto-dust ratio of 150 (Kramer et al. 2010), this implies typical hydrogen column densities of 2.4·10 21 cm −2 per beam, implying moderate [O i] opacities of 0.7. Hydrogen column densities peak at ∼ 10 22 cm −2 (A V ∼ 10 mag) (Fig. B.4), with corresponding The standard PDR models of Kaufman et al. (2006); Pound & Wolfire (2008), used below, assume homogeneous slabs of an optical extinction of A V = 10 mag. The models compute a simultaneous solution for the chemistry, the thermal balance, 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 0 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 with J ν (T ex ) = 228 (exp( 228 T ex ) − 1) −1 . The two source components are added, allowing for foreground absorption: 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 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 Sec. 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 0 (solid contours) from PDR models (Kaufman et al. 1999(Kaufman et al. , 2006. We note that for many ratios there are two solutions for n, G 0 . Upper panel: Solutions with moderate n, G 0 . Lower panel: Low-FUV solutions. (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 ± 1.8. The median corresponds to a total column density of ∼ 1.17·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 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).

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 exinction 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 indi-cate 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 Figure 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 < ∼ log(n/cm −3 ) < ∼ 4.25. The average ratios of all four regions can be explained by FUV-fields of n ∼ 2 · 10 2 cm −3 and G 0 ∼ 60 (Fig. B.6, Table 5). The observed FUV fields G 0,obs range be- 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 < ∼ log G 0 < ∼ 0.75 and 3.5 < ∼ log(n/cm −3 ) < ∼ 4.5. The average of all ratios is best fit by n ∼ 10 4 cm −3 , G 0 ∼ 1.5 (Fig.B.6, Table 5). Given the range of observed G 0,obs , this would indicate beam filling factors Φ 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 ∼ 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 · 10 3 cm −3 , G 0 = 3 (Fig.B.6). For this solution, the ratio of observed to modeled [C ii] intensity is Φ CII ∼ 10, which again seems difficult to reconcile with the models. For the moderate solution at the [C ii] peak position, on the other hand, Φ CII is ∼ 1.
The degeneracy of the PDR model solutions has been discussed, for example by Hughes et al. (2015); Parkin et al. (2013b); Kramer et al. (2005); Higdon et al. (2003), andMalhotra 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 FUVfields, while the less dense, colder, outer, and more remote, diffuse parts of the GMCs are illuminated by much lower FUVfields. 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 reagions 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 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.

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. 2015Kapala et al. , 2017.
M 33 exhibits about half solar metallicities together with shallow gradients of dropping metallicities with increasing galacto-centric radius (Magrini et al. 2010 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.
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 2 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,9).
For optically thin emission, the TIR equals with the total dust mass surface density Σ d , the line-of-sight (l.o.s.) weighted averaged dust temperature T d , the Planck function B ν , the opacity τ d , the average grain cross-section per gram κ 0 at frequency ν 0 , and the l.o.s. averaged dust emissivity index β, assuming optically thin emission. Figure 9 shows the observed log[C ii]/TIR vs. logTIR 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 · 10 −5 erg s −1 cm −2 sr −1 , 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 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 /beam (Fig. B.5), a variation which is a factor 3.5 larger than the median of the estimated errors, which is 170 M /beam. The observed steep drop of [C ii]/TIR with TIR in a region of about constant metallicity is hence naturally explained.

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 2 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,obs ∼ 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. 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  (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 2 cm −3 , G 0 ∼ 60, and a low-FUV solution with n ∼ 10 4 cm −3 , G 0 ∼ 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.  (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 −1 ), 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.