MID-J CO EMISSION FROM NGC 891: MICROTURBULENT MOLECULAR SHOCKS IN NORMAL STAR-FORMING GALAXIES

Turbulence and shocks play a profound role in the process of star formation and nuclear gas accretion. Shocks compress the interstellar medium (ISM) and can trigger star formation, while turbulence prevents the collapse of the ISM. These mechanisms can also remove angular momentum from the gas, leading to gas infall toward galactic nuclei, a process that is important in interacting galaxies.

Both mechanisms are very important in converting mechanical energy into thermal energy and can dominate the heating budget of molecular clouds. The main coolants of the shock-heated molecular gas are H₂ and CO lines. With their improved sensitivity, recent space observatories (Infrared Space Observatory (ISO); Spitzer) were able to trace extragalactic rotational H₂ emission and revealed the importance of this emission and the prevalence of shock-excited molecular gas. For example, H₂ emission associated with shock- or microturbulence-excited molecular gas was found in the group Stephan's Quintet (Appleton et al. 2006; Cluver et al. 2010), the radio galaxy 3C 326 (Ogle et al. 2007), nearby edge-on galaxies NGC 4565, NGC 5907 (Laine et al. 2010), and NGC 891 (Stacey et al. 2010, hereafter SCB10), as well as in the extended cold neutral medium in the Milky Way (Falgarone et al. 2005). Ground-based mid-J CO observations are also good probes of shock- or turbulence-excited molecular gas. Bradford et al. (2005) have invoked dissipation of magnetohydrodynamic (MHD) turbulence to explain their mid-J CO observations of the circumnuclear ring in the Galactic center. In the nuclear region of NGC 253, the strong mid-J¹²CO (e.g., Bradford et al. 2003) and ¹³CO(6–5) (Hailey-Dunsheath et al. 2008) emission is consistent with heating of the warm molecular gas by either decay of supersonic turbulence through shocks or by an enhanced cosmic-ray density.

Due to its proximity (9.5 Mpc), and edge-on presentation (i ⩾ 89°; Sancisi & Allen 1979), NGC 891 is the primary target for studies of the effects of energetic processes associated with star formation on the ISM in normal spiral galaxies. The edge-on view of NGC 891 results in greatly enhanced gas columns so that relatively weak lines such as the low-lying quadrupole rotational transitions of H₂ become observable. The mid-infrared (mid-IR) H₂ rotational transitions in NGC 891 were studied for the first time with ISO (Valentijn & van der Werf 1999) and then more recently with Spitzer (SCB10). The latter study combined the H₂ lines observed with Spitzer (S(2), S(1), and S(0)) together with [C ii], [O i], and IR continuum observations from ISO to constrain gas excitation mechanisms. They found that most of the S(0) line emission arises in photodissociation regions (PDRs), and most of the S(2) line emission arises in low-velocity microturbulent shocks, with the S(1) line arising from both types of sources with the ratio 70% PDRs and 30% shocks.

Here, we expand on the analysis of SCB10 by including submillimeter observations of the CO(6–5), CO(7–6), and [C i] 370 μm lines. The mid-J CO lines can arise from gas heated by a variety of mechanisms, including far-ultraviolet (FUV) photons (PDR), X-rays (X-ray dominated region), cosmic rays, and shocks. The [C i] emission arises mainly from PDRs. Since the Spitzer H₂ study showed that the main heating sources of the warm molecular gas in NGC 891 are turbulence and PDRs we will focus on these two heating mechanisms and use our new submillimeter lines to both test and constrain the models of SCB10. As in their analysis, we apply the PDR model of Kaufman et al. (1999) and the shock model of Flower & Pineau Des Forêts (2010).

We observed the ¹²CO J = 6 → 5 (433.56 μm) and J = 7 → 6 (371.65 μm) rotational transitions and the [C i] ³P₂ → ³P₁ (370.41 μm) fine-structure line with our submillimeter grating spectrometer ZEUS (Hailey-Dunsheath 2009) on the Caltech Submillimeter Observatory (CSO) on Mauna Kea, Hawaii, in 2006 December.

