{'fraction_non_alphanumeric': 0.04257680922731449, 'fraction_numerical': 0.03922119016025314, 'mean_word_length': 4.006835313657851, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Herschel-PACS Observations of Far-IR CO Line Emission in NGC 1068: Highly Excited Molecular Gas in the Circumnuclear Disk

We report the detection of far-IR CO rotational emission from the prototypical Seyfert 2 galaxy NGC 1068. Using Herschel-PACS, we have detected 11 transitions in the J_upper=14-30 (E_upper/k_B = 580-2565 K) range, all of which are consistent with arising from within the central 10"(700 pc). The detected transitions are modeled as arising from 2 different components: a moderate excitation (ME) component close to the galaxy systemic velocity, and a high excitation (HE) component that is blueshifted by ~80 km s^{-1}. We employ a large velocity gradient (LVG) model and derive n_H2~10^{5.6} cm^{-3}, T_kin~170 K, and M_H2~10^{6.7} M_sun for the ME component, and n_H2~10^{6.4} cm^{-3}, T_kin~570 K, and M_H2~10^{5.6} M_sun for the HE component, although for both components the uncertainties in the density and mass are plus/minus (0.6-0.9) dex. We compare the CO line profiles with those of other molecular tracers observed at higher spatial and spectral resolution, and find that the ME transitions are consistent with these lines arising in the ~200 pc diameter ring of material traced by H_2 1-0 S(1) observations. The blueshift of the HE lines may also be consistent with the bluest regions of this H_2 ring, but a better kinematic match is found with a clump of infalling gas ~40 pc north of the AGN. We discuss the prospects of placing the HE component near the AGN, and conclude that while the moderate thermal pressure precludes an association with the ~1 pc radius H_2O maser disk, the HE component could potentially be located only a few parsecs more distant from the AGN, and might then provide the N_H~10^{25} cm^{-2} column obscuring the nuclear hard X-rays. Finally, we also report sensitive upper limits extending up to J_upper=50, which place constraints on a previous model prediction for the CO emission from the X-ray obscuring torus. [Abridged]


INTRODUCTION
The excited molecular gas in the centers of Seyfert galaxies offers a sensitive probe of the nature of active galactic nucleus (AGN) feedback on the surrounding ISM. Observational studies of this material in Seyfert nuclei have typically used the H 2 rotational (Lutz et al. 2000;Rigopoulou et al. 2002;Roussel et al. 2007) and the well-studied ro-vibrational (e.g., Thompson et al. 1978;Mouri 1994;Maloney 1997;Davies et al. 2005;Rodríguez-Ardila et al. 2005) transitions. The pure H 2 rotational lines (E upper /k B ∼ > 500 K) are easily thermalized at moderate (n H2 ∼ > 10 3 cm −3 ) densities, while the ro-vibrational lines (E upper /k B ∼ > 7000 K) may be excited through collisions in dense (n H2 ∼ > 10 5 cm −3 ) gas or through UV fluorescence (Sternberg & Dalgarno 1989). Observations of these tracers in Seyferts have identified a number of potentially important excitation mechanisms, including X-rays from the AGN (Maloney 1997;Rodríguez-Ardila et al. 2005), shocks associated with supernova remnants, radio jets, and gravitational instabilities (Roussel et al. 2007), and stellar far-UV (FUV) radiation (Davies et al. 2005), with no clear consensus on a single dominant excitation source. The far-IR (FIR) CO rotational transitions (CO[J upper → J upper − 1], with J upper ≈ 13 − 50) arise from states 500 − 7, 000 K above ground and have critical densities of ∼ 10 6 − 10 8 cm −3 , and complement the H 2 transitions for studies of warm and dense material. Compared with H 2 , the FIR CO lines trace similar energy levels, but have higher critical densities, and are less sensitive to extinction. Additionally, the smaller energy gaps between levels leads to a finer sampling of density-temperature phase space. These lines have been proposed as potential tools for studying the obscuring medium of type 2 AGN (Krolik & Lepp 1989), determining the energy budgets of composite starburst/AGN systems (Meijerink et al. 2007), and identifying accreting black holes in the early universe Schleicher et al. 2010), but previous facilities were unable to detect this line emission from extragalactic sources. Here we take advantage of the superb sensitiv-ity of Herschel-PACS to conduct the first extragalactic study of FIR CO emission, from the prototypical Seyfert 2 galaxy NGC 1068.
NGC 1068 is one of the brightest and best studied Seyfert 2 galaxies. The paradigm of an optically and geometrically thick molecular torus accounting for the Seyfert type 1 and 2 dichotomy followed the detection of scattered broad line emission from this source (Miller & Antonucci 1983), and NGC 1068 has been at the center of subsequent studies of the ISM in Seyfert nuclei. The molecular gas in the central 1 ′ of NGC 1068 has been well studied, and here we review some of the key results. Interferometric observations of CO(1-0) have identified a pair of ≈15 ′′ (≈1.1 kpc) radius spiral arms (Planesas et al. 1991;Helfer & Blitz 1995;Schinnerer et al. 2000), which may be modeled as forming in response to a ≈17 kpc bar (Schinnerer et al. 2000). These arms are bright in Brγ (Davies et al. 1998), PAH emission (Le Floc'h et al. 2001, and submillimeter continuum (Papadopoulos & Seaquist 1999a), and contain most of the star formation in the central region.
Centered on the AGN is the ∼5 ′′ (∼350 pc) circumnuclear disk (CND), which is visible in CO and H 2 1-0 S(1), but becomes particularly prominent in images of HCN  and other high density tracers .
Strong emission in CO(4-3) and HCN(1-0) indicate the gas in the CND is both warm and dense Sternberg et al. 1994;Israel 2009) (although see Krips et al. (2011) for a lower density model).
The high abundances of HCN, CN, H 3 O + , and other molecules in the CND suggest an X-ray−driven chemistry (Usero et al. 2004;García-Burillo et al. 2010;Aalto et al. 2011), and X-ray heating has also been invoked to explain the strong H 2 1-0 S(1) and [Fe II] emission (Rotaciuc et al. 1991;Maloney 1997;Galliano & Alloin 2002). At ∼0.3 ′′ resolution the H 2 1-0 S(1) observations resolve the CND into a ∼1 ′′ ring-like structure (Galliano & Alloin 2002), while the line spectral profiles show evidence for rotation, expansion, and more complex kinematics (Galliano & Alloin 2002;Galliano et al. 2003;Davies et al. 2008). Shocks following this non-circular motion, and possibly associated with jet-ISM interactions, may also be important in heating the molecular CND (Krips et al. 2011). At ∼0.1 ′′ resolution the H 2 1-0 S(1) images reveal two clumps of infalling molecular material at ∼0.1−0.4 ′′ scales that likely play an important role in both fueling and obscuring the AGN, with an estimated infall rate of ∼15 M ⊙ yr −1 to within a few parsecs of the nucleus (Müller Sánchez et al. 2009). Finally, milli-arcsec resolution radio observations identify a series of H 2 O maser spots that trace out the inner surface of a ∼0.5 pc radius molecular disk, centered on the AGN (Gallimore et al. 2004, and references therein).
Hard X-ray observations of NGC 1068 indicate large obscuration to the nucleus (Iwasawa et al. 1997;Matt et al. 1997;Colbert et al. 2002) by an intervening medium with a column density possibly exceeding N H > 10 25 cm −2 . Interferometric mid-IR observations of NGC 1068 identify a parsec scale structure of hot dust that, along with the H 2 O maser disk, may represent the dusty molecular torus responsible for the X-ray obscuration (Jaffe et al. 2004; Raban et al. 2009). However, other investigators have found evidence that at least some of the nuclear obscuration occurs on few to ten parsec scales from the AGN (Cameron et al. 1993;Hönig et al. 2006;Müller Sánchez et al. 2009).
The organization of this paper is as follows. In section 2 we describe the Herschel-PACS observations of the FIR CO lines in NGC 1068. In section 3 we analyze the gas excitation and estimate physical parameters, and in section 4 we compare the physical parameters and line profiles with those of other molecular gas tracers. In sections 5 and 6 we discuss potential heating mechanisms, in section 7 we discuss non-detections, and in section 8 we summarize our findings. Throughout this paper we adopt a distance to NGC 1068 of 14.4 Mpc (Bland-Hawthorn et al. 1997), and a systemic velocity V LSR = 1125 km s −1 .

