© ESO, 2014

1. Introduction

Starburst galaxies are nearby laboratories that allow us to study intense star formation. Studying the heating and cooling mechanisms in these galaxies gives us insight into which excitation or feedback mechanisms are dominant in fueling starbursts. Molecular emission lines play a crucial role in the cooling of excited gas, and with the Herschel Space Observatory we are able to observe the rich molecular spectrum of the nucleus of NGC 253.

Carbon monoxide (CO) is one of the most abundant molecules after H₂, and a good and easily observable tracer of the condition of the molecular gas in the interstellar medium (ISM) of these galaxies. CO in the ISM is mostly located in molecular clouds and photon dominated regions (PDRs), where the radiation can penetrate the cloud and excite the gas. Tielens & Hollenbach (1985) created models of PDRs that predict the intensities of atoms and molecules in the PDR as a function of density, radiation environment, and column density. These have been expanded to include models for X-ray dominated regions (XDRs; Meijerink & Spaans 2005; Maloney et al. 1996), PDRs with enhanced cosmic ray ionization rates (Meijerink et al. 2006), and PDRs with additional mechanical heating taken into account (Kazandjian et al. 2013).

The source NGC 253 is a nearby, D_L = 2.5 Mpc (Davidge & Pritchet (1990), 12 pc/”), edge-on barred spiral galaxy (Scoville et al. 1985). The central Kiloparsec of NGC 253 is considered an archetypal starburst nucleus (L_IR ~ 2 × 10¹⁰L_⊙), which is heavily obscured at optical wavelengths by dust lanes (Prada et al. 1996). However, in the far-infrared and submillimeter wavelength regimes, NGC 253 exhibits extremely bright molecular line transitions (Henkel et al. 1991), originating from large molecular clouds in the nuclear region (Israel et al. 1995; Mauersberger et al. 1996; Houghton et al. 1997; Bradford et al. 2003; Martín et al. 2009). This gas also appears to be highly excited. Observations of HCO⁺ and HCN suggest that at least some of the gas has densities greater than 10⁴ cm^-3 and temperatures over 100 K (Paglione et al. 1995, 1997).

There have been many studies that have attempted to derive the excitation mechanism in NGC 253. The warm, excited molecular gas phase, often associated with PDRs, excites the mid- to high-J CO transitions (J> 4). The near-infrared H₂ emission lines shows that PDRs are an important excitation mechanism (Rosenberg et al. 2013). However hot H₂ gas only traces the very edges of molecular clouds. In order to study the excitation of the bulk of the molecular gas, we can use CO as a probe. Bradford et al. (2003) have observed ¹²CO up to J = 6 − 5 along with ¹³CO up to J = 3 − 2 and derive a kinetic temperature of 120 K and an H₂ density of 4.5 × 10⁴ cm^-3 for the warm phase. However, they suggest that the CO is excited by cosmic rays, and not only by PDRs. Martín et al. (2006) found that the chemistry and heating of NGC 253 is dominated by large-scale, low-velocity shocks. The presence of shocked molecular material is indeed evident through the presence of widespread SiO emission throughout the nuclear region (García-Burillo et al. 2000). In addition, Hailey-Dunsheath et al. (2008) detected the first extragalactic ¹³CO J = 6 − 5 transition, and using this determined that shocks are the dominant excitation mechanism in the nuclear region of NGC 253. Martín et al. (2009) suggest that, although NGC 253 is dominated by shock chemistry, PDRs also play a crucial role in the chemistry, since there are very high abundances of PDR tracing molecules, namely HCO⁺, CO⁺.

To better constrain the dominant excitation mechanism in NGC 253, we present the full ¹²CO ladder up to J = 13 − 12 and the ¹³CO up to J = 6 − 5 as observed with the Herschel Space Observatory. We combine this with observations of HNC and HCN transitions and apply these observations to models of PDRs, XDRs, enhanced cosmic ray PDRs, and enhanced mechanical heating PDRs to model the excitation directly. These observations were taken as part of the Guaranteed Key Program Herschel EXtra GALactic (HEXGAL, PI: R. Güsten). In Sect. 2, we will describe the observations and data reduction techniques. In Sect. 3, we will present our spectra and line fluxes. Using models of CO emission, in Sect. 4 we will constrain the parameters of the molecular gas phases and introduce a methodology for understanding the degeneracies of the models. To constrain the densest phase of the ISM, we use the HCO⁺ and HCN in Sect. 5 to determine the excitation mechanism. In Sect. 6 we analyse the implications of our results and in Sect. 7 we summarize our main conclusions.

2. Observations and data reduction

Observations of NGC 253 were taken on 5 December, 2010 with the Herschel Spectral and Photometric Imaging Receiver (SPIRE) in staring mode centered on the nucleus of NGC 253 (Obs ID: 1342210847). The SPIRE instrument is an imaging Fourier Transform Spectrometer (FTS; Griffin et al. 2010). The high spectral resolution mode was used with a resolution of 1.2 GHz over both observing bands. The low-frequency band covers ν = 447 − 989 GHz (λ = 671 − 303μm) and the high-frequency band covers ν-958-1545 GHz (λ = 313 − 194μm). A reference measurement was used to subtract the emission from the telescope and instrument.

The data were reduced using version 9.0 of the Herschel interactive processing environment (HIPE). Since NGC 253 is an extended source, a beam correction factor is necessary to compensate for the wavelength dependent beam size. Using an archival Submillimetre Common-User Bolometer Array (SCUBA) 450 μm map of NGC 253, we convolve the map with a 2D Gaussian with FWHM the same size as our beam sizes. The archival 450 μm SCUBA map of NGC 253 is shown in Fig. 1; the largest, smallest, and normalized beam sizes are shown in blue, green, and red, respectively. All beam sizes include the brightest part of the nucleus.

Fig. 1 Three SPIRE beam sizes are shown overplotted on a SCUBA 450 μm archival image, 40.5′′ (blue), 17.5′′ (green), and 32.5′′ (red). The SPIRE beam size changes as a function of wavelength, thus the blue circle represents the largest beam size and the green represents the smallest beam size. The red beam size (32.5′′) represents the beam size of the CO J = 5 − 4 transition and all other lines are convolved to this beam size.

Fig. 2 SPIRE spectra of NGC 253: short wavelength and long wavelength bands are separated by a dotted black line. 12CO transitions are marked in red, 13CO transitions in light blue, atomic transitions in green and H2O and HCO+ lines in dark blue.