ZEUS employs a 1 × 32 thermistor sensed detector array that is optimized for observations in the submillimeter wavelength regime. It provides 32 spectral elements at one spatial position. For the 2006 December observing run, we mounted two bandpass filter appropriate for the 350 and 450 μm telluric window adjacent to each other at the entrance slit of the detector housing, each filter covering 16 spectral pixels. This allowed us to observe in the 350 and 450 μm telluric windows simultaneously and it still provided a full bandwidth coverage of about 5000 km s⁻¹ in each band. The spectral resolution of ZEUS is λ/Δλ ∼ 1000 across the 350 μm and 450 μm telluric windows. The CO(7–6) and [C i] lines, which are separated in velocity by only 1001 km s⁻¹, were observed simultaneously within a single spectrum. We verified the spectral calibration through observation of CO gas absorption lines from a calibration unit mounted in front of the ZEUS dewar.

The CO lines were observed on two different nights. On the first night, we observed only the CO(6–5) emission at the position ${\rm R.A.}={\rm 02^h22^m33\mbox{$.\!\!^{\mathrm s}$}10}$ and decl. = +42°20'551 (J2000). On the second night, we observed the CO(6–5), (7–6), and [C i] lines, but centered ∼2'' further to the north. This difference is, however, within the pointing accuracy of the telescope. All observations were made in chop-nod mode with the secondary chopping at 2 Hz and a chop throw of 30''.

We established a pointing model from observations of Saturn and Uranus at various elevation and azimuth positions. The derived pointing accuracy is ≲ 5''. We estimated the coupling of the planets to the beam by assuming the planets to be disk-like sources with uniform emissivity and assuming the ZEUS beam to be Gaussian. During the observations, the beam size was 10'' in the 350 μm band and 11'' in the 450 μm band. At the distance of NGC 891, 10'' corresponds to 460 pc.

All observations were corrected for atmospheric absorption through use of the CSO radiometers. Each spectrum was corrected for the response function of the grating and detector array, then flux calibrated with respect to Saturn, which was assumed to emit like a blackbody with a temperature of T = 113 K at 370 μm and 434 μm (Hildebrand et al. 1985). The correction for telluric transmission is the largest source of error. Overall, we estimate a line flux uncertainty of ∼25%.

The CO(6–5), CO(7–6), and [C i] 370 μm lines were all strongly detected (Figure 1), and their observed line intensities are listed in Table 1.

Zoom In Zoom Out Reset image size

The warm molecular gas traced by mid-J CO emission can be excited by several mechanisms. For spiral galaxies with quiescent star formation, the most relevant are shocks in microturbulent gas or FUV radiation in PDRs. Spitzer observations of the H₂ mid-IR rotational transitions in NGC 891 have revealed the presence of warm, shock-excited molecular gas. The line intensities can only be explained if part of the H₂ emission arises from this warm, microturbulent gas phase (SCB10). While the H₂ S(0) likely originates entirely from PDRs, 30% of the H₂ S(1) and the majority (80%) of the H₂ (S2) emission likely emerges from gas excited by C-type shocks with a pre-shock gas density of ∼2 × 10⁴ cm⁻³, shock velocities between 20 and 30 km s⁻¹, and a shock filling factor of 2.8 (SCB10). The properties of this shock-excited gas should also be imprinted in the line ratio of our CO observations.

We are taking a slightly different approach in our shock analysis than SCB10. For their analysis, they considered H₂ emission from within large areas of the galaxy since they compared their Spitzer observations with far-IR lines obtained with ISO, which had a beam size of 75''. Here, we focus on the central 11'' of NGC 891 and with our new ZEUS mid-J CO observations we can refine the shock modeling. We combine the Spitzer H₂ observations (Table 2) with our ZEUS observations and apply the shock models of Flower & Pineau Des Forêts (2010). The spatial resolution of the Spitzer H₂ observations is very close to our ZEUS observations. The slit orientation of the long-high Spitzer InfraRed Spectrograph (IRS) H₂ S(0) observation was such that the slit length (223) is along the major axis and the slit width (111) is perpendicular to the major axis. For the short-high Spitzer/IRS observations of the H₂ S(1) and S(2) lines, the slit is rotated, with the slit length (113) oriented perpendicular to the major axis of NGC 891 and the slit width (47) along the plane of the galaxy. Thus, perpendicular to the plane of the galaxy the slits cover about the same area. This width is also nearly identical to the FWHM of the ZEUS beam. However, there is a gap of 30'' between the individual positions of the Spitzer observations along the major axis. Hence, to convert the Spitzer H₂ observations to the ZEUS beam size we have to assume a distribution of the H₂ emission.

