A DEEP HERSCHEL/PACS OBSERVATION OF CO(40-39) IN NGC 1068: A SEARCH FOR THE MOLECULAR TORUS

Low-J rotational lines (up to J $\approx$ 7) of CO have long been used as a tracer for molecular gas. With the Herschel Space Telescope⁹ (Pilbratt et al. 2010) it was possible to observe CO lines of higher J levels, extending the CO spectral line energy distribution (SLED) to $J=50$ (for $J\to J-1$). The PACS instrument (Poglitsch et al. 2010) on board Herschel samples the highest frequencies, with $13\;\lt \;J\;\lt \;50$, which we will call the high-J CO lines. These lines with wavelengths $186.0\;\mu {\rm{m}}\lt \lambda \lt 52.85\;\mu {\rm{m}}$ arise from states ∼500–7000 K above ground, with critical densities of $\sim {10}^{6}-{10}^{8}\;{\mathrm{cm}}^{-3}$, and are thus good tracers of the photon dominated regions (PDR) in the warm and dense interstellar matter (ISM) around stars, in shocks, or in the X-ray dominated regions (XDRs) of active galactic nuclei (AGNs).

AGN unification models require a structure that obscures the UV and optical emission from the broad line region of Type-2 AGN from our view (Antonucci & Miller 1985). Further, Compton-thick AGN like NGC 1068 (Elvis & Lawrence 1988) require an X-ray absorbing medium with a column density of at least ${N}_{H}={10}^{24}\;{\mathrm{cm}}^{-2}$, potentially in the form of molecular gas. These two effects are often assumed to be caused by a dusty torus, although the gas and the dust are not necessarily distributed in the same way. In this paper, we investigate whether high-J CO lines can trace the X-ray absorbing gas (from now on called "the torus"), be it in the form of a classical torus (Krolik & Begelman 1988; Krolik & Lepp 1989) or of a more clumpy and more extended cloud structure. The CO emission from a parsec-scale torus was first modeled by Krolik & Lepp (1989, hereafter KL89), who found that the line fluxes should scale as J³ up to $J\;\approx \;58$ (in this particular model for NGC 1068), and with absolute fluxes proportional to the total absorbed hard X-ray luminosity. More recent torus models have envisioned a dynamical structure of small clouds distributed over several parsecs surrounding the AGN (Nenkova et al. 2002; Hönig et al. 2006). However, these newer clumpy torus models do not yet incorporate a gas phase and have been used so far only to predict the dust emission from an AGN.

In Hailey-Dunsheath et al. (2012; hereafter HD12) we carried out a detailed analysis of the full CO SLED of NGC 1068 up to $J=30$ (for higher J only upper limits were found). We concluded that two components are responsible for the CO SLED between $10\lt J\lt 31$: the ∼200 pc diameter ring and the northern streamer, both traced by H₂ 1–0 S(1) observations (Müller Sánchez et al. 2009). Would the molecular torus, if it exists, be observable in the $J\gt 30$ lines? In order to answer this question, we carried out an additional, deeper observation of CO(40-39).

The paper is organized as follows: a description of the existing multiwavelength observations of NGC 1068 is given in Section 2. Section 3 covers the PACS observations and data reduction and presents the full (55–200 $\mu {\rm{m}}$) spectrum of the central spatial pixel (spaxel). Section 4 describes updates on the PDR, XDR, and shock fits presented in HD12, including the new upper limit on CO(40-39) and a small correction to the CO SLED. A brief description of the XDR code used to model the torus emission is given in Section 5. A discussion of the results follows in Section 6.

NGC 1068 is a Seyfert 2 galaxy at a distance of 14.4 Mpc, implying a 72 pc arcsec⁻¹ scale (Bland-Hawthorn et al. 1997). Because of its proximity and luminosity, NGC 1068 is a favorable object for testing hypotheses on Seyfert galaxies.

Early X-ray observations showed emission from the AGN nucleus, but only in scattered light (Elvis & Lawrence 1988) because the AGN is completely obscured from our view by a column of ${N}_{H}={10}^{25}\;{\mathrm{cm}}^{-2}$ (Bauer et al. 2014). The intrinsic 1–100 keV luminosity, estimated from the observed scattered X-rays, lies between ${10}^{43}$ and ${10}^{44}\;\mathrm{erg}\;{{\rm{s}}}^{-1}$. The largest uncertainty in the X-ray luminosity is the fraction of X-rays that is scattered into our line of sight. Pier et al. (1994) give a comparison of different methods to determine this fraction, from which we derive an intrinsic luminosity ${L}_{1-100\;\mathrm{keV}}={10}^{43.5}\;\mathrm{erg}\;{{\rm{s}}}^{-1}$, and a spectral slope ${F}_{\nu }\sim {\nu }^{-1}$ (or photon index ${\rm{\Gamma }}=2$).

Further out from the obscuring torus, the AGN is surrounded by a ring of molecular gas with a radius of $\sim 1^{\prime \prime} -5^{\prime \prime}$ (72–360 pc), also called the circumnuclear disk (CND). The CND is observed in CO and H₂ (Schinnerer et al. 2000; Müller Sánchez et al. 2009), and observations of SiO, HCN, and CN suggest an X-ray driven chemistry (Lepp & Dalgarno 1996; Usero et al. 2004; García-Burillo et al. 2010).

Inside this shell, the near-infrared integral field spectrograph SINFONI provided 0075 resolution observations of H₂ 1-0 S(1), which reveal two streamers within 04 (30 pc; Müller Sánchez et al. 2009). The southern streamer lies between us and the AGN's center. The SINFONI observations furthermore reveal the presence of the northern streamer (04 or 30 pc north of the AGN's center), which may fuel the central black hole.