OBSERVATIONS
The observations were made with the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) on board the Herschel Space Observatory (Pilbratt et al. 2010), as part of the SHINING guaranteed time key program. Ten high resolution range scans were concatenated to cover the 52 − 98 µm and 104 − 196 µm ranges, and we additionally obtained deeper integrations of , , and . These observations amounted to a total of 11.5 hours of integration time. The data reduction was done using the standard PACS reduction and calibration pipeline (ipipe) included in HIPE 5.0 975. However, for the final calibration we normalized the spectra to the telescope flux and recalibrated it with a reference telescope spectrum obtained from dedicated Neptune observations.
In Figure 1 we show the spectra from the central spatial pixel (spaxel) centered on 11 of the 12 CO transitions falling in the 104 − 196 µm range. The CO(25-24) line at λ rest = 104.44 µm lies in a noisy region at the edge of this range, and is not included. The spatial distribution over the PACS array of each detected CO line is consistent with that of an unresolved source. All fluxes and upper limits presented here were therefore extracted from the central spaxel, and referenced to a point source by dividing by the recommended beam size correction factors. The CO line fluxes were measured by fitting a Gaussian profile to the baseline-subtracted spectra. In most cases the three free parameters of the Gaussian were allowed to vary, while baseline uncertainties necessitated fixing the widths of the CO(15-14) and  lines to typical values derived from other line fits (corresponding to a deconvolved FWHM of 250 km s −1 ; see Figure 2). The CO(16-15) line is blended with the 163 µm OH doublet, and CO(17-16) with a pair of flanking OH + lines. In both cases we estimate the CO flux by simultaneously fitting all features. CO(23-22) is blended with a strong H 2 O line with a rest wavelength 209 km s −1 to the red. If the combined feature were attributed solely to H 2 O it would be both broader and more blueshifted than any of the other 6 H 2 O lines detected in the PACS scans, and we interpret this as evidence for significant contamination by . A comparison with the average of the CO(22-21) and CO(24-23) profiles suggests the CO(23-22) line is weaker and/or more redshifted than these lines (Figure 2). However, due to the uncertain- a Total uncertainties combine a 30% calibration error with statistical errors in line fits. Upper limits are 3σ, and refer to the flux density integrated over a 600 km s −1 bin. Some lines are blended with a strong feature, and it is not possible to obtain a flux measurement or useful upper limit.  was not covered in the PACS scans. b Relative to V LSR = 1125 km s −1 . Total uncertainties combine a spectral calibration accuracy of 10% of the spectral resolution (20 − 30 km s −1 ) with statistical errors in line fits.
c Fixed to an intrinsic FWHM of 250 km s −1 (see text).
ties involved in deconvolving the CO and H 2 O lines, we simply exclude CO(23-22) from our analysis. All other transitions from CO(14-13) through CO(24-23) are well detected (Table 1). The 52 − 98 µm range includes the CO(27-26) through CO(50-49) transitions, and none of these lines are detected. We estimate upper limits by first binning the data to 600 km s −1 bins, and then estimating the 3σ noise levels (Table 1). We also detect no emission from 13 CO. The 13 CO(14-13) transition at λ rest = 194.55 µm lies at the noisy edge of a scan, while the J upper ∼ > 21 transitions at λ ≥ 104 µm are blended with 12 CO lines. For 13 CO(15-14) through 13 CO(20-19) we estimate 3σ upper limits of (2 − 4) × 10 −17 W m −2 , also for a 600 km s −1 bin size.  beams that contain a mixture of the CND and the more extended emission (Israel 2009). The second set of Jupper = 1 − 3 points are interferometric measurements integrated over the central 4 ′′ from Krips et al. (2011), and the fainter CO(1-0) flux is the same from Schinnerer et al. (2000). The shaded region shows the range of good-fitting LVG models (χ 2 ≤ 3/2; section 3.2.3), and the solid black curve is the single best fitting model, with the red and blue curves showing the individual contributions from the ME and HE components. Middle: Central velocities of the FIR CO lines, with average values of the ME and HE line centers indicated with horizontal red and blue lines. Bottom: Measured FWHM of the FIR CO lines, with tracks indicating the expected line widths for sources with intrinsic widths of 0, 100, 200, and 400 km s −1 .

Evidence for Two Components
In the top panel of Figure 2 we show the line fluxes and upper limits measured here, along with lower-J measurements obtained from the literature. The middle and bottom panels show the central velocities and widths obtained from the Gaussian fitting. The inflection point seen in the FIR CO line SED at J upper ≈ 19 suggests the presence of multiple components, as does the shift in central velocities between the lowest-and highest-J transitions. For simplicity we assume the FIR CO lines are produced by 2 discrete components: a moderate excitation (ME) component near the systematic velocity, and a blueshifted high excitation (HE) component. Our excitation analysis described below indicates that the J upper ≤ 17 and J upper ≥ 20 transitions are dominated by the ME and HE components, respectively. Separately averaging the central velocities of these two sets of lines gives V ME = 17 ± 12 km s −1 and V HE = −54 ± 21 km s −1 (Figure 2), with a difference of V ME − V HE = 71 ± 25 km s −1 .

LVG Modeling
To quantitatively analyze the FIR CO line SED we employ a large velocity gradient (LVG) model. We use the LVG calculation described in Hailey-Dunsheath et al. (2008), with updated CO-H 2 collisional coefficients from Yang et al. (2010), and a thermalized H 2 ortho/para ratio. In this model the shape of the CO line SED is determined by the gas density (n H2 ), kinetic tem-perature (T kin ), and CO abundance per velocity gradient ([CO/H 2 ]/(dv/dr)). The source is assumed to consist of a number of unresolved clouds, and the absolute line luminosities scale with the total CO mass (M CO ). In the following analysis we use a CO abundance of [CO/H 2 ] = 10 −4 to reparameterize [CO/H 2 ]/(dv/dr) as dv/dr, and M CO as the H 2 mass (M H 2 ). This results in an 8 parameter model, with 4 parameters for each of the ME and HE components.

Background Radiation
The CO excitation may be affected by the background radiation, and we must therefore estimate the local radiation density. At millimeter wavelengths the background is dominated by the cosmic microwave background (CMB) (Kamenetzky et al. 2011), but the FIR background arises from within the galaxy. The PACS integral field spectra provide a continuum map of the central 47 ′′ × 47 ′′ with 9.4 ′′ pixels, and demonstrate that the continuum emission in the central pixel is dominated by sources within the central ≈10 ′′ . The flux density in the central pixel may be modeled as an optically thin modified blackbody with β = 1.5, T dust = 48 K, and normalized to F ν (λ = 100 µm) = 49 Jy, and we adopt this as an estimate of the continuum brightness of the ∼5 ′′ CND. If the flux in the central spaxel is indeed due solely to the CND than this is a moderate (by a factor of 2) underestimate, while emission from within the central ≈10 ′′ but outside of the CND may also be contributing. The measured continuum is in reasonable agreement with previous calculations and observations. For comparison, this measured F ν is a factor of 1 → 4 times larger in the λ = 50 → 200 µm range than obtained from the radiative transfer modeling in Spinoglio et al. (2005). Additionally, extrapolating our modeled SED to λ = 450 µm yields F ν (λ = 450 µm) = 1.2 Jy, comparable to the peak value of ∼1.5 Jy beam −1 measured in a ∼9 ′′ beam by Papadopoulos & Seaquist (1999a). The 60 µm/100 µm ratio in this model is 1.27, and the FIR flux 10 is 2.7×10 −12 W m −2 , giving L FIR = 4πd 2 F FIR = 1.7×10 10 L ⊙ .
We include the effects of background radiation on the equations of statistical equilibrium following Poelman & Spaans (2005). The important parameter in this approach is the mean specific intensity of the external radiation field at the cloud surface, which we define as J ν,ext . We estimate J ν,ext using a simple geometrical model in which the gas clouds are uniformly distributed in a sphere with an observed angular size Ω, and are evenly mixed with the FIR-emitting dust grains. For optically thin continuum, the mean value of J ν,ext is then related to the observed continuum flux density F ν,obs as J ν,ext = I ν,CB + 9 16 where I ν,CB is the sum of the CMB and cosmic IR background (CIB), and we take Ω to correspond to a circular diameter of 4 ′′ . We have run calculations both including and ignoring the local contributions to the background, and see negligible difference in the results. In part this is due to the fact that the background radiation temperatures are only T rad = 13 − 21 K at the wavelengths of the detected FIR lines, while the typical excitation temperatures achieved through collisions alone for the best-fitting models are T ex ≈ 100 K and T ex ≈ 400 K for the ME and HE transitions, respectively. In addition, most of the lines are optically thick for the best-fitting models, and hence the CO is insulated from the external radiation field.

Parameter Limits
We explore two component fits to the FIR CO emission over a large volume of 8-dimensional parameter space, applying physical limits to the model parameters. The most important prior restrictions are placed on the velocity gradient. For self-gravitating clouds in virial equilibrium, we can approximate (dv/dr) vir ≈ 10 km s −1 pc −1 (n H2 /10 5 cm −3 ) 1/2 (Goldsmith 2001). The actual velocity gradient may be larger due to additional sources of gravitational potential, a high pressure inter-cloud medium, or non-virialized motion (Bryant & Scoville 1996), but smaller values are unlikely. Defining K vir as the ratio between dv/dr and (dv/dr) vir (Papadopoulos et al. 2007): we restrict parameter space to K vir ≥ 1. The largest measured velocity gradient in NGC 1068 is in the H 2 O maser disk associated with the AGN. The line of sight velocities of the maser spots shift by ∼600 km s −1 over a ∼2 pc linear range (Gallimore et al. 2001), corresponding to an effective dv/dr ∼ 300 km s −1 pc −1 . To accommodate the maser disk and other high dispersion regions in our models, we extend our calculations up to dv/dr = 1000 km s −1 pc −1 . We note that restricting K vir ≥ 1 and dv/dr ≤ 1000 km s −1 pc −1 combine to limit the density to n H2 ≤ 10 9 cm −3 , but as we discuss below, such high densities are ruled out by other considerations. We calculate the total gas mass in our models as M gas = 1.36 × M H 2 , including the contribution from helium. Schinnerer et al. (2000) estimate a dynamical mass of M dyn = 9 × 10 8 M ⊙ for the CND, and we discard any of our models in which the total M gas of the two components exceeds M dyn . The range of parameter space allowed by the CO data can be further reduced by considering the H 2 pure rotational lines, which arise from states with similar upper energy levels as the FIR CO transitions. We calculate the H 2 rotational spectrum for each model, under the simplifying assumption that the lines are optically thin and thermalized, and the H 2 ortho/para ratio is thermalized. We then rule out any model that overpredicts the flux in any of the lines measured in the large (14 ′′ − 27 ′′ ) ISO-SWS apertures (Lutz et al. 2000). These prior constraints on the LVG model parameters are summarized in Table 2.