In our analysis, we also include the CO(3–2) line, which is available in the literature (Mauersberger et al. 1999; Dumke et al. 2001; Bayet et al. 2006). Unfortunately, the line intensities and flux densities are not consistent. Despite very similar beam sizes, the reported values vary by a factor of about 6.5. For our analysis, we use the intensity given by Mauersberger et al. (1999; Table 2), which are about a factor of three less than the largest reported intensity.

The CO(1–0) (e.g., Scoville et al. 1993; Sofue & Nakai 1993) and CO(3–2) (Mauersberger et al. 1999; Dumke et al. 2001) emission shows a plateau of extended emission with strong central peak and local maxima roughly 1' to the northeast and southwest. The H₂ Spitzer observations also show extended emission with a central peak, a trough at the next sampling positions 30'' away from the center, local maxima at the second sampling position at 60'' from the center on either side, and a decrease in the following sampling positions. Thus, within the limited sampling of the Spitzer observations the distribution of the H₂ and CO emission appears similar. Dumke et al. (2001) estimated a source size of about 17'' × 8'' (deconvolved) for the central gas concentration from their CO(3–2) map. The profile of the CO(1–0) emission along the major axis is similar, but slightly larger, showing an FWHM of about 25'' above the plateau (Figures 2 and 3 in Scoville et al. 1993; Sofue & Nakai 1993, respectively). It is plausible that the higher excitation CO emission is more concentrated toward the center than the CO(1–0) emission. Here, we use the source size estimated from the CO(3–2) emission to convert the H₂ Spitzer observations and the CO(7–6) and [C i] intensities to the 11'' beam size of our ZEUS CO(6–5) observations, resulting in conversion factors of 1.4 for the H₂ S(0) line and 0.9 for the H₂ S(1) and S(2) lines. The converted intensities are given in Table 5.

To model the shock parameters, we use a different method than SCB10. They compared the H₂ emission with the far-IR ISO lines on the scale of the large ISO beam, determined the PDR model that best fits the data (biased toward the ISO far-IR lines), and attributed the remaining H₂ emission to shock-excited, microturbulent gas. Since the ZEUS observations cover a much smaller region we cannot use the ISO observations to constrain the line emission fraction. Scaling the ISO observations to the much smaller beam size of ZEUS is subject to large uncertainties without knowing the morphology of the gas-phase locked in PDRs. While a reasonable scaling factor could be inferred for some far-IR lines, it would be questionable for others. Several observations with the Kuiper Airborne Observatory (KAO) have found that the distribution of the [C ii] and CO(1–0) emission is correlated in galaxies (cf. Crawford et al. 1985; Stacey et al. 1991). It is therefore reasonable to use the observed CO(1–0) morphology in NGC 891 to scale the [C ii] ISO observation to the smaller ZEUS beam. However, no such correlation has been shown for the [O i] lines. In our approach, we therefore started by only using the higher spatial resolution observations, the CO(7–6), (6–5), and CO(3–2) transition together with the H₂ S(0), S(1), and S(2) lines, and applied a constrained least-squares algorithm to simultaneously fit the individual line emission fraction, beam filling factor, and the shock parameters for only a C-shock and a J-shock model. Unfortunately, trying to fit all these parameter only to a shock model results in an underconstrained problem and no unique solution is possible. This is apparent in unreasonably small χ² values for a variety of shock parameters, line fractions, and beam filling factors. In addition, the beam filling factor and the line fractions are degenerate.

To obtain useful model fits, it is necessary to constrain the fraction of the line emission from shocks. This can be achieved by introducing and simultaneously fit a second gas component for the remaining line emission. An obvious choice for this component is PDRs, which have been traced in NGC 891 through the far-IR line emission on large scales. In our fitting routine for the PDR phase, we include our [C i] 370 μm line and the IR continuum.

To scale the far-IR continuum to our 11'' beam, we use the observed morphology of the submillimeter continuum. Serabyn et al. (1999) obtained a map of NGC 891 in the 350 μm continuum using SHARC on the CSO. They estimated that about 4.8% of the continuum flux density from the entire galaxy emerges from within the central 12''. Here, we assume that the total IR luminosity of NGC 891 of L(8–1000 μm) = 1.89 × 10¹⁰ L_☉ (RBGS; Sanders et al. 2003; converted to a distance of 9.5 Mpc) scales the same. We convert the IR luminosity from the central region into an intensity and multiply that value by 1.09 to scale to our 11'' FWHM beam. A large fraction of the IR continuum emission in NGC 891 likely arises from a diffuse dust component and is probably not associated directly with PDRs. Considering the entire galaxy as much as 69% of the IR continuum might arise from this component (Popescu et al. 2004). Given that our observation is focused on the center of NGC 891 the percentage might be lower. Here, we follow SCB10 and assume a value of 50% as the IR-fraction arising from cirrus in NGC 891, and allow an uncertainty of 50%.

