A Thorough View of the Nuclear Region of NGC 253: Combined Herschel, SOFIA, and APEX Data Set

Several results have been published in the last few years, reporting observations with the SPIRE Fourier Transform Spectrometer (FTS; Griffin et al. 2010) on board the Herschel Space Observatory⁷   toward various galaxies, e.g., M82, NGC 1068, Mrk 231, Arp 220, NGC 6240, and NGC 253 (Panuzzo et al. 2010; van der Werf et al. 2010; Rangwala et al. 2011; Kamenetzky et al. 2012; Spinoglio et al. 2012; Meijerink et al. 2013; Rosenberg et al. 2014), and toward the Galactic Center (Goicoechea et al. 2013). The availability of a comprehensive number of transitions in the ¹²CO ladder (up to the J = 13 $\to$ 12 transition) provided by the SPIRE-FTS has made possible the analysis of the ¹²CO line spectral energy distribution (LSED) using a variety of models of photon-dominated regions (PDRs), hard X-ray-dominated regions (XDRs), cosmic-ray (CR)-dominated regions, and shocks. These models were used to estimate the source of excitation and the ambient conditions of the molecular gas in all these galaxies while constraining their parameters using not only the ¹²CO spectral lines, but also the ¹³CO lines (e.g., Kamenetzky et al. 2012) as well as complementary ground-based observations of other molecules, like HCN (e.g., Rosenberg et al. 2014).

From these galaxies, so far only NGC 1068 has been studied by combining the ¹²CO (and other molecule) lines from SPIRE-FTS and PACS spectrometry data (Hailey-Dunsheath et al. 2012; Spinoglio et al. 2012). The PACS ¹²CO lines were shown to trace very different ambient conditions driven by X-rays in the nuclear region of NGC 1068 (Spinoglio et al. 2012). Since the PACS data for NGC 253 are also available, we can now revisit and extend the analysis and interpretation of the ¹²CO LSED done by Rosenberg et al. (2014) based on SPIRE-FTS data only.

The Sculptor galaxy NGC 253 is a nearby (D ≈ 3.5 Mpc; e.g., Mouhcine et al. 2005; Rekola et al. 2005), nearly edge-on (i ∼ 72°–78°; Pence 1981; Puche et al. 1991), isolated spiral galaxy of type SAB(s)c and is one of the closest galaxies outside the Local Group. Its angular size in the visible range is 275 × 68. Together with M82, it is the best example of a nuclear starburst (Rieke et al. 1988). Although an active galactic nucleus (AGN) has also been suggested to coexist with the nuclear starburst (e.g., Weaver et al. 2002; Müller-Sánchez et al. 2010), the corresponding low-luminosity AGN is not energetically dominant (Forbes et al. 2000; Weaver et al. 2002), and its IR/radio luminosity ratio indicates that the nature of its AGN candidate is more similar to that of the low-accretion-rate supermassive black hole Sgr A* at the center of the Milky Way galaxy than to that of an actual AGN driven by a more luminous central supermassive black hole (Fernández-Ontiveros et al. 2009).

The bulk of the NGC 253 starburst is confined to a ∼60 pc region centered southwest of the dynamical nucleus, according to the distribution of the 10–30 μm continuum (Telesco et al. 1993). The estimated star formation rate in the nuclear region is ∼2–3 M_☉ yr⁻¹ (Radovich et al. 2001; Ott et al. 2005). The 1–300 μm luminosity of NGC 253 is, within ∼30'', 1.6 ×10¹⁰ L_☉ (Telesco & Harper 1980). The total IR luminosity detected by IRAS is L_IR ≈ 2 × 10¹⁰ L_☉ (Rice et al. 1988). Four luminous super star clusters were discovered with the Hubble Space Telescope, and a bolometric luminosity of 1.3 × 10⁹ L_☉ was estimated for the brightest cluster (Watson et al. 1996). A bar can be seen in the Two Micron All Sky Survey image of NGC 253, and the kinematics show evidence for orbits in a bar potential (Scoville et al. 1985; Das et al. 2001, and references therein). Thus, the starburst activity is thought to be supported by the material brought to the nucleus by this bar (e.g., Engelbracht et al. 1998; Jarrett et al. 2003, and references therein).

Current estimates of the neutral gas mass in the central 20''–50'' range from 2.5 × 10⁷ M_☉ (Harrison et al. 1999) to 4.8 × 10⁸ M_☉ (Houghton et al. 1997). Studies of near-infrared and mid-infrared lines showed that the properties of the initial mass function (IMF) in the starburst region are similar to those of a Miller–Scalo IMF that has a deficiency in low-mass stars (Engelbracht et al. 1998).

A supernova rate of ≤0.3 yr⁻¹ has been inferred for the entire galaxy from radio (Antonucci & Ulvestad 1988; Ulvestad & Antonucci 1997) and infrared (van Buren & Greenhouse 1994) observations. The rate is most pronounced in the central starburst region, where a conservative estimate yields a supernova rate of ∼0.03 yr⁻¹, which is comparable to that in our Galaxy (Engelbracht et al. 1998). This suggests a very high local CR energy density. The mean density of the interstellar gas in the central starburst region is n ≈ 600 protons cm⁻³ (Sorai et al. 2000), which is about three orders of magnitude higher than the average density of the gas in the Milky Way. This extraordinary combination of high-density gas and enhanced local CR energy density was predicted to produce gamma-rays at a detectable level (Paglione et al. 1996; Domingo-Santamaría & Torres 2005; Rephaeli et al. 2010). In fact, very high-energy (VHE; >100 GeV) gamma-rays were effectively detected later in the nuclear region of NGC 253 with the High Energy Stereoscopic System array of imaging atmospheric Cherenkov telescopes (Acero et al. 2009) and with the Large Area Telescope on board the Fermi Gamma-ray Space Telescope (Abdo et al. 2010). The integral gamma-ray flux of the source above 220 GeV is ∼5.5 × 10⁻¹³ cm⁻² s⁻¹, which corresponds to ∼0.3% of the VHE gamma-ray flux observed in the Crab Nebula (Aharonian et al. 2006b) and is consistent with the original prediction (Paglione et al. 1996). The detection of VHE gamma-rays in NGC 253 implies a high energy density of CRs in this galaxy. Based on the observed gamma-ray flux, the CR density was estimated to be 4.9 × 10⁻¹² cm⁻³, with a corresponding energy density of CRs of ∼6.4 eV cm⁻³ (Acero et al. 2009). Thus, the CR density in NGC 253 is about 1400 times the value at the center of the Milky Way (Aharonian et al. 2006a).

The observations of Hα emission together with the earlier Einstein and ROSAT X-ray data (Ptak et al. 1997 and references therein) revealed a starburst-driven wind emanating from the nucleus along the minor axis of the galaxy. This wind was also detected later on by Chandra (Strickland et al. 2000). Extraplanar outflowing molecular gas was also mapped in the ¹²CO J = 1 $\to$ 0 line with the high spatial resolution provided by the Atacama Large Millimeter/submillimeter Array (ALMA), and it was found to follow closely the Hα filaments (Bolatto et al. 2013). The estimated molecular outflow rate is 3–9 M_☉ yr⁻¹, implying a ratio of mass-outflow rate to star formation rate of about 1–3. This ratio is indicative of the suppression of the star formation activity in NGC 253 by the starburst-driven wind (Bolatto et al. 2013).

NGC 253 is also the brightest extragalactic source in the submillimeter range, so its nucleus has been observed in various lines of CO and C (Harris et al. 1991; Wall et al. 1991; Guesten et al. 1993; Israel et al. 1995; Harrison et al. 1999; Sorai et al. 2000; Israel & Baas 2002; Bradford et al. 2003; Bayet et al. 2004; Güsten & Philipp 2006b), and it was the selected target for the first unbiased molecular line survey of an extragalactic source (Martín et al. 2006). This survey showed that NGC 253 has a very rich molecular chemistry, with strong similarities to that of the Galactic central molecular zone, and even molecules like H₃O⁺ have been detected in emission in the nuclear region of NGC 253 (Aalto et al. 2011). High-resolution SiO observations show bright emission resulting from large-scale shocks, as well as gas entrained in the nuclear outflow (García-Burillo et al. 2000). High-resolution observations of HCN and HCO⁺ $J=1\to 0$ show strongly centrally concentrated emission (Knudsen et al. 2007).

After the [C ii] 58 μm and [O i] 63 μm lines, the ¹²CO transitions are the most important cooling lines of the molecular gas in the interstellar medium (ISM). Therefore, several studies of the LSED of ¹²CO have been done to estimate the ambient conditions and excitation mechanisms of the molecular gas in the nuclear region of NGC 253 (e.g., Bradford et al. 2003; Bayet et al. 2004). Gas temperatures of T_kin ≈ 60 K and H₂ density of ∼10⁴ cm⁻³ were found to be sufficient to explain the velocity-integrated intensities of the lower-J (up to $J=7\to 6$) ¹²CO transitions, using a single-component (temperature/density) Large Velocity Gradient (LVG) model (Güsten & Philipp 2006b). More recent observations using SPIRE on Herschel (Pilbratt et al. 2010) provided a more extended ¹²CO LSED, including transitions up to $J=13\to 12$, which was reproduced using three gas components (Rosenberg et al. 2014). Three possible combinations of excitation mechanisms were explored, and all were found to be plausible explanations of the observed molecular emission. However, mechanical heating was found to be more dominant than CR heating in some of the models by Rosenberg et al. (2014).

We present a combined set of updated and extended Herschel SPIRE-FTS, PACS (Poglitsch et al. 2010), and HIFI (de Graauw et al. 2010) observations of the nuclear region of NGC 253. These observations are part of the Herschel EXtra GALactic (HEXGAL; PI: R. Güsten) Guaranteed Key Program. We also present data obtained with the Stratospheric Observatory For Infrared Astronomy (SOFIA; Young et al. 2012), as well as from the ground-based Atacama Pathfinder EXperiment (APEX⁸ ; Güsten et al. 2006a) telescope. Together, this set of data advances previous results on NGC 253 in the (sub-)millimeter and far- and mid-IR ranges, providing much more information about the excitation processes at play, and their effects on the molecular and atomic line emissions, than previously reported from ground-based observations and SPIRE spectra alone.

The organization of this article is as follows. In Section 2, we describe the observations and the procedure followed to reduce and extract the spectral and photometry data. The results (maps, spectra, and fluxes of detected and identified lines) are presented in Section 3. The analysis and modeling of the combined SPIRE and PACS continuum emission are discussed in Section 4. In Section 6, we propose a new non-LTE radiative transfer model for the ¹²CO SLED, and present the analysis and discussions of the model results and the interpretation of the line ratios with other molecular lines. A summary and final remarks are presented in Section 7.

A summary of all the observations that provided data used in this work and that are available in the Herschel Science Archive (HSA⁹ ) is shown in Table 1. All of the Herschel data presented here were obtained as part of the key program guarantee time HEXGAL: KPGT_rguesten_1 (P.I. Rolf Güsten). Herschel data were processed and reduced using the Herschel Interactive Processing Environment (HIPE¹⁰ ) v.14 and v.15 (Ott 2010). The difference in calibration obtained for the SPIRE line fluxes is less than 1% between these two versions, but it can be as large as 15%–20% with respect to older versions (e.g., HIPE v.11).

2.1. The APEX SABOCA and CHAMP⁺ maps

2.2. Correcting the SPIRE Intensity Levels

The SPIRE spectral cubes calibrated with the point-like source pipeline were used because they lead to lower noise levels, and the match between the two spectral bands is better. This is consistent with the large SPIRE beam compared with the estimated size of the emitting region, as discussed below.

Most works found in the literature present SPIRE spectra corrected using a photometric flux at one or more wavelengths, integrated in an aperture equivalent to the largest beam size (at the lowest frequency) of the SPIRE-FTS. Instead, we used the semiExtendedCorrector task (available in HIPE since v.11), which stitches together the long- and short-wavelength FTS bands (SLW and SSW, respectively). This tool corrects the SPIRE spectra by simulating a source size convolved to a ∼40'' beam. Details of this task and a description of the method used were reported by Wu et al. (2013), and our application of it is described in Appendix A. There are two main reasons to prefer this method: first, we can obtain an estimate for the size of the emitting region, and second, we identify frequency ranges in the SPIRE spectra that may still be affected by some calibration inaccuracies (Swinyard et al. 2014).