A ring of H₂O masers, with a radius of ∼20 mas (1.5 pc), is centered at the AGN (Gallimore et al. 1996). The water masers trace densities ${n}_{H}\gt {10}^{7}\;{\mathrm{cm}}^{-3}$, and could well locate the inner edge of a molecular torus. The dust in the torus, on the other hand, is best observed in the IR continuum. The NGC 1068 nucleus has been observed several times with the mid-infrared interferometer (MIDI) in order to determine the dust distribution. López-Gonzaga et al. (2014) presented the most recent observations, which resolve the central 1–10 pc of the AGN. The observations are best fit with three gray bodies with a Gaussian intensity distribution: a $\sim 1.4\;\mathrm{pc}\times 0.5\;\mathrm{pc}$ component of 700 K, a $\sim 3\;\mathrm{pc}\times 2\;\mathrm{pc}$ component of $\sim 300\;{\rm{K}}$, and a component offset to the west (by $1.7\;\mathrm{pc}$) and north (by $5.6\;\mathrm{pc}$) with a projected size of $\sim 14\;\mathrm{pc}\times 3.5\;\mathrm{pc}$ and a temperature of $\sim 360\;{\rm{K}}$.

The CO ladder of the NGC 1068 nucleus has been observed up to $J=30$ with the Photo detector Array Camera and Spectrometer (PACS) on board the Herschel Space Observatory, and is described in HD12. Here we present a deeper observation taken around $\lambda =65.94\;\mu {\rm{m}}$, where we expect the CO(40-39) line with a rest wavelength $\lambda =65.69\;\mu {\rm{m}}$ to be at a system velocity of ${v}_{\mathrm{LSR}}=1125\;\mathrm{km}\;{{\rm{s}}}^{-1}$. The observations were taken as part of our GT Key Project SHINING; the OBSIDs are listed in Table 1. Among the high-J CO lines, CO(40-39) is particularly suitable because the spectra of bright galaxies have little contamination from other lines around $\lambda =65.69\;\mu {\rm{m}}$. Moreover, the line traces gas with densities $\gt {10}^{6}\;{\mathrm{cm}}^{-3}$ and temperatures $\gt 500\;{\rm{K}}$ (Table 2), and distinguishes a highly excited (HE) component from the gas in the CND and streamers.

Table 1.  Observed Wavelengths

OBSID	Blue Band	Red Band
	($\mu {\rm{m}}$)	($\mu {\rm{m}}$)
1342203129	65.0–70.5	130.2–141.2
1342203128	58.5–66.0	117.4–132.0
1342203127	52.0–59.5	103.8–119.4
1342203126	94.6–98.0	188.6–196.0
1342203125	91.0–95.1	181.4–190.4
1342203124	87.3–91.5	174.0–183.0
1342203123	83.2–87.7	166.0–176.0
1342203122	78.9–83.7	157.0–167.7
1342203121	74.4–79.4	148.0–159.0
1342203120	70.0–75.0	139.6–150.0
1342259623	65.4–66.5	⋯
1342259624	65.4–66.5	⋯

Notes. OBSIDs and covered wavelength of the full Herschel/PACS scan of NGC 1068. All observations for the full spectrum have been performed on observational day (OD) 460. The last two OBSIDs are from the deeper observations on CO(40-39).

Download table as:  ASCIITypeset image

Table 2.  Details of the CO(40-39) Line

Transition	Rest Wavelength	Upper Limit	ncrit	Eupper/KB
	($\mu {\rm{m}}$)	(${10}^{-17}\;{\rm{W}}\;{{\rm{m}}}^{-2}$)	(${\mathrm{cm}}^{-3}$)	(K)
CO(40-39)	65.69	2	107	$4\times {10}^{3}$

Download table as:  ASCIITypeset image

We used HIPE 11 for the data reduction. The previous data were analyzed with HIPE 6 (HD12), so we re-reduced some of the $J\leqslant 30$ lines with HIPE 11. The high-J lines detected in NGC 1068 are all consistent with the scenario that they arose in the central spatial pixel that covers the central $\sim {(9\buildrel{\prime\prime}\over{.} 4)}^{2}\;\mathrm{or}\;{(700\;\mathrm{pc})}^{2}$ of NGC 1068. Therefore, the data reduction discussed here applies to the central spaxel only. The spectra were reduced with the standard pipeline for a "chopped line scan and short range scan with background normalization," including the point source correction factor (Poglitsch et al. 2010), and setting the up-sample parameter to 1 and the over-sample parameter to 2. We found that the CO(30-29) flux should be slightly modified, so we now measure $(2\pm 1)\times {10}^{-17}\;{\rm{W}}\;{{\rm{m}}}^{-2}$ instead of $(4.2\pm 1.9)\times {10}^{-17}\;{\rm{W}}\;{{\rm{m}}}^{-2}$. There is one other small modification we make to the error bars published in HD12: the error in CO(22-21) is larger than 30% of the total flux. This line is observed at the edges of two adjacent spectral scans, and therefore the uncertainty in the flux is larger than in other lines. The modified line flux is $({4}_{-1.4}^{+4})\times {10}^{-17}\;{\rm{W}}\;{{\rm{m}}}^{-2}$.