We exclude the [C i] 609 μm observation (Table 2) since the intensity appears unreasonably small compared to the [C i] 370 μm line intensity. The ratio of the (converted to 11'') [C i] 370/[C i] 609 intensity is 17.3. In a PDR this ratio increases with density and FUV field, and for a density of n = 10⁶ cm⁻³ and an FUV field of G_o = 10⁶ (in units of the Habing field: 1.6 × 10³ erg s⁻¹ cm⁻²) this line ratio is only 7.8 (Kaufman et al. 1999). For a PDR density of n = 10³ cm⁻³ and an FUV field of G_o = 10², the expected line ratio would be 2.4. The CO(3–2) intensity quoted by Bayet et al. (2006) is also the smallest value compared with the values reported by Mauersberger et al. (1999) and Dumke et al. (2001). It could be possible that our ZEUS [C i] 370 μm intensity is too high. However, we have measured the [C i] 370 μm line simultaneously with the CO(7–6) line. If the [C i] 370 μm intensity is too high would also mean the CO(7–6) is too high. Both the [C i] lines and the cirrus-subtracted IR continuum are expected to originate in PDRs. The [C i] 370/IR and [C i] 609/IR intensity ratios are 2.7 × 10⁻⁴ and 1.6 × 10⁻⁵, respectively. For a PDR density of n = 10³ cm⁻³ and an FUV field of G_o = 10², the expected ratios are 3.2 × 10⁻⁴ and 1.3 × 10⁻⁴ for the [C i] 370/IR and [C i] 609/IR intensity ratios, respectively. Although we cannot rule out errors in estimating the IR continuum arising from PDRs in the central region of NGC 891 or a systematic error in our calibration, it appears that the ZEUS [C i] 370 μm intensity is consistent with other measurement and PDR predictions. A possibility for the smaller-than-expected [C i] 609 μm intensity is self-absorption by cold foreground gas. The [C i] 609 μm transition connects to the ground, has an upper level energy equivalent to 24 K, and atomic carbon is the fourth-most-abundant element. Thus, in an edge-on galaxy self-absorption of fine-structure line transitions between the lowest energy levels of abundant elements is a plausible scenario. For example, [O i] 63 μm line also connects to the ground, and self-absorption by colder foreground gas has also been invoked to explain the relatively weak [O i] 63 μm emission from galaxies (e.g., Vasta et al. 2010 and references therein).

We compare the observed to the predicted line intensities from the combination of a PDR and shock model and apply a least-squares algorithm to fit the individual line fractions and the beam filling factors. The model grid for shocks (Flower & Pineau Des Forêts 2010) is rather coarse, with just two values for the pre-shock densities (2 × 10⁴ cm⁻³ and 2 × 10⁵ cm⁻³) and velocity ranges from 10 to 30 km s⁻¹ for J-shocks and 10–40 km s⁻¹ for C-shocks, with a resolution of 10 km s⁻¹. The model grid for the PDR parameters covers a large range, with densities, log (n/cm⁻³), and FUV fields, log (G_o), between 1 and 5, with a step size of 0.25. In the fitting algorithm, we consider the model intensities from a single PDR surface. We also take into account that the CO emission is usually optically thick (e.g., Bradford et al. 2005) and that 50% of the IR continuum originates from cirrus in NGC 891. For an externally FUV-irradiated molecular cloud, PDRs are created on the outer surface. Thus, along the line of sight a beam intercepts two projected PDR surfaces from a single cloud, one at the front of the cloud and one at the back of the cloud. An optically thin line would be observed from both PDRs, while an optically thick line would only be detected from the front side. To compare the predicted line intensities from a single PDR surface with the observations, we therefore multiply the predicted intensities of the optically thin lines by a factor of two. The beam filling factor derived for PDRs is then with respect to single PDR surfaces and not for entire clouds. We apply this fitting routine to two combinations: a PDR and C-type shock combination and a PDR and a J-type shock combination. Note that the beam filling factors can be larger than unity if the beam intersects many PDR surfaces or shocks along the line of sight, with each contributing to the observed intensities. We do not include the far-IR ISO lines in the fitting algorithm because of the uncertain scaling factors. Instead, we only list the predicted intensities for the [O i], [C ii], and [C i] 609 μm lines from the shock and PDR models. The results are shown in Tables 3–5, where the predicted intensities from shocks and PDRs are already multiplied by the appropriate beam filling factor, and for the PDRs the optically thin lines are multiplied by an additional factor of two.