To cross-check this method, choose a source size and hence, the final spectra to use; we compare the continuum level of the corrected spectra with the actual dust continuum emission at given wavelengths, as observed with an equivalent beam size (or aperture). We first extracted the fluxes of the SPIRE photometric maps of NGC 253 (obs. ID 1342199387), obtained from the HSA, as well as the integrated flux at 350 μm obtained with the Submillimeter APEX Bolometer Camera (SABOCA) on the APEX telescope (Siringo et al. 2010). The SPIRE photometric maps were reprocessed with the pipeline for the large photometer maps provided in HIPE v.15. Details of the reprocessing can be found in Appendix A. Figure 1 shows the SPIRE (left) and SABOCA (right) 350 μm flux density maps. The SABOCA bolometer, with a higher spatial resolution (HPBW ∼ 8''), was used to map only the nuclear region of the galaxy.

Zoom In Zoom Out Reset image size

The 40'' annular sky aperture integrated fluxes are 278.7 ±19.0 Jy, 93.6 ± 10.8 Jy, and 20.1 ± 5.1 Jy for the 250, 350, and 500 μm images, respectively. The integrated APEX/SABOCA flux for an aperture of 40'' is 131.4 ± 11.5 Jy, i.e., ∼30% higher than the corresponding SPIRE flux. Note, however, that using the DAOphot algorithm with automatic aperture correction, the SPIRE flux at 350 μm would be 117.0 ± 13.5 Jy instead, more consistent, to within their (1σ) uncertainties, with the SABOCA flux. The difference in absolute values may be due to the different calibration schemes, the fact that an atmospheric contribution affects the SABOCA map, or that the SPIRE photometric map at 350 μm is also affected by small calibration uncertainties.

Fitting a 2D Gaussian distribution to the SABOCA map yields a beam-deconvolved source size (FWHM) of 173 × 92 (about 210 × 112 pc at an assumed distance of 2.5 Mpc; Houghton et al. 1997), which corresponds to an elliptical shape (shown in Figure 1, right) with eccentricity 0.85. This source size is about 24% smaller than the size (30'' × 16'') found by Weiß et al. (2008) from the APEX/LABOCA 870 μm map, with a larger beam size (192). However, the eccentricity of the latter case is the same (0.85), which indicates that the 2D Gaussian intensity distribution of the continuum emission is consistent at these two wavelengths with the two different beam sizes. The SPIRE-FTS spectra corrected assuming the source size estimated from the SABOCA map are shown in Figure 2.

Zoom In Zoom Out Reset image size

2.3. Extracting the PACS Spectra

2.4. The HIFI Spectra

2.5. The SOFIA GREAT and upGREAT Spectra

More than 60 lines (in absorption and emission) were detected and identified in the wavelength range (57–671 μm) covered by both SPIRE and PACS. The corrected SPIRE apodized spectrum of NGC 253 is shown in Figure 3. All of the line fitting and analysis, however, were done using the unapodized spectrum due to its higher spectral resolution, the less blended lines, and the more accurate fluxes obtained from fitting sinc functions compared to the about 5% less flux obtained when fitting Gaussians to the apodized spectrum (cf. SPIRE data reduction guide, Section 7.10.6 in version 3). We detected 35 lines in the SPIRE spectra, including a few unidentified lines not reported here. We detected eight H₂O lines in emission, and CH⁺ $J=1\to 0$, three OH⁺, and two o-H₂O⁺ lines in absorption, among others. The emission part of the OH⁺ P-Cygni feature is more evident at 907 GHz than the N = 1−0 line observed at 971 GHz, which is also observed and velocity resolved with HIFI (van der Tak et al. 2016). The fluxes (in units of W m⁻²) and the equivalent luminosities are summarized in Tables 2 and 3.

Zoom In Zoom Out Reset image size

Table 2.  Fluxes from the SPIRE-FTS Emission Lines Extracted from the 40'' Aperture Observations of NGC 253

Line	${\nu }_{\mathrm{rest}}$ a	Fluxb	Luminosityc
	(GHz)	(${10}^{-16}$ W m−2)	(104 L☉)
12CO $J=4\to 3$	461.041	13.48 ± 1.37	51.2 ± 6.5
12CO $J=5\to 4$	576.268	18.44 ± 1.86	70.0 ± 8.8
12CO $J=6\to 5$	691.473	18.72 ± 1.89	71.0 ± 9.0
12CO $J=7\to 6$	806.652	19.19 ± 1.94	72.8 ± 9.2
12CO $J=8\to 7$	921.800	17.67 ± 1.79	67.1 ± 8.5
12CO $J=9\to 8$	1036.912	16.29 ± 1.72	61.8 ± 8.0
12CO $J=10\to 9$	1151.985	13.41 ± 1.45	50.9 ± 6.7
12CO $J=11\to 10$	1267.014	10.72 ± 1.20	40.7 ± 5.5
12CO $J=12\to 11$	1381.995	7.79 ± 0.95	29.6 ± 4.2
12CO $J=13\to 12$	1496.923	5.69 ± 0.79	21.6 ± 3.4
13CO $J=5\to 4$	550.926	0.92 ± 0.25	3.5 ± 1.0
13CO $J=6\to 5$	661.067	0.74 ± 0.24	2.8 ± 0.9
13CO $J=7\to 6$	771.184	0.66 ± 0.24	2.5 ± 0.9
13CO $J=8\to 7$	881.273	0.58 ± 0.24	2.2 ± 0.9
C18O $J=5\to 4$	548.831	0.34 ± 0.22	1.3 ± 0.8
[C I] ${}^{3}{P}_{1}-{}^{3}{P}_{0}$	492.161	4.93 ± 0.56	18.7 ± 2.5
[C I] ${}^{3}{P}_{2}-{}^{3}{P}_{1}$	809.342	12.25 ± 1.25	46.5 ± 5.9
[N II] ${}^{3}{P}_{1}-{}^{3}{P}_{0}$	1460.977	20.63 ± 2.13	78.3 ± 10.0
o-H2O ${1}_{10}-{1}_{01}$	556.936	0.37 ± 0.22	1.4 ± 0.8
p-H2O ${2}_{11}-{2}_{02}$	752.033	2.97 ± 0.37	11.3 ± 1.6
p-H2O ${2}_{02}-{1}_{11}$	987.927	6.17 ± 0.76	23.4 ± 3.4
o-H2O ${3}_{12}-{3}_{03}$	1097.365	4.29 ± 0.62	16.3 ± 2.7
o-H2O ${3}_{12}-{2}_{21}$	1153.127	3.03 ± 0.54	11.5 ± 2.2
o-H2O ${3}_{21}-{3}_{12}$	1162.912	7.44 ± 0.87	28.2 ± 3.9
p-H2O ${4}_{22}-{4}_{13}$	1207.639	1.55 ± 0.47	5.9 ± 1.8
p-H2O ${2}_{20}-{2}_{11}$	1228.789	6.45 ± 0.78	24.5 ± 3.5
CH ${}^{2}{{\rm{\Pi }}}_{1/2}$ $J=3/2-1/2$ d	532.730	0.99 ± 0.25	3.8 ± 1.0
CH ${}^{2}{{\rm{\Pi }}}_{1/2}$ $J=3/2-1/2$ d	536.760	0.95 ± 0.25	3.6 ± 1.0
OH+ ${1}_{01}-{0}_{12}$ e	907.500	1.49 ± 0.45	5.7 ± 1.8

Notes.

^aObtained from the LAMDA, CDMS, and NASA/JPL databases. ^bThe flux errors include the statistical uncertainty of the instrument, 6% of the calibration uncertainty (Swinyard et al. 2014), and 8% of the uncertainty in the source size used to correct the continuum levels (Section 2.1). ^cLuminosity estimated assuming a flat space cosmology (H₀ = 70 $\mathrm{km}\,{{\rm{s}}}^{-1}$ Mpc⁻¹, ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$=0.73, ${{\rm{\Omega }}}_{M}$=0.27) and a distance of 3.5 ± 0.2 Mpc for NGC 253 (Rekola et al. 2005). The luminosity errors include the relative uncertainty of the respective fluxes and the distance of the galaxy, as well as a 5% uncertainty for the assumed cosmology model. ^dThese lines are blended with the HCN and HCO⁺ $J=6\to 5$ lines at 531.716 GHz and 535.062 GHz, respectively, as detected in the HIFI spectra by Rangwala et al. (2014). ^eEmission part of the OH⁺ P-Cygni feature. The flux in absorption centered at 909 GHz is shown in Table 3.

Download table as:  ASCIITypeset image

Table 3.  Fluxes from the SPIRE-FTS Absorption Lines Extracted from the 40'' Aperture Observations of NGC 253

Line	${\nu }_{\mathrm{rest}}$ a	Flux	Luminosityb
	(GHz)	(10−16 W m−2)	(104 L☉)
CH+ $J=1\to 0$	835.138	−1.71 ± 0.25	−6.5 ± 1.1
o-NH2 111−000	952.578	−1.60 ± 0.24	−6.1 ± 1.0
OH+ 101−012	909.159	−2.48 ± 0.49	−9.4 ± 2.0
OH+ 122−011	971.805	−5.63 ± 0.69	−21.4 ± 3.1
OH+ 112−012c	1033.119	−6.53 ± 0.76	−24.8 ± 3.4
p-H2O ${1}_{11}-{0}_{00}$	1113.343	−2.86 ± 0.54	−10.8 ± 2.2
o-H2O+ 111−000 ${J}_{\tfrac{3}{2}-\tfrac{1}{2}}$	1115.204	−4.57 ± 0.65	−17.4 ± 2.8
o-H2O+ ${1}_{11}-{0}_{00}$ ${J}_{\tfrac{1}{2}-\tfrac{1}{2}}$	1139.654	−3.16 ± 0.55	−12.0 ± 2.3
HF $J=1\to 0$	1232.476	−4.99 ± 0.63	−18.9 ± 2.8

Notes.

^aObtained from the LAMDA, CDMS, and NASA/JPL databases. ^bLuminosity estimated assuming a flat space cosmology (H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_λ = 0.73, Ω_M = 0.27) and a distance of 3.5 ± 0.2 Mpc for NGC 253 (Rekola et al. 2005). The luminosity errors include the relative uncertainty of the respective fluxes and the distance of the galaxy, as well as a 5% uncertainty for the assumed cosmology model. ^cThis line is likely blended with the OH⁺ 1₁₁−0₁₁ line at 1032.998 GHz. Both lines have comparable Einstein A coefficients (1.41 × 10⁻² s⁻¹).

Download table as:  ASCIITypeset image