The CO(40-39) line observations have been reduced with the same pipeline and parameter setting as described above. Figure 1 shows the region of the CO(40-39) line in the rest frame with a spectral resolution of $0.04\;\mu {\rm{m}}$. There is no clear emission at $\lambda =65.69\;\mu {\rm{m}}$, but there could be an absorption feature around $\lambda =65.6\;\mu {\rm{m}}$. Potential absorbers are ¹⁸OH ${{\rm{\Pi }}}_{3/2}-{{\rm{\Pi }}}_{3/2}\frac{{9}^{+}}{{2}}-\frac{{7}^{-}}{{2}}$ at $65.69\;\mu {\rm{m}}$, the CO(40-39) line itself, ${{\rm{H}}}_{2}{{\rm{O}}}^{+}$ ${3}_{31}-{2}_{20}\;7/2-5/2$ at $65.61\;\mu {\rm{m}}$, and NH₂ at $65.57\;\mu {\rm{m}}$. From our PACS observations of other galaxies at both 65 and $120\;\mu {\rm{m}}$ the line flux ratio of ¹⁶OH to ¹⁸OH is always found to be equal to or greater than 10 (González-Alfonso et al. 2012, 2014). Because there is no ¹⁶OH absorption visible in Figure 1 at $65.28\;\mu {\rm{m}}$, ¹⁸OH cannot be responsible for the absorption feature at $\lambda =65.6\;\mu {\rm{m}}$. Moreover, if identified with ¹⁸OH or CO(40-39), the line would be blueshifted by $400\;\mathrm{km}\;{{\rm{s}}}^{-1}$ with respect to the host galaxy, while the highest velocity features in the AGN nucleus are the water masers, with $\sim 300\;\mathrm{km}\;{{\rm{s}}}^{-1}$ (Gallimore et al. 1996).

Zoom In Zoom Out Reset image size

Absorption by NH₂ is unlikely because it should have several features of comparable intensity as indicated in Figure 1, and these are not observed. Moreover, the other NH₂ lines in NGC 1068 are faint. See Figure 2 for the full (55–200 $\mu {\rm{m}}$) PACS spectrum from the central spaxel of NGC 1068. ${{\rm{H}}}_{2}{{\rm{O}}}^{+}$ absorption cannot be excluded completely, since the $65.61\;\mu {\rm{m}}$ line is the strongest of the three lines in Figure 1, but it would be an exception in an otherwise emission dominated spectrum. A detailed analysis of the PACS spectra of NGC 1068 and other objects will be given in J. Fischer et al. (2015, in preparation).

Zoom In Zoom Out Reset image size

Therefore, we assume that the feature around $\lambda =65.6\;\mu {\rm{m}}$ is an artifact, and treat it as noise. This results in a noise level of 0.12 Jy, and a 3σ upper limit on CO(40-39) of $2\times {10}^{-17}\;{\rm{W}}\;{{\rm{m}}}^{-2}$, for a line width of $250\;\mathrm{km}\;{{\rm{s}}}^{-1}$ (HD12). Figure 1 shows the upper limit superposed on the spectrum. If the feature is masked out during the fit, the 3σ upper limit becomes $1\times {10}^{-17}\;{\rm{W}}\;{{\rm{m}}}^{-2}$, and the absorbed flux is $4\times {10}^{-17}\;{\rm{W}}\;{{\rm{m}}}^{-2}$. The feature could have absorbed the blue part of the CO(40-39) line, but it is hard to quantify by how much. Because there is no red wing visible, the artifact would have to be coincidentally just equal to that of a possible CO(40-39) emission line in order to mask it. We find this scenario unlikely and use the aforementioned upper limit of $2\times {10}^{-17}\;{\rm{W}}\;{{\rm{m}}}^{-2}$ for our analysis. The CO SLED with the upper limit and the modifications to HD12 is given in Figure 3.

Zoom In Zoom Out Reset image size

With the slightly modified CO(30-29) value and the new CO(40-39) upper limit we repeated our PDR, XDR, and shock fits with the same models as in HD12. For comparison we briefly summarize the main results from HD12 here. The observed CO SLED was analyzed with  large velocity gradient (LVG), shock, PDR, and XDR models. Only the PACS data were fitted; the $3\;\lt \;J\;\lt \;14$ CO lines that trace the colder gas were analyzed by Spinoglio et al. (2012). The PACS SLED is best explained with two components: a moderately excited (ME) component with ${T}_{\mathrm{kin}}\approx 169\;$ K, ${n}_{{{\rm{H}}}_{2}}\approx {10}^{5.6}\;{\mathrm{cm}}^{-3}$ and a total mass of ${M}_{{{\rm{H}}}_{2}}\approx {10}^{6.7}\;{M}_{\odot }$, and an HE component with ${T}_{\mathrm{kin}}\approx 571$ K, ${n}_{{{\rm{H}}}_{2}}\approx {10}^{6.4}\;{\mathrm{cm}}^{-3}$ and a total mass of ${M}_{{{\rm{H}}}_{2}}\approx {10}^{5.7}\;{M}_{\odot }$. The ME component can arise from X-ray- or shock-heated gas in the CND, while the HE component can arise from either the CND or the northern streamer. In order to determine whether the HE component arises from the CND or the northern streamer, we followed H12 and compared the PACS CO line shifts with the SINFONI observations of H₂ 1-0 S(1). The PACS HE lines are blueshifted by $59\pm 20\;\mathrm{km}\;{{\rm{s}}}^{-1}$, which best matches the $\sim 25\;\mathrm{km}\;{{\rm{s}}}^{-1}$ blueshift of the northern streamer.

To fit our new results, we use the XDR grid from Meijerink et al. (2007). This grid contains the CO line emission for gas with densities of ${n}_{H}={10}^{4}-{10}^{6.5}$ cm⁻³ and X-ray fluxes of ${F}_{{\rm{X}}}=1.6-160\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$. In HD12, the HE component was best fit with ${F}_{{\rm{X}}}=160\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$, i.e., the maximum value in the grid. This time we used an extended grid, with ${F}_{{\rm{X}}}=9-510\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ (R. Meijerink 2015, private communication).