The fit of combined PDR and shock models to the data is slightly better for J-type than for C-type shocks. However, most of the resulting model parameters are similar and the difference in χ² is not enough to be conclusive. In both cases, the best fit is for a pre-shock gas density of n = 2 × 10⁴ cm⁻³, a PDR density of n_PDR = 3.2 × 10³ cm⁻³, and an FUV field of G_o = 100. For the J-shock/PDR combination, the remaining parameters that best fit the data are a shock velocity of 20 km s⁻¹, a shock beam filling factor of 0.87, and a PDR filling factor of 1.34. For the C-shock/PDR combination, the best shock velocity is 10 km s⁻¹ and the beam filling factors are 0.25 for the C-shock and 1.28 for the PDR component. The fraction of the individual line emissions is also very similar for both types of shocks. In both cases, 100% of the CO(7–6) and 99% of the CO(6–5) intensity would arise from shock-excited gas. Of the CO(3–2) emission, 47% and 32% would arise from shocks in the C-shock/PDR model and J-shock/PDR model, respectively. For the H₂ emission, 94% and 93% of the S(2) intensity and 45% and 44% of the S(1) intensity would arise from shocks in the C-shock/PDR and J-shock/PDR model, respectively. For both model combinations only about 3% of the H₂ S(0) intensity would originate in shocks. The main difference between the results of the model combination is in the beam filling factor of the shocks.

The parameter solution for the shock models is very sensitive to a change in observed intensities, and more observations would be very helpful to provide better constraints. This is mainly due to the coarse grid of the shock model. Marginalization of the probability density distribution shows that the PDR solutions are strongly peaked with a 1σ range that is less than a step size in the model grid for density and FUV parameters. While the marginalized probability density distribution for the shock velocities also strongly peaks at the best-fit values with a 1σ range that is less than a step size, running the fitting routine with the CO(3–2) intensity increased by about 10% changes the densities of the C-shock/PDR model to a pre-shock density of n = 2 × 10⁵ cm⁻³, a shock filling factor of 0.11, and a PDR density of n_PDR = 5.6 × 10³ cm⁻³. In contrast, J-shock/PDR solutions are hardly affected. The fractions of the line emission originating from shocks are quite robust. Given the uncertainty in the estimate of the IR continuum and the observed CO(3–2) intensity, for which we assumed a relative error of 50% in the fitting routine, the predicted intensities are not better than a factor of two in most cases. In fact, a change in shock velocity by 10 km s⁻¹ can especially change the predicted [O i] intensities originating from shocks by an order of magnitude. Nevertheless, except for high-velocity J-shocks the [O i] intensities from shocks would always be smaller compared to the predicted PDR intensities.

We note that the observed [C i] 609 μm intensity is too weak by roughly an order of magnitude compared to the J-shock/PDR model result and by about a factor of about three compared to the C-shock/PDR model result.

Since there are no high spatial resolution observations available for the far-IR lines, we can only compare the model predictions among each other. Including the beam filling factors, the models predict a fraction of <1% of the [O i] 63 μm and 146 μm intensity from C-shocks or J-shocks. For our best-fit shock/PDR model combinations, the entire [O i] emission is predicted to originate from PDRs. It is nevertheless interesting to note that the predicted intensity of the [O i] 63 μm line is similar to the intensity in the H₂ S(0) line in shocks. However, both intensities are small compared to the intensity in the higher H₂ transitions.

We now compare the predicted far-IR line emission that likely emerges from the central 11'' with the ISO observations. According to the CO(1–0) cut along the major axis of NGC 891 (e.g., Sofue & Nakai 1993), about 20% of the line emission within the 75'' (FWHM) central region originates from the central 11'' (FWHM). If the ISO far-IR fine-structure lines scale similarly, as would be a reasonable assumption for the [C ii] line, then the predicted line intensities from our PDR/shock model are only about a factor of two larger than the scaled ISO far-IR lines.

We can also just fit the PDR model grid to the ISO observations of the center of NGC 891. If the IR continuum would scale similar to the [C ii] emission as observed by ISO then about 23% of the IR continuum from the entire galaxy would arise from the inner 75'' (FWHM). As before, we also take into account that probably only 50% of the IR continuum might arise from PDRs (with the rest originating from cirrus) and that the [O i] 63 μm line is optically thick. The best PDR fit to the observed [C ii], and [O i] 63 μm and 146 μm ISO observations and the IR continuum is obtained for a gas density of n_PDR = 1 × 10³ cm⁻³, an FUV field of log (G_o) = 2.2, and a PDR beam filling factor of 0.07. The observed intensities and the intensities predicted by the PDR model are shown in Table 6. Both the estimated density and the FUV field are similar to the predictions we obtained for the combined shock/PDR model within the much smaller ZEUS beam. In contrast, the beam filling factor is much smaller than the one we obtained previously for the ZEUS observation, as would be expected for the large ISO beam, whose diameter is much larger than the thickness of the edge-on galactic disc. This small PDR beam filling factor within the ISO beam would correspond to a filling factor of 3.26 within the ZEUS beam. It thus appears that the PDR conditions are similar within a factor of a few (<10) within small region covered by ZEUS and the large region covered by ISO. If we apply the beam filling factor corresponding to the ZEUS beam (3.26) instead of the ISO beam filling factor, which translates into multiplying the predicted values in Table 6 by 3.26/0.07, then the predicted intensities from the PDR are within factors of two to three for the [C i] 370 μm intensity observed with ZEUS and the CO(3–2) and H₂ S(0) intensities. The PDR model, however, fails to predict the mid-J CO intensities by more than two orders of magnitude and the H₂ S(2) emission by about one order of magnitude. In fact, the large observed CO(7–6)/CO(6–5) intensity ratio can hardly be modeled by a PDR at all. There is no single PDR model that can fit both the CO(7–6)/CO(6–5) and CO(7–6)/[C i] 370 μm intensity ratios (Figure 2), clearly indicating the presence of an additional heating mechanism.

Zoom In Zoom Out Reset image size

We have shown that the warm ISM in the nuclear region of NGC 891 is likely heated by a combination of FUV radiation and shocks, with the mid-J CO transitions probably entirely due to shocks. Since NGC 891 is an isolated galaxy with only a moderate star formation rate (SFR) and no active galactic nucleus, we do not consider shock excitation through powerful mechanism that are encountered in ultraluminous infrared galaxies or galaxy interactions, where strong X-ray emission from an active nucleus, nuclear jets, or galaxy-wide gas flows likely drive internal or external shock fronts that are suggested by their strong H₂ emission. Instead, we have considered moderate internal energy sources as heating mechanisms for the PDRs, shocks, and turbulence.

NGC 891 has a modest SFR of 3.8 M_☉ yr⁻¹ (Popescu et al. 2004), which could be responsible for the PDR excitation of the far-IR and submillimeter fine-structure lines. Sources of shocks and turbulence that could excite the H₂ and CO rotational line emission include cloud–cloud collisions, outflows from young stellar objects (YSOs), or supernovae (SNe) driven outflows. For example, the warm molecular gas that strongly emits mid-J CO lines in the circumnuclear ring in the Milky Way is likely heated by dissipation of MHD turbulence (Bradford et al. 2005). Cloud–cloud collisions in the material that falls into the central region is possibly the source for the turbulence there. However, a 3 pc scale circumnuclear ring, like the one in the Milky Way, would not contribute significantly within our 500 pc scale ZEUS beam on NGC 891.

In the center of NGC 253, the mid-J CO emission is probably either due to the decay of turbulent shocks or due to cosmic-ray heating (Hailey-Dunsheath et al. 2008). Both mechanisms could provide enough energy to heat the warm molecular gas in NGC 253. The enhanced cosmic-ray flux in the center of NGC 253 is likely due to the high supernova rate (SNRa) of 0.1 yr⁻¹ in the central ∼100 pc (Bradford et al. 2003). SN outflows could also drive turbulence.