The ratio of the flux within each convolved map and the flux within the beam size of the CO J = 5 − 4 transition (beam size = 32.5′′) is the beam correction factor (κ_S), where ${\mathit{F}}_{\mathrm{corr}}\mathrm{=}\frac{{\mathit{F}}_{\mathrm{obs}}}{{\mathit{\kappa }}_{\mathit{S}}}\mathrm{·}$(1)Thus, all fluxes are normalized to a beam size of 32.5′′ (i.e., 394 pc). In addition, we use the ground-based ¹²CO, ¹³CO, HCN, and HNC fluxes from Israel et al. (1995) and Israel (priv. comm.), all normalized to a beam size of 32.5′′. We add the 27′′ aperture integrated ground-based observations of HCO⁺ 1−0 and 4−3 transitions from Knudsen et al. (2007), also normalized to a beam size of 32.5′′. Fluxes were first extracted using FTFitter¹, a program specifically created to extract line fluxes from SPIRE FTS spectra. This is an interactive data language (IDL) based graphical user interface, that allows the user to fit lines, choose line profiles, fix any line parameter, and extract the flux. We define a polynomial baseline to fit the continuum and derive the flux from the baseline subtracted spectrum. In order to more accurately determine the amplitude of the line, we fix the FWHM for transitions higher than J = 8 − 7 to the expected line width of ¹²CO at each source, using the velocity widths measured by Israel et al. (1995).

In the case of very narrow linewidths (~780 km s^-1 at 650 μm), narrower than the instrumental resolution (FWHM_max = 1000 km s^-1 at 650 μm), we do not fix the FWHM, but fit the lines as an unresolved profile, which is the case for CO J = 4 − 3 through J = 8 − 7. We use an error of 30% for our fluxes, which encompasses our dominant source of error, specifically the uncertainty of the beam size correction using the 450 μm SCUBA map (15%) and line flux extraction (10%). We also have some uncertainty in the SPIRE calibration error of ~5% for extended sources.

3. Results

The spectra of NGC 253 are presented in Fig. 2. The ¹²CO transitions are visible from J = 4 − 3 to J = 13 − 12 and labeled in red and the two ¹³CO transitions detected are marked in light blue. There is also a strong detection of [NII] at 1461 GHz and [CI] at 492 GHz and 809 GHz in the rest frame, shown in green. We detect five strong water emission lines and H₂O 1₁₁ − 0₀₀ in absorption, marked in blue. We also detect HCO⁺ 7−6, also marked in blue. In addition, we also find CH⁺, CH, OH⁺, H₂O⁺, and HF, but do not point them out since they will not be used in our analysis. As seen in Fig. 2, there is a discontinuity between the high- and low-frequency bands of the spectrograph, marked with a dotted line that occurs because the beam size for the high-frequency band is much smaller than for the low-frequency band. A beam correction factor (κ_S) for each wavelength is calculated using the method described in Sect. 2, and displayed in Table 1. In addition, the baseline ripple seen in the spectrum, specifically in the inset of Fig. 2, is due to the sinc profile of the strong CO transitions and does not represent noise.

4. Dissecting the CO excitation ladder

In order to analyse the excitation conditions in NGC 253, we can create a spectral line energy distribution or what is referred to as a CO ladder, which plots the flux of each CO transition as a function of the upper J number. Since we have multiple transitions of both ¹²CO and ¹³CO, we can use both in our analysis. Since we take a beam size of 32.5′′ (~400 parsec diameter), we get emission from multiple phases of the ISM in one spectrum and cannot spatially separate the distinct ISM environments. Thus, when analysing the properties of each ISM phase, it is important to realize that these properties are average representative values of the dominant ISM environments that are responsible for the particular emission.

4.1. LVG analysis

To better understand the ISM excitation in the center of NGC 253, we start by modeling the ground-based fluxes, i.e., the first four ¹²CO and first three ¹³CO line intensities and ratios with the RADEX large velocity gradient (LVG) radiative transfer models (van der Tak et al. 2007). These codes provide model line intensities as a function of three input parameters per molecular gas phase: molecular gas kinetic temperature T_k, density n(H₂), and the CO velocity gradient N(CO)/dV. By comparing model to observed line ratios, we identify the physical parameters that best describe the actual conditions.

In the modeling, we assume a constant CO isotopical abundance ¹²CO/¹³CO = 40 throughout. This is close to values generally found in starburst galaxy centers (Henkel et al. 1991). We identify acceptable fits by exploring a large grid of model parameter combinations (T_k = 10 − 150 K, n(_H2) = 10² − 10⁵ cm^-2, and N(CO)/dV = 6 × 10¹⁵ − 3 × 10¹⁸ cm^-2/km s^-1) for (combined) line ratios matching those observed. No single-phase gas provides a good fit to the observed line intensities. Consequently, we have also modeled the ¹²CO and ¹³CO lines simultaneously with two molecular gas phases, the relative contribution of the two phases being a free parameter. Because we lack a detailed physical model for the distribution of clouds and their sources of excitation, two phases is the most we can fruitfully explore for J< 5, especially because only the high signal-to-noise lower J¹³CO transitions have observations that allow us to break the inherent T_k − n(H₂) degeneracy. This degeneracy is present because the intensity ratios of optically thick ¹²CO lines exhibit very little change going from hot, low-density gas to cool, high-density gas (Jansen 1995).

No unique solution is obtained; however, in all combinations the parameters of the first molecular gas phase are well-established and narrowly constrained to a low kinetic temperature of about 60 K (± 10) and a density of about log(n_H) = 3.5 (± 0.5). In all combinations the second phase has a reasonably well-determined higher density typically of log (n_H) = 5.0, but the kinetic temperature is not at all constrained by the lower J transitions included, nor is the proportion of molecular gas in either phase.

We now use this result, especially the firmly established parameters of the coldest and least dense gas from phase 1, to optimize a three-phase LVG model that includes the higher transitions as well, $\mathrm{Model}\mathrm{=}{\mathrm{\Omega }}_{\mathit{I}}\mathit{LV}{\mathit{G}}_{\mathit{I}}\mathrm{+}{\mathrm{\Omega }}_{\mathit{II}}\mathit{LV}{\mathit{G}}_{\mathit{II}}\mathrm{+}{\mathrm{\Omega }}_{\mathit{III}}\mathit{LV}{\mathit{G}}_{\mathit{III}}\mathit{,}$(2)where LVG_I, LVG_II, and LVG_III are three LVG models of specific density, temperature, and column density in units of W m^-2, and Ω_I, Ω_II, and Ω_III represent the respective filling factors of each ISM phase. Filling factors traditionally represent how much of the beam is filled, so they only range from 0 to 1. However, in the case of our filling factors, we model one cloud with a δv = 1 km s^-1 and allow for multiple clouds in our line of sight, such that Ω is not only a beam filling factor, but also a volume filling factor, which accounts for it being greater than one. In the case of Ω < 1, only a fraction of the beam is filled by clouds, whereas Ω > 1 represents the number of clouds in our beam volume, each with a slightly different velocity. We also attempted the same procedure for cloud models with δv = 5 and 10 km s^-1 without any significant difference. We create a composite model for both ¹²CO and ¹³CO where each model ISM phase and filling factor is the same for ¹²CO and ¹³CO. Since we already have an idea of the first LVG phase (LVG I), we only vary the second and third LVG phases and all three filling factors.