General Features of Good-Fitting Models
We proceed by generating model SEDs over a regular 8-dimensional grid, and for each model calculating χ 2 in the normal manner. With 10 data points and 8 free parameters, our modeling has 2 degrees of freedom (dof), and here we discuss the general properties of the set of solutions satisfying χ 2 /dof ≤ 3/2. In Figure 2 we show the range of SEDs covered by this set of good solutions, and for the single best fit model we show the decomposition into the ME and HE components. The  and  lines typically receive comparable contributions from the ME and HE components, while the lower-and higher-J transitions are dominated by the ME and HE components, respectively. The shape of the ME SED is relatively well constrained, and peaks in the J upper = 13 − 16 range. The single dish measurements of CO(1-0), CO(2-1), CO(3-2), and CO(4-3) we show in Figure 2 were obtained with 11 ′′ − 21 ′′ beams, in some cases comparable to the Herschel-PACS resolution (Israel 2009). However, these low-J lines receive strong contributions from the lower excitation gas in the ≈15 ′′ radius starburst ring, and do not constrain our models. The interferometric measurements of the CO(1-0) flux in the CND range from 20 − 120 Jy km s −1 (Schinnerer et al. 2000;Krips et al. 2011), larger than the median ME model flux of 8 Jy km s −1 , suggesting the ME component contributes no more than a minor fraction of the observed low-J emission. The HE component is much less well constrained, with the best fitting SEDs that peak at J upper ≈ 23 only providing a marginally lower χ 2 than those that turn over at J upper ∼ > 30. For these latter models, our upper limits to CO(30-29) and CO(34-33) become useful constraints.
The FIR CO emission is an important coolant of the nuclear molecular ISM, but does not dominate. The total luminosity emitted in the 10 transitions detected here is L CO,FIR = 3.1 × 10 6 L ⊙ , and summing the modeled ME and HE emission over all transitions yields L CO,ME + L CO,HE = (4.3 − 10.5) × 10 6 L ⊙ . For the highest excitation models some additional cooling may arise from the J upper > 40 transitions not included in our LVG calculation, and significant emission is also expected in the J upper ≤ 13 submillimeter transitions. The total emission in the H 2 0-0 S(1), S(3), S(4), S(5), and S(7) rotational lines detected by ISO-SWS is L H2 = 1.7 × 10 7 L ⊙ (Lutz et al. 2000). Treating the upper limits to S(0), S(2), and S(9) as detections increases this by a factor of 1.5, while at the same time some fraction of the lowest-J emission measured with the largest apertures may arise from the starburst ring. Our PACS scans have also detected a number of OH and H 2 O transitions. The bulk of the emission in these molecules detected in the central spaxel likely arises from the unresolved CND, and with this assumption we estimate nuclear luminosities of L OH,FIR ≈ 1.5 × 10 7 L ⊙ and L H2O,FIR ≈ 3.0 × 10 6 L ⊙ . The FIR range includes the strongest OH lines at 79 µm, 119 µm, and 163 µm, while for H 2 O (as with CO) the longer wavelength emission should be strong. In the FIR the CO and H 2 O cooling is therefore comparable, while the FIR CO luminosity is weaker by a moderate factor than the OH and H 2 rotational emission.

Bayesian Analysis
We follow the Bayesian approach outlined by Ward et al. (2003) to quantify the probable values of the model parameters. We consider an 8-dimensional array of bins centered on the grid points for which we have generated a model SED, and calculated a χ 2 value. The probability that the actual parameter set falls within a given bin is proportional to the product of bin size, likelihood L ∝ exp(−χ 2 ), and an assumed prior probability. We choose priors that are flat in the logarithm of each parameter, and that go to zero for any model that violates one of the restrictions listed in Table 2.
In Figure 3 we show the joint density-temperature probability density functions for both components, with contours at 68%, 95%, and 99% of the enclosed probability. The close similarity between the 95% and 99% (and in some regions also the 68%) contours is due to a truncation of the density function following the violation of one of our model restrictions. Starting with the ME component, the behaviour of the contours may be understood as follows. For either a fixed K vir or dv/dr, the shape of the CO SED is approximately conserved if an increase in density is matched by an approriate decrease in temperature. The pair of blue curves in Figure 3 show the density-temperature relation best fitting the data for K vir = 1 and dv/dr = 1000 km s −1 pc −1 . As the shape of the ME SED is well constrained, the region of acceptable parameter space for either set of models corresponds to a narrow band centered on the best fit curve. These two bands intersect at n H2 = 10 9 cm −3 and diverge at lower densities, thereby bracketing a region of high probability filled by models with intermediate velocity gradients. For lower temperatures a decreasing fraction of the CO is excited to the higher-J states, and a larger total mass is needed to reproduce the absolute line fluxes. For n H2 ∼ > 10 8 cm −3 the gas mass exceeds the dynamical mass, and consequently a small region of parameter space is excluded. For higher temperatures the model emission does not fall off with increasing J in the J upper ≈ 18 region as rapidly as the data require, and the lower quality fits limit the 68% confidence interval to T kin 230 K. For T kin ∼ > 400 and n H 2 10 5 cm −3 the models begin to overproduce the detected S(1) and upper limit to S(2) H 2 rotational lines, and this accounts for the sharp cutoff in the probability density in the upper left region.
The shape of the HE SED is not as well constrained as the ME SED, and the H 2 rotational lines provide a stronger constraint to the extent of the probability density contours in Figure 3. The best fit HE SEDs peak at J upper ≈ 23, and while lower excitation solutions are much less probable, higher pressure models peaking at J upper ∼ > 30 also provide reasonable fits. As a result, the dv/dr ≤ 1000 km s −1 pc −1 restriction no longer forms the upper right boundary of the density function, and models are allowed to extend into this higher density and temperature region until the modeled H 2 emission exceeds the observations. The measured S(0) upper limit provides the most important constraint at high densities -Top: Joint density-temperature probability density function for the ME component. Contours are drawn at 68%, 95%, and 99% enclosed probability. The mean of log(M H 2 /M ⊙ ) needed to reproduce the absolute line fluxes is shown by the red curves, and the logarithm of the thermal pressure (log(n H 2 × T kin ) in K cm −3 ) is shown by the dotted black lines. Blue curves show the density-temperature relation giving the best solution for K vir = 1 and dv/dr = 10 3 km s −1 pc −1 (see text). Bottom: Same as the top, but for the HE component.
(n H 2 ≈ 10 8.5 cm −3 ) and low temperatures (T kin ≈ 200 K). The higher excitation H 2 lines become more important at higher temperatures, and the S(7) transition limits the upper left region at T kin ≈ 1000 K.
In Figure 4 we show the fully marginalized distribution functions for each of the 4 primary model parameters, as well as for K vir and the thermal pressure P/k B = n H 2 × T kin . The distribution functions for the HE T kin , n H 2 , and P/k B are shifted to higher values than for the corresponding ME parameters, although the density distributions for the two components contain significant overlap. The lower limit imposed on K vir translates into a lower limit to dv/dr that increases with density (see equation 2). As such, lower density solutions are found over a broad range of velocity gradients, while higher density models are limited to larger values of dv/dr. This accounts for the positive slope of the dv/dr distribution. Similarly, the upper limit placed on dv/dr generates a negative slope in the K vir distribution. As the distribution functions for these two parameters are heavily influenced by our prior restrictions, our modeling produces no meaningful constraint on the dynamical state of the gas.
We take the median of each distribution as the single best estimate of the parameter value, and indicate this number with a vertical line in each panel of Figure 4. Most of the distribution functions are relatively symmetric, and we calculate the equivalent 1σ uncertainties in the parameter estimation by finding the range of values symmetrically enclosing 68% of the total probability. This range is shown as the hatched area under each distribution. The results of this analysis are summarized in Table 3. The asymmetry in the ME T kin distribution results from the limit to low temperature, high density parameter space following the K vir ≥ 1 and M gas ≤ M dyn . For this parameter, as well as for dv/dr and K vir , we extend the acceptable range listed in Table 3 to the truncated edge of the distribution.