Table 3 and Figure 3 show the parameters of the best XDR, PDR, and shock fits. Only detected lines were used during the fit, but the best fits were in all cases consistent with the upper limits of the high-J CO lines. The XDR fit of the HE component is slightly different from HD12, primarily due to the lower value for CO(30-29) and to the extension of the grid: the best fit now has an input X-ray flux ${F}_{{\rm{X}}}=510\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$, which is again the high endpoint of the grid. We tested whether an even higher X-ray input flux would better fit the HE component, using a different XDR grid (described in Section 5.1) that contains X-ray fluxes up to ${F}_{{\rm{X}}}={10}^{6}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$. For a thermal X-ray input spectrum like that in the grid from Meijerink, the CO SLEDs from the two codes agree within 20% for J $\lt 30$. In this model an input flux of ${F}_{{\rm{X}}}={10}^{3}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ could still fit the HE component, but is no improvement over the fit with ${F}_{{\rm{X}}}=510\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$.

Table 3.  Heating Mechanisms

	ME	HE	Full	${\chi }^{2}$
XDR	${n}_{H}={10}^{5.25}\;{\mathrm{cm}}^{-3}$	${n}_{H}={10}^{6.5}\;{\mathrm{cm}}^{-3}$	⋯
	${F}_{{\rm{X}}}=9\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$	${F}_{{\rm{X}}}=510\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$	⋯	1.23
	$A\approx {(108\;\mathrm{pc})}^{2}$	$A\approx {(20\;\mathrm{pc})}^{2}$	⋯
PDR	⋯	${n}_{H}={10}^{6.5}\;{\mathrm{cm}}^{-3}$	${n}_{H}={10}^{6}\;{\mathrm{cm}}^{-3}$
	⋯	${G}_{0}={10}^{5}$	${G}_{0}={10}^{5}$	2.16
	⋯	${L}_{\mathrm{FUV}}\approx 3\times {10}^{9}{L}_{\odot }$	${L}_{\mathrm{FUV}}\approx {10}^{10}{L}_{\odot }$
C-shock	${n}_{0}=2\times {10}^{5}\;{\mathrm{cm}}^{-3}$	${n}_{0}={10}^{6}\;{\mathrm{cm}}^{-3}$	⋯
	$v=20\;\mathrm{km}\;{{\rm{s}}}^{-1}$	$v=30\;\mathrm{km}\;{{\rm{s}}}^{-1}$	⋯	1.88
	$A\approx {(150\;\mathrm{pc})}^{2}$	$A\approx {(19\;\mathrm{pc})}^{2}$	⋯

Note. Parameters for the best fits with XDR, PDR, and shock models. The XDR and PDR models are from Meijerink et al. (2007), with the XDR grid extended to ${F}_{{\rm{X}}}=510\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$. Two models were employed for the shock fit: the ME component is modeled with Flower & Pineau Des Forêts (2010), and the HE component with Kaufman & Neufeld (1996; like in HD12). A is the illuminated area. The ${\chi }^{2}$ values apply to the two-component XDR and shock fit, and to the full PDR fit.

Download table as:  ASCIITypeset image

Our shock modeling uses a two-component fit: the ME C-shock model is from Flower & Pineau Des Forêts (2010), and the HE C-shock model is from Kaufman & Neufeld (1996). The lower value of the CO(30-29) line has a small effect on the best fit of a shock model to the HE component. The HE component can now be fit with a 30 km s⁻¹ velocity, instead of 40 km s⁻¹. The effect on the PDR modeling is small (using the PDR model from Meijerink et al. 2007). A PDR cannot be ruled out, but it is hard to match the $J=30$ line without overproducing the lower J lines. A single PDR component only marginally fits the entire SLED (see Figure 3). Because the kinematics of the lines suggest two components, we discard this option. The PDR fit to the HE component ($J\geqslant 19$) alone is slightly improved with the lower CO(30-29) flux, but is not as good as an XDR or shock fit.

The upper panel in Figure 3 also shows the expected CO emission from a parsec-scale molecular torus with ${N}_{H}={10}^{24}\;{\mathrm{cm}}^{-2}$ and ${f}_{\mathrm{cov}}{f}_{x44}=0.08$, as predicted by KL89, where ${f}_{\mathrm{cov}}=0.25$ is the surface covering factor of the torus and ${F}_{{\rm{X}}}={F}_{x44}\times {10}^{44}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$. This is further discussed in Section 5.2.

We conclude that the revised value for CO(30-29) and the new upper limit on CO(40-39) do not significantly change the results from HD12. They found that the ME component can arise from X-ray or shock-heated gas in the CND, while X-ray heated gas in the northern streamer is the preferred explanation for the HE component, based on the fits and the kinematics.

In the previous section we attributed the HE component of the high-J CO SLED to the northern streamer, at about 30 pc from the nucleus. If this is correct then we have not detected CO emission from the putative torus that is required by AGN unification models. In this section we want to address the following question: what does the upper limit on CO(40-39) tell us about the existence or non-existence of a torus? If there is such a torus, is it reconcilable with the upper limit, and what can be learned about its physical parameters (like size, density, distance from the center)?

In the following sections we model the molecular torus as a semi-infinite slab or as a system of such slabs, with strong irradiation from X-rays. The intrinsic X-ray luminosity ${L}_{1-100\;\mathrm{keV}}={10}^{43.5}\;\mathrm{erg}\;{{\rm{s}}}^{-1}$ indicates that input fluxes up to ${F}_{{\rm{X}}}={10}^{6}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ are needed to account for gas as close as ∼0.5 pc from the AGN. Differing from the previous section (and from HD12) we here use our own XDR code as described in Bruderer et al. (2012; including benchmark tests), with updates in Bruderer (2013). This code, originally developed for the analysis of protoplanetary disks, was extended to the higher X-ray input fluxes we need here.