We perform a modified Pearson’s χ² minimized fit for ¹²CO and ¹³CO simultaneously, where the modified Pearson’s χ² is ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\frac{{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{data}}}\mathrm{(}\frac{{\mathrm{obs}}_{\mathit{i}}\mathrm{-}{\mathrm{model}}_{\mathit{i}}}{{\mathrm{model}}_{\mathit{i}}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{N}}_{\mathrm{data}}}$(3)with obs_i as the observed flux of a particular transition, model_i is the composite model flux (Eq. (2)) of a particular transition, and N_data is the total number of transitions. Thus, the χ² represents a χ² per transition for each molecule. We calculate the χ² value as the sum of the χ² for ¹²CO and ¹³CO for every possible combination of models in our grid. In doing this, we are able to plot the full parameter space for each phase, and show the χ² value in gray scale for each combination of density and temperature to see where the degeneracies lie. Figure 3 shows these degeneracy plots for the second and third LVG phases (LVG II and LVG III).

Fig. 3 Degeneracy plots of LVG II (right) and LVG III (left). Each square represents a model with specific density (x-axis) and temperature (y-axis) for the best fitting N(CO)/δv. The color scale represents the χ2 value of that model fit to the CO data. Since we use a modified Pearson’s χ2, the actual values of the gray scale are not meaningful, only their relative differences. The large black asterisk represents the best fitting model (minimum χ2).

Although there is a best-fit model, denoted by the black asterisk in Fig. 3, there are many models for each phase that have similarly low χ² values, especially for LVG II. In both LVG phases we see the general diagonal trend displaying the temperature-density degeneracy, trading temperature for density. For LVG II, there is a large range of models that fit well, however the best-fits tend to cluster at the low-temperature, high-density range (lower right corner). Instead, LVG III fills a much smaller part of parameter space. The diagonal degeneracy still exists, but this phase requires both high temperature and high density to be fit well. It is also important to note that LVG II and LVG III occupy completely independent diagonal regions of parameter space, and since the diagonal trend trades temperature for density, we can think of these separate regions as unique pressure phases. In Fig. 3, we only compare the degeneracies of temperature and density, yet we also analyse the degeneracies of N(CO)/δv and find that within a factor of ± 5 the column density is well constrained.

Although we cannot fully resolve the degeneracies, we will select a range of model parameters with a low χ² value to continue our analysis. We will define the density range of LVG II as 4 < log ₁₀(n_H/ cm^-3) < 6 and the range of LVG III as 4.5 < log ₁₀(n_H/ cm^-3) < 6. Figure 4 shows our three best-fit ISM LVG phases. Here we do not include the dense gas tracers to constrain the density of our LVG models since the abundances may vary between ISM phases; however, the densities of our best-fit LVG models are consistent with the HCN and HNC ratios shown in Fig. 7 and Table 4.

Fig. 4 12CO and 13CO excitation ladders of NGC 253 with the flux of each transition plotted as black asterisks with red error bars. In dotted blue, green, and yellow lines we plot the LVG I, II, and III ISM phases with their filling factors. The composite model is plotted with a solid red line. The model density, temperature, and column density are shown in the legend along with the relative contribution of each phase in terms of emission and column density.

The model parameters of our best three models (Fig. 4) are summarized in Table 2. The column density stated in this table is the model gradient, log(N(CO)/δv). We show a relative contribution of the emission by summing the flux over all transitions of each model phase and comparing that to the flux of the composite model (solid red line). We compare the total emission (summed flux of CO transitions) to the integrated emission of each ISM phase. We find a relative contribution in terms of emission (luminosity) to the total CO flux from LVG I, LVG II, and LVG III of 0.2:0.2:0.6, respectively. Next, we compare the relative scaled column densities which take the filling factors into account such that N_i = Ω_i × N_CO,i. We see a relative contribution in terms of column density of 0.6:0.4:0.02 for LVG I, II, and III respectively, which is also proportional to the relative mass contributions of each phase. Although LVG III is responsible for the bulk of the line emission, and all of the emission for the high-J CO lines, it represents a very small fraction (~2%) of the mass. It is also interesting to note that LVG III contributes almost nothing to the ¹³CO emission. This is due mainly to the low optical depth of LVG III, restricting the ¹²CO/¹³CO to much lower level than in LVG I and II. The fact that LVG III has such a high density, but such a low column density implies that the path length is very small or that the CO abundance is very low. The low CO column density combined with a high temperature suggests the physical environment of a PDR where the high UV flux keeps most of the CO dissociated; we explore this further in Sect. 4.2. We estimate a total gas mass of 3 × 10⁷M_⊙, which is only a lower limit for the mass since we assume a [C]/[H] abundance ratio of 1.4 × 10^-4 and that all of the [C] is in CO. However, in the case of a high-density, high-temperature LVG phase (LVG III), it is likely that the radiation environment is strong, and thus most of the CO would be dissociated, which would suggest a lower H₂ mass than there is in reality. Because of these uncertainties, we will perform a more detailed mass calculation in the following sections.

4.2. PDR analysis

Now that we have constrained the density and temperature ranges of the three dominant molecular ISM phases, we can determine which physical processes are exciting the molecular gas. Since the nucleus of NGC 253 experiences high star formation rates, there is an increase in number density of OB stars and thus a high UV energy density. Through photoelectric heating and FUV pumping of H₂, the FUV photons heat the outer layers (A_V< 5) of molecular clouds. This area of the molecular cloud is the PDR, and is responsible for warm molecular gas emission. The thermal state of PDRs is determined by processes such as photo-electric heating; heating by pumping of H₂ followed by collisional de-excitation; heating by cosmic rays; [OI] and [CII] fine-structure line cooling; and CO, H₂O, H₂, and OH molecular cooling. The ionization degree of the gas is driven by FUV photo-ionization, and counteracted by recombination and charge transfer reactions with metals and PAHs. The ionization degree is at most x_e ~ 10^-4 outside of the fully ionized zone. The chemistry exhibits two fundamental transition, H to H₂, and C⁺ to C to CO. Using PDR models from Kazandjian et al. (2012) that solve for chemistry and thermal balance throughout the layers of the PDR, we use the predictions of the ¹²CO and ¹³CO emission as a function of density, radiation environment (G, the Habing radiation field), and column density. The density profile is constant, as in the LVG models, and the Habing field is defined as 1 G = G₀ = 1.6 × 10^-3 erg cm^-2 s^-1 from photons between 6 eV and 13.6 eV. The fact that these models use G instead of temperature provides a direct diagnostic of the radiation source responsible for the gas excitation. The model is a semi-infinite slab model which is evaluated at A_V = 3,4,5,10,20, and 30 mags, and we calculate the column density from the A_V and density that gives us the best-fit emission.