Physical Parameters and Mass Fraction
The bulk of the emission traced by high resolution millimeter and near-IR molecular gas maps in the central 10 ′′ arises from the ∼5 ′′ (∼350 pc) CND (Schinnerer et al. 2000;Galliano & Alloin 2002;García-Burillo et al. 2010). Our LVG modeling indicates that no more than half of the CO(1-0) emission in the CND is generated by our ME and HE components, indicating lower excitation material is present. The physical conditions of this low excitation component have been studied by many groups. Tacconi et al. (1994) detected strong CO(4-3) emission toward the center of NGC 1068, and combined this with an interferometric CO(1-0) measurement to show that the gas was both warm (T kin ≥ 70 K) and dense (n H 2 ≥ 2 × 10 4 cm −3 ). Subsequent modeling of J upper ≤ 4 transitions of 12 CO and 13 CO have typically adopted a fixed T kin = 50 K, and derived n H 2 ∼ 10 4 − 10 5 cm −3 Helfer & Blitz 1995;Usero et al. 2004). Krips et al. (2011) have ob-tained fluxes of the lowest 3 transitions of both 12 CO and 13 CO with 2 ′′ resolution. For a subsection of the CND in which the gas appears perturbed by the radio jet, and hence possibly shock-heated, Krips et al. (2011) derive T kin ≥ 200 K and n H 2 = 10 3.5−4.5 cm −3 . While this section may not be representative of the CND as a whole, this analysis does suggest that the global temperatures may be higher than assumed by previous authors, and similar to that derived for our ME component. We also note that Kamenetzky et al. (2011) have similarly derived a globally high temperature (T > 100 K) for the CND based on an analysis of CS and other high density tracers. Given this range of temperature estimates, the clearest difference between the ME CO component and the lower excitation material traced by the low-J CO lines is the higher density (n H2 ∼ 10 5.7 cm −3 vs n H2 10 5 cm −3 ). This suggests a scenario in which the FIR lines trace denser material in the CND, which coexists with a more diffuse medium that generates the millimeter CO emission.
The mass fraction of the high excitation gas may be estimated by comparing the ME and HE masses derived here with the mass traced by CO(1-0). With a standard Galactic conversion factor N (H 2 )/I CO = 2 × 10 20 cm −2 (K km s −1 ) −1 , Schinnerer et al. (2000) estimate a mass of M H 2 = 5 × 10 7 M ⊙ for the CND. This is similar to previous estimates from Planesas et al. (1991, af-ter correcting their mass to the d = 14.4 Mpc used here) and Helfer & Blitz (1995), both of which used similar conversion factors. However, Usero et al. (2004) have derived a much lower N (H 2 )/I CO = 0.3 × 10 20 cm −2 (K km s −1 ) −1 (for the [CO/H 2 ] = 10 −4 used here) from their excitation modeling of the low-J CO emission from the CND. Papadopoulos & Seaquist (1999b) and Israel (2009) have also derived conversion factors significantly lower than the Galactic value by analyzing the low-J CO emission averaged over beam sizes of ∼1 ′ and ∼21 ′′ , respectively. These authors have suggested the lower conversion factor arises from the gas not being virialized. At the same time, Krips et al. (2011) have reported a CO(1-0) flux from the central 4 ′′ of 120 Jy km s −1 , much larger than the 20 Jy km s −1 value reported by Schinnerer et al. (2000). For N (H 2 )/I CO = 0.3 × 10 20 cm −2 (K km s −1 ) −1 , and F CO(1−0) = 20 − 120 Jy km s −1 , we estimate M H 2 = (0.5 − 3) × 10 7 M ⊙ . This is comparable to the M H 2 = (0.1 − 8) × 10 7 M ⊙ range we associate with our ME component. While the uncertainties of both numbers are high, this suggests that the ME component makes a non-negligible contribution to the total mass budget. The HE emission traces a lower mass of warmer gas, albeit still cooler than the small amount (M ∼ 10 3 M ⊙ ) of hot (T ∼ 2000 K) gas detected in the near-IR H 2 lines (Rotaciuc et al. 1991;Blietz et al. 1994;Galliano & Alloin 2002).

PACS Lines
Further insight into the nature of the ME and HE components may be obtained by comparing the FIR CO line profiles with those of other molecular tracers. In particular, we seek to understand the physical origins of the velocity shift between the two components. To bet- 10 5.7 10 4.9 − 10 6.8 10 6.3 10 5.7 − 10 7.3 P/k B [K cm −3 ] 10 7.9 10 7.4 − 10 8.8 10 9.0 10 8.5 − 10 9.
ter demonstrate the spectral shift, in Figure 5  With each line having a calibration error of σ i , we estimate the error (σ) in the composite spectra as σ 2 =< σ 2 i > /N . This gives σ ME = 11 km s −1 and σ HE = 14 km s −1 , and we show these as horizontal error bars in Figure 5.
In addition to CO, the PACS scans detect strong molecular emission from several transitions of OH and H 2 O, and weaker emission from OH + , H 2 O + , and other molecules. Some of the OH and H 2 O lines show extended emission over the PACS array, but for each line the flux in the central spaxel is dominated by a compact source that we associate with the CND. The OH + and H 2 O + emission is also unresolved, although an association of these lines with the same molecular gas is less certain. In Figure 6 we compare the central velocities and widths of these lines. For OH, each point represents the average properties of a doublet, while for OH + and H 2 O + each point is an average of multiple fine-structure transitions. On the right hand side of the top panel we show the mean central velocity of each molecule, with the ME and HE CO components separately. The HE CO is the only tracer systematically offset from the systemic velocity. In the bottom panel we show the widths of the ME and HE composite profiles as open points at λ = 109 µm, where the smoothed 311 km s −1 resolution is equal to the instrumental resolution. Overplotted are a set of curves showing the expected measured linewidths for intrinsic linewidths of 0, 250, and 400 km s −1 . The individual CO lines, as well as the composite profiles, are consistent with an intrinsic FWHM ∼ 250 km s −1 for both the ME and HE components. The other molecular lines are somewhat narrower.
The ∼50 − 70 km s −1 offset between the HE lines and the other molecular tracers in Figure 6 is 1/3 of the PACS spectral resolution, and at this level instrumental effects must be considered. The PACS spectral response depends on the illumination of the slit, such that a physical translation of a source on the sky in the dispersion direction will produce a wavelength shift in the spectral profile (Poglitsch et al. 2010). In the top panel of Figure 6 we show the equivalent velocity shift resulting from moving a point source ±2 ′′ from the center of the slit. A comparison with the CO velocities shows that a ∼2 − 3 ′′ offset between the centers of the slit and the HE CO emission could be partially responsible for the measured wavelength shifts. With the exception of CO(17-16) and CO(24-23), each point in Figure 6 corresponds to a line measured in the same set of concatenated range scans. Therefore if all of the molecular emission is presumed to arise from the same region on the sky, a wavelength shift in the HE CO lines induced by a pointing offset should be matched by a comparable shift in the λ ≈ 100 − 150 µm OH and H 2 O lines, and this is not observed. Alternatively, the HE CO lines may be modeled as arising from a region centered ∼2 − 3 ′′ away from the source of the OH and H 2 O emission. However, the CO(24-23) line at λ = 109 µm was observed 6 months prior to the series of concatenated scans, and the slit position angles of these two epochs differ by ≈180 • . Any pointing-induced shift to the CO(24-23) line should then have a comparable magnitude but be in the opposite direction of the shifts in the other HE CO lines, and this is also not the case. We therefore argue that the observed wavelength shift represents a real velocity offset, with instrumental effects playing no more than a minor role.