We detected more than 30 lines in the PACS long-range SED spectrum. The individual emission lines (including their corresponding Gaussian fit, rms, and signal-to-noise ratio (S/N) are shown in Figure 4. The absorption lines detected are shown in Figure 5. In addition to several ionized species ([C ii], [N ii], [N iii] and [O iii]), we also detected five more ¹²CO lines, extending the ladder observed with SPIRE-FTS. We only detected an upper limit for the ¹²CO $J=19\to 18$ transition, and the $J=17\to 16$ transition is in a spectral region with very low S/N (with an uncertainty of 45%), so we do not trust the flux obtained for this transition. We also detected two OH doublet lines in emission and two doublets in absorption, as well as H¹⁸O in absorption. The fluxes and luminosities of all the detected (and identified) PACS lines are listed in Tables 4 and 5. There are 30 lines currently identified in the PACS spectra of NGC 253.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Fluxes of the Emission Lines Extracted from the PACS 5 × 5 Spectrum of NGC 253

Line	${\lambda }_{\mathrm{rest}}$ a	${\nu }_{\mathrm{rest}}$ a	Fluxb	Luminosityc
	(μm)	(GHz)	(10−16 W m−2)	(104 L☉)
[N III] 57 μm	57.320	5230.155	54.75 ± 9.00	207.8 ± 37.6
[N ii] 122 μm	121.888	2459.566	95.40 ± 14.36	362.1 ± 61.1
[C ii]	157.741	1900.537	479.87 ± 72.08	1821.2 ± 306.5
[O III] 88 μm	88.356	3392.991	73.32 ± 11.05	278.3 ± 47.0
[O i] ${}^{3}{P}_{0}-{}^{3}{P}_{1}$	145.525	2060.070	42.51 ± 6.41	161.3 ± 27.2
[O i] ${}^{3}{P}_{1}-{}^{3}{P}_{2}$	63.184	4744.775	347.60 ± 52.53	1319.2 ± 223.1
12CO $J=14\to 13$	186.000	1611.787	5.64 ± 1.01	21.4 ± 4.2
12CO $J=15\to 14$	173.630	1726.617	3.95 ± 0.75	15.0 ± 3.1
12CO $J=16\to 15$	162.812	1841.346	2.64 ± 0.58	10.0 ± 2.3
12CO $J=17\to 16$	153.267	1956.018	1.06 ± 0.49	4.0 ± 1.9
12CO $J=18\to 17$	144.780	2070.676	2.00 ± 0.48	7.6 ± 1.9
12CO $J=19\to 18$	137.200	2185.076	1.45 ± 0.58	5.5 ± 2.2
p-H2O 322−313	156.190	1919.409	2.09 ± 0.48	7.9 ± 1.9
p-H2O 313−202	138.530	2164.098	1.07 ± 0.30	4.0 ± 1.2
p-H2O 331−322	126.710	2365.973	3.21 ± 1.00	12.2 ± 3.9
o-H2O 303−212	174.630	1716.729	5.37 ± 1.00	20.4 ± 4.1
o-H2O 423−414	132.410	2264.122	1.24 ± 0.36	4.7 ± 1.4
o-H2O 441−432	94.705	3165.533	3.85 ± 1.02	14.6 ± 4.0
OH ${}_{\tfrac{1}{2}\tfrac{3}{2}-\tfrac{1}{2}\tfrac{1}{2}}$	163.400	1834.715	13.30 ± 2.14	50.5 ± 9.0
OH ${}_{\tfrac{1}{2}\tfrac{3}{2}-\tfrac{1}{2}\tfrac{1}{2}}$	163.120	1837.865	12.51 ± 2.03	47.5 ± 8.5
OH ${}_{\tfrac{1}{2}\tfrac{1}{2}-\tfrac{3}{2}\tfrac{3}{2}}$	79.179	3786.257	14.57 ± 2.87	55.3 ± 11.7
OH ${}_{\tfrac{1}{2}\tfrac{1}{2}-\tfrac{3}{2}\tfrac{3}{2}}$	79.115	3789.301	6.68 ± 1.75	25.3 ± 6.9

Notes.

^aObtained from the LAMDA, CDMS, and NASA/JPL databases. ^bThe flux errors include the statistical uncertainties of the instrument and 15% of the calibration uncertainty, adopted for the 5 × 5 spaxel spectra. ^cLuminosity estimated assuming a flat space cosmology (H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_λ = 0.73, Ω_M = 0.27) and a distance of 3.5 ± 0.2 Mpc for NGC 253 (Rekola et al. 2005). The luminosity errors include the relative uncertainty of the respective fluxes and the distance of the galaxy, as well as a 5% uncertainty for the assumed cosmology model.

Download table as:  ASCIITypeset image

Table 5.  Fluxes of the Absorption Lines Extracted from the PACS 5 × 5 Spectrum of NGC 253

Line	${\lambda }_{\mathrm{rest}}$ a	${\nu }_{\mathrm{rest}}$ a	Fluxb	Luminosityc
	(μm)	(GHz)	(10−16 W m−2)	(104 L☉)
p-H2O 322−211	89.988	3331.479	−6.04 ± 1.19	−22.9 ± 4.8
p-H2O 331−220	67.090	4468.512	−4.20 ± 0.92	−15.9 ± 3.7
o-H2O 212−101	179.526	1669.906	−15.40 ± 2.54	−58.4 ± 10.6
o-H2O 221−110	108.070	2774.058	−4.85 ± 1.00	−18.4 ± 4.1
o-H2O 321−212	75.380	3977.082	−15.05 ± 2.47	−57.1 ± 10.3
o-H2O 330−221	66.440	4512.228	−6.87 ± 1.43	−26.1 ± 5.8
o-H2O 432−321	58.700	5107.197	−2.72 ± 0.92	−10.3 ± 3.6
OH ${}_{\tfrac{3}{2}\tfrac{5}{2}-\tfrac{3}{2}\tfrac{3}{2}}$	119.440	2509.984	−60.74 ± 9.69	−230.5 ± 40.7
OH ${}_{\tfrac{3}{2}\tfrac{5}{2}-\tfrac{3}{2}\tfrac{3}{2}}$	119.230	2514.405	−67.72 ± 10.84	−257.0 ± 45.5
OH ${}_{\tfrac{3}{2}\tfrac{7}{2}-\tfrac{3}{2}\tfrac{5}{2}}$	84.600	3543.646	−10.24 ± 3.60	−38.9 ± 14.0
OH ${}_{\tfrac{3}{2}\tfrac{7}{2}-\tfrac{3}{2}\tfrac{5}{2}}$	84.420	3551.202	−10.57 ± 2.01	−40.1 ± 8.2
H18O ${}_{\tfrac{1}{2}\tfrac{1}{2}-\tfrac{3}{2}\tfrac{3}{2}}$	79.080	3791.002	−3.56 ± 2.90	−13.5 ± 11.1
CH ${2}_{\tfrac{3}{2}2-}-{1}_{\tfrac{1}{2}1+}$	149.390	2006.771	−9.66 ± 1.52	−36.7 ± 6.4
CH ${2}_{\tfrac{3}{2}2+}-{1}_{\tfrac{1}{2}1-}$	149.092	2010.787	−9.82 ± 1.54	−37.3 ± 6.5

Notes.

^aObtained from the LAMDA, CDMS, and NASA/JPL databases. ^bThe flux errors include the statistical uncertainties of the instrument and 15% of the calibration uncertainty, adopted for the 5 × 5 spaxel spectra. ^cLuminosity estimated assuming a flat space cosmology (H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_λ = 0.73, Ω_M = 0.27) and a distance of 3.5 ± 0.2 Mpc for NGC 253 (Rekola et al. 2005). The luminosity errors include the relative uncertainty of the respective fluxes and the distance of the galaxy, as well as a 5% uncertainty for the assumed cosmology model.

Download table as:  ASCIITypeset image

The HIFI maps of ¹²CO $J=9\to 8$, [C I] ${}^{3}{P}_{2}-{}^{3}{P}_{1}$, and [C ii] are shown in Figure 6. The spectra obtained by convolving the maps with an equivalent (HPBW) 40'' beam is also shown in order to compare them with the corresponding SPIRE and PACS data. Although the horizontal and vertical polarization spectra were used independently to create the final maps, the grid maps of Figure 6 show the average spectrum of the two polarizations for clarity. The spectra of the single-pointing observations, at their respective beam resolutions, can be found in Appendix C.

Zoom In Zoom Out Reset image size

Even though we clearly see more than one component in the velocity-resolved HIFI lines, we fit a single Gaussian component to the spectra, since this fit is good enough to extract the total flux and the width (FWHM) of the lines. In Table 6, we list the FWHM and velocity-integrated temperatures (in T_mb scale). For the HIFI maps, we also list the FWHM and total flux (W m⁻²) of the spectra convolved with an equivalent (HPBW) 40'' beam, as well as the flux ratio for the lines that were observed with other instruments.

Table 6.  Line Widths and Fluxes^a from HIFI Spectra of NGC 253

Line	FWHMb	FWHMc	Intensityd	Fluxc	Instruments
	(km s−1)	(km s−1)	(K km s−1)	(10−16 W m−2)	Ratio
12CO $J=5\to 4$	195.8 ± 24.6		215.9 ± 27.1
12CO $J=6\to 5$	188.3 ± 23.7		201.4 ± 25.3
12CO $J=9\to 8$	157.8 ± 19.8	164.3 ± 20.7	88.4 ± 11.1	14.6 ± 1.8	0.94e
13CO $J=5\to 4$	171.9 ± 21.6		12.0 ± 1.5
13CO $J=6\to 5$	191.8 ± 24.3		12.9 ± 1.6
13CO $J=9\to 8$ f	⋯		⋯
[C i] ${}^{3}{P}_{1}-{}^{3}{P}_{0}$	193.7 ± 24.4		71.7 ± 9.0
[C I] ${}^{3}{P}_{2}-{}^{3}{P}_{1}$	174.0 ± 21.9	184.5 ± 23.2	98.9 ± 12.4	10.7 ± 1.3	0.86e
[C ii]	195.7 ± 24.6	205.7 ± 25.9	637.5 ± 80.2	406.7 ± 51.2	0.79g
[C ii]			264.9 ± 33.0	56.6 ± 7.0	1.45h

Notes.

^aThe errors quoted include the rms obtained from the baseline subtraction and uncertainties of 6% in the sideband ratio, 3% in the planetary model, 10% in the beam efficiency, 2% in the pointing, and 3% in the correction for standing waves. ^bFWHM obtained from a single-component Gaussian fit of the corresponding HIFI spectrum. ^cFWHM and flux convolved to a 40'' HPBW from the ¹²CO J = 9–8 and [C I] ${}^{3}{P}_{2}-{}^{3}{P}_{1}$ HIFI maps. ^dVelocity-integrated temperatures obtained from the Gaussian fit of the HIFI spectra at their respective beam sizes. ^eRatio between the HIFI and SPIRE fluxes for the 40'' HPBW spectra. ^fThe ¹³CO $J=9\to 8$ line was observed but not clearly detected since the spectrum is severely affected by standing waves. ^gRatio between the HIFI and PACS fluxes for the 40'' HPBW spectra. ^hRatio between the HIFI and upGREAT fluxes observed at about the offset position (−115, −82) for the 151 HPBW spectra.

Download table as:  ASCIITypeset image

For the [C ii] map, we obtained a source size of about 279 × 183 (average of ∼231) from a 2D Gaussian fit. Unfortunately, the [C i] and ¹²CO maps are too small (only 5 × 5 pixels) to fit a 2D Gaussian (the fitting procedure does not converge). So, we assumed that the [C i] emitting region has the same size as that of [C ii]. In the case of line ¹²CO $J=9\to 8$, instead, we used the ¹²CO $J=6\to 5$ map, with a resolution (HPBW) of 94, obtained with CHAMP⁺ (Kasemann et al. 2006; Güsten et al. 2008) on APEX, assuming the emitting regions of these two lines have similar sizes (a discussion on the sizes can be found in the next section). A 2D Gaussian fit of the $J=6\to 5$ map gives a CO source size of 208 × 125 (or ∼167 on average), which is consistent with the source size found from the SABOCA map (cf. Figure 1) when considering a 10% uncertainty in the estimates.

The fluxes reported in Table 6 were obtained using the average source sizes estimated for ¹²CO $J=9\to 8$, [C i], and [C ii]. These fluxes are, respectively, 6%, 14%, and 21% smaller than the corresponding fluxes obtained with SPIRE and PACS. The larger difference in the PACS line may be because the PACS calibration could be affected by the bright [C ii] line and the continuum level around it see Section 4. On the other hand, the HIFI ¹²CO $J=6\to 5$ integrated intensity is only about 10 K km s⁻¹ (5%) higher than the intensity obtained from the APEX/CHAMP⁺ map convolved with the same equivalent beam (HPBW = 367) of HIFI at the ¹²CO $J=6\to 5$ frequency. Considering the uncertainties of the data from all instruments, there is no significant difference between the fluxes of, for instance, ¹²CO $J=6\to 5$ from APEX/CHAMP⁺ and Heschel/SPIRE (Tables 8 and 2, respectively). Similarly, the fluxes of ¹²CO $J=9\to 8$ obtained with SPIRE and HIFI (Tables 2 and 6, respectively) are practically the same, given the uncertainties. On the other hand, the [C ii] flux obtained with PACS is 21% more than that obtained with HIFI. Since the uncertainties of both instruments are similar (15% and 13% for PACS and HIFI, respectively), and given that the HIFI spectra of [C ii] do not have a fully covered baseline (due to the relatively short bandwidth of HIFI at ∼1.9 THz), we conclude that the [C ii] flux obtained with PACS is more reliable. Unfortunately, we cannot yet compare the PACS flux with the SOFIA/upGREAT flux because we did not manage to map the full central region of NGC 253 in the [C ii] line with upGREAT. But we can compare the spectrum of the central pixel of the latter with the associated HIFI spectrum, as discussed below.

Table 7.  Line Fluxes from SOFIA/upGREAT and APEX/CHAMP⁺ Toward the Southwest Position in the Nuclear Region of NGC 253

Line	Intensitya	Fluxb
	(K km s−1)	(10−16 W m−2)
12CO $J=6\to 5$	398.91 ± 39.89	4.10 ± 0.41
12CO $J=11\to 10$	24.29 ± 2.55	3.47 ± 0.36
[C ii] 158 μm	182.35 ± 20.18	38.95 ± 4.31

Notes.

^aThe errors quoted include the rms obtained from the baseline subtraction and 10% accounting for calibration and pointing uncertainties. ^bFluxes were estimated using the corresponding beam sizes of 151 for [C ii] and ¹²CO $J=6\to 5$, and 227 for ¹²CO $J=11\to 10$. ^cThe ¹²CO $J=11\to 10$ flux would be about 56% smaller if 151 is considered instead.

Download table as:  ASCIITypeset image

The footprint of the SOFIA/GREAT/upGREAT observations and the spectra obtained are shown in Figure 7, and the fluxes observed with the central pixel are reported in Table 7. The central pixel of the upGREAT array correspond to the [C ii] line observed at the offset position (−115, −82) SW from the nuclear region of NGC 253. This location encloses part of the densest gas observed in the HCN $J=1\to 0$ high-resolution map reported by Paglione et al. (2004). The rotation of the gas can be seen in the (shifted to higher velocities) shape of the velocity-resolved [C ii], CO $J=11\to 10$, and CO $J=6\to 5$ lines shown in Figure 8. The CO $J=6\to 5$ line observed with APEX/CHAMP⁺ was convolved to the same 151 resolution of the [C ii] line. We convolved the [C ii] map of HIFI to the 151 HPBW of SOFIA/upGREAT for comparison. The obtained HIFI flux of [C ii] is about 45% brighter than the value obtained from upGREAT. Such a large difference may be due to the different calibration schemes and relative pointing errors (after regridding, the HIFI spectrum is about (15, 06) closer to the central region than the upGREAT spectrum), but is mostly due to the difficulty and uncertainty of fitting a good baseline to the HIFI spectra due to its narrower instant bandwidth. We also note that the beam coupling efficiencies of these instruments differ significantly. While for SOFIA/upGREAT we estimated a 70% efficiency, Herschel/HIFI achieved only 57% efficiency when the [C ii] map was observed.

The CHAMP⁺ fluxes, convolved with different beam sizes, are listed in Table 8. The effect of the convolution on the line profiles is shown in the bottom panel of Figure 9. Contrary to what could be expected, the dynamical range (or full width at zero intensity) of the line remains the same (∼420 km s⁻¹), while the FWHM widens with larger beams, not because the beams cover emitting regions with exceeding kinematical components not seen at the central region (as shown in Figure 9, middle panel), but just because the peak temperature of the line decreases (due to the beam smearing effect). This is an effect that needs to be taken into account when interpreting the line shapes from extragalactic observations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dust properties in NGC 253 have been investigated by Radovich et al. (2001), Melo et al. (2002), and Weiß et al. (2008) using the mid- and far-IR data from ISOPHOT and IRAS, and the submillimeter data from LABOCA on APEX (Siringo et al. 2009), respectively. We have reanalyzed the dust temperatures, mass, optical depths, and column densities using the SPIRE and PACS photometry fluxes, complemented at the shorter wavelengths by archival data from Spitzer/MIPS (24 μm, AOR: 22610432) and MSX (21 μm, 15 μm, 12 μm and 8 μm, bands E, D, C and A, respectively; only these four images are available for NGC 253 in the MSX data archive¹⁵ ).

Following Weiß et al. (2008) and Vlahakis et al. (2005), the dust emission was modeled using the gray body formulation

$$\begin{eqnarray}&&{S}_{\nu }=(1-{e}^{-{\tau }_{\nu }})({B}_{\nu }({T}_{i})-{B}_{\nu }({T}_{\mathrm{cmb}})){{\rm{\Omega }}}_{s}{{\rm{\Phi }}}_{i},\end{eqnarray} \tag{ 1 }$$

where B_ν is the Planck function, τ_ν the dust optical depth, Ω_s the source solid angle, T_cmb = 2.73 K the cosmic microwave background temperature, and T_i and Φ_i the dust temperature and beam area filling factor of each component. The dust optical depth was computed as

$$\begin{eqnarray}&&{\tau }_{\nu }={k}_{d}(\nu ){M}_{\mathrm{dust}}/({D}^{2}{{\rm{\Omega }}}_{s}{{\rm{\Phi }}}_{c}),\end{eqnarray} \tag{ 2 }$$

where M_dust is the dust mass, D the distance to the source, Φ_c the filling factor of the coldest component, and following Weiß et al. (2008) and Kruegel & Siebenmorgen (1994), the adopted dust absorption coefficient ${k}_{d}(\nu )$ was

$$\begin{eqnarray}&&{k}_{d}{(\nu )=0.04(\nu /250\mathrm{GHz})}^{\beta }\ [{{\rm{m}}}^{2}\,{\mathrm{kg}}^{-1}],\end{eqnarray} \tag{ 3 }$$

with ν in GHz and β = 2 (Priddey & McMahon 2001). We used the flux observed at 500 μm (the most optically thin emission in our data set) to compute the dust mass for each component

$$\begin{eqnarray}&&{M}_{\mathrm{dust},i}=\displaystyle \frac{{F}_{500}{D}^{2}}{{k}_{d,500}}{({B}_{500}({T}_{i})-{B}_{500}({T}_{\mathrm{cmb}}))}^{-1}.\end{eqnarray} \tag{ 4 }$$

The total dust mass was estimated to be M_dust = 3.0 ±0.9 × 10⁶ M_☉, while the total gas mass is M_gas = 4.5 ± 1.3 ×10⁸, assuming the same gas-to-dust mass ratio of 150 used by Weiß et al. (2008) and assuming a shorter distance of 2.5 Mpc.

The dust temperatures and masses depend on the underlying source solid angle and the area filling factor of each component. We used the source size of 173 × 92 (deconvolved from the SABOCA map). This is about half the size (30'' × 17'') derived by Weiß et al. (2008) from a 80'' beam. Note that Weiss et al. adopted a distance of 2.5 Mpc, estimated assuming that the observed stars in NGC 253 were similar to the asymptotic giant branch stars in Galactic globular clusters (Davidge & Pritchet 1990; Davidge et al. 1991; Houghton et al. 1997). Instead, we used a more recent estimate of 3.5 ± 0.2 Mpc based on models of planetary nebula accounting for dust (Rekola et al. 2005), which is consistent with estimates based on measurements of the magnitude of the tip of the red giant branch (Mouhcine et al. 2005).

The dust SED fit is shown in Figure 10, and the parameters are summarized in Table 9. The uncertainties of the dust temperatures consider a total of 10% error for the source size and filling factors adopted for each component. Given the uncertainties, the temperature ∼37 K of the cold component is practically the same as that (30–35 K) found by Weiß et al. (2008). Our second component, on the other hand, is considerably (∼16 K) higher than the one found before. This may be due to the fact that we include fluxes at shorter wavelengths that were not used by Weiß et al. (2008) in their SED fit. Our third component, however, has the highest flux uncertainties of 30% from the MSX data. This is because the flux at 21 μm does not follow the trend of the MIPS 24 μm flux (indicating a different flux scale between these instruments) and because the flux at 8 μm and (to a lesser extend at 12 μm) contain emission from PAHs (e.g., Povich et al. 2007 and references therein), which we did not correct for since they are difficult to assess in an unresolved source. Hence, the dust temperature of the third component should be considered an upper limit. Besides, a spectral index β = 2, as assumed above, may not be appropriate for the warmest dust. Leaving it as a free parameter only for the third component would lead to a poorly constrained value anyway because of the contamination by PAHs. Hence, we did not investigate further on this matter.

Zoom In Zoom Out Reset image size

Table 9.  Dust Continuum SED Fit Parameters of NGC 253

	Component Parameters
Quantity	Cold	Warm	Hot
Φa	5 × 10−1	1 × 10−2	1 × 10−4
Tdust [K]	36.6 ± 3.7	70.0 ± 7.0	187.7 ± 37.5
G0	$3.4\times {10}^{5}$	8.7 × 106	1.2 × 109
Mdust [106 M⊙]b	1.0 ± 0.3	0.4 ± 0.1	0.14 ± 0.04

Notes.

^aUncertainties in the filling factors are of the order of 10%. ^bDust mass obtained using a source size solid angle of ${{\rm{\Omega }}}_{s}=17.3\times 9.2$ arcsec² as obtained from a 2D Gaussian intensity distribution fit of the SABOCA map.

Download table as:  ASCIITypeset image

Assuming FUV heating, we also estimated the FUV flux G₀, in units of the equivalent Habing flux (1.6 × 10⁻³ erg cm⁻² s⁻¹), from Hollenbach et al. (1991, their Equation (7)), as

$$\begin{eqnarray}&&{G}_{0}=3.7\times {10}^{-3}{\tau }_{100\mu {\rm{m}}}{T}_{d}^{5},\end{eqnarray} \tag{ 5 }$$

using the dust opacity at 100 μm (τ_{100 μm} ≈ 1.4) estimated from the dust SED fit of each component and assuming that the dust temperatures are similar to the actual equilibrium dust temperatures (T₀) at the surface of their respective emitting regions.

Assuming most hydrogen is in molecular form, we can also estimate the molecular hydrogen column density from the dust opacity and absorption coefficient following the formulation by Kauffmann et al. (2008, their Equation A.9) as

$$\begin{eqnarray}&&N({{\rm{H}}}_{2})={10}^{-4}\times \displaystyle \frac{{\tau }_{\nu }}{{\mu }_{{{\rm{H}}}_{2}}{m}_{{\rm{H}}}{k}_{\nu }}[{\mathrm{cm}}^{-2}],\end{eqnarray} \tag{ 6 }$$

where m_H is the hydrogen atom mass (in kg) and μ_H2 is the molecular weight per hydrogen molecule. We use a value of 2.8 for the latter, which is the value needed to compute particle column densities, while the classical value of 2.33, used sometimes in the literature, actually corresponds to the mean molecular weight per free particle (μ_p), which is used to estimate other quantities, like thermal gas pressure. The factor 10⁻⁴ is used to convert the dust absorption coefficient k_ν from units of m² kg⁻¹ (from Equation (3)) to cm² kg⁻¹. Combining Equation (6) with Equations (2) and (3), we obtained N(H₂) = (5.2 ± 2.3) × 10²¹ cm⁻².

We can also estimate the visual extinction A_V (mag) from the standard conversion factor $N({{\rm{H}}}_{2})=9.4\,\times \,{10}^{20}{A}_{V}$ from which the atomic hydrogen column density can be estimated using the relation N(H) = 2.21 × 10²¹ A_V found for the Milky Way (Güver & Özel 2009). We obtained A_V = 5.5 ± 2.5 mag and $N({\rm{H}})=(1.2\pm 0.5)\times {10}^{22}$ cm⁻². All observed fluxes, dust absorption coefficients, and optical depths are summarized in Table 10.

Table 10.  Flux and Dust Properties at the Observed Wavelengths

Wavelength	Observed Flux	${k}_{\nu }$ a	${\tau }_{\nu }$ b
(μm)	(Jy)	(m2 kg−1)
500	20.1 ± 5.1	0.23	0.06 ± 0.02
350	93.6 ± 10.8	0.47	0.11 ± 0.03
250	278.7 ± 19.0	0.92	0.22 ± 0.06
160	914.1 ± 69.5	2.25	0.54 ± 0.15
100	1383.4 ± 102.7	5.75	1.39 ± 0.39
70	1271.0 ± 94.9	11.74	2.84 ± 0.80
24	45.0 ± 5.3	99.86	24.19 ± 7.11
21	53.7 ± 7.3	130.43	31.60 ± 9.54
15	21.7 ± 4.7	255.65	61.93 ± 21.35
12	17.6 ± 4.2	399.45	96.77 ± 34.83
8	7.7 ± 2.8	898.76	217.72 ± 98.03

Notes.

^aThe uncertainties in the absorption coefficients are assumed to be of the order of 10%. ^bThe errors in the optical depths consider the uncertainties of the temperatures of the three components and 10% uncertainties in the source size solid angle and the corresponding uncertainty of the distance to NGC 253, d = 3.52 ± 0.18 Mpc (Rekola et al. 2005).

Download table as:  ASCIITypeset image

The formation of hydrogen fluoride (HF) is dominated by a reaction of F with H₂, making the HF/H₂ abundance ratio more reliably constant than ¹²CO/H₂, specially for clouds of small extinction A_v (Neufeld et al. 2005). Therefore, HF has been proposed as a potentially sensitive probe of the total column density of the diffuse molecular gas (e.g., Neufeld et al. 2005; Monje et al. 2011).

Because of its very large A coefficient (A₁₀ = 2.42 × 10⁻² s⁻¹), this transition is generally observed in absorption (e.g., Neufeld et al. 1997, 2005, 2010; Phillips et al. 2010; Sonnentrucker et al. 2010; Monje et al. 2011; Rangwala et al. 2011; Kamenetzky et al. 2012; Pereira-Santaella et al. 2013). This high A coefficient translates into a simple excitation scenario, where most HF molecules are expected to be in the ground J = 0 state from where they can be excited into the J = 1 state by absorbing a photon at 1232.5 GHz under ambient conditions common to the diffuse and even dense ISM. Only an extremely dense region ($n({{\rm{H}}}_{2})\,\gt {10}^{9}$ cm⁻³, at ∼50 K), with a strong radiation field, could excite HF and generate a $J=1\to 0$ feature in emission (e.g., Neufeld et al. 1997, 2005; van der Werf et al. 2010; Spinoglio et al. 2012; Pereira-Santaella et al. 2013). For a more extended reference list, see van der Wiel et al. (2016, their Section 1).

From Equation (3) in Neufeld et al. (2010), and assuming all HF molecules are in the ground state, we can estimate the total HF column density from the absorption optical depth as

$$\begin{eqnarray}&&\int \tau {dv}=\displaystyle \frac{{A}_{{ul}}{g}_{u}{\lambda }^{3}}{8\pi {g}_{l}}N(\mathrm{HF})\end{eqnarray} \tag{ 7 }$$

where g_u = 3 and g_l = 1, which yields $\int \tau d\nu =4.16\,\times {10}^{-13}N(\mathrm{HF})$ cm² km s⁻¹. The optical depth of HF can be estimated from a double sideband receiver as τ =ln(2F_l/F_c − 1), with F_l/F_c the line-to-continuum ratio (Neufeld et al. 2010). In the case of SPIRE (a single-sideband spectrometer), τ can simply be estimated as τ = −ln(F_l/F_c) (cf. Kamenetzky et al. 2012; van der Wiel et al. 2016, their Section 4.1; who discuss the caveat of line smearing by a spectrometer that does not resolve the spectral profile). Figure 11 shows the estimated optical depth of HF (bottom panel) at each frequency element. In order to reduce the uncertainties and noise (ringing effect) introduced by the sinc convolution of the SPIRE-FTS, we first fit all of the prominent ¹²CO, [N ii], and H₂O lines (including the near by p-H₂O 2₂₀−2₁₁ at 1228.8 GHz) and then subtract their combined fluxes from the SSW band to produce the residual spectrum (normalized by the continuum) of HF $J=1\to 0$ shown in Figure 11 (top panel).

Zoom In Zoom Out Reset image size

Integrating the optical depth (of the normalized absorption feature below unity), we find a 40'' beam averaged column density N(HF) ≈ (1.07 ± 0.11) × 10¹⁴ cm⁻². The uncertainty for this column was estimated as the fraction (∼0.105) between the rms value (∼5.07 Jy), computed for the residual spectrum (around the HF line) between 1220 GHz and 1244 GHz, and the peak flux (∼48.23 Jy) of the HF absorption feature. This column density is a factor of ∼2.3 lower than the HF column density derived from the velocity-resolved HIFI spectrum of HF (Monje et al. 2014). However, the latter shows blueshifted absorption and redshifted emission (i.e., a P-Cygni profile), which are unresolved in the SPIRE spectrum. The P-Cygni profile of HF is suggestive of an outflow of molecular gas with a mass of ∼10⁷ M_☉ and an outflow rate of ∼6.4 M_☉ yr⁻¹ (Monje et al. 2014), which is in agreement with the outflow rate derived from the ¹²CO $J=1\to 0$ high-resolution map obtained with ALMA (Bolatto et al. 2013).

Because of the unresolved line profile, we quote our estimated HF column density as a lower limit. From the N(HF)/N(H₂) = 2.94 × 10⁻⁸ abundance ratio observationally determined by Indriolo et al. (2013; for the warm component of AFGL 2136 IRS 1, their Table 3), which is similar to the value 3.6 × 10⁻⁸ predicted by Neufeld & Wolfire (2009), we obtain a molecular hydrogen column density of $(3.64\pm 0.37)\,\times$ 10²¹ cm⁻² for the 40'' beam. This hydrogen column density is comparable to the column obtained in Section 4 from the dust emission at 100 μm (cf. Table 10). Besides the unresolved line profile of HF, there are other uncertainties to be considered in this calculation. First, the column density we derive is also a lower limit of the total column density, since we only observe the HF gas in front of the continuum emission. Second, the molecular abundance, and whether or not all HF molecules are truly in the ground state, is a debatable assumption since non-equilibrium chemistry could be at play in the environment with enhanced CR density of the nuclear region of NGC 253.

6.1. Note on the CO Line Widths

6.2. Non-LTE Excitation Analysis

Following previous works in the literature, we used the radiative transfer code RADEX¹⁶ (van der Tak et al. 2007) to explore a wide range of possible excitation conditions that can lead to the observed line fluxes of a particular molecule. Those line intensities are sensitive to the kinetic temperature (T_k), the volume density of the collision partner (n(H₂)), and the column density per line width (N/ΔV). For our analysis, we use only H₂ as collision partner since it is the most abundant molecule and has the largest contribution to the excitation of the CO lines. The code uses a uniform temperature and density of the collision partner to model a homogeneous sphere. Therefore, our analysis is not depth dependent. RADEX assumes the LVG (large velocity gradient/expanding sphere) formalism for the escape probability calculations. Hence, these models can only reproduce a clump that represents the average physical conditions of the gas from which the CO emission emerges. This is a well-fitted model for single-dish observations of unresolved emissions convolved with the telescope beams. The physical conditions were modeled using the collisional data available in the LAMDA¹⁷ database (Schöier et al. 2005). The collisional rate coefficients for ¹²CO and H₂O are adopted from Yang et al. (2010) and Daniel et al. (2011), respectively.

For the volume density, we explored ranges between 10² and 10⁷ cm⁻³, the kinetic temperature varies from 4 to 300 K, and the column density per line width lies between 10¹⁰ and 10²⁰ cm⁻² km⁻¹ s. In order to obtain the actual column density, the values reported must be multiplied by the local velocity dispersion (line width) of a single cloud. For comparison, a ${\rm{\Delta }}V=23$ km s⁻¹ was derived for the nuclear region of the active galactic nucleus (AGN)-driven galaxy NGC 1068 from high-resolution maps (Schinnerer et al. 2000). Since we do not have a good estimate for NGC 253, a conservative value of ΔV = 10 km s⁻¹ was adopted.

Since the optical depths obtained from the dust SED fit (Section 4) are not negligible (i.e., the dust emission is not optically thin in the whole frequency/wavelength range), and considering that the gas and dust must be well mixed in the emitting region, we modified the original RADEX code in order to include a more representative background emission I_bg(ν) as a diluted blackbody radiation field, in a similar way to that done by Poelman & Spaans (2005) and Pérez-Beaupuits et al. (2009). We considered the first two dust components at 37 K and 70 K (as estimated in Section 4) as well as the contribution from the cosmic microwave background at T_CMB = 2.73 K, according to the following equation

$$\begin{eqnarray}{I}_{\mathrm{bg}}(\nu ) & = & {B}_{\nu }({T}_{\mathrm{CMB}})\\ & & +(1-{e}^{-{\tau }_{\nu }})[{B}_{\nu }({T}_{c})+{f}_{w}{B}_{\nu }({T}_{w})],\end{eqnarray} \tag{ 8 }$$

where ${B}_{\nu }$ is the Planck function and the dust optical depth τ_ν is computed for each transition line using Equation (2), with a fixed dust mass M_dust = 3 × 10⁶ M_☉ estimated from the 500 μm photometric flux (Section 4). The factor f_w corresponds to the relative contribution of the warm component with respect to the cold component and is defined as the ratio between the corresponding area filling factors, f_w = Φ_w/Φ_c = 0.02 (cf. Table 9). The contribution factor f_w is needed in order to mimic the observed dust continuum emission in the spherical clump. Otherwise, the warm dust component at T_w = 70 K would dominate the background radiation field in the radiative transfer calculations, which would not be realistic. To be very rigorous, the second term of Equation (8) should be multiplied by a geometrical dilution factor η_d, which indicates the fraction of the dust emission actually seen by the molecules. However, we do not have a way to constrain this parameter from the convolved (unresolved) emission of the entire nuclear region of NGC 253, which was collected by the single dish of Herschel. Hence, for simplicity we assume η_d = 1, which is equivalent to assuming that the dust and the gas arise from the same volume.

The spectrum in the millimeter regime is usually dominated by the cosmic background blackbody radiation field at 2.73 K, which peaks at 1.871 mm. Therefore, this component of the radiation field of Equation (8) is generally considered (in the literature) to dominate the radiative excitation of the lower-J levels of heavy molecular rotors, such as CO, CS, HCN, HCO⁺, and H₂CO. Hence, there seems to be a general agreement in the (sub-)millimeter astronomy community that knowing the specific background radiation field of a single molecular cloud (or an ensemble of clouds, as in the case of extragalactic astronomy) is not really needed.

On the other hand, the far- and mid-infrared radiation fields (mainly from dust emission, especially in circumstellar material or in star-forming regions) are important for molecules with widely spaced rotational energy levels (e.g., the lighter hydrides OH, H₂O, H₃O⁺, and NH₂), as well as for the higher-J levels of the heavy rotors mentioned before. Since the dust is usually at higher temperatures than 2.73 K, its diluted blackbody radiation field will peak at shorter wavelengths (cf. Figure 10), increasing the radiative excitation of the higher-J levels and, hence, leaving fewer available molecules to populate the lower-J levels. This effect is particularly important for Herschel observations, in which several of the higher-J levels in the far- and mid-infrared regimes have become available for a number of molecules.

The actual effect of a background radiation field (including dust emission) on the redistribution among rotational levels depends on the local ambient conditions of the emitting gas. That is because at high densities (or temperatures), the collisions are expected to dominate the excitation of the mid- and high-J levels of molecules, such as CO, while at lower densities (or temperatures), the radiative excitation, as well as the spontaneous decay from higher-J levels, is expected to be the dominant component driving the redistribution of the level populations.

To demonstrate this, Figure 12 shows the fluxes (W m⁻²) of several transitions of the ¹²CO and ortho-H₂O molecules, considering (F₁) and not considering (F₀) the dust emission in the background radiation field (cf. Equation (8)), for different volume densities and kinetic temperatures. The ratio between these two fluxes is shown in the inset. Column densities (per line width) N/ΔV = 10¹⁷ and N/ΔV = 10¹⁶ cm⁻² km⁻¹ s were used for ¹²CO and o-H₂O, respectively. The original RADEX fluxes were scaled up assuming an average line width ΔV = 190 km s⁻¹ (from the average FWHM of the ¹²CO and ¹³CO $J=6\to 5$ HIFI lines; cf. Table 6) for each transition.

Zoom In Zoom Out Reset image size

Although the difference in fluxes of the ¹²CO transitions is barely noticed in the logarithmic scale, the absolute fluxes obtained without using the dust emission in the background field are more than 40% brighter than the fluxes obtained when the dust emission is included in the background field, for low densities (10³ cm⁻³) and moderate temperatures (100 K). On the other hand, at high densities (10⁶ cm⁻³) and relatively low temperatures (50 K), the fluxes F₀ of the lower-J levels (J_up < 5) are just a few percent brighter than the fluxes F₁. The difference in higher-J levels (J_up ≥ 5) varies up to J_up = 18, where the relation between the two fluxes is inverted.

In the case of o-H₂O, the relation between the few fluxes F₀ and F₁ up to the energy level E_up ∼ 200 K varies depending on the ambient conditions. Above that level, the fluxes obtained with the dust emission in the background radiation field are always brighter (for the ambient conditions explored) by factors of a few and up to three orders of magnitude. Since the H₂O lines observed with SPIRE are not spectrally resolved, and knowing (from HIFI spectra) that some of them are blended with other lines, the excitation analysis and abundance estimates of H₂O are not addressed here. Instead, the analysis and more sophisticated models of the H₂O lines in NGC 253 (and other galaxies) were presented in a parallel work based on HIFI velocity-resolved spectra by Liu et al. (2017). In the next sections, we present the excitation analysis of ¹²CO, ¹³CO, and HCN.

6.3. Excitation of the CO Lines

From the SPIRE and PACS spectra, we have ¹²CO transitions from $J=4\to 3$ to $J=19\to 18$, although the $J=17\to 16$ transition was not detected with PACS because it is found in a very noisy spectral range. The lower-J transitions are taken from the values reported for a 43'' beam by Israel et al. (1995) and Wall et al. (1991), and they were corrected for a 40'' beam, assuming the average source size of 167 as found from the ¹²CO $J=6\to 5$ map (Section 2.1).

First we tried to fit the full ¹²CO LSED using two components. The low-J (J ≤ 7) lines can be fit with one component, but the second component can fit either the mid-J (J ≤ 12) or the high-J (J ≥ 13) transitions, but not both simultaneously. So, we need three components to fit the full ¹²CO LSED simultaneously. The model we use is described by

$$\begin{eqnarray}&&{F}_{\mathrm{tot}}(\nu )={{\rm{\Phi }}}_{1}{F}_{1}+{{\rm{\Phi }}}_{2}{F}_{2}+{{\rm{\Phi }}}_{3}{F}_{3}\end{eqnarray} \tag{ 9 }$$

where the Φ_i are the beam area filling factors and the F_i are the estimated fluxes for each component in units of W m⁻². The estimated fluxes are a function of three parameters: the density of the collision partner (usually H₂) n(H₂) (cm⁻³), the kinetic temperature of the gas T_k (K), and the column density per line width N/ΔV (cm⁻² km⁻¹ s) of the molecule being studied. These are the input parameters for the modified RADEX code that uses the background radiation field as described in Section 6.2.

In contrast with previous works in the literature, we prefer to fit all of the ¹²CO LSED simultaneously so we do not have to guess or speculate about up to which transition we should fit first and then subtract the modeled fluxes from the remaining higher transitions. This is also because the latter method considers the effect of the first component on the higher-J lines, but it does not take into account the effect of the second component on the lower- and mid-J lines, which we note is not negligible. We also try the excitation conditions (temperature, volume, and column densities) obtained by Rosenberg et al. (2014). In their models, gas densities of up to ∼3 × 10⁵ cm⁻³ were found for their third component. The beam filling factors they reported are larger than 1, which we find nonphysical for a galaxy with an unresolved source size (see the discussion in Section 6.3.1), so we have to scale our fluxes estimated with RADEX by the appropriate filling factors. We found that the ¹²CO $J=14\to 13$ and $J=15\to 14$ lines observed with PACS are underestimated by factors of ∼2.5 and ∼4, respectively, using the excitation conditions from Rosenberg et al. (2014), while the higher-J lines are underestimated by more than one order of magnitude.

Considering all of the transitions would require 12 parameters to fit the ¹²CO LSED alone, so methods like the Bayesian likelihood analysis used in the literature (e.g., Ward et al. 2003; Kamenetzky et al. 2012) become impractical due to the large number of combinations of input parameters that need to be explored. Instead, we use the simplex method (e.g., Nelder & Mead 1965; Kolda et al. 2003) to minimize the error between the observed and estimated fluxes, using sensible initial values and constraints of the input parameters as described below. Following Rosenberg et al. (2014), we also included all of the available ¹³CO fluxes (from SPIRE and ground-based telescopes; e.g., Israel et al. 1995) to constrain the column density of the lower-J lines (up to $J=8\to 7$), as well as the HCN fluxes (from Paglione et al. 1997; Knudsen et al. 2007), to break the dichotomy between density and temperature for the high-J transitions. The RADEX fluxes of the ¹²CO and ¹³CO lines were corrected by the FWHM (estimated from a Gaussian fit) of ΔV = 190 km s⁻¹ (Wall et al. 1991, consistent with the average FWHM of the ¹²CO and ¹³CO $J=6\to 5$ HIFI lines; cf. Table 6) and by ΔV = 120 km s⁻¹ for the HCN lines (Paglione et al. 1997, their Table 2). All fluxes from ground-based telescopes were corrected to our 40'' beam.

6.3.1. Constraints of the Model Parameters

From the high spatial resolution maps by Sakamoto et al. (2006, 2011) and the two-dimensional Gaussian fit of the continuum and ¹²CO emission (Sections 2 and 3), we know that the size of the ¹²CO-emitting region (∼167) is smaller than the beam size (40''), so the beam area filling factors Φ_i must be strictly lower than unity, irrespective of the number of clouds or clumps found along the line of sight. Also, high-resolution maps (Sakamoto et al. 2011) and SOFIA/GREAT observations of the ¹²CO $J=16\to 15$ toward Galactic molecular clouds (e.g., Pérez-Beaupuits et al. 2015a) indicate that the size of the CO-emitting region decreases with J-transition. Therefore, the beam area filling factors of the three components should also decrease. Hence, the following condition was imposed in the fitting procedure:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{3}\leqslant {{\rm{\Phi }}}_{2}\leqslant {{\rm{\Phi }}}_{1}\lt 1\end{eqnarray} \tag{ 10 }$$

Following Ward et al. (2003) and Kamenetzky et al. (2012), we also restricted the density n(H₂) and column density N(CO) to physically plausible values. That is, the total molecular mass of the emitting region (M_region) cannot be larger than the dynamical mass 2.4 × 10⁹ M_☉ of the galaxy (Houghton et al. 1997), and the column lengths cannot be larger than the size of the emitting region. These restrictions eliminate models with very large column density and too low volume densities. The molecular gas mass contained in the beam is estimated as

$$\begin{eqnarray}&&{M}_{\mathrm{mol}}=\displaystyle \frac{{A}_{\mathrm{beam}}1.5{m}_{{{\rm{H}}}_{2}}{\displaystyle \sum }_{i=1}^{3}{{\rm{\Phi }}}_{i}{N}_{i}}{{X}_{\max }}\end{eqnarray} \tag{ 11 }$$

where A_beam is the area (in cm²) subtended by the beam size, Φ_i and N_i are the beam area filling factors and column densities of the three components, and the factor 1.5 multiplying the molecular hydrogen mass ${m}_{{{\rm{H}}}_{2}}$ accounts for helium and other heavy elements (Kamenetzky et al. 2012). Following Ward et al. (2003), we assumed a conservative value X_max = 5 × 10⁻⁴ for the [¹²CO]/[H₂] fractional abundance, since the average value found in starburst galaxies may be even higher than the values (e.g., 2.7 × 10⁻⁴) measured in warm star-forming molecular clouds like NGC 2024 (Lacy et al. 1994).

The circumnuclear gas layer extends about 680 × 255 pc (∼40$^{\prime\prime} \,\times$ 15'') at position angle 58° as estimated from the CO $J=2\to 1$ map by Sakamoto et al. (2006). A smaller extension, however, is expected for the higher-excitation gas. From the 2D Gaussian fit of the ¹²CO $J=6\to 5$ map, we estimate a CO-emitting gas extension of about 350 × 210 pc (∼208 × 125), assuming a distance D = 3.5 Mpc (Rekola et al. 2005). So, we used the smallest extent of 210 pc across, to constraint the equivalent length of the ¹²CO column density. The latter can be approximated from the area filling factor Φ_i, assuming a circular (Gaussian) homogeneous emitting region of size 210 pc and a circular homogeneous cloud of size S_cloud ≈ N_i/(n(H₂)X_max). In the same way, the beam filling factor can be estimated as (Ω_source/Ω_beam)², assuming a homogeneous source size Ω_source and a Gaussian beam size Ω_beam, and the area filling factor of our models can be estimated as the area of the cloud size over the area of the emitting region. So, the cloud size can be constrained using the smallest extension of 210 pc across as the upper limit using the following expression:

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{i}}{n({{\rm{H}}}_{2}){X}_{\max }}\leqslant \sqrt{{{\rm{\Phi }}}_{i}}(210\,\mathrm{pc}).\end{eqnarray} \tag{ 12 }$$