We narrow the PDR model parameter space by allowing only the densities based on Fig. 3 (4 < log ₁₀(n_H/ cm^-3) < 6) for the second phase and (4.5 < log ₁₀(n_H/ cm^-3) < 6) for the third phase. In order to further constrain our fits, we include ¹²CO, ¹³CO, HCN, and HNC to our χ² calculation. Since HCN and HNC are typically dense gas tracers (both have n_crit> 10⁵ cm^-3), they provide an important additional constraint on the second and third phases, since according to our LVG analysis both phases are high density. We also have multiple HCO⁺ transition observations, but we do not use them in our analysis because of the complicated nature of the emission of HCO⁺. Its abundances are very sensitive to the interplay between heating and ionization rate and the resulting chemical state of the gas (Meijerink et al. 2011) and is thus an untrustworthy diagnostic.

The parameters of the best fitting models are listed in Table 3. In Fig. 5 the PDR fit is shown for ¹²CO, ¹³CO, HCN, and HNC. Each PDR phase is represented by a dotted blue, green, and yellow line, while the composite model is shown as a solid red line. We can achieve a relatively good fit for all phases, yet it is important to check if this fit makes physical sense.

Fig. 5 Left: 12CO and 13CO excitation ladders, right: HCN and HNC excitation ladders of NGC 253 with flux of each transition plotted as black asterisks with red error bars. In dotted blue, green and yellow lines we plot the PDR I, II and III ISM phases with their filling factors. The composite model is plotted with a solid red line. The model density, G0 and column density are shown in the legend along with the relative contribution of each phase in terms of emission and column density.

Using these three model PDR phases, we can estimate the molecular gas mass in each phase by ${\mathit{M}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{=}\sum _{\mathit{i}}^{\mathit{n}}\frac{{\mathrm{\Omega }}_{\mathit{i}}{\mathit{N}}_{{\mathrm{H}}_{\mathrm{2}}\mathit{,i}}{\mathit{A}}_{\mathrm{beam}}{\mathit{m}}_{{\mathrm{H}}_{\mathrm{2}}}}{{\mathit{M}}_{\mathrm{\odot }}}\mathit{,}$(4)

where N_H2 is the H₂ column density in cm^-2 which is consistently calculated in the PDR models, A_beam is the beam area in cm², and m_H2 is the mass of a hydrogen molecule. The mass of each PDR phase is shown in Table 3 and the total mass of the system is 1.9 × 10⁸M_⊙, which is almost equal to the gas mass measured for a beam size of 80′′, 3 × 10⁸M_⊙ (Weiß et al. 2008). Harrison et al. (1999) and Bradford et al. (2003) found a gas mass of ~3 × 10⁷M_⊙ for a beam size of 15′′, roughly half the size of our beam, yet we find a gas mass 6.5 times larger. In addition, our LVG models predict a mass of ~3 × 10⁷M_⊙, which is a reasonable (albeit a lower limit) gas mass estimate. The LVG mass is a lower limit since we can see that the CO becomes dissociated for a PDR with a high radiation field. Thus, using the standard abundance of CO to H₂ will underestimate the true H₂ column.

We test the physicality of our PDR models with greater precision by comparing the atomic gas to the molecular gas. Carral et al. (1994) uses [CII] 158 μm to trace the atomic gas mass in the inner 40′′ of NGC 253 and finds an atomic mass of 2.4 × 10⁶M_⊙. This atomic mass is 30 times smaller than the molecular gas that we observe in our 32.5′′ beam, which is similar to the findings of Bradford et al. (2003) and Hailey-Dunsheath et al. (2008). Tielens & Hollenbach (1985) predict that in Galactic PDRs, the first ~3A_Vs are irradiated and then cooled via atomic lines, whereas the next 3 A_V are cooled through molecular lines. From this they estimate that in a PDR, the masses of atomic and molecular mass should be about the same. Since our molecular mass is 30 times greater than our atomic mass, we can calculate the A_V where the atomic heating must transition to molecular heating to preserve this ratio and in which conditions this takes place.

First, we assume our PDR is a 2D finite slab with an ionized atomic outer layer and a molecular deeper layer; UV radiation penetrates the slab from one side and the slab has a constant density and a finite depth. In the case of our PDRs the depth is A_V = 5, as shown in Table 3. We can measure the mass in each layer, which is proportional to volume×density, and we assume a constant density and beam area (V = BA × r), where r is the depth into the cloud, so these both drop out. The volume of the first layer is proportional to R_x, the radius where the atomic to molecular transition occurs. The volume of the molecular layer is proportional to R_{AV = 5} − R_x, where R_{AV = 5} is the distance to the end of the slab. We can then say that the ratio of the molecular to atomic mass is proportional to (R_{AV = 5}/R_x) − 1 and we observe this ratio to be 30. Solving this equation for R_x, we find that in order to observe a molecular to atomic mass ratio of 30, the transition between atomic to molecular must occur at A_V = 0.2. Although this is a simplified calculation, it demonstrates the difficulty in producing 30 times more molecular mass than atomic. We use our PDR models to test in which environments this very low transition depth occurs. We find that the lowest transition depth in our model parameter space is for a PDR with the highest density (10⁶ cm^-3) and the lowest radiation field (10²G₀). The transition depth for our three PDR models is A_V = 3.7, 2.2, and 4.5 for PDR I, II, and III respectively. Since all of these are much above the necessary A_V = 0.2, we suggest that an additional heating source is necessary in order to account for the large amount of bright molecular gas.

In order to determine the cause of this discrepancy, we start with our most constrained ISM phase, LVG/PDR I. We have determined that this ISM phase has a density of 10^3.5 cm^-3 and a kinetic temperature of 60 K. We fit this with a PDR with the same density, the same column density, and a G = 10^2.5G₀. From G₀ we can calculate the average temperature of the PDR, which is the dominant factor in determining the emission intensity of the lines. We find that in PDR , we only have a temperature of 13.5 K, while we need a temperature of 60 K. If we increase the radiation field in order to heat the gas further, we get up to a peak temperature of 18 K and then begin to dissociate the CO. The only way in which we can reproduce a PDR with n = 10^3.5 cm^-3 and T = 60 K is by including an additional heating term.

5. Dominant molecular excitation mechanisms

Since our analysis suggests that a PDR cannot be the sole source of the high-J CO transitions, we investigate excitation from other sources for comparison.

5.1. Excitation mechanisms

In addition to PDR heating, the three dominant heating mechanism of molecular gas are X-rays, cosmic rays, and mechanical heating. X-rays heat gas in regions called X-ray dominated regions (XDRs) that are similar to PDRs except that the chemistry is driven by X-ray photons instead of FUV photons; the X-ray photons are able to penetrate farther into the cloud without efficiently heating the dust at the same time. These X-rays are mostly produced by active galactic nuclei (AGN) or in areas of extreme massive star formation. The observed X-ray flux between 0.1 keV and 2.4 keV is 5.72 ± 0.59 × 10^-12 erg s^-1 cm^-2, which translates to a luminosity of L_X = 4 × 10³⁹ ergs s^-1 at a distance of 2.5 Mpc (Pietsch et al. 2000). Since NGC 253 lacks an X-ray bright AGN and has a relatively low X-ray luminosity compared to regions where we see X-ray heating, we exclude this mechanism from farther analysis. Cosmic rays can also heat gas in cosmic ray dominated regions (CDRs). Cosmic rays are able to penetrate into the very centers of molecular clouds, where even X-rays have trouble reaching and are typically produced by supernovae. Similarly, PDRs with additional mechanical heating (mPDRs) are due to turbulence in the ISM and may be driven by supernovae, strong stellar winds, or jets.

In order to test for mechanical heating, we can add a mechanical heating term to the PDR models (mPDR models). Since turbulence can penetrate a cloud at all scales and depths, as a simple approximation the mechanical heating (Γ_mech) is introduced uniformly throughout the 1D PDR models. We parameterize the strength of Γ_mech with α, which represents the fractional contribution of mechanical heating in comparison to the total heating at the surface of a pure PDR (excluding mechanical heating). At the surface the heating budget is dominated by photoelectric heating. When α = 0, there is only photoelectric heating and cosmic ray heating and when α = 1, the mechanical heating rate is equal to the photoelectric heating at the edge of the cloud. These rates can be related to supernovae rates, which in turn can be related to star formation rates (as was done in Loenen et al. (2008)). The freedom in the choice of the values of α allows for a flexibility in modeling the source of the absorbed mechanical heating. The surface heating rates for the PDR models can be recovered from Figure 1 in Kazandjian et al. (2012). More details about relating the mechanical heating to the star formation rate can be found in the methods section of that paper and in Loenen et al. (2008). In addition to the mechanical heating, there is a heating rate due to cosmic rays applied throughout the PDR. The heating by cosmic rays is only dependent on H₂ density and ionization rate, in this case the galactic cosmic ray rate (0.2−3.5 × 10^-16 s^-1) (Goldsmith & Langer 1978; van Dishoeck & Black 1986; Indriolo & McCall 2012) and the H₂ density of the selected model.

Since both cosmic rays and mechanical heating derive from supernovae, it is likely that if we have excitation from one mechanism, we should expect excitation from both. Thus, in addition to the mPDR models, we have run mCDR models. These models not only have photoelectric heating from the PDR but they also have mechanical heating and an increased cosmic ray ionization rate of 750 ζ_gal or 3.75 × 10^-14 s^-1, as suggested by Bradford et al. (2003). The enhanced cosmic ray ionization rate is paired with a varying mechanical heating rate from α = 0, or no mechanical heating, up to α = 1.0 which translates to a heating rate of Γ_mech = 1.3 × 10^-15 erg s^-1 cm^-3. This allows for pure cosmic ray excitation, where α = 0, and for the case where both mechanical heating and cosmic rays play a role in the gas heating.

5.2. Application to NGC 253

In order to determine the alternate heating mechanism, we first start with our most constrained ISM phase, PDR I/LVG I. With only photoelectric heating, we cannot reach the temperature measured in the LVG model (60 K), and thus must include an additional heating source. In order to compare various heating sources, we plot the LVG model (Fig. 6, left panel) in black and overplot other similar models. Like the pure PDR models, pure CDR with no mechanical heating dissociates the CO before the temperature can get high enough to produce the observed emission. First, we overplot the PDR model that we used in the PDR fit (Sect. 4.2) in green. The emission from the pure PDR model is approximately a factor of 10 lower than the LVG model of the same density and T = 60 K. This explains the high filling factors necessary to fit the model emission to our observations, and thus the high mass estimates. We can then use the same parameters as the pure PDR fit, but add additional heating to see if we can reproduce the LVG I emission profile. First we try to fit the emission with an CDR, identical to the PDR except with an enhanced cosmic ray rate of 750 times the galactic CR rate. This curve is overplotted in orange. The ¹²CO is ten times less luminous, the ¹²CO/¹³CO ratio is not preserved and the ¹³CO is underproduced in comparison to the LVG I model, which makes this an unlikely model. We can then test a mechanically heated PDR (mPDR) with an α of only 1%. This matches the emission of LVG I for both ¹²CO and ¹³CO and requires only a small amount of mechanical heating. Therefore, we choose mechanical heating as the necessary additional heating mechanism for the first ISM phase.

Fig. 6 Comparison of PDR (green), mPDR (blue), and mCDR (orange) models to LVG I (black) for both 12CO and 13CO.

Next, we must test this first component to see if it is well constrained. In the right panel of Fig. 6, we have varied radiation strength, column density, and amount of mechanical heating to show that this component is not well constrained. We can match the emission of the LVG model with any model with mechanical heating between 1%−25% and radiation in the range 10^1.0−2.5G₀. However, it is important to note that while these parameters change, so does the column density in a way that preserves mass. Therefore, we cannot determine which of these models represents the true physical conditions, but we can say that mechanical heating is necessary and the mass of the component is around 5 × 10⁷M_⊙. In order to proceed with the analysis, we select the model with 1% mechanical heating (blue in Fig. 6) since it is a slightly better fit to the ¹³CO.

Now that we have determined the range of models applicable for the first ISM phase, we fit our second and third ISM phases with the mPDR and mCDR models to determine their dominant heating source. We hold the first ISM phase constant and equal to that of the 1% mechanical heating mPDR in Fig. 6, but we allow the filling factor to vary. The second and third phases are held to the same requirements as in the PDR fitting, namely the second phase density is between 10⁴ and 10⁶ cm^-3 and the third phase density is between 10^4.5 and 10⁶ cm^-3. The χ² for the ¹²CO, ¹³CO, HCN, and HNC ladders are summed to find an all around best fit. The best fitting models are shown in Fig. 7 for

• 1.

  one mPDR and two mCDRs;

• 2.

  two mPDRs and one mCDR;

• 3.

  three mPDRs.

The parameters of the fit for each case are shown in Table 4. The mCDRs can have mechanical heating or can have α = 0, representing a pure CDR.

Fig. 7 12CO and 13CO (top) with HCN and HNC (bottom) excitation ladders of NGC 253 with flux of each transition plotted as black asterisks with red error bars. In dotted blue, green and yellow lines we plot the PDR/mCDR/mCDR (left), PDR/mPDR/mCDR (center), and PDR/mPDR/mPDR(right) ISM phases with their filling factors. The composite model is plotted with a solid red line. The model H2 density [log cm-3], G0 [log Habing flux], and percentage mechanical heating (α) are shown in the legend along with the relative contribution of each phase in terms of emission and column density.

Bradford et al. (2003) estimate an ionization rate of (1.5−5.3) × 10^-14 s^-1, which is 750 times higher than the galactic ionization rate. Acero et al. (2009) measured a cosmic ray ionization rate up to 1000 times our galactic ionization rate, which is consistent with the results from Bradford et al. (2003). We consider both of these cosmic ray ionization rates to be upper limits, and suggest that the true cosmic ray ionization rate is somewhere between galactic and 750 times the galactic value. The difference between a cosmic ray ionization rate of 1000 and 750 times the galactic value is negligible, so we will use the values 750 throughout the rest of the paper. Martín et al. (2006) finds evidence that large-scale, low-velocity shocks drive the chemistry and heating of NGC 253. In support of this, diffuse SiO emission has been detected throughout the nuclear region (García-Burillo et al. 2000). In addition, Martín et al. (2005) have measured an overabundance of OCS which is caused by low velocity shocks injecting these molecules into the gas phase. Mechanical heating can also be introduced into the ISM through supernova remnants (Loenen et al. 2008) or strong jets and we test its efficiency at heating molecular gas using the mPDR and mCDR models. Strong molecular outflows have been observed by Bolatto et al. (2013), and they estimate an outflow rate of 9 M_⊙ yr^-1, implying the outflows are a strong source of turbulence in the galactic nucleus.

Since Bradford et al. (2003) suggest that cosmic rays are responsible for most of the molecular gas excitation in NGC 253, we first use case 1, two mCDR models, to fit our observations and see if they are consistent with cosmic ray excitation. The mCDR models allow for mechanical heating, but can also have α = 0, representing a pure CDR. We define cosmic ray excitation as ionization by cosmic rays, which produces slow electrons that excite atomic and molecular hydrogen, and in turn produces photons that heat the gas (Glassgold et al. 2012). As seen in the left panel of Fig. 7, the best-fitting model requires no mechanical heating for the second ISM phase, thus it is a pure CDR and requires only 1% mechanical heating for the third ISM phase. The fit of case 1 to CO, HCN, and HNC is very good, yet this suggests that the first component could also have an enhanced CR rate. We tested a fit with all three ISM phases including an enhanced cosmic ray rate and similarly to Fig. 6 the ¹³CO is not fit at all. If we tweak the χ² calculation to make CO more important than HCN and HNC, we can fit ¹²CO and ¹³CO but completely destroy the fits of HCN and HNC. However, the fact that we cannot fit three components with a high CR rate does not mean that the fit of case 1 is not possible. Instead, it could represent the physical situation where the diffuse ISM, traced by the first ISM phase, has a diluted cosmic ray rate that can be assumed to be equivalent to the galactic rate.

Since cosmic rays are predominately from supernovae, they may only affect a small percentage of the total gas mass. Therefore, we attempt to recreate our observed line fluxes in case 2 with one PDR, one mPDR, and one mCDR. This scenario represents the case where mechanical heating is present throughout the ISM, and in addition, cosmic rays excite a small amount of gas around supernova remnants. The best-fit models for this scenario are shown in the center panel of Fig. 7. The models fit the ¹²CO, ¹³CO, HCN, and HNC reasonably well, if not within the error bars of all transitions. Similarly to case 1, the HNC emission is off by a factor of ~2, yet since our models are not very advanced in their treatment of the chemistry, a factor of 2 is not enough to rule out this case.

Finally, we test the case where mechanical heating is the dominant heating mechanism and cosmic rays are not important to the heating of the molecular gas. The right panel of Fig. 7 shows the best-fit for case 3, three mPDR models; mPDR II has an α = 0.1 which translates into a Γ_mech = 1 × 10^-20 erg s^-1 cm^-3 and mPDR III has an α = 0.1 which is Γ_mech = 1.6 × 10^-18 erg s^-1 cm^-3. This model fits all the ¹²CO and ¹³CO transitions within the error bars. In addition, HCN and HNC are better fit than in both case 1 and case 2. Since all three cases reproduce very similar, and very good fits of the molecular transitions, we cannot determine the precise mechanisms. Even in case 1 mechanical heating is needed in the first ISM phase and in the third ISM phase needing an α = 0.01, which is a non-negligible heat rate (1 × 10^-18 erg s^-1 cm^-3) as shown in Fig. 8. Thus, it is clear that mechanical heating is necessary to reproduce the observed emission, regardless of the added heating of cosmic rays. This finding is similar to that of Nikola et al. (2011) in NGC 891 and Hailey-Dunsheath et al. (2008) for NGC 253 itself, where they use a combination of shock and PDR models to determine that microturbulence and shocks are responsible for heating the mid-J CO gas.

6. Discussion

Since we have concluded that in NGC 253 mechanical heating plays an important role in the gas excitation, we can compare the total mechanical heating with the other heating sources to see which type of heating is dominant. In order to understand the evolution of heating throughout the cloud, we plot the heating rates as a function of A_V, or depth into the cloud. We do this for each ISM phase separately as well as for the total heating of the composite system. We compare all three of our model scenarios: case 1 with two (m)CDRs, case 2 with one mPDR and one mCDR, and case 3 with two mPDRs. The photoelectric heating includes all related heating mechanisms such as H₂ pumping, CII ionization, viscous heating, etc. The heating breakdown is shown in Fig. 8.

Fig. 8 Heating rate compared with AV, or depth into cloud. The total heating is shown in black, the cosmic ray heating in purple, the mechanical heating in blue, and the photoelectric heating in orange. The top row is for case 1 with one mPDR, one CDR, and one mCDR. The middle row is case 2 with two mPDRs and one mCDR and the bottom row is case 3 with three mPDRs. The first column is for the first component in each case, the mPDR, and it is the same in all three cases. The second and third columns represent the second and third ISM phases for each case. The vertical dotted line in these plots represents the edge of the PDR, as some of the models we use have AV< 5. Finally, the last column is the integrated heating rates for all ISM phases.

6.1. Case 1

For the first ISM phase (the mPDR) the total heating (in black) is produced mostly by photoelectric heating from the PDR since there is only 1% mechanical heating. However, at A_V = 3.5 the mechanical heating overtakes the photoelectric heating and becomes the dominant heating source. In all cases, the first PDR phase is kept constant, thus for the next sections we will not discuss it further.

For each ISM phase and for the integrated heating, at A_V = 0, or the edge of the cloud, photoelectric heating is the dominant heating source, as expected in a PDR. However, in the second ISM phase (the pure CDR) the cosmic ray heating quickly becomes equally, and slightly more important than photoelectric heating at A_V = 1. The third phase, the mCDR has an α = 0.01, and despite a very low mechanical heating rate, the mechanical heating rate is at least ten times higher than the cosmic ray heating and the photoelectric heating rate dominates the heating through A_V = 5. Looking at the integrated heating, we see that this case (and all cases) is dominated by photoelectric heating. Even though this case includes two ISM phases with an enhanced cosmic ray ionization rate, the CR heating never exceeds that of the mechanical heating.

6.2. Case 2

In contrast, the integrated heating for case 2, with two mPDRs and one mCDR, shows that although photoelectric heating is the dominant mechanism, mechanical heating still contributes a significant fraction of the heating and CR heating is insignificant. We see that for the second ISM phase, the emission is dominated by photoelectric heating up to A_V = 2 and then mechanical heating becomes equally important. Here cosmic rays are relatively unimportant. The mCDR model has an enhanced cosmic ray ionization rate, and is reflected in the heightened CR heating rate. However, in order to reproduced the observed HCN and HNC fluxes, mechanical heating is still necessary and is consistently more important than the cosmic ray heating rate by a factor of at least 50. Thus, in case 2, although we included an enhanced cosmic ray ionization rate in one of the ISM phases, it is unimportant to the total heating of the gas.

6.3. Case 3

Finally we study case 3, where there are three mPDRs. In this case, both the second and the third phase have an α = 0.1 and neither has an enhanced cosmic ray ionization rate. In both mPDR phases, the photoelectric heating drops off sharply, and is replaced by mechanical heating at A_V = 1 and A_V = 3 for phase 2 and 3, respectively. Studying the integrated heating, photoelectric heating still dominates, but the mechanical heating is even more important here, and the CR heating is about 5 orders of magnitude less important.

From these plots, we can conclude that in a PDR with any amount of mechanical heating, mechanical heating is more important the cosmic ray rate, regardless of whether it is enhanced. However, mechanical heating being the dominant heating source at the edge of a particular cloud does not mean that it is the overall most efficient source of heating (see the integrated heating plots). For example, in case 3 with three mPDRs, mechanical heating becomes the dominant heating source in each individual PDR by at least A_V = 3.5, yet if we integrate the heating of the composite model we find that mechanical heating is only responsible for 20% of the total heating, while photoelectric heating heats 80%. Similarly, if we integrate the contributions from each heating source in the composite model for case 2, we find 84% photoelectric heating and 16% mechanical heating. In both cases 1 and 2, cosmic rays play a negligible role in heating the gas, which may be expected since only the third phase in case 2 has an enhanced cosmic ray ionization rate. However, when we integrate the heating phases of case 1, with two phases with an increased ionization rate, we find 98% of the heating from photoelectric heating and the remaining 2% from mechanical heating. The cosmic ray heating contributes only 0.17% of the total heating. Therefore, we suggest that although we cannot rule out cosmic ray contributions to heating the gas, the dominant heating mechanisms are photoelectric heating and mechanical heating. In addition, although mechanical heating may dominate the heating after A_V = 3 in individual clouds, it still only contributes a maximum of 20% of the heating in the cloud, in agreement with the upper limit prediction from Rosenberg et al. (2013). This suggests that the source of feedback exciting the high-density gas is not directly from star formation, but instead from the molecular outflows which provide mechanical feedback without the production of cosmic rays Bolatto et al. (2013).

We estimate the mass of each ISM phase and find that regardless of the excitation mechanisms, the mass is relatively well constrained for each case. The first phase, an mPDR with 1% mechanical heating, has a mass of 6 × 10⁷M_⊙. The second phase (either CDR or mPDR) has a mass of 2 − 5 × 10⁷M_⊙ and the third phase (either mPDR or mCDR) has a mass of 2 − 4 × 10⁶M_⊙. With the exception of the third phase, the masses of each ISM phase are all approximately the same. In addition, the total mass of each system is well conserved. We find a mass of 8 − 11 × 10⁷M_⊙. This mass is in agreement with 3 × 10⁸M_⊙ in an 80′′ beam (Weiß et al. 2008) since if the gas were evenly distributed, we would expect a mass of 5 × 10⁷M_⊙. The gas mass is concentrated in the center of the galaxy, thus it is reasonable to expect a mass a few times higher than the evenly distributed mass. In all cases the third phase contributes at least 60% of the CO emission but represents a maximum of 15% of the mass. In addition, in the third ISM phase of all cases, the CO column density is less than 1% of the total CO column. This suggests that the CO is mainly dissociated in this phase, but still responsible for ~ 60% of the C.

The source NGC 253 is an ideal case in which we have any observation desired. It is a nearby, well-studied, star-forming galaxy that gives us a representative view of environments undergoing rapid star formation. We have identified that in these galaxies, although mechanical heating is important, heating by UV photons is still the dominant heating source. We have also determined that in this environment, cosmic rays at best play a minor role in heating the molecular gas. In this sense, any star forming galaxy within a similar luminosity class to NGC 253, would have similar contributions of heating mechanisms. Yet in (ultra) luminous infrared galaxies ((U)LIRGs), there are many other factors in play, including AGN which produce X-rays, mergers, and powerful outflows, all of which increase turbulent heating. It is conceivable that mechanical heating is the domninant source of excitation is these galaxies, or that X-rays are much more important.

For example, studies of Mrk 231 show that the ISM is being processed both by star-formation and supernovae (González-Alfonso et al. 2014; Fischer et al. 2010), while NGC 6240 shows clear evidence of shocks and mechanical heating (Meijerink et al. 2013). In NGC 1068, Hailey-Dunsheath et al. (2012) find that the gas is either excited by X-rays or shocks, but cannot differentiate between the two. These galaxies are some examples of (U)LIRGs that show the complexity of disentangling heating mechanisms in these luminous systems. They also demonstrate that in the case of (U)LIRGs it is rarely only UV photons that contribute to heating the molecular gas. Although NGC 253 is not a LIRG, the starburst in the center represents a similar environment, and already demonstrates the need to include additional heating mechanisms.

7. Conclusions

We observe the starburst nucleus of NGC 253 with Herschel SPIRE in order to understand the heating and excitation of the gas.

• 1.

  We extract the fluxes of all the lines in our spectra using FTFitter,and correct them for the changing beam size of SPIRE. We do abasic LVG analysis, using both ¹²CO and ¹³CO in order to constrainthe parameter space. By mapping the χ² values for the full parameterspace, we can visualize the degeneracies and limit the parameterspace further.

• 2.

  Using the limits set by the LVG analysis, we perform a PDR analysis. We use not only the ¹² and ¹³CO, but also include the HCN and HNC ladders, to constrain the dense gas. Although we get a good fit with three PDR models, we can see that the the solutions are not reasonable and an additional heating source is required. The additional heating source is needed because we cannot reproduce a high enough temperature with just a PDR model to match the observed CO emission.

• 3.

  Using PDR models with a constant mechanical heating rate (mPDRs) and PDR models with both a constant mechanical heating rate and a cosmic ray ionization rate equal to 750 times our galactic value (mCDRs), we can fit ¹²CO, ¹³CO, HCN, and HNC. We fit the four molecular ladders for three different cases, case 1 is one mPDR with two mCDRs, case 2 is two mPDRs and one mCDR, and case 3 is three mPDRs.

• 4.

  We find that all cases provide a good fit. We cannot fit all four molecules with only enhanced cosmic ray models, we require the first model to be mechanically heated. Although we require mechanical heating to fit the observed transitions, we cannot rule out the existence of cosmic ray heating in addition to the mechanical heating.

• 5.

  Taking these three different cases of fits, we study the heating balance as a function of depth into the cloud. We find that in case 1, with an α of just 1% mechanical heating dominates cosmic ray heating at all depths of the cloud. In both cases 2 and 3, mechanical heating quickly becomes the dominant heating mechanism as the cloud depth increases for individual models

• 6.

  Regardless of mechanical heating dominating the heating sources, in our most heavily mechanically heated case (case 3), it only contributes about 20% of the total heating, while photoelectric heating is still the dominant overall heating source. In all cases, including case 1 with an enhanced ionization rate for two of the ISM phases, cosmic ray heating is negligible to the total heating of the cloud, as seen in the integrated heating rates of Fig. 8.

• 7.

  We estimate the mass of each ISM phase and regardless of the case, the mass is well constrained. The first phase, excited by predominantly photoelectric heating, has a mass of 6 × 10⁷M_⊙. The second phase has a mass of 3 × 10⁷M_⊙ and the third phase has a mass of 9 × 10⁶M_⊙. Although the phases have similar masses, the third phase is responsible for ~ 60% of the CO emission.

Acknowledgments

We would like to thank Edward Polehampton for his help preparing the SPIRE observations. SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including: Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). The Herschel spacecraft was designed, built, tested, and launched under a contract to ESA managed by the Herschel/Planck Project team by an industrial consortium under the overall responsibility of the prime contractor Thales Alenia Space (Cannes), and including Astrium (Friedrichshafen) responsible for the payload module and for system testing at spacecraft level, Thales Alenia Space (Turin) responsible for the service module, and Astrium (Toulouse) responsible for the telescope, with in excess of a hundred subcontractors. HIPE is a joint development by the Herschel Science Ground Segment Consortium, consisting of ESA, the NASA Herschel Science Center, and the HIFI, PACS and SPIRE consortia.