High Resolution Observations
To search for kinematic substructure in the CND that might explain the velocity shifts among the FIR CO lines, we compare the FIR line profiles with those of high resolution observations of CO(2-1) (Schinnerer et al. 2000) and H 2 1-0 S(1) (Müller Sánchez et al. 2009). In Figure 7a we show maps of these two tracers in the central 4 ′′ × 6 ′′ . At 0.7 ′′ resolution the CO(2-1) emission from the CND is dominated by a pair of knots centered ∼1 ′′ east and ∼1.5 ′′ west of the AGN, connected by lower surface brightness emission. Similar morphology is also seen in the ∼1 ′′ resolution aperture synthesis maps of CN and SiO (García-Burillo et al. 2010) and of HCN and HCO + (Krips et al. 2011). Observations of the H 2 1-0 S(1) line at 0.075 ′′ resolution have resolved the CND into a ∼1.2 ′′ radius ring (hereafter the "H 2 ring") centered ∼0.6 ′′ southwest of the AGN. The eastern knot is also prominent in the H 2 map, while the western knot is much fainter. The H 2 map also shows strong emission from a clump ∼1 ′′ to the north of the AGN that is not prominent in the CO(2-1) map.
In panel b we show separate maps of the CO(2-1) emission integrated over blue and red velocities (see also Figure 6 in Krips et al. (2011)). The strong emission to the east is largely to the blue of the systemic velocity, while the emission to the west separates into a blueshifted SW component and a redshifted NW component. The H 2 velocity field is generally similar, with the exception of the northern clump as discussed below. Schinnerer et al. (2000) model the CO(2-1) kinematics as a warped disk, while subsequent study of the H 2 and millimeter tracers has shown evidence for additional noncircular motion, possibly indicating a radially expanding component (Galliano & Alloin 2002;Davies et al. 2008;García-Burillo et al. 2010;Krips et al. 2011). In panels e through i we compare the composite ME and HE spectral profiles with those of CO(2-1) and H 2 extracted from selected apertures. For each panel we have smoothed the CO(2-1) and H 2 spectra to the same resolution (311 km s −1 ) and resampled to the same grid as used for the ME and HE composite spectra.
The most important result of this comparison is the mismatch between the FIR line profiles and the broad and highly blueshifted H 2 emission from the bright knot ∼1 ′′ to the north of the AGN, shown in panel f. In the high spectral resolution H 2 map presented in Galliano & Alloin (2002) this is the one region in the H 2 ring that displays a double-peaked profile (the two peaks become blended following the smoothing done here), with an extra emission component to the blue that is a clear outlier to their simple rotation plus expansion model. This region is also the site of the strongest Brγ emission in the H 2 ring, which displays a broad (FWHM ∼ 1000 km s −1 ) line. The radio jet exits the nucleus to the north (Gallimore et al. 1996), and may interact with material in the narrow-line region (Axon et al. 1998). Galliano & Alloin (2002); Galliano et al. (2003) interpret the broad and complex Brγ and H 2 line profiles as arising from gas perturbed by the jet, and additionally suggest that the ionization resulting from the jet-ISM interaction may be responsible for the strength of the Brγ emission in this region. The smoothed H 2 profile in panel f is inconsistent with either the ME or HE composite spectrum, and we conclude that this region of the CND does not dominate the PACS CO emission. In sections 5 and 6 we discuss possible excitation mechanisms for the FIR CO, including shock heating. The fact that we can rule out an origin in the region of the H 2 ring with the best evidence for molecular gas perturbed by the jet leads us to suggest that such jet-driven shocks in the H 2 ring do not excite the high-J CO.
Aside from the bright H 2 knot in the north, the ME profile is generally consistent with much of the rest of the H 2 ring. The best fit is with the profile of the strong H 2 emission from the east. The profiles from the NW and SW are moderately red and blueshifted with respect to the ME composite, but a combination of these regions, and indeed of the H 2 emission integrated over the entire CND (excluding the northern knot), would also generate a reasonable fit. We note that the H 2 linewidths (FWHM = 390 − 440 km s −1 ) are larger than the CO(2-1) linewidths (FWHM = 330−350 km s −1 ) in each panel, and better match the widths of the composite ME and HE profiles (FWHM = 400 − 410 km s −1 ). This may  1) (black solid) and CO(2-1) (black dashed) profiles from selected apertures in the H 2 ring compared with the ME (red) and HE (blue) composite profiles, and ±40 km s −1 horizontal error bar indicating the calibration uncertainty between the FIR CO and H 2 spectra, d,j) smoothed H 2 1-0 S(1) profiles from the northern and southern streamer compared with the ME and HE composite profiles, and ±40 km s −1 horizontal error bar.
indicate that the hot gas probed by the H 2 is a better tracer of the material producing the high-J CO, although due to extinction of the H 2 line (see, e.g., Galliano et al. 2003), the morphology of the FIR CO emission may be different than the H 2 image.
The blueshift of the HE emission is more challenging to match. The H 2 spectra from the blue regions in the E and SW are ≈60 and ≈40 km s −1 to the red, respectively, while the CO(2-1) profiles from the same regions are too narrow. However, the uncertainty in the mean centroid of the HE lines is ≈20 km s −1 , while the H 2 calibration error is smaller. Registration errors between the velocity frames of the PACS spectra and the ground-based H 2 spectra are also likely to be important at the level of ∼10 − 20 km s −1 , although more difficult to quantify. With a conservative estimate of 40 km s −1 uncertainty in the relative calibration of the PACS CO and the H 2 spectra, the HE CO emission only differs by 1 − 1.5σ with the H 2 emission from the E or SW regions. We conclude that while the HE emission profile is not natu-rally matched to the H 2 or CO(2-1) emission from the H 2 ring, an association with the bluest material to the E or SW may be within the measurement errors, and should not be excluded.
In panel c we show the 0.025 ′′ resolution H 2 1-0 S(1) map, which identifies two gas clouds streaming toward the nucleus on highly elliptical orbits from the north and south (hereafter the "northern" and "southern" streamers; Müller Sánchez et al. 2009). The northern streamer is connected to the H 2 ring in both H 2 emission and mid-IR continuum, and has been proposed as a means by which material is transported to the AGN (Tomono et al. 2006;Müller Sánchez et al. 2009). The southern streamer is detected to within ∼10 pc of the AGN, and may play a role in obscuring the nuclear emission. The southern streamer is modeled to lie in front of the AGN in the plane of the galaxy, and has a redshifted velocity that increases from ∼50 km s −1 to ∼80 km s −1 as the gas moves to within a projected distance of 0.1 ′′ . The northern streamer approaches the AGN from behind, and the brightest emission occurs 0.4 ′′ from the center at a blueshifted velocity of ∼ −30 km s −1 . In panels d and j we compare the smoothed profiles of these two regions with the PACS lines. The emission from the southern streamer is detected with lower S/N, but displays a profile reasonably consistent with that of the ME CO composite. The line profile of the northern streamer produces an excellent match to the HE composite. The centroids of the southern streamer and the HE CO composite differ by ≈100 km s −1 , which we argue is too large to be plausibly explained by calibration uncertainties between the two datasets. We conclude that in addition to an origin in the H 2 ring as discussed above, an origin of the ME and HE components with the southern and northern streamers, respectively, would be consistent with the line profiles.

HEATING THE ME COMPONENT
The ME CO emission arises from a warm (T kin ∼ 146 K) and dense (n H2 ∼ 10 5.7 cm −3 ) component, and with a total mass of M H2 ∼ 10 6.9 M ⊙ , represents an important fraction of the ISM in the CND. The kinematic analysis presented in section 4.2 shows that this component may be readily attributed to the H 2 ring, or possibly to the southern streamer. Here we consider potential heating mechanisms, and conclude that X-ray and shock heating are both plausible, while heating by far-UV photons is less likely. We further argue that no plausible heating mechanism is consistent with an origin in the southern streamer, and hence the emission is likely associated with the H 2 ring.

X-Ray Heating
The AGN in NGC 1068 emits a hard X-ray luminosity of L 2−10 keV = 10 43−44 erg s −1 (Iwasawa et al. 1997;Colbert et al. 2002). Our view of the AGN is obscured by a Compton-thick medium, but the extended emission detected by Chandra in the 6 − 8 keV band demonstrates that the nuclear X-rays irradiate the ISM over the central ∼kpc (Ogle et al. 2003;García-Burillo et al. 2010). Hard X-rays penetrate deeply into clouds, and efficiently heat large columns of molecular gas through photoionization heating (Maloney et al. 1996). The bright H 2 1-0 S(1) emission in the CND has been attributed to X-ray heated gas (Rotaciuc et al. 1991;Maloney 1997;Galliano & Alloin 2002), and X-ray heating should also produce strong emission in the FIR CO lines (Krolik & Lepp 1989). Meijerink & Spaans (2005) and Meijerink et al. (2007) present a detailed photochemical modeling of X-ray dominated regions (XDRs) that includes predictions for the emergent CO line intensities as a function of the gas density (n H ) and incident hard X-ray flux (F X = F 2−10 keV ). We use their type A models, which calculate the emission from a parsec thick cloud over a grid covering n H = 10 4 − 10 6.5 cm −3 and F X = 1.6 − 160 erg cm −2 s −1 . These models generate CO line SEDs with similar shapes as the isothermal models used in our LVG analysis, and an analogous two component fit reproduces the FIR CO line fluxes. In Figure 8a we show a model that uses n H = 10 5.75 cm −3 and F X = 9 erg cm −2 s −1 for the ME emission (see Table 4). For an AGN luminosity of L 2−10 keV = 10 43−44 erg s −1 , geometric dilution of the radiation field at the d ∼ 100 pc distance of the H 2 ring yields a flux of F X = 8.4 − 84 erg cm −2 s −1 , broadly consistent with this modeled flux. For the plane-parallel geometry employed by Meijerink et al. (2007) the absolute line luminosities scale with the total XDR surface area, and the normalization of the ME component of the model shown in Figure 8a requires A ∼ (120 pc) 2 . Galliano et al. (2003) model the H 2 ring as a section of a 40 pc thick disk. For a radius of 100 pc the inner surface area exposed to the AGN is then ∼(160 pc) 2 , similar to the XDR model requirement. In sum, we conclude that if a substantial fraction of the H 2 ring is exposed to nuclear hard X-rays, then both the shape of the ME segment of the CO SED and the absolute line fluxes are naturally reproduced in this XDR model.