Roussel et al. (2007) have studied H₂ emission in the SINGS sample, which is comprised mainly of "normal" galaxies and some LINER/Seyfert galaxies. For the galaxies powered by star formation, they found a good correlation between the combined power in the H₂ S(0), S(1), and S(2) lines and the strength of the 7.7 μm polycyclic aromatic hydrocarbon (PAH) emission. In contrast, the H₂/PAH ratio in the LINER/Seyfert galaxies shows a much larger scatter and the ratio is usually elevated for stronger H₂ emission. This suggests that if the combined power in the H₂ S(0), S(1), and S(2) lines is above a certain threshold, in relationship to the strength of the 7.7 μm PAH feature, then the H₂ lines are too bright to be excited exclusively within PDRs. For H₂ sources above this threshold, they suggest alternative sources of heat ranging from the X-rays from a nuclear engine to interactions with nuclear outflows via shocks. SCB10 have compared the combined H₂ rotational lines with the 7.7 μm PAH feature along the plane of NGC 891 and found that the H₂/PAH ratios are near the top (threshold) values for star-forming galaxies in the SINGS sample for regions within ∼5 kpc of the nucleus. These inner regions contain most (>80%) of the observed H₂ line luminosity for the galaxy so that most of the line emission is at the threshold between PDR and shock-excited origins. However, for the outer regions of NGC 891 the H₂/PAH ratio appears to rise by a factor of two, which places these regions more into the shock-excited regime as defined by Roussel et al. (2007). The rise in the H₂/PAH ratio is due to a faster decline in the PAH emission compared with the H₂ emission in the outer galaxy, and might be due to a variety of effects including a lessened importance of star formation in the outer galaxy, or lowered metallicity in the outer regions of NGC 891. This later effect is a plausible cause for the enhanced H₂/PAH ratio seen for the low-metallicity dwarf galaxy NGC 6822 in the Roussel et al. (2007) sample.

For our modeled C-type shocks with a velocity of 10 km s⁻¹ and a pre-shock gas density of 2 × 10⁴ cm⁻³, the total expected CO intensity (up to J = 20) is 3.2 × 10⁻⁴ erg s⁻¹ cm⁻² sr⁻¹ (Table A1 in Flower & Pineau Des Forêts 2010). The total intensity of the H₂ rotational transitions is a factor of two higher, and together with CO they provide basically the entire radiative cooling for the C-shock. Flux energy in the neutral oxygen line is less than 1% and therefore negligible. About 54% of the mechanical flow energy is converted into line radiation. Applying the derived beam filling factor of 0.25 for C-shocks, the total combined luminosity of H₂ and CO within the 11'' beam is then 2.17 × 10⁶ L_☉, which is roughly half the mechanical heating rate.

For our J-shock solution, the total expected CO intensity (up to J = 20) is 1.4 × 10⁻³ erg s⁻¹ cm⁻² sr⁻¹ (Table A2 in Flower & Pineau Des Forêts 2010), a factor of 4.3 higher than for the C-shocks. The strongest CO transitions are around J_up = 16, so the total CO intensity given above is a lower limit. The energy in all H₂ lines is a factor of about eight higher than from all CO lines (Flower & Pineau Des Forêts 2010). This results in a combined luminosity of all the H₂ and CO lines of 9.7 × 10⁷ L_☉ within the 11'' beam and including the shock filling factor of 0.87. Both lines account for about 82% of the cooling, and the combined luminosity is more than 40 times larger than for C-shocks. H₂O would account for an additional 15% cooling, and about 97% of the flow energy of this J-shock is converted into line radiation. As for the C-type shocks cooling via the neutral oxygen lines is negligible for J-shocks with the above parameters.

In NGC 891, the total SNRa is about 0.042 yr⁻¹ (Strickland et al. 2004). Following Loenen et al. (2008) and assuming an energy output of 10⁵¹ erg per SN and a transfer efficiency of 10%, the expected mechanical heating rate is L_SN = 8.34 × 10⁹ × 0.1 × SNRa L_☉ (Equation (4) in Loenen et al. 2008), or L_SN = 3.5 × 10⁷ L_☉. This is a factor of 10 larger than the required heating rate for the C-shocks, but about a factor of 3 less than required to explain the J-shocks.