References

All Tables

All Figures

Fig. 1 Three SPIRE beam sizes are shown overplotted on a SCUBA 450 μm archival image, 40.5′′ (blue), 17.5′′ (green), and 32.5′′ (red). The SPIRE beam size changes as a function of wavelength, thus the blue circle represents the largest beam size and the green represents the smallest beam size. The red beam size (32.5′′) represents the beam size of the CO J = 5 − 4 transition and all other lines are convolved to this beam size.

In the text

	Fig. 2 SPIRE spectra of NGC 253: short wavelength and long wavelength bands are separated by a dotted black line. 12CO transitions are marked in red, 13CO transitions in light blue, atomic transitions in green and H2O and HCO+ lines in dark blue.
In the text

Fig. 3 Degeneracy plots of LVG II (right) and LVG III (left). Each square represents a model with specific density (x-axis) and temperature (y-axis) for the best fitting N(CO)/δv. The color scale represents the χ2 value of that model fit to the CO data. Since we use a modified Pearson’s χ2, the actual values of the gray scale are not meaningful, only their relative differences. The large black asterisk represents the best fitting model (minimum χ2).

In the text

Fig. 4 12CO and 13CO excitation ladders of NGC 253 with the flux of each transition plotted as black asterisks with red error bars. In dotted blue, green, and yellow lines we plot the LVG I, II, and III ISM phases with their filling factors. The composite model is plotted with a solid red line. The model density, temperature, and column density are shown in the legend along with the relative contribution of each phase in terms of emission and column density.