From the RADEX documentation and several practical tests we performed, we know that the cloud excitation temperature becomes too dependent on optical depth at high column densities. So, very high optical depths can lead to unreliable temperatures due to convergence uncertainties in RADEX. Therefore, we excluded column densities that lead to an optical depth τ ≥ 100 in any of the transition lines. For the volume densities and kinetic temperatures explored, we usually met this condition with ${\mathrm{log}}_{10}({N}_{\mathrm{CO}}/{\rm{\Delta }}V)\gtrsim 18.2$ cm⁻² km⁻¹ s.

Since ¹²CO and ¹³CO are supposed to coexist in the same emitting gas, we used the same volume density and kinetic temperature for ¹³CO as obtained in the three components of ¹²CO. For the column density of ¹³CO, we used the ¹²C/¹³C isotope ratio of ∼40 confirmed by Henkel et al. (2014). From the high-resolution maps by Sakamoto et al. (2011) and observations of Galactic molecular clouds (e.g., Pérez-Beaupuits et al. 2010, 2012), we know that the ¹³CO emission is less extended than that of ¹²CO. Therefore, we restricted the beam area filling factors of ¹³CO components to be lower than those of ¹²CO, and they are the only three free parameters used to fit the ¹³CO LSED. We found that only the first two components of ¹²CO are sufficient to fit the ¹³CO LSED as well as the lower (J_up < 9) transitions of ¹²CO, while the third component contributes significantly for J_up > 10 transitions.

On the other hand, we found that the HCN LSED can be reproduced using the same volume density and temperature of the second and third components of ¹²CO, while the HCN column densities are free parameters. Because the critical densities of the mid- and high-J CO lines are comparable to those of the low-J HCN lines, and from the extension of the HCN $J=4\to 3$ map by Sakamoto et al. (2011), we inferred that the beam area filling factor of the first HCN component should be ${{\rm{\Phi }}}_{1}(\mathrm{HCN})\leqslant {{\rm{\Phi }}}_{2}{(}^{13}\mathrm{CO})$. We set the area filling factor of the second HCN component to be equal to Φ₃(¹²CO), given that the extension and distribution of the high (J_up > 13) ¹²CO transitions are similar to those of the HCN lines, as observed in Galactic star-forming regions (e.g., M17 SW; Pérez-Beaupuits et al. 2015a, 2015b).

6.3.2. The LSED Model

The best model fit for ¹²CO is shown in Figure 13. The insets show the model fit for ¹³CO and HCN. The resulting parameters of each component are summarized in Table 11. The third component fitting the higher-J (PACS) lines is a totally new addition surpassing works previously reported. It shows that the molecular gas in the central 350 × 210 pc of NGC 253 is much more highly excited than that traced with only the lower- and mid-J lines observed with ground-based telescopes or even Herschel/SPIRE alone. The HCN fluxes help to constrain well the parameters of the third component. We found that temperatures ≳180 K do not reproduce the slope described by the J_up > 11 ¹²CO fluxes, overestimating the J_up > 14. Lower volume densities could compensate for the overestimation, but n(H₂) < 10⁵ cm⁻³ does not reproduce the three available HCN fluxes.