5.1. XDR Grid

In the model, X-rays with photon energies between 1 and 100 keV deposit energy in the gas, which has two effects: they heat the dust and drive ionization reactions. The code iterates on radiative, thermal, and chemical processes, until all are balanced. From this balanced situation, the gas temperature and CO rotational line emission is written out as a function of distance and column density into the slab.

In a semi-infinite slab model, the observer is located on the same side as the X-ray source. X-rays enter the semi-infinite slab from one side and deposit their energy in the gas. CO is emitted in all directions, but only the CO emission toward the observer is calculated. CO can become optically thick between the point of emission and the observer (we did not assume the gas to have a large velocity gradient); in this case the attenuated CO emission is given. However, NGC 1068 is a Type-2 AGN, so the gas is located between the observer and the X-ray source. To check whether the location of the observer introduces a large uncertainty in our results, we extracted the optical depths as a function of column density from the XDR code, for all CO lines. With these optical depths, the unattenuated CO emission was calculated as a function of column density, as well as the attenuated CO emission toward an observer on the opposite side of the slab. Indeed, for some combinations of density and X-ray flux, the CO becomes very optically thick and the observed CO emission differs by more than an order of magnitude from one side of the slab to the other. For the parameters that are relevant for this paper however, the difference is at most 30% in CO line flux at column densities $\lt {10}^{25}\;{\mathrm{cm}}^{-2}$. Another potential issue is geometrical dilution; this is discussed further in Section 5.3.

The expected CO emission is calculated for a grid of semi-infinite slabs with densities of $n\in$ [10², ${10}^{2.5}$ ... 10⁹] cm⁻³ and F_X $\in$ [${10}^{-4}$, ${10}^{-3}$ ... 10⁶] $\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$. The input SED is a power law with the shape ${F}_{\nu }\propto {\nu }^{-1}$, which best matches the observed XDR spectrum of NGC 1068 (Pier et al. 1994). We also ran the grid for a softer intrinsic spectrum (spectral index of −2 instead of −1). In this case, the X-ray absorption and CO emission take place within smaller column densities than for the hard X-ray spectrum.

Line absorption by dust particles is not taken into account because this depends strongly on the geometry and viewing angle. Since we do not know the geometry of gas structures in the central few parsecs of NGC 1068, we would rather assume no line absorption by dust. This can be done safely up to a column ${N}_{H}=8.5\times {10}^{23}\;{\mathrm{cm}}^{-2}$, where the dust becomes optically thick at 66 $\mu {\rm{m}}$ (Draine 2003). At higher column densities, the results from our model are thus upper limits to the line fluxes (see Section 5.4).

In all grid points, the line emission is calculated up to a column of ${N}_{H}={10}^{26}\;{\mathrm{cm}}^{-2}$ into the slab. High-J CO emission arises from a certain combination of density n, X-ray flux ${F}_{x}$, and column density, best represented by the effective ionization parameter ${\xi }_{\mathrm{eff}}$. ${\xi }_{\mathrm{eff}}$ is calculated from formula 12 in Maloney et al. (1996); ${\xi }_{\mathrm{eff}}=4\pi \frac{{F}_{x}}{n\;{N}_{22}^{\phi +1}}$, with ${N}_{22}$ as the column density in units of ${10}^{22}\;{\mathrm{cm}}^{-2}$, $\phi =(\alpha -1)/\gamma$, in which $\alpha =-1$ is the spectral index, and $\gamma =8/3$ is a parameter that reflects that we are interested in the high-energy X-rays ($E\gt 1$ keV). The high-J CO lines (around J = 40) are well-produced for an effective ionization parameter ${\xi }_{\mathrm{eff}}\approx 0.02$. ${\xi }_{\mathrm{eff}}$ is calculated for three grid points with $n={10}^{7}\;{\mathrm{cm}}^{-3}$, and ${F}_{{\rm{X}}}={10}^{4},{10}^{5}\;\mathrm{or}\;{10}^{6}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$, and plotted in Figure 4, together with the gas temperature and the H and CO- abundances. The central grid point produces most CO(40-39) emission: at higher F_X, the gas remains hot and the CO abundance remains low, while at lower F_X the gas quickly becomes molecular and cools rapidly.

Zoom In Zoom Out Reset image size

Although H₂ is the main collision partner of CO in a molecular gas, atomic hydrogen also contributes to the collisional excitation of CO if the gas is partly atomic or ionized. Because collision rates between H and CO are currently still uncertain, they are approximated by the H₂ collision rate, scaled with the reduced mass of the colliding particles. The next section treats the CO SLED that results from our model.

5.2. An X-Ray Irradiated Slab

A model of a torus with a smooth plane-parallel slab that is intensively radiated with X-rays was already given by KL89, who concluded that most of the torus cooling happens in (very) high-J CO lines (with J up to ∼57). However, a certain column density is required to produce observable CO fluxes because the gas has a stratified structure: the first layer, through which the X-rays enter the gas, is fully ionized with temperatures of a few thousand Kelvin. The CO forms only further into the slab, and cools the gas rapidly.