In the text

Fig. 5 Left: 12CO and 13CO excitation ladders, right: HCN and HNC excitation ladders of NGC 253 with flux of each transition plotted as black asterisks with red error bars. In dotted blue, green and yellow lines we plot the PDR I, II and III ISM phases with their filling factors. The composite model is plotted with a solid red line. The model density, G0 and column density are shown in the legend along with the relative contribution of each phase in terms of emission and column density.

In the text

	Fig. 6 Comparison of PDR (green), mPDR (blue), and mCDR (orange) models to LVG I (black) for both 12CO and 13CO.
In the text

Fig. 7 12CO and 13CO (top) with HCN and HNC (bottom) excitation ladders of NGC 253 with flux of each transition plotted as black asterisks with red error bars. In dotted blue, green and yellow lines we plot the PDR/mCDR/mCDR (left), PDR/mPDR/mCDR (center), and PDR/mPDR/mPDR(right) ISM phases with their filling factors. The composite model is plotted with a solid red line. The model H2 density [log cm-3], G0 [log Habing flux], and percentage mechanical heating (α) are shown in the legend along with the relative contribution of each phase in terms of emission and column density.

In the text

Fig. 8 Heating rate compared with AV, or depth into cloud. The total heating is shown in black, the cosmic ray heating in purple, the mechanical heating in blue, and the photoelectric heating in orange. The top row is for case 1 with one mPDR, one CDR, and one mCDR. The middle row is case 2 with two mPDRs and one mCDR and the bottom row is case 3 with three mPDRs. The first column is for the first component in each case, the mPDR, and it is the same in all three cases. The second and third columns represent the second and third ISM phases for each case. The vertical dotted line in these plots represents the edge of the PDR, as some of the models we use have AV< 5. Finally, the last column is the integrated heating rates for all ISM phases.

In the text