Another possibility for the heating source are outflows from YSOs. Using Equation (2) of Loenen et al. (2008), the ejected energy rate from YSOs can be estimated as (dE_flow/dt) ≈ 3 × 10²⁹r²_pcv³_flown erg s⁻¹, where v_flow is the outflow velocity in km s⁻¹, n is the density of the surrounding medium in cm⁻³, and r_pc is the size of the bubble (in parsecs) that is affected by the outflow. YSOs only affect the molecular ISM locally, so for a cloud with a radius of 0.1 pc (Loenen et al. 2008) the resulting luminosity is ∼16 L_☉, and thus more than 10⁵ YSOs would be required within an 11'' beam (∼500 pc) to account for the line luminosity in the (C-type) shock-excited gas. Given the SFR of 3.8 M_☉ yr⁻¹ (Popescu et al. 2004) and a lifetime of the YSO phase of 10⁵–10⁶ yr (Loenen et al. 2008) the approximate luminosity ranges from 6.1 to 61 × 10⁶ L_☉, or 2.4 to 24 times the line luminosity from C-type shocks. Lada et al. (2010) suggest a relationship between the star-forming rate and the number of YSO based on observations of Galactic star-forming clouds: N(YSO) = 4 × 10⁶ × SFR M_☉ yr⁻¹. Applying this relationship to the SFR of the entire galaxy NGC 891 we would expect about 1.5 × 10⁷ YSOs. NGC 891 is often thought of as an edge-on analog to the Milky Way. Within the central ≈400 pc of the Milky Way the SFR was estimated to be ∼0.07 M_☉ yr⁻¹ (An et al. 2011), which would correspond to about 3 × 10⁵ YSOs. If the energy transfer efficiency is similar to SN, YSO outflows could provide enough energy to drive C-type or J-type shock excitation.

We have observed the CO(6–5), CO(7–6), and [C i] 370 μm lines from the nuclear region of NGC 891 using our submillimeter grating spectrometer ZEUS on the CSO. To analyze the data, we have include observations of the CO(3–2) (Mauersberger et al. 1999) and H₂ S(0), S(1), S(2) (SCB10), and IR continuum from the literature. We find the emission is best explained by a combination of PDRs and shocks. While J-shocks provide a slightly better fit than C-shocks the type of shocks is not conclusive.

The CO(7–6)/CO(6–5) intensity ratios are large for a PDR paradigm and would require extremely high FUV and density values. In addition, the CO ratio and the CO(7–6)/[C i] 370 μm ratio cannot be fitted by a single PDR model. However, the mid-J CO ratios are consistent with shock-excited warm molecular gas, and our ZEUS observations together with Spitzer H₂ observations can be explained by a combined shock/PDR model. This result confirms recent studies of the H₂ S(0), S(1), and S(2) transitions with Spitzer (SCB10), which required C-shocks to explain the higher H₂ rotational transitions. The best-fitting shock models require a pre-shock gas density of 2 × 10⁴ cm⁻³ and shock velocities of 10 km s⁻¹ and 20 km s⁻¹ and shock filling factors of 0.25 and 0.87 for C-shocks and J-shocks, respectively. The most likely energy sources for the shocks in the center of NGC 891 are SNe and outflows from YSOs. While both types of shocks fit our data about equally well, from energy considerations it is significantly easier to provide the requisite mechanical energy for C-shock excitation of the lines than it is for J-shock excitation.

We have estimated the fraction of line emission originating from shocks and PDRs and find that for both types of shocks all of the CO(7–6) intensity and 99% of the CO(6–5) intensity arise from shock-excited gas. For the CO(3–2) emission, about 47% and 32% would originate from C-shocks or J-shocks, respectively. Our model fits also suggest a larger fraction of the H₂ S(2) (94%) and S(1) (45%) emission from shocks than estimated by SCB10. Only about 3% of the H₂ S(0) emission would emerge from shock-excited gas.

The [C i] 370 μm emission arises from PDRs. When fitting the combined PDR/shock model to the ZEUS and Spitzer data we find a best solution for a PDR density of 3.2 × 10³ cm⁻³ and an FUV field of G_o = 100 for both types of shocks. The PDR filling factors are 1.28 and 1.34 for a combination with C-shocks and J-shocks, respectively. Our best PDR solution is consistent with [O i] 63 μm and 146 μm and [C ii] observations on the scale of the larger ISO beam. Comparing the line predictions from our best-fitting shock and PDR models shows that the [O i] emission arises almost entirely from PDRs.

We thank the CSO staff for their excellent support of the ZEUS observations. We also thank the anonymous referee whose comments led to a great improvement in the data analysis. This work was supported by NSF grants AST-0096881, AST-0352855, AST-0705256, and AST-0722220, and by NASA grants NGT5-50470 and NNG05GK70H.