In the following discussion the illuminated slab has a covering factor of 0.25 since 20%–25% of all AGNs are Compton-thick (Malizia et al. 2009; Burlon et al. 2011) and its illuminating X-ray flux depends solely on the intrinsic luminosity of NGC 1068 (${10}^{43.5}\;\mathrm{erg}\;{{\rm{s}}}^{-1}$, assuming an isotropically radiating source) and the distance toward the AGN. The absolute emitted CO line fluxes of these scenarios are subject to large uncertainties (Röllig et al. 2007), especially at high energy levels. One of the reasons for the uncertainties is that the main collision partner for CO, H₂, forms on dust grains. The H₂ abundance thus depends strongly on the dust abundance, its composition, and its temperature, which are all uncertain parameters. Moreover, the dust and gas temperatures decouple in the extreme conditions considered here, with ${T}_{\mathrm{dust}}\;\sim$ a few 100 K, and ${T}_{\mathrm{gas}}\;\sim$ a few 1000 K. So unlike in the study of protoplanetary disks, where absolute line fluxes can be given within a factor of a few (Bruderer et al. 2012), we adopt an uncertainty of 1 order of magnitude in the modeled line fluxes.

A general discussion on what part of parameter space can be constrained by the upper limit is given in Section 6. But first we analyze a simple slab model with parameters based on the water masers described in Gallimore et al. (1996). In this model, the gas has a density $n={10}^{7}\;{\mathrm{cm}}^{-3}$ and an impinging X-ray flux of ${10}^{5}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ because the water masers imply gas densities of $\geqslant {10}^{7}\;{\mathrm{cm}}^{-3}$ at 1.5 pc from the nucleus. These are the same parameters as used by KL89, and with an intrinsic X-ray flux ${F}_{{\rm{X}}}={10}^{43.5}\;\mathrm{erg}\;{{\rm{s}}}^{-1}$, they represent a torus with inner radius of 1.6 pc (our intrinsic luminosity and spectral range differs slightly from KL89). Our model prediction of the CO(40-39) flux is a factor of ∼8 lower than presented in KL89—with ${f}_{\mathrm{cov}}{f}_{x44}=0.08$ to account for a covering factor 0.25—as can be seen in the right panel of Figure 5. The reason is that KL89 assume the high-J CO lines to be optically thick up to J ∼ 57, and therefore the luminosity of the line ${L}_{J}\propto {\nu }^{3}\propto {J}^{3}$. Instead, we calculate the optical depth as function of column density into the slab, and find that the CO(40-39) line does not become optically thick within ${N}_{H}={10}^{25}\;{\mathrm{cm}}^{-2}$. For larger column densities our model converges toward the KL89 prediction (but needs ${N}_{H}\gt {10}^{26}\;{\mathrm{cm}}^{-2}$ to match it). At ${N}_{H}={10}^{24}\;{\mathrm{cm}}^{-2}$, as used in KL89, the CO(40-39) flux differs by a factor 100 between the 2 models. It clearly matters at what column density the high-J CO lines become optically thick.

Zoom In Zoom Out Reset image size

The upper limit corresponds to a column density of ${N}_{H}=1.5\times {10}^{25}\;{\mathrm{cm}}^{-2}$, which is still similar to the expected obscuring gas column density in NGC 1068 (Bauer et al. 2014). So although the CO(40-39) line traces the gas we are interested in ($n\sim {10}^{7}\;{\mathrm{cm}}^{-3}$, ${N}_{H}\sim {10}^{25}\;{\mathrm{cm}}^{-2}$, located within a few pc of the AGN), the observation did not reach the required sensitivity to constrain the gas distribution.

5.3. A Clumpy Cloud System

In the above scenario the gas is distributed smoothly. However, there are indications that the dust in a torus is distributed in clouds, and most recent torus models treat a clumpy, rather than a smooth, dust distribution (see Hönig et al. 2006 for such a model of NGC 1068).

Assuming that the gas could also be distributed in clouds, we approximate a clumpy torus with 10 plane-parallel slabs. The slabs have the same density as in our previous example and the layer closest to the source receives an X-ray input flux of ${10}^{5}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$. All other layers lie at larger distance from the source, and are obscured by the intermediate layers, as shown in Figure 6. The total column density along the line of sight to the source is ${10}^{25}\;{\mathrm{cm}}^{-2}$, and the illuminated surface of the slabs increases with distance to the source, to keep a covering factor of 0.25. The total amount of gas has therefore increased with respect to the previous example.

Zoom In Zoom Out Reset image size

The first layer receives the highest X-ray flux, while layers further away from the center receive lower X-ray fluxes due to geometrical dilution and due to X-ray absorption in the intermediate layers. The geometrical dilution is accounted for by taking the grid point with the X-ray flux that corresponds to the geometrically diluted—but unattenuated—X-rays from the central source. In case this flux lies between two grid points, the CO SLED is interpolated between these grid points. Each layer thus has a cumulative line distribution attached to it (as shown in the left panel in Figure 5), from which the line emission is extracted between two column densities. For example, the emission of the first layer could thus be extracted from ${N}_{H}=0\to {10}^{24}\;{\mathrm{cm}}^{-2}$, and that of the second from ${N}_{H}={10}^{24}\to 2\times {10}^{24}\;{\mathrm{cm}}^{-2}$. This accounts for the X-ray attenuation by intermediate layers because the XDR model calculates the X-ray attenuation as function of column density.

The 10 slabs have column densities of ${10}^{24}\;{\mathrm{cm}}^{-2}$. The distance between the slabs is determined by the volume filling factor ${{\rm{\Phi }}}_{{\rm{V}}}$: ${{\rm{\Phi }}}_{{\rm{V}}}=0.5$ implies that a 0.03 pc thick layer is followed by a 0.03 pc thick gap. Figure 7 shows the CO SLED with ${{\rm{\Phi }}}_{{\rm{V}}}\;=$ 1, 0.5, and 0.2. The three plots look very similar, there is no sign of a lower output flux or change in SLED for lower ${{\rm{\Phi }}}_{{\rm{V}}}$. The reason is that the increasing layer surface compensates the geometrical dilution (both $\propto {r}^{2}$). Only when ${{\rm{\Phi }}}_{{\rm{V}}}$ becomes very small (${{\rm{\Phi }}}_{{\rm{V}}}\sim 0.05$), the outermost layers are located so far from the X-ray source that the CO SLED changes significantly. We therefore conclude that the volume filling factor is of little importance in our case.