Zoom In Zoom Out Reset image size

Table 11.  LVG Model Results for NGC 253

	Component Parameters
Quantity	First component	Second component	Third component
N(H2) [cm−2]	(8.0 ± 1.8) × 1021	(1.6 ± 0.4) × 1022	(3.2 ± 0.7) × 1021
Scloud [pc]	1.6 ± 0.4	(1.5 ± 0.3) × 10−2	(2.6 ± 0.6) × 10−4
Mmol [M☉]	(1.9 ± 0.4) × 107	(7.6 ± 1.7) × 106	(6.1 ± 1.3) × 104
12CO
Φ	(2.8 ± 0.6) × 10−1	(5.5 ± 1.2) × 10−2	(2.2 ± 0.5) × 10−3
TK [K]	90 ± 10	50 ± 6	160 ± 12
n(H2) [cm−3]	(1.6 ± 0.3) × 103	(3.2 ± 0.8) × 105	(3.9 ± 0.8) × 106
N(12CO) [cm−2]	(4.0 ± 1.5) × 1018	(7.9 ± 3.5) × 1018	(1.6 ± 0.4) × 1018
13CO
Φ	(2.5 ± 0.6) × 10−1	(1.2 ± 0.3) × 10−2
TK [K]	90 ± 10	50 ± 6
$n({{\rm{H}}}_{2})$ [cm−3]	(1.6 ± 0.3) × 103	(3.2 ± 0.8) × 105
N(13CO) [cm−2]	(1.0 ± 0.3) × 1017	(2.0 ± 0.8) × 1017
HCN
Φ		(1.2 ± 0.3) × 10−2	(2.2 ± 0.5) × 10−3
TK [K]		50 ± 6	160 ± 12
n(H2) [cm−3]		(3.2 ± 0.8) × 105	(3.9 ± 0.8) × 106
N(HCN) [cm−2]		(1.3 ± 0.5) × 1014	(1.6 ± 0.6) × 1015

Download table as:  ASCIITypeset image

Table 12.  Line Flux (W m⁻²) Ratios from the 40'' Aperture (SPIRE and PACS) and Toward the Southwest (SW) Position (SOFIA/upGREAT and APEX/CHAMP+)

Line	40'' Central	151 SW
Ratio	Region	Position
[C ii]/CO(11–10)	5.31 ± 1.12	1.23 ± 0.22a
[C ii]/CO(6–5)	16.29 ± 3.35	5.61 ± 1.01
CO(11–10)/CO(6–5)	3.01 ± 0.55	4.45 ± 0.78b

Notes.

^aIf a 151 beam is considered instead of a 227 beam to compute the flux of ¹²CO $J=11\to 10$, this ratio would be a factor of ∼2.25 larger. ^bIf a 151 beam is considered for ¹²CO $J=11\to 10$, this ratio would be about 56% smaller.

Download table as:  ASCIITypeset image

From the N(CO) columns, we derive the column density of molecular hydrogen for each component. As discussed above, we assumed a ¹²CO abundance relative to H₂ of 5 × 10⁻⁴, which leads to values of N(H₂) for the first and third components that are similar (within the uncertainties) to the H₂ columns estimated from the dust continuum emission (Section 4). Our assumed [¹²CO]/[H₂] value is a factor of ∼2.3 larger than the relative abundance of 2.2 × 10⁻⁴ derived by Harrison et al. (1999) based on an assumed carbon gas-to-dust ratio and the measured fractions of gaseous carbon-bearing species. On the other hand, our assumed [¹²CO]/[H₂] value is a factor of ∼6 larger than the value of 8 × 10⁻⁵ assumed by Bradford et al. (2003) for a much smaller 15'' beam. If we use the latter [¹²CO]/[H₂] value instead, we would obtain molecular hydrogen column densities that are even larger than the N(H) column density derived from dust continuum emission.

6.3.3. Gas Mass Traced by ¹²CO

The gas mass of the cloud associated with each component of the model can be estimated using Equation (11). Using our assumed relative abundance value of [¹²CO]/[H₂] = 5 × 10⁻⁴ and the adopted local velocity dispersion of 10 km s⁻¹, we find molecular gas masses of about 2 × 10⁷, 8 × 10⁶, and 6 × 10⁴ M_☉, for the first, second, and third components, respectively (cf. Table 11). If we add up the masses of the three components, we obtain a total gas mass of ∼2.7 × 10⁷ M_☉, which is similar (within the uncertainties) to the range of masses (1–5 × 10⁷ M_☉) found by Harrison et al. (1999) and Bradford et al. (2003), based on ground-based observations of low- and mid-J¹²CO transitions. This gas mass, however, is about one order of magnitude lower than the gas mass derived from the 870 μm and 500 μm dust continuum emission of APEX/LABOCA (Weiß et al. 2008) and in our Section 4, as well as the gas mass derived by Houghton et al. (1997) from the ¹²CO $J=1\to 0$ line intensity alone.

The total mass derived from our ¹²CO LSED model is in agreement with the mass found by Bradford et al. (2003) and the LVG models by Rosenberg et al. (2014). As noted by Rosenberg et al., this mass value should be considered a lower limit, since CO becomes dissociated in the presence of high radiation fields, and thus, our assumed [¹²CO]/[H₂] abundance ratio may underestimate the actual column density of H₂ gas. In order to match the gas mass obtained from the dust continuum emission at 500 μm (Section 4), a ¹²CO to H₂ relative abundance of 2.2 × 10⁻⁵ would be needed, which is a factor of ∼3.6 smaller than the value assumed by Bradford et al. (2003).

6.3.4. Cloud Sizes and Their Relation with Star-forming Regions

From Equation (12) we obtain a characteristic cloud size (including only the molecular gas) of about 1.6 ± 0.4 pc for the first ¹²CO component. This is similar to the cloud size of 2 pc found by Bradford et al. (2003, their Section 4.2) based on visual extinction arguments and including the atomic gas. The size of this component is comparable to the size of diffuse clouds or individual dark clouds in the Milky Way (e.g., ζ Ophiuchi; Stahler & Palla 2005).

The characteristic sizes of the second and third components are much smaller than 1 pc (cf. Table 11). The second component has the largest column density, as well as the lowest gas temperature of the three components, and it has a high volume density of ∼3 × 10⁵ cm⁻³. So this component can be associated with starless cores, or dense and relatively cold cores where the star formation process may be at play.

From the third component, instead, we derive an equivalent size of about 3 × 10⁻⁴ pc (∼9.3 × 10⁹ km or ∼62 AU). This is just about two orders of magnitude larger than the size of the supergiant star Rigel in the Orion constellation and is about two orders of magnitude smaller than the ∼0.5 pc size of the small clouds (SCs) found in some SNRs like IC 443, although these SCs are also expected to be clumped and have small filling factors in 45'' and 55'' FWHM beams (e.g., Lee et al. 2012). This estimated size is also much smaller than its estimated Jeans length (∼12 pc) derived from its gas density and temperature (assuming all the gas is molecular), so the objects/clouds this component is associated with are not likely to have a purely gravitational origin.

6.3.5. Energetics and Excitation

The observed CO LSED measures the luminosity of the molecular gas in the central 40'' of NCG 253. The total cumulative flux of all available CO lines (cf. Figure 13) is 1.65 × 10⁻¹⁴ W m⁻². This corresponds to ∼34% of the [C ii] 158 μm intensity and ∼42% of the combined [O i] 63 μm and [O i] 145 μm intensities (Table 4).

At a distance of 3.5 Mpc, the molecular (CO) gas flux corresponds to a luminosity of 6.3 × 10⁶ L_☉, giving a luminosity-to-mass ratio of ∼0.23 L_☉/M_☉ considering the total gas mass contained in our three CO components. This ratio is about a factor of 2 larger than the ratio found by Bradford et al. (2003), who consider a distance of 2.5 Mpc instead.

The total CO flux in the inner 40'' region represents a 3.2 × 10⁻⁴ fraction of the total IR luminosity of the galaxy observed with IRAS (Rice et al. 1988). This fraction is considerably large because of the large fraction of gas mass (represented by the first and second CO components in our model) that is highly excited. The L_CO/L_IR observed in NGC 253 is only about a factor of 2 lower than the luminosity ratio found in the luminous infrared galaxy NGC 6240 (Meijerink et al. 2013), and almost one order of magnitude higher than the ratio found for Mrk 231 (van der Werf et al. 2010) and Arp 220 (Rangwala et al. 2011). Although a large line-to-continuum ratio can be explained by gas compressed and heated by shocks, like in the case of NGC 6240 (Meijerink et al. 2013), we do not think that this is the case for NGC 253. That is because the bulk of the ¹²CO luminosity is contained in the low- and mid-J lines that have either low density (first component) or low temperature (second component), which do not match a shock- or turbulence-dominated scenario. Besides, as mentioned by Bradford et al. (2003), the near-IR H₂ emission observed in the nuclear region of NGC 253 is more characteristic of UV fluorescence rather than the thermal spectrum produced by shock heating, as found in Galactic outflow sources (cf. Engelbracht et al. 1998). This was confirmed by high-resolution VLT/SINFONI maps of the near-IR H₂ and [Fe ii] emissions, where no strong correlation between H₂ and [Fe ii] (a strong near-IR shock tracer) was found, while a good match between H₂ and the ISAAC PAH 3.21 μm (a tracer of the fluorescently excited gas, including excitation by both O and B stars) maps was observed (Rosenberg et al. 2013). We believe the above is enough observational evidence to conclude that shocks or turbulence heating is not likely to be the main excitation source of the bulk of the CO line intensities, described by the first and second components in our LSED model. This may seem in contradiction with (mechanically and CR-heated) PDR models previously reported by Rosenberg et al. (2014), where they concluded that mechanical heating is necessary to reproduce the observed CO emission from ground-based and SPIRE data alone. However, it is sufficient to recall Rosenberg et al. (2014, their Section 6), where they discuss that the relative contribution of mechanical heating is dominant over CR heating in their models, but where they also state that the main source of heating in all of the models they tested is actually photoelectric heating. This is then in agreement with our statement above, that shock (or turbulent/mechanical) heating is not the main source of excitation for the low- and mid-J CO emission.

On the other hand, the third component in our model describing the ¹²CO J_up > 13 lines does have high density and temperature, matching best a shock- or turbulence-dominated scenario. The emission of the high-J¹²CO lines detected with PACS may originate from shocked clumps in SN remnants as well as in the turbulence-dominated clumps along the molecular outflow traced by the ¹²CO $J=1\to 0$ high-resolution ALMA observations (Bolatto et al. 2013). However, this molecular outflow is observed neither in the ¹²CO $J=2\to 1$ (neither its isotopes) nor in the ¹²CO $J=3\to 2$ high-resolution observations done with the Submillimeter Array (SMA) by Sakamoto et al. (2011). Although it can be argued that the emission of the $J=2\to 1$ and $J=3\to 2$ ¹²CO transitions may indeed be present in the molecular outflow identified by Bolatto et al. (2013), and they may just be below the sensitivity level of SMA (in comparison with ALMA), it is a fact that the intensity of the J_up > 13 lines of ¹²CO are actually even fainter than the $J=3\to 2$ line, as observed in the PACS fluxes (cf. Figure 13). That means the fluxes of the ¹²CO J_up > 13 lines observed in our 40'' beam must be dominated by the emission arising from the starburst ring. This scenario is well supported if we consider that the high-J CO lines present a remarkable spatial correlation with the HCN and HCO⁺ lines, as observed with SOFIA/GREAT in some Galactic molecular clouds (e.g., M17 SW; Pérez-Beaupuits et al. 2015b). And the HCN $J=4\to 3$ high-resolution map by Sakamoto et al. (2011) does not trace the molecular outflow observed in ¹²CO $J=1\to 0$. In fact, the actual outflow, originally observed in Hα, may lack sufficiently dense gas to excite any HCN emission (as well as the high-J CO lines), as shown with the HCN $J=1\to 0$ OVRO map by Knudsen et al. (2007, their Figure 5).

Another plausible scenario could be an internal source of heating within the dense gas environment of the starburst ring, i.e., hot cores, which have comparable sizes to the one derived above for the third CO component. The detection of H₂S in NGC 253 by Martín et al. (2005) can be considered to arise from the massive star-forming cores in the nuclear starburst. That is because sulfur is largely depleted (by a factor of 100–1000) in the ISM, and the major gas-phase formation routes to H₂S are mainly endothermic (≥7000 K; Forets et al. 1993; Rodgers & Charnley 2003). Therefore, the H₂S emission is generally associated with sputtering on dust grains due to either intense UV radiation from star-forming regions or shocks generated by young stellar objects, as is being observed in the Orion KL outflow (Minh et al. 1990). Likewise, García-Burillo et al. (2000) reported enhanced abundance of the shock tracer SiO from high-resolution observations, arguing that the SiO emission may arise in bipolar outflows powered by young massive stars associated with the nuclear starburst and/or due to large-scale shocks induced by the nuclear bar. The SiO emission in NGC 253 is located in two regions, between 10'' and 20'' away from the center and opening out in a spiral-like structure, as well as in an inner ring of radius ∼4'' in the center of NGC 253, interpreted as the inner Inner Lindblad Resonance (iILR). The latter SiO emission coincides with the high-resolution H₂S-emitting area reported by Minh et al. (2007). Thus, like SiO, the H₂S emission could originate from shock waves. However, the high rotation temperature derived for H₂S is considered a signature of hot core chemistry, where H₂S is released from dust mantles by the heating of massive star-forming regions (Rodgers & Charnley 2003). Besides, the detection of the H₂ S 2_2,0–2_1,1 transition, which has an upper-state energy level of 84 K, indicates the presence of hot gas. Therefore, Minh et al. (2007) favor the hot core chemistry scenario for the H₂S emission, which is supposed to trace the ongoing star formation through hot core activity. A rough estimate indicates that several thousands of Orion KL–like cores may exist toward the H₂S-emitting area in the inner 20'' × 20'' nuclear region of NGC 253 (Minh et al. 2007). Therefore, the high-J¹²CO lines may as well be associated with these hot cores.