Far-UV Heating
Far-UV (FUV; 6 eV < hν < 13.6 eV) photons offer another means of heating the molecular gas. Photodissociation regions (PDRs) powered by the FUV radiation from OB stars in the Galaxy have indeed been identified as prominent sources of FIR CO emission (e.g., Kramer et al. 2004), as have the FUV-irradiated cavities of protostellar outflows (van Kempen et al. 2010). However, a comparison with the high-J CO emission from the prototypical starburst galaxy M82 suggests that the CO emission from NGC 1068 is too strong to be powered by stellar FUV. Herschel observations of M82 have detected the submillimeter CO transitions (J upper = 4−13), which trace emission from gas heated by shocks (Panuzzo et al. 2010) and/or in PDRs (Loenen et al. 2010). For an L FIR ∼ 3 × 10 10 L ⊙ (Telesco & Harper 1980;Joy et al. 1987), the CO(13-12)/FIR ratio integrated over the central 43.4 ′′ (∼821 pc) is ∼6×10 −6 . The FIR continuum in M82 is produced by the FUV output of young stars, and we interpret this CO(13-12)/FIR ratio as a benchmark estimate of the fraction of FUV radiation that may be converted to FIR CO emission in an extragalactic starburst. The CO(14-13)/FIR ratio we estimate for NGC 1068 is ∼3 × 10 −5 , a factor of ∼5 larger than the CO(13-12)/FIR ratio in M82. This discrepancy is further increased by noting that only a minor fraction of the FIR continuum from the CND in NGC 1068 originates in ... n = 10 6.5 cm −3 n = 10 6 cm −3 ... G 0 = 10 4.75 G 0 = 10 5 ... L FUV ∼ 2 × 10 9 L ⊙ L FUV ∼ 10 10 L ⊙ shock C-shock C-shock ... n 0 = 2 × 10 5 cm −3 n 0 = 10 6 cm −3 Note. -Details for the models used in Figure 8. XDR and PDR models are from Meijerink et al. (2007), ME C-shock model is from Flower & Pineau Des Forêts (2010), and HE C-shock model is from Kaufman & Neufeld (1996).  Table 4. young stars. The nuclear stellar cluster in NGC 1068 has a characteristic age of 200 − 300 Myr, and hence L bol /L K 70 (for a starburst with exponential decay timescale of τ SF = 100 Myr; Davies et al. 2007). For the L K measured by Davies et al. (2007) this corresponds to L bol 3×10 9 L ⊙ , more than 5 times lower than the L FIR measured here. PDRs heated by the output of this cluster would therefore have to be ∼25 times more efficient in converting the incident FUV radiation to FIR CO emission than the PDRs in the nucleus of M82. Loenen et al. (2010) model the J upper = 1 − 13 CO emission from M82 with a 3 component PDR model, in which the CO(13-12) emission is almost entirely generated by the highest excitation component. The CO(13-12)/FIR ratio of this component is comparable to the measured CO(14-13)/FIR ratio in NGC 1068, but the total modeled ratio is a factor of ∼14 lower due to the FIR emission from the lower excitation material. Increasing the net ratio by a factor of ∼25 would require the total FIR be dominated by the highest excitation component, and hence a significant suppression of the lower excitation PDR emission. We consider that such an extreme repartitioning of the stellar FUV between M82 and NGC 1068 less likely than the XDR scenario discussed in section 5.1.
The CO emission is more plausibly attributed to PDRs powered by the AGN FUV radiation, and in many ways this is an attractive option. The intrinsic AGN luminosity in NGC 1068 of L FUV ∼ 7.4 × 10 43 erg s −1 (Pier et al. 1994) is similar to the observed L FIR , supporting the idea that the FIR continuum is emission from dust grains heated by the AGN. The measured 60 µm/100 µm flux ratio of ∼1.3 implies that the incident FUV flux on these grains corresponds to G 0 ∼ 10 5 (where (Abel et al. 2009), which would also be expected if the AGN L FUV is absorbed at the d ∼ 100 pc distance of the H 2 ring. To estimate the CO emission produced in these PDRs we employ the type A PDR models of Meijerink et al. (2007), which adopt the same geometry and gas densities as the XDR models, and span a FUV intensity range of G 0 = 10 2 − 10 5 . In dense PDRs most of the CO forms at cloud depths greater than A V ≈ 2 at gas temperatures T 100 K, but a secondary CO abundance peak at A V ≈ 0.6 with T ∼ 800 K follows from a local enhancement of OH (Sternberg & Dalgarno 1995). The CO line SEDs generated in the Meijerink et al. (2007) models are characterized by two distinct peaks, which we associate with these warm and hot CO phases. In Figure 8b we show a model with n H = 10 6 cm −3 and G 0 = 10 5 that produces an intriguing match to the entire FIR CO line SED (solid black line), with the warm and hot phases reproducing the ME and HE components, respectively. Assuming the emitted FIR continuum is equal to the incident FUV, this model predicts a CO(14-13)/FIR ratio of 4 × 10 −5 . This is similar to the observed ratio in high-G 0 Galactic PDRs (Kramer et al. 2004), and the value measured here for NGC 1068.
However, there are at least 2 reasons to reject this scenario. First, the AGN UV/optical emission escapes to d ∼ 100 pc distances only in the ionization cone (PA ≈ 15 • Macchetto et al. 1994), while the brightest emission from the H 2 ring is to the east. If the warm molecular gas traced by the ME CO emission is assumed to be cospatial with the hot gas traced by the H 2 1-0 S(1) line, then we would expect little interaction between this warm gas and the nuclear FUV. Additionally, any single PDR model that simultaneously matches the ME and HE emission is inconsistent with the kinematic evidence that these two sets of lines are tracing physically distinct components. In general, we require that any model reproducing the J upper 17 lines must underpredict the fluxes in the J upper ∼ > 20 transitions. In contrast, all of the models produced by Meijerink et al. (2007) that fit the ME section of the CO SED either match or exceed the observed HE lines fluxes. Consequently we exclude PDRs as a potential source of the ME emission. We stress that this latter argument depends sensitively on an accurate modeling of the J upper ∼ > 20 CO emission from the hot surfaces of PDRs. Future Herschel-PACS studies of the FIR CO line SEDs in Galactic PDR templates will be useful in further evaluating this PDR scenario.

Shock Heating
Shock heating offers a simple means of exciting the high-J CO emission, and is generally the preferred mechanism for producing the highest excitation molecular gas in Galactic sources (e.g., Sempere et al. 2000). Using the C-shock models of Flower & Pineau Des Forêts (2010), we find that the ME component of the CO SED can be fit with a preshock density of n H = 2 × 10 5 cm −3 and shock velocities of v s = 10 − 20 km s −1 . The contributions of H 2 O and H 2 to the total cooling budget increase rapidly with shock velocity in these models, and keeping v s ≤ 20 km s −1 is necessary to prevent the predicted H 2 O and H 2 line fluxes from exceeding the measured values. Decreasing the shock velocity from v s = 20 → 10 km s −1 generates weaker CO lines, and requires increasing the total cross-section from A ∼ (130 pc) 2 → (200 pc) 2 to match the absolute line fluxes. In Figure 8c we combine the v s = 20 km s −1 model with a separate C-shock fit to the HE CO emission.
What shock mechanism could produce the ME CO line emission? Mechanical heating from stellar feedback has been proposed as the energy source behind the submillimeter CO emission in the M82 starburst (Panuzzo et al. 2010) and in other galactic nuclei (Hailey-Dunsheath et al. 2008;Nikola et al. 2011), but the contrast in CO/FIR ratios between NGC 1068 and M82 discussed in section 5.2 argues against a similar mechanism here. Furthermore, we can use the results of Davies et al. (2007) to estimate the mechanical power injected into the ISM by supernovae (SN) in the nuclear cluster of NGC 1068. These authors model the supernovae rate (SNR) to L K ratio as 6 × 10 −11 yr −1 L −1 ⊙ in this cluster, which combined with their measured L K yields SNR 2.4 × 10 −3 . Following Loenen et al. (2008) and estimating a mechanical energy release of 10 51 erg per SN, with 10% dissipated in molecular gas, yields a heating rate of 2 × 10 6 L ⊙ . This is at least a factor of ∼8 − 28 times lower than the mechanical luminosity L = 1/2ρv 3 s A required by the shock models discussed above. Jet-ISM interactions are another potential source of shock heating in Seyfert nuclei, but as discussed in section 4.2, the mismatch in line profiles between the FIR CO lines and the disturbed H 2 1-0 S(1) profile ∼1 ′′ north of the AGN suggests jet-driven shocks in the H 2 ring may not be important.
Alternatively, we note that the cross-sections required to normalize these plane-parallel shock models are similar to the estimated cross section of the H 2 ring as viewed from the center. The kinematics drawn from the high resolution H 2 1-0 S(1) and millimeter-wave molecular gas maps require a radial expansion component to the H 2 ring (Galliano & Alloin 2002;Davies et al. 2008;Krips et al. 2011). Galliano & Alloin (2002) model the H 2 dynamics by combining a rotational component with a v = 140 km s −1 radially expanding component that generates 1/3 of the total line emission, and Krips et al. (2011) construct a similar model to explain the CO dynamics. We suggest that the interaction of the outflowing gas with non-outflowing material in the H 2 ring offers the most plausible source of shock heating. We recall that our excitation model requires the dense gas associated with the ME CO emission to be mixed with lower density (n H 2 10 5 cm −3 ) material responsible for the millimeter-wave CO emission (section 4.1). Assuming the nv 2 product is conserved for shocks propagating through an inhomogeneous medium (Klein et al. 1994), the moderate (v 20 km s −1 ) velocities we require could ultimately be produced in v ∼ 140 km s −1 shocks triggered in lower density gas.

Southern Streamer
In section 4.2 we showed that the ME CO line profile may be consistent with that of the H 2 1-0 S(1) emission from the southern streamer. However, such an association is inconsistent with either the X-ray or shock heating scenarios outlined above. In both of these models the required cross-section of A ∼ (120 pc) 2 − (130 pc) 2 is significantly larger than the ∼20 pc size of the southern streamer. Equivalently, the absolute intensity of the CO emission from this cloud would have to be more than an order of magnitude larger than predicted by the models. As such, we argue that the ME CO emission does not arise from the southern streamer, but most likely arises from the H 2 ring.

Summary and Discussion
In summary, we conclude that the ME CO emission arises from either X-ray or shock-heated gas in the H 2 ring.
The challenge of unambiguously determining the heat source for this gas is similar to the situation for the warm and hot gas traced by the H 2 rotational and ro-vibrational lines in NGC 1068 (Rotaciuc et al. 1991;Lutz et al. 2000), and in larger samples of Seyferts (Davies et al. 2005;Rodríguez-Ardila et al. 2005;Roussel et al. 2007). The link between the ME CO and the H 2 emission in NGC 1068 is not immediately clear, although at least to within the PACS resolution the similarity of the ME CO and H 2 1-0 S(1) line profiles suggests a connection (Figure 7). Additionally, we note that the H 2 1-0 S(1) brightness is quantitatively reproduced with a similar XDR model as used here for the ME CO (although using a lower density of n = 10 5 cm −3 ; Maloney 1997; Galliano et al. 2003), while the shock model shown in Figure 8c also produces an H 2 1-0 S(1) flux within a factor of ∼2 of the total CND emission. Further joint modeling of the CO, H 2 , and other molecular emission in NGC 1068, and FIR CO data on a larger sample of comparison sources, will be useful in better understanding the nature of the ME CO component.
Here we offer two reasons for preferring the X-ray heating scenario. First, the hard X-ray luminosity of NGC 1068 is reasonably well established through either an analysis of the directly observed (scattered) emission, or by applying scaling relations established for type 1 systems (Iwasawa et al. 1997;Colbert et al. 2002). As discussed above, combining this luminosity with reasonable estimates of the gas density and the CND geometry naturally produces the ME CO line SED. In contrast, it is not clear whether the shocks in the CND are sufficient to dissipate enough mechanical energy at the low velocities needed to power the CO. As we argued above, jet-driven shocks in the H 2 ring do not provide a good fit, while the energetics of the possible shocks arising from the radial expansion of the H 2 ring have yet to be demonstrated. Secondly, a number of authors have noted that the chemical composition of the CND is best described through a model of X-ray-driven chemistry (Usero et al. 2004;Krips et al. 2008;García-Burillo et al. 2010;Krips et al. 2011). The OH + and H 2 O + emission detected in our PACS scans offers further evidence (van der Werf et al. 2010;Rangwala et al. 2011), and will be discussed in a future paper. If the nuclear X-rays are responsible for the anomalous molecular abundances in the CND, it is likely they also play an important role in the energetics.

HEATING THE HE COMPONENT
The HE CO emission arises from a small mass (M H2 ∼ 10 5.8 M ⊙ ) of warm (T kin ∼ 435 K) and dense (n H2 ∼ 10 6.3 cm −3 ) material that represents only a minor fraction ( 10%) of the total gas in the CND. The kinematic analysis presented in section 4.2 suggests that this component may potentially be associated with the most blueshifted emission in the east or west of the H 2 ring, or with the clump of infalling gas ∼0.4 ′′ north of the AGN. Here we consider potential heating mechanisms, and argue that the HE CO emission arises from either Xray-heated gas in the northern streamer, or shock-heated material in either the northern streamer or H 2 ring.
6.1. X-ray Heating In section 5.1 we used the Meijerink et al. (2007) XDR models to generate a two component fit to the FIR CO emission, and overploted a sample solution in Figure 8a. For the n H = 10 5 → 10 6.5 cm −3 density range the ME component requires F X = 5.1 → 16 erg cm −2 s −1 , while for the same densities the HE component requires the maximum (or higher) F X = 160 erg cm −2 s −1 used in the Meijerink et al. (2007) model grid. An important outcome of this XDR modeling is therefore that the HE component requires irradiation by an X-ray field at least an order of magnitude stronger than the ME component. This is difficult to achieve if both components are presumed to arise from the H 2 ring. The eastern segment of the ring lies no more than a factor of ∼1.5 closer to the AGN than the western segment, so a variation in the distance to the AGN across the ring appears unlikely to produce a factor of ∼10 variation in the radiation field strength. García-Burillo et al. (2010) have used the 6 − 8 keV emission detected with Chandra to trace the penetration of nuclear hard X-rays into the cold neutral ISM. Their image suggests a relatively isotropic illumination of the central few hundred parsecs (see also Ogle et al. 2003). Specifically, while the hard X-rays may be expected to escape to ∼100 pc scales more easily in the ionization cone, the 6 − 8 keV band image is not significantly enhanced along this direction within the H 2 ring.
The most straightforward way of generating a higher X-ray intensity is to attribute the HE emission to the northern streamer. In projection, the northern streamer lies a factor of ∼3 closer to the AGN than does the H 2 ring (∼0.4 ′′ vs ∼1.2 ′′ ; Müller Sánchez et al. 2009), and is therefore expected to see a ∼9 times stronger F X . A fit to the HE component with n H = 10 5.75 cm −3 and F X = 160 erg cm −2 s −1 (Figure 8a) requires a surface area of ∼(25 pc) 2 , only a moderate factor larger than the ∼14 pc (∼0.2 ′′ ) projected lateral size of this clump. We therefore argue that if X-ray heating accounts for both the ME and HE components, the HE component is most plausibly associated with this infalling gas.

Far-UV Heating
Using the Meijerink et al. (2007) PDR models we find that for a proper choice of high density (n H ∼ > 10 6 cm −3 ) and strong FUV field (G 0 ∼ > 10 4 ), the models reproduce the J upper ∼ > 20 lines while underpredicting the lower-J fluxes. An example model with n H = 10 6.5 cm −3 and G 0 = 10 4.75 is shown in Figure 8b (thin blue line). The FUV continuum associated with this component is L FUV ∼ 2 × 10 9 L ⊙ , a factor of ∼10 less than the observed L FIR . Thus we may consider a scenario in which the ME CO emission arises from X-ray or shock-heated gas in the H 2 ring, while the HE emission originates in a trace amount of PDR material that makes only a minor contribution to the total continuum emission.
Any gas in the CND exposed to the FUV from the AGN must also be irradiated by hard X-rays. For a FUV origin of the HE CO lines the physical properties and radiative environment of the clouds must therefore be such that the FUV is more effective than the X-rays in generating J upper ∼ > 20 CO emission. Under what conditions would this occur? A first approach to addressing this question is to separately consider the CO emission from a PDR with that of a properly selected XDR. In Figure 9 we compare the CO(19-18) (the highest-J transition calculated over the full model grid) flux from the Meijerink et al. (2007) PDR models with that of a Meijerink et al. (2007) model XDR with ≈1/6 Fig. 9.-Ratio of CO(19-18) emission from a PDR to an XDR (black ), using the models of Meijerink et al. (2007). The incident X-ray flux in the XDR model is 6 times weaker than the FUV flux in the PDR model. The fraction of the measured L FIR that must be attributed to this PDR is shown in red. The region of high density and high-G 0 that can match the HE CO line SED without overpredicting the lower-J fluxes is indicated with the dotted lines.
of the incident flux. This flux ratio is chosen to match the intrinsic L FUV /L X ≈ 6 ratio of the AGN in NGC 1068 (Pier et al. 1994), and is therefore appropriate for a cloud seeing the unattenuated AGN emission. The red curves in Figure 9 show the fraction of the observed FIR continuum (assuming L FUV = L FIR ) that must arise from the model PDR in order to account for the absolute CO(20-19) line flux. The region in the upper right quadrant enclosed by the dotted lines indicates the subset of high excitation PDR models that reproduce the HE CO lines while underpredicting the lower-J fluxes -a requirement for this solution given the velocity shifts between the ME and HE lines. For clouds with n H ∼ > 10 6 cm −3 , G 0 ∼ > 10 4 , and G 0 /n H 10 −1.5 cm 3 , this approach suggests that the FUV is more important than the X-rays in generating J upper ∼ > 19 CO emission, while only a minor fraction of the AGN FUV luminosity would be required to power such a PDR.
Where in the nuclear region might these conditions be satisfied? The UV/optical emission from the AGN escapes to large distances within the ionization cone, which runs roughly north-south in projection. If the HE CO emission arises from the H 2 ring, however, the line blueshifts would suggest an origin to the east or west (section 4.2), where little nuclear FUV penetrates. The northern streamer may have a direct view of the AGN, which at a distance of d ≈ 40 pc would correspond to G 0 ∼ 10 5.4 . However, the covering factor of this material (for a size of ∼14 pc) is only ∼0.01. The contents of Figure 9 indicate that achieving the absolute HE CO line fluxes with such a small covering factor would require densities far larger than the n H 2 ∼ 10 6.3 cm −3 indicated by our LVG modeling (section 3.2). We therefore argue that the HE CO emission is unlikely to be FUV-powered, although more detailed modeling of the conversion of FUV radiation to high-J CO at high densities (n H > 10 6.5 cm −3 ) and FUV intensities (G 0 > 10 5 ), and in the presence of an additional X-ray field, would be useful for a more quantitative evaluation of this scenario.

Shock Heating
The HE CO component is fit by a C-shock model with preshock density n H 2 = 10 6 cm −3 , velocity v s = 30 km s −1 , and cross-section A ∼ (21 pc) 2 (Kaufman & Neufeld 1996) (Figure 8c). As discussed in sections 4.2 and 5.3, jet-ISM interactions in the H 2 ring are unlikely to produce the FIR CO emission, but a jet-induced shock in the northern streamer is a more plausible source of the HE CO lines. The radio jet changes direction in the vicinity of the brightest clump in the northern streamer, evidence that the fast moving, low density material in the jet is colliding with and being diverted by the dense molecular gas (Müller Sánchez et al. 2009). This interaction should produce shocks in the molecular material, and the presence of H 2 O masers in this cloud supports this picture (Gallimore et al. 2004). The modeled A ∼ (21 pc) 2 cross-section matches the ∼14 pc size of this clump, consistent with this scenario. The jet mechanical power estimated in the bow shock model of Wilson & Ulvestad (1987) is ∼2 × 10 8 L ⊙ , a factor of ∼3 larger than the mechanical luminosity L = 1/2ρv 3 s A of the shock model discussed here. Jet-induced shocks are therefore energetically feasible, although would require an efficient mechanism for converting kinetic energy at high velocities into slow molecular shocks. Alternatively, shocks in the bluest regions of the H 2 ring, possibly associated with the ring expansion, could also generate the HE CO.
An important constraint on this shock modeling is the need to match the CO brightness while not overproducing the ro-vibrational H 2 emission. For the C-shock model discussed above the predicted H 2 1-0 S(1) flux is a factor of ∼70 larger than observed in the northern streamer (Müller Sánchez et al. 2009), and ∼7 larger than in the bright eastern clump in the H 2 ring (Galliano & Alloin 2002). In this region of the Kaufman & Neufeld (1996) model parameter space the H 2 /CO intensity ratio is a steep function of shock velocity, while the CO emission is sensitive to both the velocity and preshock density. Finely-tuning v s and n 0 (going to lower v s and higher n 0 ) may therefore yield a solution that reduces the H 2 emission while conserving the CO line SED. Extinction of the 2 µm H 2 emission may also be important, particularly in the northern streamer. Müller Sánchez et al. (2009) estimate a column of N H ∼ 8 × 10 24 cm −2 in the southern streamer. Assuming a similar value for the northern streamer, and for A K /N H = (2×10 22 cm −2 ) −1 , the resulting A K ∼ 400 is more than sufficient to provide the necessary attenuation in a mixed dust model. A J-shock model for the HE CO would also generate much weaker H 2 emission (Hollenbach & McKee 1989), although the required cross-section would then be a factor of ∼5 higher than for a C-shock, and no longer match the northern streamer size. In sum, a number of scenarios may be envisioned to attribute the HE CO lines to a shock in the H 2 ring or the northern streamer, while also satisfying the H 2 1-0 S(1) brightness.

Summary and Discussion
We conclude that the HE component most plausibly arises from X-ray or shock heated gas in the northern streamer, or shock heated gas in the H 2 ring. For the ME component, we argued in section 5.5 that energetic considerations and the ancillary evidence for an XDR chemistry favored X-ray over shock heating, but the situation for the HE component is less clear. The amount of mechanical power available in the jet or other sources is indeed uncertain, but the strong evidence for a jet interaction with the northern streamer suggests that jetinduced shocks in this cloud be given full consideration. Additionally, the OH + and H 2 O + emission lines that pinpoint an XDR chemistry are emitted close to the galaxy systemic velocity, with no detected emission at the blueshifted velocity of the HE CO lines (Figure 6). Thus while nuclear X-rays may dominate the energetics and chemistry in the CND, an origin of the HE CO lines in a small amount of shock-heated gas must also be considered.

NON-DETECTIONS, AND CONSTRAINTS ON THE NUCLEAR OBSCURATION
In addition to the detected J upper = 14−24 transitions, our upper limits to the J upper ≤ 50 lines provide a useful constraint on the nuclear molecular component. One of the earliest discussions of high-J CO in extragalactic sources was presented in Krolik & Lepp (1989), who modeled the emission expected from the Compton-thick material obscuring type 2 AGN. These authors calculated the CO emission from an N H = 10 24 cm −2 cloud located ∼1 pc from a luminous (L X ∼ 10 44 erg s −1 ) hard X-ray source, in which the heating is dominated by Compton scattering of 10 − 100 keV photons. In this model the CO SED scales as J 3 upper up to J upper ≈ 58, and the absolute line luminosities are proportional to the total absorbed 10 − 100 keV luminosity (f abs L x44 ; in units of 10 44 erg s −1 ). Several of our upper limits to CO transitions with J upper ≈ 30 − 50 independently place an upper limit to this component corresponding to f abs L x44 0.1, while the non-detection of the CO(44-43) line achieves the most stringent limit of f abs L x44 0.09 ( Figure 8a). This limit may be reduced by stacking the individual non-detections. We have obtained such a stack by first scaling the spectrum of each undetected line by (44/J upper ) 2 , which references each spectrum to that of  under the assumption that the CO fluxes scale as J 3 upper . We then calculated a weighted average of 8 lines with low noise that fall in clean spectral regions, and binned to 600 km s −1 . The result is shown in Figure 10. The upper limit to this stacked line pushes the limit of the Krolik & Lepp (1989) model to f abs L x44 0.038 (Figure 8a).
Given this range of estimates we adopt L 2−10 keV = 10 43.5 erg s −1 , and increase this by 1.43 to correct to the 10 − 100 keV range (for L ν ∝ ν −1 ). Assuming 25% of this power is absorbed in the phase modeled by Krolik & Lepp (1989) gives an expected f abs L x44 = 0.11, a modest factor of ∼2.9 higher than our measured upper limit. Given the uncertainties in the AGN luminosity and covering factor of the Comptionthick material in NGC 1068, as well as uncertainties in the model itself, this factor of ∼2.9 overprediction is too small to reject the model of a parsec-scale torus envisioned by Krolik & Lepp (1989). This non-detection may, however, provide a useful constraint on the more recent set of detailed torus models (e.g., Hönig et al. 2006;Nenkova et al. 2008).

SUMMARY AND FUTURE WORK
We have detected 10 CO transitions in the J upper = 14 − 24 range from the central 10 ′′ (700 pc) of NGC 1068, and obtained sensitive upper limits to most other transitions up to J upper ≤ 50. These are the first extragalactic detections of FIR CO, which represent a new probe of excited molecular gas in Seyfert nuclei. The detected transitions are modeled as arising from 2 different components: a moderate excitation (ME) component close to the galaxy systemic velocity, and a high excitation (HE) component that is blueshifted by ∼70 km s −1 . Our main results are as follows: 1) Using a two component LVG model we derive n H2 ∼ 10 5.7 cm −3 , T kin ∼ 146 K, and M H 2 ∼ 10 6.9 M ⊙ for the ME component, and n H2 ∼ 10 6.3 cm −3 , T kin ∼ 435 K, and M H 2 ∼ 10 5.8 M ⊙ for the HE component. The 1σ uncertainties on the derived temperatures are ∼±0.25 dex, while for density and mass this is ∼±1 dex. Extending the measured CO line SED to lower-and higher-J lines would likely reduce these uncertainties, as would a joint analysis with H 2 , OH, H 2 O, and other complementary molecular tracers. We will follow both approaches in forthcoming papers.
2) These two components are denser than the gas traced with millimeter CO observations, and the HE (and possibly the ME) component is also warmer. The ME component makes a non-negligible contribution to the nuclear mass budget, although large uncertainties in the masses estimated from both the FIR CO lines and from CO(1-0) prevent a more quantitative statement.
3) Comparing the FIR CO line profiles with those of high spatial and spectral resolution observations of CO(2-1) and H 2 1-0 S(1) allows a first estimate of the origins of the ME and HE components within the central 10 ′′ . Good matches are found with H 2 1-0 S(1), which for the ME component suggests an origin in the ∼200 pc diameter H 2 ring. The blueshifted HE lines may also be consistent with an origin in the bluest regions of the H 2 ring, but are better matched to the clump of infalling molecular gas ∼40 pc north of the AGN.
4) The ME component is nicely consistent with models of X-ray heating of gas in the CND. A shock model is also possible, although due to the uncertainties in the amount of mechanical power available for dissipation in slow shocks, and the evidence for X-ray driven chemistry in the CND, we favor the X-ray heating scenario. Far-UV heating is unlikely.
5) The HE component is also consistent with either X-ray or shock heating. X-ray heating would best fit with an origin in the cloud ∼40 pc north of the AGN, supporting the results of the line profile matching (point 3). Shocks triggered by the interaction of the radio jet with this clump, or arising from the H 2 ring, are also plausible. Far-UV heating is unlikely. 6) Our non-detections of J upper = 30 − 50 lines place a constraint on any gas associated with a putative parsecscale torus. Our upper limit to the Krolik & Lepp (1989) prediction is a factor of ∼2.9 lower than the expected signal, although due to the uncertainties involved in applying this model, this non-detection is insufficient to rule out the torus paradigm. The inclusion of a gas phase in the current set of clumpy torus models, and a comparison of the predicted CO line SED with our detections and upper limits, would be a useful next step.
We thank Eva Schinnerer for generously providing her IRAM PdBI CO(1-0) and CO(2-1) data, and Francisco Müller Sánchez for providing his SINFONI/VLT H 2 1-0 S(1) map. Basic research in IR astronomy at NRL is funded by the US ONR; J.F. also acknowledges support from the NHSC. S.V. acknowledges support from a Senior NPP Award from NASA and thanks the host institution, the Goddard Space Flight Center. PACS has been developed by a consortium of institutes led by MPE (Germany)