Zoom In Zoom Out Reset image size

It is also worth noting that the layered model with ${{\rm{\Phi }}}_{{\rm{V}}}=1$ looks very similar to the single slab model that was presented in Section 5.2. Geometrical dilution is thus not important for the grid point with $n={10}^{7}$ cm⁻³, ${F}_{{\rm{X}}}={10}^{5}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$, and ${N}_{H}={10}^{25}\;{\mathrm{cm}}^{-2}$. The same applies to the other grid points that are used in this paper.

5.4. Line Absorption by Dust

Figure 5 shows that column densities of $\sim {10}^{25}\;{\mathrm{cm}}^{-2}$ are needed to match the CO(40-39) upper limit in NGC 1068. Within such high column densities, dust absorbs significant fractions of the line fluxes, and especially affects the high-J CO lines.

For simplicity, we assume that the gas and dust follow the same distribution, and that the dust resembles galactic dust. The dust cross-sections are taken from Draine (2003) and follow $\sigma [{\mathrm{cm}}^{-2}\mathrm{per}\;{\rm{H}}-\mathrm{atom}]\propto {\lambda }^{-2}$ at the relevant wavelengths.¹⁰ In this case the position of the observer with respect to the slab does matter for the final SLED. The right panel in Figure 8 shows the CO SLED when the torus is located between the X-ray source and the observer, taking line absorption by dust into account. The modeled CO(40-39) emission decreased by a factor 10. Because we do not know the real dust distribution in NGC 1068, we do not include line absorption by dust in the remainder of this paper, but keep in mind that it can attenuate the CO(40-39) flux by a factor of 10.

Zoom In Zoom Out Reset image size

A smooth molecular torus (with $n={10}^{7}\;{\mathrm{cm}}^{-3}$ at 1.6 pc from a ${L}_{1-100\;\mathrm{keV}}={10}^{43.5}\;\mathrm{erg}\;{{\rm{s}}}^{-1}$ source and with a column density of ${N}_{H}=1.5\times {10}^{25}\;{\mathrm{cm}}^{-2}$) could have been observable in CO(40-39) within the sensitivity of the PACS observation. But the uncertainty in the model at these high-J lines is an order of magnitude. Moreover the dust could attenuate the CO(40-39) line flux by a factor of 10, so this scenario is not strongly constrained by the upper limit.

Can we exclude anything else? Is there any grid point where the CO emission exceeds the upper limit by more than an order of magnitude? Figure 9 shows the cumulative CO(40-39) emission as a function of column density and $\frac{{F}_{{\rm{X}}}}{n}$. At ${N}_{H}\gt {10}^{25}\;{\mathrm{cm}}^{-2}$ the uncertainty rises above 30% for a gas located between the X-ray source and the observer. Line absorption by dust is not taken into account. The white broken line illustrates where the observed upper limit is matched for an illuminated area of $4\pi {r}^{2}$, with r = 5, 1.6 or 0.5 pc, as given in Table 4.

Zoom In Zoom Out Reset image size

Table 4.  Matching the CO(40-39) Upper Limit

FX	n	Inner Radius	Outer Radius	Covering Factor	Required Emission
($\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ )	(cm−3)	(pc)	(pc)	(0–1)	($\mathrm{erg}\;{{\rm{s}}}^{-1}\;{\mathrm{cm}}^{-2}\;{\mathrm{sr}}^{-1}$)
106	108	0.5	0.53	4.5	1.18
105	107	1.6	1.9	0.4	0.118
104	106	5	8.2	0.06	0.0118

Note. Surface covering factors required to fit the $J=40$ upper limit face value (line absorption by dust is not taken into account). The third column in the table lists the inner radius of the torus, assuming that ${L}_{1-100\;\mathrm{keV}}={10}^{43.5}\;\mathrm{erg}\;{{\rm{s}}}^{-1}$. The outer radius follows from a column of 10²⁵ cm⁻². A 10⁸ cm⁻³ density slab located at 0.5 pc from the AGN, needs a covering factor $\gt 1$ to be observable, which is not physical. The scenarios with 10⁷ and 10⁶ cm⁻³ could be observable in CO(40-39).

Download table as:  ASCIITypeset image

In all cases, the maximum emission arises at column densities $\geqslant {10}^{25}\;{\mathrm{cm}}^{-2}$, and at $-2\lt \mathrm{log}(\frac{{F}_{{\rm{X}}}}{n})\lt -1.5$. So our initial analysis with $n={10}^{7}\;{\mathrm{cm}}^{-3}$ and ${F}_{{\rm{X}}}={10}^{5}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ already produced CO(40-39) very efficiently, and most grid points will therefore not be constrained by the upper limit. Sticking to a total column density of ${10}^{25}\;{\mathrm{cm}}^{-2}$, there are four grid points that produce observable amounts of CO(40-39) emission: those with $n={10}^{8}$ and ${F}_{{\rm{X}}}={10}^{6}$, $n={10}^{7}$ and ${F}_{{\rm{X}}}={10}^{5}$, $n={10}^{6}$ and ${F}_{{\rm{X}}}={10}^{4}$, and $n={10}^{5}\;{\mathrm{cm}}^{-3}$ and ${F}_{{\rm{X}}}={10}^{3}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$. We investigate here what sizes the gas structures with these parameters must have to match the CO(40-39) upper limit. However, the slab with $n={10}^{5}\;{\mathrm{cm}}^{-3}$ and ${F}_{{\rm{X}}}={10}^{3}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ is excluded from the analysis because the CO lines become optically thick, and the observed CO emission depends strongly on the position of the observer with respect to the slab.