Nevertheless, as discussed by Rosenberg et al. (2014), the effects of CRs (although perhaps not dominant) cannot be ruled out as an external source of heating in hot cores, or in addition to shock generated by YSOs, bipolar outflows, or mechanical heating. The high star formation activity, along with the relatively high SN rate in NGC 253, allowed an enhanced CR density to be estimated, which, in turn, allowed VHE (>100 GeV) gamma-rays from the nuclear region of NGC 253 (e.g., Paglione et al. 1996; Domingo-Santamaría & Torres 2005; Acero et al. 2009; Rephaeli et al. 2010, and references therein) to be predicted (and detected). The gamma-rays are basically the product of the enhanced CR rates interacting with dense gas (e.g., Paglione et al. 1996; Hewitt et al. 2009). The observed gamma-ray flux in NGC 253 indicates a CR density that is three orders of magnitude higher than that found in the center of the Milky Way, and therefore, it is expected to play a significant role in the excitation of the high-J¹²CO lines as well as in the abundances and line intensities of other species.

6.4. Molecular Lines as Diagnostic of Enhanced Cosmic Rays

6.5. Line Ratios as Diagnostic of Ionization, Density, and Temperature of the ISM

We presented a large set of molecular and atomic lines detected in NGC 253, using the three instruments SPIRE, PACS, and HIFI on board the Herschel Space Observatory. About 35 lines were detected and identified in the SPIRE spectra (while few lines still remain unidentified), and 30 lines were identified in the four PACS spectral ranges. A significant number of lines is still unidentified in the PACS spectra, which will be reported in a follow-up work.

Because NGC 253 is a very rich molecular laboratory outside the Milky Way, we were able to detect exotic molecules such CH⁺ J = 1 → 0 in absorption (with SPIRE) and CH⁺ $J=2\to 1$ in emission (with PACS). Other molecules like OH, OH⁺, and H¹⁸O were also detected, both in emission and in absorption.

The APEX/SABOCA high-resolution map of the dust continuum emission at 350 μm allowed us to estimate an average angular size of 167 for the emitting region covered in the 40'' equivalent beam size derived for the SPIRE spectra. The average source size derived for the continuum emission is in agreement with the source size estimated from the 2D Gaussian fit of the APEX/CHAMP⁺ map of the ¹²CO $J=6\to 5$ line.

The velocity-resolved spectra of the [C ii] 158 μm fine-structure line and the high-J CO $J=11\to 10$ transition obtained with SOFIA/upGREAT toward an SW position in the nuclear region were compared with the corresponding spectra of the CO $J=6\to 5$ spectra from APEX/CHAMP⁺. We found that the corresponding [C ii] obtained from the HIFI map matches neither the line shape nor the intensity obtained with SOFIA/upGREAT. We believe that the SOFIA/upGREAT spectrum is correct and that the HIFI spectrum shows excess spectral baseline structure. The CO(11–10)/CO(6–5) line ratio observed toward the SW position is indicative of high densities, in agreement with the position of the brightest HCN $J\,=1\to 0$ emission in the nuclear region of NGC 253.

A thorough data reduction and careful combination of the full set of available SPIRE and PACS spectral and photometric data allowed us to merge both data products in order to study the dust continuum SED of NGC 253, as seen with Herschel. We found a cold dust component with temperature ∼37 K (in agreement with previous results quoted in the literature) and a warm dust component at ∼70 K, about 20 K higher than previously estimated, which is significant when considering a ±10 K error in the temperatures reported. A third component with higher dust temperature of ∼188 K was also identified from the continuum flux observed at shorter (<50 μm) wavelengths. A first-order estimate of the incident FUV fluxes that heat the three dust components yielded values of G_0,c ∼ 3.4 × 10⁵, G_0,w ∼ 8.7 × 10⁶ and G_0,h ∼ 1.2 × 10⁹ (in units of the Habing flux), respectively. Total gas mass ∼4.5 × 10⁸ M_☉ and column density of ∼1.2 × 10²² M_☉ were estimated from the SPIRE continuum flux observed at 500 μm.

Combining the SPIRE and PACS data, we obtained the most extended ¹²CO ladder of the 40'' nuclear region of NGC 253, including the $J=4\to 3$ up to $J=13\to 12$ transitions from the SPIRE-FTS, and the $J=14\to 13$ up to $J=19\to 18$ transitions from the PACS spectral ranges. A non-LTE excitation analysis showed that at least three components are needed in order to fit the ¹²CO LSED. The ¹³CO (from J_up = 5 to J_up = 8) fluxes detected with SPIRE, as well as ground-based observations of the lower-J¹³CO and HCN lines, were used to constrain the parameters of the models. The total molecular gas mass derived from the three ¹²CO components is in agreement with the gas mass derived from previous CO observations found in the literature. A diffuse and rather warm component, with density ∼2 × 10³ cm⁻³ and temperature ∼90 K, was found to fit mostly the lower (J_up < 7) ¹²CO transitions. The second component corresponds to gas with high density of ∼3 × 10⁵ cm⁻³ and a relatively low temperature of ∼50 K. Because of their densities and temperatures, none of these components can be associated with shocks or mechanical heating. The third component, however, which fits mostly the higher-J¹²CO lines detected with PACS, has both high-density (∼4 × 10⁶ cm⁻³) and high-temperature (∼160 K) gas, which makes it a better candidate for shock-/mechanical-heating-driven gas. We also argue that hot cores are another plausible scenario for the excitation of the HCN and PACS ¹²CO lines, given the detection and spatial distribution of H₂S (likely probing hot core chemistry) and HCN, as observed in high spatial resolution maps. However, the effect and role of the enhanced CRs present in the circumnuclear starburst ring (as derived from SN rates and gamma-ray observations) cannot yet be ruled out, nor accounted for, at this stage. The OH lines detected with PACS show fluxes of about one order of magnitude higher than the H₂O lines, which can be a signature of the enhanced CR ionization rates in the nuclear region of NGC 253. Similarly, the detection of the ionic species OH⁺ and H₂O⁺ is also indicative of high ionization fractions due to the enhanced CRs.

SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including University of Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, University of Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, University of Sussex (UK); and Caltech, JPL, NHSC, University of 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).

SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NAS2-97001 and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 and 50 OK 1301 to the University of Stuttgart. We thank G. Sandell, E. Chambers, and the SOFIA operations crew for their outstanding work during the flight and observing campaigns in Palmdale (USA) and Christchurch (NZ) and for delivering high-quality calibrated data.

J.P.P.-B. (ESO/MPIfR) was sponsored by the Alexander von Humboldt Foundation between January 2010 and December 2012, the period during which part of this work was done. Molecular databases that have been helpful include the CDMS (Müller et al. 2005), NASA/JPL (Pickett et al. 1998), LAMDA (Schöier et al. 2005), and NIST (Lovas 2004).

The authors thank the referee for constructive and insightful remarks that helped to improve this work.

Facilities: Herschel (SPIRE - , PACS, HIFI) - , APEX (SABOCA & CHAMP⁺) - , SOFIA (GREAT, upGREAT) - .

Software: HIPE v14,15 (Ott 2010), KOSMA atmospheric calibrator (Guan et al. 2012), APEX calibration software (Muders et al. 2006), BoA (Schuller 2012), GILDAS/CLASS (Pety 2005), RADEX (van der Tak et al. 2007).

We applied the correction for relative gains for bolometers between the product levels 0.5 and 1.0 to obtain a better extended map. In order to extract the integrated fluxes for a given aperture (or beam size), the SPIRE photometry images must be converted first from units of Jy/beam to Jy/pixel using the beam areas (in arcsecs²) corresponding to a 1'' beam as assumed in the pipeline, i.e., this process is independent of the spectral index of the source. The color corrections instead (needed to extract the fluxes) depend on the spectral index of the source. This was estimated as α = −0.7 for NGC 253 from high-resolution radio observations of SNRs (e.g., Ulvestad & Antonucci 1997; Ulvestad 2000; Tingay 2004). In Table 6.15, Chap. 6.9, Recipe for SPIRE Photometry (in the HIPE v.15 Help System), there are tabulated values for color corrections of the three wavelengths (250 μm, 350 μm, and 500 μm) of the SPIRE photometric images for a number of spectral indexes. We interpolated the tabulated color correction values for alpha = −0.5 and alpha = −1.0 to arrive at one that is applicable for NGC 253 at alpha = −0.7. The 40'' aperture integrated fluxes were obtained using the annularSkyAperturePhotometry task in HIPE.

Using a Gaussian shape in the optimization method of the semiExtendedCorrector, and the eccentricity (0.8463) found from the 2D Gaussian fit of the SABOCA map, we found that a semimajor axis (FWHM) of 155 would be the optimal fit to match the spectra in the overlap region (between about 960 GHz and 990 GHz). This is shown in Figure 14 (left) with the original and corrected spectra using the optimized source size, and a zoom in to the overlap region showed in the inset. On the other hand, using the semimajor axis (FWHM) of 173 from the Gaussian fit on the SABOCA map, the corrected spectra show a better match between the OH⁺ line detected in absorption in both bands (right panel in Figure 14).

Zoom In Zoom Out Reset image size

Note that the change in the line fluxes introduced by the semiExtendedCorrector tool is between 4% and 8% weaker than without correction in the SLW band, and 22%–30% brighter than without correction in the SSW band. Different source sizes and shapes (or methods) used to correct and match the SPIRE-FTS bands will have different effects on the line fluxes.

Since the difference between the SLW continuum level and the photometric flux at 250 μm is the largest (∼27 Jy) of the three wavelengths when using a source size of 173 × 92 (while the continuum level at the other two wavelengths are within the errors), we interpret this as an indication that the spectral continuum level at frequencies above 1000 GHz in the SLW band may still be affected by calibration uncertainties. Nevertheless, the spectrum in Figure 14, right panel, is the best calibrated and corrected SPIRE spectrum reported so far for NGC 253, and it will be the spectrum used throughout this work.

The PACS spectral data were reprocessed using the default pipeline for Chopped Large Range Scan SED observations including the telescope background normalization. The leakage regions of the spectra were excluded from the specFlatFieldRange procedure of the pipeline, and a polynomial of order three was used instead of the default order five (used in HIPE versions <v.14) or wavelet (in the latest versions of HIPE) in order to improve the normalization of the population of spectra used to obtain a single spectrum for each spaxel.

For consistency with the 40'' beam-corrected SPIRE spectra (Section 2.1), the PACS spectral ranges were obtained as the total cumulative spectra from all spaxels, corrected by the 3 × 3 point-source correction factor included in the PACS SPG v14.2.0 calibration tree. Using the 5 × 5 point-source correction factor overestimates the continuum flux compared with the corresponding photometric fluxes at 70 μm, 100 μm, and 160 μm.

We compared the continuum level of the PACS spectra with the corresponding 40'' aperture fluxes of the PACS photometry maps (from the HSA; obs. IDs 1342221744 and 1342221745) at 70, 100, and 160 μm. The aperture fluxes were obtained using the annularSkyAperturePhotometry task, with inner and outer radii of 750'' and 800'', respectively. The PACS flux uncertainties include the errors estimated with the annular sky aperture and the 7% absolute point-source flux calibration for scan maps (Balog et al. 2014).

The 5 × 5 corrected PACS spectra, obtained with the modified Chopped Large Range Scan SED pipeline, are shown in Figure 15. The continuum levels of some wavelength ranges do not match the PACS photometry fluxes, especially the wave range between 100 and 150 μm. This mismatch is due mainly to the fact that NGC 253 is a bright source; the default pipeline used in earlier versions of HIPE was not optimized for bright sources. The mismatch is corrected when using the same modified pipeline mentioned above, but using the background normalization procedure. This is now included in the standard SPG of HIPE v.15.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The final 5 × 5 corrected, and background-normalized, PACS spectra used in this work are shown in Figure 2 (right). The section of the spectrum affected by spectral leakage is shown by the gray filled bands, and they were not used in our analysis.

Some of the key lines ([C II], [C I], ¹²CO, ¹³CO) detected in the SPIRE and PACS spectra were also observed with HIFI as single-pointing spectra. These observations can be found in Figure 16. The spectra shown, at their respective beam resolutions, correspond to the averages between the horizontal and vertical polarizations.