Assuming a toroidal gas distribution, with ${N}_{H}={10}^{25}\;{\mathrm{cm}}^{-2}$, we calculate the surface covering factor (1 for a $4\pi \;{{\rm{r}}}^{2}$ shell covered with gas) that is needed to match the CO(40-39) upper limit. The results are given in Table 4 and are interpreted as follows.

• 1.  
  A gas layer with $n={10}^{8}\;{\mathrm{cm}}^{-3}$ at 0.5 pc from the AGN needs a covering factor of 4.5 to be observable, which is not physical. Therefore, such a compact component could not be observed in CO rotational lines at the current sensitivity, even if it has a column high enough to completely block the X-rays from our view.
• 2.  
  The slabs placed at 1.6 and 5 pc do give reasonable covering factors of 0.4 and 0.06, respectively.

MIDI observations of the hot dust in the NGC 1068 nucleus are resolved in two central components with scales of 0.5–3 pc (López-Gonzaga et al. 2014), and a third component that is offset. If the gas is distributed at scales similar to these central components, the scenario with a radius of 5 pc is not valid. This leaves only one scenario left under which the CO(40-39) could be matched: the one that was presented in Section 5.2.

If the CO(40-39) line were detected at the upper limit value, it would trace gas between 1.6 and 5 pc from the AGN, with densities $\sim {10}^{6}-{10}^{7}\;{\mathrm{cm}}^{-3}$, which covers at least 6% to 40% of a shell around the AGN. Column densities of ${10}^{25}\;{\mathrm{cm}}^{-2}$ are required. Only a few AGN satisfy this condition, and the detection of CO(40-39) in an XDR would point to very high column densities (although this line can be observed at lower column densities in shocks; Herczeg et al. 2012). Currently, there are no instruments that can perform an even deeper observation of this line, but it may be possible with future observatories such as the Space Infrared Telescope for Cosmology and Astrophysics (SPICA).

The deep observation on the CO(40-39) line in the NGC 1068 central region resulted in an upper limit of $2\times {10}^{-17}\;{\rm{W}}\;{{\rm{m}}}^{-2}$, more than 10 times lower than the previous upper limit (HD12). The upper limit does not change the previous conclusions, as to what mechanism heats the highly excited component in the CO SLED (HD12): all models (shocks, PDR, XDR) that fit the highest excited component in the observed SLED, are consistent with the new $J=40$ upper limit.

We extended an XDR code, used for protoplanetary disks, with densities up to ${10}^{9}\;{\mathrm{cm}}^{-3}$ and input fluxes up to ${F}_{{\rm{X}}}={10}^{6}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$ to investigate how far the upper limit constrains the gas distribution around the AGN. The spectral index of the modeled X-ray spectrum is −1, and the slabs have column densities of ${10}^{26}\;{\mathrm{cm}}^{-2}$. In the semi-infinite slab model, the X-ray source and the observer are located on the same side of the slab. The observed CO(40-39) flux differs by at most 30% if the observer is placed at the opposite side of the slab, as long as ${N}_{H}\leqslant {10}^{25}\;{\mathrm{cm}}^{-2}$. However, the difference increases rapidly for higher column densities. We investigated the CO emission from a gas with a density of ${10}^{7}\;{\mathrm{cm}}^{-3}$, ${F}_{{\rm{X}}}={10}^{5}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1}$, which peaks around $J=40$. The modeled CO(40-39) flux matches the upper limit for ${N}_{H}=1.5\times {10}^{25}\;{\mathrm{cm}}^{-2}$.

A line detected at an upper limit value would trace gas with densities of $\sim {10}^{6}-{10}^{7}\;{\mathrm{cm}}^{-3}$ at 1.6–5 pc from the AGN, with column densities of at least ${10}^{25}\;{\mathrm{cm}}^{-2}$, and surface covering factors of at least 6%–40%. Especially because of the high column density required to be observable within the PACS sensitivity, we expect to see only few AGNs with CO(40-39) emission.

Because tori are thought to be clumpy structures rather than a smooth medium, we approximated a clumpy distribution with 10 slabs and volume filling factors 1, 0.5, and 0.2. This did not change the CO SLED at CO(40-39) because the increasing layer surface compensates the geometrical dilution. A difference is only seen when ${{\rm{\Phi }}}_{{\rm{V}}}=0.05$.

On the other hand, line absorption by dust is a very important effect. Dust attenuates the high-J CO lines, which makes a difference of a factor of 10 in the modeled CO(40-39) line flux at a column density  of ${N}_{H}={10}^{25}\;{\mathrm{cm}}^{-2}$.

In HD12 we speculated that there might be no molecular torus, based on the non-detection of several high-J CO lines. But we did not have an XDR code that addresses all of the caveats and uncertainties that are described above. With our new model, which is adapted to the high densities and X-ray fluxes in an AGN torus, considering the position of the observer and line absorption by dust, we conclude that the non-detection of CO(40-39) is still consistent with the existence of a molecular torus.

We thank the anonymous referee for helpful comments that improved the paper. We thank Rowin Meijerink for providing us with an extended XDR grid. Basic research in infrared astronomy at NRL is funded by the US ONR; J.F. also acknowledges support from the NHSC. E.G.A is a Research Associate at the Harvard-Smithsonian Center for Astrophysics. A.S. thanks the DFG for support via German-Israeli Project Cooperation grant STE1869/1-1.GE625/15-1. S.V. acknowledges support by NASA through Herschel contracts 1427277 and 1454738. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); and IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain).