A CENSUS OF X-RAY GAS IN NGC 1068: RESULTS FROM 450 ks of CHANDRA HIGH ENERGY TRANSMISSION GRATING OBSERVATIONS

X-ray spectra demonstrate that many compact sources are viewed through partially ionized gas. This gas manifests itself as a rich array of lines and bound–free features in the 0.1–10 keV energy range. The existence of this "warm absorber" gas and the fact that it often shows Doppler shifts indicating outflows, has potential implications for the mass and energy budgets of these sources. These features have been used to infer mass outflow rates for many accreting sources, notably those in bright, nearby active galactic nuclei (AGNs; Turner & Miller 2009; Crenshaw et al. 2003; Miller et al. 2009; Krongold et al. 2003; Kaastra et al. 2011; Kaspi et al. 2002).

Intrinsic to the problem is the fact that the gas is not spherically distributed around the compact object. If it were, resonance lines would be absent or would have a P-Cygni character. P-Cygni profiles are observed in X-ray resonance lines from the X-ray binary Cir X-1 (Schulz et al. 2008); most AGNs do not show this behavior clearly. Understanding the quantity of warm absorber gas, its origin and fate, are key challenges. Interpretation of absorption spectra is often based on the assumption that the residual flux in the line trough is solely due to finite optical depth or spectral resolution. That is, the effect of filling in by emission is neglected. Understanding the possible influence of emission is desirable, both from this point of view and since emission provides complementary information about a more extended region. Comparison between warm absorbers and the spectra of objects in which the gas is viewed in reflection rather than transmission provides an added test for our understanding of the geometry of the absorbing/emitting gas. Emission spectra also are possibly less affected by systematic errors associated with the observation, such as internal background or calibration errors which could affect the residual flux in the core of deep absorption features.

Emission spectra are observed from compact objects when the direct line of sight to the central compact object is blocked. This can occur in X-ray binaries or in Seyfert galaxies where the obscuration comes from an opaque torus (Antonucci & Miller 1985). A notable example is in Seyfert 2 galaxies, and the brightest such object is NGC 1068 (Bland-Hawthorn et al. 1997). This object has been observed by every X-ray observatory with sufficient sensitivity, most recently by the grating instruments on Chandra and XMM-Newton. These reveal a rich emission line spectrum, including lines from highly ionized medium-Z elements, fluorescence from near neutral material, and radiative recombination continua (RRCs; Liedahl et al. 1990) which are indicative of recombination following photoionization.

Apparent emission can be caused by various physical mechanisms. These include cascades following recombination, electron impact excitation, inner shell fluorescence, and radiative excitation by absorption of the continuum from the central object in resonance lines. Radiative excitation produces apparent emission as a direct complement to the line features seen in warm absorbers; rates can exceed that due to other processes. It is also often referred to as resonance scattering, since it is associated with scattering of continuum photons in resonance lines. However, past treatments in the context of Seyfert galaxy X-ray spectra have been limited to consideration of a single scattering of an incident continuum photon by a resonance line, after which the photon is assumed to be lost. Since there is a considerable literature on the topic of transfer of photons in resonance lines, also sometimes referred to as resonance scattering, here and in what follows we will not use this term.

Since radiative excitation preferentially affects lines with large oscillator strengths arising from the ground term of the parent ion, it affects line ratios such as the n = 1–2 He-like lines, causing them to resemble the ratios from coronal plasmas. This was pointed out by Kinkhabwala et al. (2002), who presented high resolution spectra of NGC 1068 obtained using the reflection grating spectrometer (RGS) on the XMM-Newton satellite. The RGS is most sensitive at wavelengths greater than approximately 10 Å, and Kinkhabwala et al. (2002) showed that the He-like lines from N and O in NGC 1068 are affected by radiative excitation. They also showed that, since radiative excitation depends on pumping by continuum radiation from the central compact source, this process is affected by the attenuation of the continuum. This in turn depends on the column density of the line emitting gas; radiative excitation is suppressed when the continuum must traverse a high column density. Kinkhabwala et al. (2002) adopted a simple picture in which the gas is assumed to be of uniform ionization and opacity and were able to then constrain the column density in the NGC 1068 line emitting region from the observed line ratios.

The High Energy Transmission Grating (HETG) on the Chandra satellite is sensitive to wavelengths between 1.6 and 30 Å, allowing study of the He-like ratios from elements Ne, Mg, Si and heavier, in addition to those from O. Ogle et al. (2003) used an HETG observation of NGC 1068 to show that the He-like ratios from Mg and Si show a stronger signature of radiative excitation than do O and Ne. They interpret this as being due to the fact that all the lines come from the same emitting region with a high column, and that attenuation of the continuum is stronger for Mg and Si owing to ionization effects. That is, that O and Ne are more highly ionized than the heavier elements, and point out that this is consistent with the results of photoionization models.

The challenge of modeling X-ray emission spectra divides into parts: (1) modeling the excitation process which leads to emission, (2) modeling the ionization of the gas, and (3) summing over the spatial extent of the emitting gas leading to the observed spectrum. Kinkhabwala et al. (2002) include a detailed treatment of the population kinetics for emitting ions which takes into account the effect of the continuum radiation from the compact object on the populations, along with the attenuation of this radiation in a spatially extended region. Ogle et al. (2003) add to this a treatment of the ionization balance. Matt et al. (2004) have used the intensities of lines formed by fluorescence to infer the covering faction of low ionization, high column density gas and shown that this is consistent with the properties of the obscuring torus.

A shortcoming of efforts so far is that none self-consistently treats the spatial dependence of the absorption effects affecting radiative excitation together with the corresponding effects on the ionization balance. That is, models such as Ogle et al. (2003) adopt temperature and ionization balance associated with optically thin photoionized gas. They then assume that this ionization is constant throughout the cloud and use this to calculate continuum attenuation and its effects on the population kinetics. They did not consider the spatial dependence of the ionization of the gas associated with attenuation of the incident continuum, along with the associated suppression of radiative excitation.

Since the work of Kinkhabwala et al. (2002) and Ogle et al. (2003), there have been numerous studies of absorption dominated warm absorber spectra. In addition, a large campaign with the Chandra HETG was carried out on NGC 1068 in 2008, resulting in a data set with unprecedented statistical accuracy (Evans et al. 2010). For these reasons, and the reasons given above, we have carried out new modeling of the photoionized emission spectrum of NGC 1068. The questions we consider include: To what extent can photoionized emission models fit the observations? Which structures, known from other studies, can account for the observed X-ray line emission? Are there patterns in the elemental abundances in the photoionized emission spectrum which provide hints about the origin or fate of the X-ray gas? What mass of gas is associated with the X-ray emission, and what flow rate does this imply? In Section 2 we present the observed spectrum and in Section 3 we describe model fitting. A discussion is in Section 4.

The data set that we use in this paper consists of an approximately 400 ks observing campaign on NGC 1068 during 2008 using the Chandra HETG. A description of these observations was reported briefly in Evans et al. (2010), but no detailed description of the data has been published until now. In addition, we also incorporate the earlier 46 ks observation with the same instrument which was carried out in 2000 and was published by Ogle et al. (2003). A log of the observations is given in Table 1.

The satellite roll angles during the observations were in the range 308°–323°. This corresponds to the dispersion direction being approximately perpendicular to the extended emission seen in both optical and X-rays (Young et al. 2001). The HETG disperses the spectrum into two arms, referred to as High Energy Grating (HEG) and Medium Energy Grating (MEG). The standard HEG and MEG extraction region size is 4.8 arcsec on the sky. This is comparable to the extent of the brightest part of the X-ray image, which is approximately 6 arcsec. It is possible to extract the spectra from regions which are narrower or wider in the direction perpendicular to the dispersion, and we will discuss this below.

The spectrum was extracted using the standard Ciao tools, including detection of zero order, assigning grating orders, applying standard grade filters, gti filters, and making response files. The dispersion axes for all the observations were nearly perpendicular to the axis of the narrow line region (NLR; position angle ≃40°). As pointed out by Ogle et al. (2003) the width of the nuclear emission region in the dispersion direction for the zeroth-order image is 0.81–0.66 arcsec, which corresponds to a smearing of FWHM = 0.015–0.018 Å over the 6–22 Å range in addition to the instrumental profile (FWHM = 0.01, 0.02 Å for HEG, MEG). Events were filtered by grade according to standard filters, streak events were removed. First-order HEG and MEG spectra of the entire data set were extracted using CIAO 4.4. Positive and negative grating orders were added.

In this paper we analyze data extracted using the standard size extraction region from the HETG, which is 1.3 × 10⁻³ deg or 4.8 arcsec. Significant line emission does originate from the "NE cloud," located 3.2 arcsec from the nuclear region, and was discussed by Ogle et al. (2003). This emission is a factor 2–3 times weaker than the emission centered on the nucleus. In Section 2.2 we briefly discuss analysis of spatially resolved data, i.e., data from smaller or larger extraction regions. A more complete examination of the spatial dependence of the spectrum will be carried out in a subsequent paper.

We fit the spectrum with a model consisting of an absorbed power law continuum plus Gaussian lines. Table 2 lists 86 lines, of which 67 are distinct features detected using the criteria described below. Gaussian line fitting was carried out using an automated procedure which fits the spectrum in 2 Å intervals. These are chosen for convenience: narrow enough that the continuum can be approximated by power law and containing a small number of features. Adjacent intervals are chosen to overlap in order to avoid artifacts associated with boundaries. Within each interval we fit a local power law continuum plus Gaussians. The Gaussians are added to the model at wavelengths corresponding to known lines. For each chosen centroid wavelength the width and normalization are varied using a χ² minimization procedure until a best fit is obtained. The centroid wavelength is also allowed to vary by a moderate amount, twice the value of a fiducial thermal Doppler width, during this procedure. The trial line is considered to be detected if the χ² improves by 10, corresponding to approximately 99% confidence for three interesting parameters (Avni 1976). The fiducial thermal Doppler width is a free parameter in this procedure (but note it does not influence the final fitted width, only the searching procedure for the line center). For the results shown here we adopt a value of 200 km s⁻¹ for this parameter. Most lines are free of blending so the results of the line search are independent of the value of this parameter.

Table 2. Lines Found

Index	Wavelength	Width	Flux	Lab	Ion	Lower Level	Upper Level	voff
	(Å)	(km s−1)	(cm−2s−1)	(Å)				(km s−1)
1	1.780	⩽2400	1.3$^{+0.1}_{-0.3} \times 10^{-5}$	1.780	Fe xxvi	1s1(2S)	1s02p1(2P)	1100$^{+ 840}_{- 40}$
2	1.855	⩽2000	2.9$^{+0.2}_{-0.2} \times 10^{-5}$	1.869	Fe xxv	1s2(1S)	1s12p1(1P)	3400$^{+ 40}_{- 400}$
3	1.945	1700$^{+ 17}_{- 17}$	5.6$^{+0.2}_{-0.3} \times 10^{-5}$	1.941	Fe i	Kα		500$^{+ 40}_{- 800}$
4	2.550	⩽2500	2.7$^{+0.1}_{-1.6} \times 10^{-6}$	2.549	Ca xx	1s1(2S)	1s03p1(2P)	1100$^{+ 2400}_{- 40}$
5	3.030	1900$^{+ 550}_{- 440}$	3.9$^{+0.1}_{-1.6} \times 10^{-6}$	3.020	Ca xx	1s1(2S)	1s02p1(2P)	150$^{+ 40}_{- 1000}$
6	3.175	⩽3600	4.7$^{+0.4}_{-0.4} \times 10^{-6}$	3.150	Ar xviii	1s1(2S)	1s03p1(2P)	...
7	3.195	⩽2600	2.3$^{+0.4}_{-0.4} \times 10^{-6}$	3.180	Ca xix	1s2(1S)	1s12p1(1P)	...
8	3.195	⩽2600	2.3$^{+0.4}_{-0.4} \times 10^{-6}$	3.190	Ca xix	1s2(1S)	1s12p1(3P)	670$^{+ 3300}_{- 40}$
9	3.215	⩽2000	2.7$^{+0.4}_{-0.4} \times 10^{-6}$	3.210	Ca xix	1s2(1S)	1s12s1(3S)	670$^{+ 5200}_{- 40}$
10	3.380	⩽2200	5.5$^{+0.4}_{-0.5} \times 10^{-6}$	3.365	Ar xvii	1s2(1S)	1s13p1(1P)	...
11	3.735	⩽2200	4.0$^{+0.1}_{-1.4} \times 10^{-6}$	3.739	Ar xviii	1s1(2S)	1s02p1(2P)	1500$^{+ 840}_{- 400}$
12	3.945	⩽2000	5.1$^{+0.2}_{-1.5} \times 10^{-6}$	3.950	Ar xvii	1s2(1S)	1s12p1(1P)	1500$^{+ 1500}_{- 40}$
13	4.000	⩽2600	5.2$^{+0.5}_{-0.4} \times 10^{-6}$	3.990	Ar xvii	1s2(1S)	1s12p1(3P)	380$^{+ 1300}_{- 40}$
14	4.000	⩽2600	5.2$^{+0.5}_{-0.4} \times 10^{-6}$	3.990	S xvi	1s1(2S)	1s03p1(2P)	380$^{+ 1300}_{- 40}$
15	4.000	⩽2600	5.2$^{+0.5}_{-0.4} \times 10^{-6}$	4.010	Ar xvii	1s2.1S	1s1.2s1.3S	1900$^{+ 1300}_{- 40}$
16	4.755	⩽2100	7.6$^{+0.4}_{-0.9} \times 10^{-6}$	4.730	S xvi	1s1.2S	1s0.2p1.2P	...
17	5.040	⩽2100	1.4$^{+0.1}_{-0.1} \times 10^{-5}$	5.040	S xv	1s2(1S)	1s12p1(1P)	1100$^{+ 1500}_{- 40}$
18	5.100	⩽2000	8.6$^{+0.8}_{-0.7} \times 10^{-6}$	5.070	S xv	1s2(1S)	1s12p1(3P)	...
19	5.100	⩽2000	8.6$^{+0.8}_{-0.7} \times 10^{-6}$	5.084	Si xiii	rrc		180$^{+ 2400}_{- 160}$
20	5.225	⩽2000	8.6$^{+0.7}_{-0.6} \times 10^{-6}$	5.212	Si xiv	1s1(2S)	1s03p1(2P)	390$^{+ 40}_{- 320}$
21	5.690	1100$^{+ 120}_{- 130}$	7.0$^{+0.3}_{-1.1} \times 10^{-6}$	5.680	Si xiii	1s2(1S)	1s13p1(1P)	600$^{+ 160}_{- 40}$
22	6.195	1300$^{+ 150}_{- 27}$	2.2$^{+0.0}_{-0.2} \times 10^{-5}$	6.180	Si xiv	1s1(2S)	1s02p1(2P)	410$^{+ 280}_{- 120}$
23	6.660	1500$^{+ 120}_{- 47}$	3.1$^{+0.0}_{-0.2} \times 10^{-5}$	6.640	Si xiii	1s2(1S)	1s12p1(1P)	230$^{+ 480}_{- 40}$
24	6.750	1700$^{+ 1300}_{- 2300}$	2.4$^{+0.1}_{-0.1} \times 10^{-5}$	6.690	Si xiii	1s2(1S)	1s12p1(3P)	...
25	6.750	1700$^{+ 1300}_{- 2300}$	2.4$^{+0.1}_{-0.1} \times 10^{-5}$	6.744	Si xiii	1s2(1S)	1s12s1(3S)	870$^{+ 120}_{- 120}$
26	6.870	⩽2000	7.4$^{+0.4}_{-0.3} \times 10^{-6}$	6.860	Si x	Kα		700$^{+ 40}_{- 120}$
27	6.995	⩽3600	7.3$^{+0.1}_{-2.2} \times 10^{-6}$	6.930	Si ix	Kα		...
28	7.005	450	2.9$^{+2.5}_{-1.6} \times 10^{-7}$	7.001	Si viii	Kα		970$^{+ 5600}_{- 40}$
29	7.005	450	2.9$^{+2.5}_{-1.6} \times 10^{-7}$	7.040	Mg xi	rrc		2600$^{+ 5600}_{- 40}$
30	7.150	⩽2700	2.2$^{+0.1}_{-0.1} \times 10^{-5}$	7.110	Mg xii	1s1(2S)	1s03p1(2P)	...
31	7.320	⩽2000	6.1$^{+0.4}_{-0.3} \times 10^{-6}$	7.310	Mg xi	1s2(1S)	1s15p1(1P)	740$^{+ 640}_{- 40}$
32	7.485	1400$^{+ 14}_{- 14}$	8.1$^{+0.2}_{-1.3} \times 10^{-6}$	7.470	Mg xi	1s2(1S)	1s14p1(1P)	530$^{+ 40}_{- 240}$
33	7.770	⩽2000	9.6$^{+0.4}_{-0.5} \times 10^{-6}$	7.758	Al xii	1s2(1S)	1s13p1(1P)	640$^{+ 4400}_{- 40}$
34	7.860	1200$^{+ 200}_{- 12}$	1.2$^{+0.1}_{-0.0} \times 10^{-5}$	7.850	Mg xi	1s2(1S)	1s13p1(1P)	750$^{+ 1200}_{- 40}$
35	8.435	1400$^{+ 96}_{- 43}$	3.3$^{+0.0}_{-0.3} \times 10^{-5}$	8.420	Mg xii	1s1(2S)	1s02p1(2P)	600$^{+ 360}_{- 40}$
36	9.190	1500$^{+ 230}_{- 30}$	4.2$^{+0.0}_{-0.4} \times 10^{-5}$	9.170	Mg xi	1s2(1S)	1s12p1(1P)	480$^{+ 200}_{- 120}$
37	9.340	⩽2500	3.5$^{+0.1}_{-0.1} \times 10^{-5}$	9.230	Mg xi	1s2(1S)	1s12p1(3P)	...
38	9.340	⩽2500	3.5$^{+0.1}_{-0.1} \times 10^{-5}$	9.320	Mg xi	1s2(1S)	1s12s1(3S)	490$^{+ 40}_{- 760}$
39	9.465	⩽2000	1.9$^{+0.0}_{-0.3} \times 10^{-5}$	9.475	Fe xxi	2p2(3P0)	2p14d1(3D1)	1500$^{+ 1900}_{- 80}$
40	9.720	⩽2000	1.9$^{+0.1}_{-0.1} \times 10^{-5}$	9.679	Fe xix	2p4(1D2)	2p35d1(3D$_3\#3$)	...
41	9.720	⩽2000	1.9$^{+0.1}_{-0.1} \times 10^{-5}$	9.800	Fe xvii	rrc		3600$^{+ 160}_{- 79}$
42	10.045	⩽2400	2.5$^{+0.0}_{-0.7} \times 10^{-5}$	10.014	Fe xix	2p4(1D2)	2p35d1(3D3)	230$^{+ 40}_{- 400}$
43	10.260	⩽2000	3.9$^{+0.1}_{-0.2} \times 10^{-5}$	10.240	Ne x	1s1(2S)	1s03p1(2P)	550$^{+ 1600}_{- 40}$
44	10.365	1600$^{+ 95}_{- 120}$	3.4$^{+0.0}_{-0.7} \times 10^{-5}$	10.388	Ne ix	rrc		1800$^{+ 320}_{- 40}$
45	10.670	⩽2200	3.6$^{+0.2}_{-0.1} \times 10^{-5}$	10.641	Fe xix	2p4(3P2)	2p34d1(3P2)	320$^{+ 40}_{- 1300}$
46	10.670	⩽2200	3.6$^{+0.2}_{-0.1} \times 10^{-5}$	10.660	Fe xvii			860$^{+ 40}_{- 1300}$
47	10.830	⩽2000	2.3$^{+0.1}_{-0.1} \times 10^{-5}$	10.770	Fe xvii	2p6(1S0)	2p56d1(1P1)	...
48	10.830	⩽2000	2.3$^{+0.1}_{-0.1} \times 10^{-5}$	10.820	Fe xix			860$^{+ 240}_{- 40}$
49	11.025	1300$^{+ 300}_{- 13}$	3.5$^{+0.1}_{-0.4} \times 10^{-5}$	11.000	Ne ix	1s2(1S)	1s14p1(1P)	460$^{+ 160}_{- 80}$
50	11.555	1400$^{+ 40}_{- 70}$	4.0$^{+0.0}_{-0.6} \times 10^{-5}$	11.500	Fe xviii			...
51	11.555	1400$^{+ 40}_{- 70}$	4.0$^{+0.0}_{-0.6} \times 10^{-5}$	11.547	Ne ix	1s2(1S)	1s13p1(1P)	930$^{+ 280}_{- 40}$
52	11.800	⩽2000	4.2$^{+0.1}_{-0.2} \times 10^{-5}$	11.762	Fe xxi	2p2(3P1)	2s12p23d1(3P2)	170$^{+ 520}_{- 80}$
53	12.165	1400$^{+ 200}_{- 14}$	9.5$^{+0.1}_{-0.8} \times 10^{-5}$	12.100	Ne x	1s1(2S)	1s02p1(2P)	...
54	12.300	⩽2400	5.7$^{+0.1}_{-1.6} \times 10^{-5}$	12.264	Fe xvii	2p6(1S0)	2p54d1(3D1)	260$^{+ 40}_{- 280}$
55	12.840	⩽2000	5.5$^{+0.2}_{-0.3} \times 10^{-5}$	12.800	Fe xx			200$^{+ 1300}_{- 40}$
56	12.840	⩽2000	5.5$^{+0.2}_{-0.3} \times 10^{-5}$	12.812	Fe xviii	2p5(2P3/2)	2s12p53p1(2D5)	480$^{+ 1300}_{- 40}$
57	13.509	⩽3600	3.2$^{+0.3}_{-0.5} \times 10^{-5}$	13.447	Ne ix	1s2(1S)	1s12p1(1P)	...
58	13.509	⩽3600	3.2$^{+0.3}_{-0.5} \times 10^{-5}$	13.500	Ne ix	1s2(1S)	1s12p1(3P)	940$^{+ 6700}_{- 40}$
59	13.725	⩽3000	1.4$^{+0.0}_{-0.1} \times 10^{-4}$	13.700	Ne ix	1s2(1S)	1s12s1(3S)	590$^{+ 360}_{- 40}$
60	14.215	⩽2000	1.1$^{+0.0}_{-0.2} \times 10^{-4}$	14.206	Fe xviii	2p5(2P3/2)	2p53d1(2D$_{5/2}\#$2)	950$^{+ 960}_{- 40}$
61	14.215	⩽2000	1.1$^{+0.0}_{-0.2} \times 10^{-4}$	14.250	O viii	rrc		1900$^{+ 960}_{- 40}$
62	14.415	⩽2100	7.2$^{+0.1}_{-2.2} \times 10^{-5}$	14.394	Fe xviii	2p5(2P3/2)	2p43d1(2D5/2)	710$^{+ 40}_{- 960}$
63	14.580	⩽2000	5.5$^{+0.4}_{-0.3} \times 10^{-5}$	14.500	Fe xx	2p3(2P1/2)	2p23s1(4P1/2)	...
64	14.835	1600$^{+ 410}_{- 100}$	3.5$^{+0.1}_{-0.8} \times 10^{-5}$	14.816	Fe xix	2p4(3P1)	2p33s1(1D2)	740$^{+ 920}_{- 40}$
65	14.835	1600$^{+ 410}_{- 100}$	3.5$^{+0.1}_{-0.8} \times 10^{-5}$	14.832	O viii	1s1.2S	1s0.5p1.2P	1100$^{+ 920}_{- 40}$
66	15.040	⩽2700	1.2$^{+0.1}_{-0.0} \times 10^{-4}$	15.015	Fe xvii	2p6(1S0)	2p53d1	640$^{+ 1600}_{- 40}$
67	15.295	⩽3600	8.5$^{+0.4}_{-0.3} \times 10^{-5}$	15.188	O viii	1s1(2S)	1s04p1(2P)	...
68	15.295	⩽3600	8.5$^{+0.4}_{-0.3} \times 10^{-5}$	15.262	Fe xvii	2p6(1S0)	2p53d1(3P1)	490$^{+ 40}_{- 990}$
69	15.295	⩽3600	8.5$^{+0.4}_{-0.3} \times 10^{-5}$	15.418	Fe xvii	2p6(1S0)	2p53d1(3P2)	3500$^{+ 39}_{- 980}$
70	15.670	⩽2800	2.5$^{+0.2}_{-0.2} \times 10^{-5}$	15.627	Fe xviii	2p5(2P3/2)	2p43s1(2P5/2)	310$^{+ 200}_{- 200}$
71	16.040	1500$^{+ 250}_{- 62}$	8.4$^{+0.2}_{-1.0} \times 10^{-5}$	16.007	Fe xviii	2p5(2P3/2)	2p43s1(2P3/2)	520$^{+ 400}_{- 120}$
72	16.785	1700$^{+ 500}_{- 17}$	1.8$^{+0.0}_{-0.3} \times 10^{-4}$	16.777	Fe xvii	2p6(1S0)	2p5(3S1)	990$^{+ 40}_{- 160}$
73	16.785	1700$^{+ 500}_{- 17}$	1.8$^{+0.0}_{-0.3} \times 10^{-4}$	16.777	O vii	rrc		990$^{+ 40}_{- 160}$
74	17.135	⩽2000	1.6$^{+0.1}_{-0.1} \times 10^{-4}$	17.050	Fe xvii	2p6(1S0)	2p5(3S1)	...
75	17.420	1100$^{+ 300}_{- 58}$	3.1$^{+0.2}_{-0.4} \times 10^{-5}$	17.396	O vii	1s2(1S)	1s15p1(1P)	720$^{+ 40}_{- 200}$
76	17.785	1400$^{+ 200}_{- 160}$	5.2$^{+0.3}_{-0.5} \times 10^{-5}$	17.768	O vii	1s2(1S)	1s14p1(1P)	850$^{+ 440}_{- 40}$
77	18.635	1700$^{+ 470}_{- 34}$	1.2$^{+0.0}_{-0.2} \times 10^{-4}$	18.627	O vii	1s2(1S)	1s13p1(1P)	1000$^{+ 40}_{- 280}$
78	19.025	1200$^{+ 56}_{- 59}$	3.1$^{+0.1}_{-0.2} \times 10^{-4}$	18.968	O viii	1s1(2S)	1s02p1(2P)	240$^{+ 80}_{- 160}$
79	20.970	1500$^{+ 360}_{- 190}$	9.0$^{+0.8}_{-1.4} \times 10^{-5}$	20.910	N vii	1s1(2S)	1s03p1(2P)	280$^{+ 40}_{- 600}$
80	21.651	1300$^{+ 230}_{- 38}$	2.4$^{+0.1}_{-0.2} \times 10^{-4}$	21.602	O vii	1s2(1S)	1s12p1(1P)	460$^{+ 440}_{- 40}$
81	21.855	⩽2900	1.7$^{+0.1}_{-0.6} \times 10^{-4}$	21.804	O vii	1s2(1S)	1s12p1(3P)	440$^{+ 40}_{- 1000}$
82	22.145	560$^{+ 27}_{- 11}$	5.8$^{+0.5}_{-0.2} \times 10^{-4}$	22.110	O vii	1s2(1S)	1s12s1(3S)	660$^{+ 80}_{- 40}$
83	23.625	⩽3600	6.5$^{+0.4}_{-2.1} \times 10^{-5}$	23.440	O i	2p4(1D2)	1s12s22p5(1P1)	...
84	23.835	⩽2600	8.6$^{+0.6}_{-2.2} \times 10^{-5}$	23.771	N vi	1s2(1S)	1s14p1(1P)	330$^{+ 160}_{- 720}$
85	24.840	1300$^{+ 210}_{- 38}$	3.1$^{+0.1}_{-0.3} \times 10^{-4}$	24.781	N vii	1s1(2S)	1s02p1(2P)	420$^{+ 200}_{- 40}$
86	24.840	1300$^{+ 210}_{- 38}$	3.1$^{+0.1}_{-0.3} \times 10^{-4}$	25.164	Ar xv	2s2(1S0)	2s13p1(3P1)	5000$^{+ 200}_{- 39}$

Download table as:  ASCIITypeset images: 1 2

Also included in Table 2 are tentative identifications for the lines, along with laboratory wavelengths. These come from the xstar database (Kallman & Bautista 2001; Bautista & Kallman 2001). We treat as potential identifications the (likely) strongest feature which lies within a wavelength range corresponding to a Doppler shift of ⩽+5000 km s⁻¹ from a known resonance line or RRC. Table 2 includes the values of the implied Doppler shift for these identifications. Note that our detection criterion does not guarantee accurate measurement of the other line parameters, i.e., centroid wavelength and width, so that weaker lines do not all have bounds on these quantities. Where they are given in Table 2, errors on these quantities and the line flux represent 90% confidence (Δχ² ≃ 3 for one parameter; Avni 1976) limits.

The NGC 1068 spectrum together with the model fits discussed in the next section are shown in Figures 1–4. These are plotted in 2 Å intervals and separated according to grating arm (HEG or MEG). We do not display wavelength regions where the grating arm has no sensitivity or where there are no features. No rebinning or grouping is performed in these plots or in our fitting procedure. We do quote values of χ² in Section 3, and these were calculated using binned data. Also, when plotting the longer wavelength regions, we bin the spectrum is binned for plotting purposes only in order to avoid the dominance of the noisy bins with low counts at long wavelengths (i.e., ⩾18 Å). In addition, the iron K region, 1.5–2.5 Å is not plotted here, and is discussed in Section 2.4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Notable line features in the spectrum include the 1s–np lines of H- and He-like ions of O, Ne, Mg, Si and S. Corresponding features are also present from Ar and Ca, though not clearly detected. The 21–23 Å wavelength region contains the He-like lines of oxygen. These are discussed in more detail in Section 2.3 below. The 17–19 Å wavelength region contains the Lα line from H-like line of O viii, the most prominent line in the spectrum. In addition, lines from Fe xvii near 16.8 and 17.1 Å, and n = 1–3, 1–4 and 1–5 lines from O vii near 17.4, 17.8 and 18.6 Å, respectively. These lines are discussed in Section 2.5. The 15–17 Å wavelength region contains the Lβ lines from O viii plus lines from ion stages of iron: Fe xvii–xix. The n = 1–2 lines from He-like Ne are at 13.5–13.8 Å. The 11–13 Å wavelength region contains the Lα lines from Ne ix at 12.16 Å, the higher order n = 1–3 and n = 1–4 lines from Ne ix at 11.6 and 11.0 Å, respectively, and additional lines from Fe xvii–xix. Also apparent between 11.15 and 11.35 Å are inner shell fluorescence lines from L shell ions of Mg: Mg v–Mg x. The 9–11 Å wavelength region contains the Lβ line from Ne x plus the H- and He-like lines from Mg near 8.4 and 9.15–9.35 Å, respectively. The RRC of Ne ix is near 10.35 Å. The 5–7 Å wavelength region contains the n = 1–2 lines from both H-like and He-like Si, near 6.2 Å and 6.7 Å, respectively. The Mg xii Lβ line is apparent at 7.15 Å. Inner shell fluorescence lines from L shell ions of Si: Si vii–Si xii are between 6.85 and 7.05. Lines due to H- and He-like S are at 4.74 and 5.05 Å, respectively.

The fluxes in Table 2 can be compared with those given by previous studies of the NGC 1068 X-ray spectrum by Kinkhabwala et al. (2002) and Ogle et al. (2003). There is incomplete overlap between the lists of detected lines by us and these authors. They detect 40 and 60 lines, respectively, which can be compared with 86 in our Table 2. Of these, some of the lines detected using the XMM-Newton RGS by Kinkhabwala et al. (2002) are outside of the spectral range of the Chandra HETG. There are additional discrepancies in the significance and identification of a small number of weak lines when compared with Ogle et al. (2003), though none of these is highly statistically significant in either our spectra or theirs. Since the Ogle et al. (2003) study was based on a much shorter observation than ours, it is not surprising that we detect a larger number of lines. We find general consistency between the fluxes in Table 2 and those of Ogle et al. (2003). There are significant discrepancies in the comparison between our line fluxes and those of Kinkhabwala et al. (2002), generally in the sense that the Kinkhabwala et al. (2002) are greater than ours by factors of several, possibly due to the different characteristics of the XMM-Newton telescope and the spatial extent of the NGC 1068 X-ray emission which allows flux from a more extended region to enter the spectrum.

2.1. Line Wavelengths

Figure 5 and Table 2 show the Doppler shifts of the line centroids for the strongest lines in Table 2. The shift is given in velocity units and is measured relative to the rest frame of the AGN. Laboratory wavelengths are taken from the xstar database (Bautista & Kallman 2001); most are from NIST.⁵ Error bars are from the errors provided by the automated Gaussian fitting procedure described in the previous section. Wavelengths are only plotted if the errors on the wavelength are less than 2000 km s⁻¹. Most wavelengths are consistent with a Doppler shift in the range 400–600 km s⁻¹. Exceptions correspond to RRCs, e.g., Ne ix 10.36 Å, for which the tabulated wavelengths are less precise, and weak lines such as O viii 1s–5p which is blended with Fe xix near 14.8 Å. This figure shows that there is evidence for wavelength shifts which are larger than the errors on the line centroid determination. Examples of this include the O viii Lα line compared with the 1s²–1s2s³S (f) line of O vii.

Zoom In Zoom Out Reset image size

Figure 5, lower panel, and Table 2 show the Gaussian widths of the lines for the strongest lines. The width is given in velocity units, i.e., σ_εc/ε where σ_ε is the Gaussian FWHM in energy units and ε is the line energy. Error bars are from the confidence errors provided by our automated fitting procedure as described above. Wavelengths are only plotted in Figure 5 if the errors on the wavelength are less than 2000 km s⁻¹. This shows that there is no single value for the width which is consistent with all the lines. Examples of broader lines include the O viii Lα line, and contrast with the forbidden line f of O vii, for which the width is smaller. The FWHM which fits to most lines is approximately 1400 km s⁻¹. This corresponds to a Gaussian sigma of ≃600 km s⁻¹ which we adopt in our numerical models discussed below.

2.2. Spatial Dependence

The discussion so far has utilized data extracted using the standard size extraction region from the HETG, which has width 1.3 × 10⁻³ deg or 4.8 arcsec. We have also explored possible spatial dependence of the spectrum in the direction perpendicular to the dispersion direction by extracting the data using regions which are half and double the angular size, i.e., 2.4 and 9.6 arcsec wide. We then carry the fitting steps described above. Figure 6 illustrates the positions of these regions superimposed on the zero-order image.

Zoom In Zoom Out Reset image size

In this way we can study the influence of adding successive regions further from the dispersion axis to the spectrum. Figure 7 shows the ratios of the line fluxes in the standard extraction region (green hexagons) and the extraction double the standard size (red diamonds) compared to the line flux in the extraction region half the standard width. This shows that the line emission outside half width region is significant; the values of the ratios of essentially all the lines are greater than unity, and many are ∼2. Plus, the emission from outside the standard 4.8 arcsec extraction region is not negligible; the values of the ratio for this region (red symbols) are on average greater than the values of the ratio for the standard region. The principal difference between standard region and the double size region is the importance of the radiative excitation emission, which is larger in the standard region. This is consistent with results from Ogle et al. (2003). For example, for the He-like lines of Ne ix, the G ratio is larger in the standard region spectrum than in the spectrum from the double size region. For the O vii lines, the G ratios are similar for the two regions, but the R ratio is lower in the standard region, implying higher density there.

Zoom In Zoom Out Reset image size

These results demonstrate that most of the line emission is contained within the standard extraction region, 4.8 arcsec wide, and most line ratios are unaffected when emission from outside the standard extraction region is included. It is worth noting that the physical length scale corresponding to this angular size is approximately 1 arcsec = 72 pc (Bland-Hawthorn et al. 1997), so the standard extraction region corresponds to a total size of 347 pc, or a maximum distance from the central source of approximately 174 pc. This is much larger than the region where the broad line region and the obscuring torus are likely to lie, which is ≃2–3 pc (Jaffe et al. 2004). For the remainder of this paper we will discuss primarily analysis of data from the standard extraction region.

2.3. He-like Lines

The He-like lines in Table 2 provide sensitive diagnostics of emission conditions: excitation mechanism, density, and ionization balance. These challenges of fitting these lines are apparent from the region containing the H and He-like lines from the elements O, Ne, Mg, Si, and S. This shows the characteristic three lines from the n = 2–n = 1 decay of these ions: the resonance (r), intercombination (i) and forbidden (f). These lines provide useful diagnostics of excitation and density. These are described in terms of the ratios R = f/i and G = (f + i)/r (Gabriel & Jordan 1969a, 1969b). R is indicative of density, since the energy splitting between the upper levels of the two lines (2³S and 2³P) is small compared with the typical gas temperature and collisions can transfer ions from the 2³S to the 2³P when the density is above a critical value. R can also be affected by radiative excitation from the 1s2s³S to the 1s2s³P levels, but this requires photon field intensities which are greater than what is anticipated for NGC 1068. In our modeling we include this process, using an extension of the global non-thermal power law continuum, and find that it is negligible. We do not consider the possibility of enhanced continuum at the energy of the 2³S to the 2³P provided by, e.g., hot stars. G is an indicator of temperature or other mechanism responsible for populating the various n = 2 levels. At high temperature, or when the upper level of the r line (2¹P) is excited in some other way, then G can have a value ⩽1; in the absence of such conditions, the 2³S and 2³P levels are populated preferentially by recombination and G has a characteristic large value, G ⩾ 4. Table 2 shows that the various elements have common values of R, in the range 1–3, indicating moderate density. More interesting is the fact that there is a diversity of values for G: O and Ne have values 2.5–3 while the other elements all have smaller values, G ⩽ 1.5.

The G ratio is most affected by the process responsible for excitation of the r line. This can be either electron impact collisions or radiative excitation. Collisions are less likely in the case of NGC 1068 due to the presence of RRCs from H-like and He-like species in the observed spectra (Kinkhabwala et al. 2002). This indicates the presence of fully stripped and H-like ions, which produce the RRCs, but at a temperature which is low enough that the RRCs appear narrow. Typical RRC widths in NGC 1068 correspond to temperatures which are ⩽10⁵ K. A coronal plasma, in which electron impact collisions are dominant, would require a temperature ⩾10⁶ K in order to produce the same ions. Radiative excitation can produce G ratios as small as 0.1. Thus, the G ratios for NGC 1068 indicate the importance of radiative excitation for Mg and Si, while O and Ne are dominated by recombination.

Radiative excitation depends on the presence of strong unattenuated continuum from the central source in order to excite the r line. A consequence of the large cross section for radiative excitation in the line is the fact that the line will saturate faster than continuum. At large column densities from the source, the radiation field is depleted in photons capable of exciting the r line, while still having photons capable of ionizing or exciting other, weaker lines. Thus, the G ratio is an indicator of column density. He-like lines formed after the continuum has traversed a large column density will have large G ratios, corresponding to primarily recombination. He-like lines formed after the continuum has traversed a small column density will have small G ratios, corresponding to radiative excitation (Kinkhabwala et al. 2002).

Figures 8 and 9 display the G ratio for the elements O, Ne, Mg, and Si from the NGC 1068 HETG observations. These are shown with error bars as the red crosses in the centers of the various frames. This shows that the ratios differ between the elements: O and Ne have G ≃ 2–3, while Si and Mg have G ≃ 1. This differing behavior suggests that the He-like O and Ne lines are emitted by recombination, while the He-like lines of Mg and Si are emitted by radiative excitation. Since radiative excitation is suppressed by large resonance line optical depths, this suggests small optical depths in Si and Mg and larger optical depths in the lines from other elements. The differences in optical depths are significant; as shown by Kinkhabwala et al. (2002) equivalent hydrogen columns of 10²³ cm⁻² or more are required to suppress radiative excitation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Also plotted in Figures 8 and 9 are the values of these ratios produced by photoionization models consisting of a single slab of gas of given ionization parameter at the illuminated face and given column density. The solid curves correspond to constant total slab column density, where black = 10^23.5, blue = 10^22.5, green = 10^21.5 and red = 10^20.5. The dashed curves correspond to constant ionization parameter in the range 1 ⩽ log (ξ) ⩽ 3. Here the ionization parameter is ξ = 4πF/n where F is the ionizing energy flux in the 1–1000 Ry range and n is the gas number density. The role of radiative excitation is apparent from the fact that the low column density models (red curves) produce smaller G values than the high column density models (black curves).

The results in Figures 8 and 9 show that no single value of the photoionization model parameters can simultaneously account for all the ratios. O and Ne require log(ξ) ≃ 1 and column ⩾10²³cm⁻², while Si requires log(ξ) ≃ 2 and column ⩽10²² cm⁻² and Mg is intermediate in both ξ and column.

Also plotted in Figures 8 and 9 are values for certain other ratios, in the same units as for the He-like lines, and also including the model values. These include the higher series allowed lines from some He-like ions, e.g., 1–4/1–2 versus 1–3/1–2 for O vii, where 1–2 includes the r, i, and f lines while 1–4 and 1–3 include only the resonance component. These ratios depend primarily on ionization parameter and only very weakly on column density. This is because recombination cascades make larger values of both ratios than radiative excitation, so recombination tends to dominate production of these lines.

2.4. Iron K Line Region

The iron K lines from NGC 1068 have been the subject of previous studies by Matt et al. (2004) and Iwasawa et al. (1997), who showed that the line consists of three components, corresponding to near neutral, and H- and He-like ion species. The HETG provides higher spectral resolution than these previous studies, but fewer total counts. Figure 10 shows the spectrum in the region between 1.5 and 2.5 Å and reveals the three components of the iron K lines. The most prominent are a feature at 6.37 keV (1.95 Å), consistent with neutral or near-neutral iron, a feature at 6.67 (1.85 Å) indicative of He-like iron, and a feature near 7 keV (1.75 Å) which may be associated with a combination of the K absorption edge from near-neutral gas and emission from H-like iron. Matt et al. (2004) also separately identify emission from Be-like iron; we are not able to make such an identification. We do find emission at other nearby wavelengths, but this is suggestive of a broad component or blended emission from various other species.

Zoom In Zoom Out Reset image size

We model the Fe K line region using "analytic" models which utilize the xstar database and subroutines. These include optically thin models which include radiative excitation (which we denote "scatemis") and models which do not include this process ("photemis"), plus power law continuum. Both these models are available for use as "analytic" models in xspec from the xstar site.⁶ Radiative excitation cannot excite the neutral line, and we model it as emission from gas at log(ξ) = −3. The emission component with radiative excitation has log(ξ) = 3, corresponding to highly ionized gas in which the H- and He-like stages are the most abundant ionization stages of Fe. In the HETG data the He-like triplet is not fully resolved. Nonetheless, the centroid of the He-like component is closer to the wavelength of the Fe xxv r line than it is to the wavelength of the f or i lines, so we find a better fit with a resonance excitation model. The fit is shown in Figure 10. The best-fit parameters are summarized in Table 3; we find χ² = 103 for 206 degrees of freedom for the 1–2.5 Å (5–12.4 keV) spectral region. The flux in the neutral-like line is 3.8$^+_-0.44 \times 10^{-5}$ cm⁻² s⁻¹.

Table 3. Model for Fe K Region Parameters

Component	log(ξ)	Norm	Width	Velocity
			(km s−1)	(km s−1)
Photemis	−3	2.47$^{+0.82}_{-0.61} \times 10^8$	1.28$^{+0.25}_{-0.23}\times 10^3$	2.41$^{+0.51}_{-0.67} \times 10^2$
Scatemis	3	33.1$^{+8.7}_{-9.7}$	1.28$^{+0.25}_{-0.23}\times 10^3$	2.41$^{+0.51}_{-0.67} \times 10^2$

Download table as:  ASCIITypeset image

We also point out what appears to be unresolved emission in the wavelength region near ≃1.9 Å, between the He-like line and the neutral-like line. We have not attempted to model this. It may be due to a broadened component of the modeled lines, or an unresolved blend from other high ionization stages of iron. These could include L shell ions between Li-like and Ne-like. The two ionization parameters which we use to fit the iron K lines are both higher and lower than the values used to fit the remainder of the emission spectrum, described in the next section. The high ionization parameter, log(ξ) = 3, may be an extension of the ionization parameter range needed to fit the lines in the energy band below ∼5 keV. The low ionization parameter, log(ξ) = −3, likely corresponds to the torus. Neither of these components produces significant emission below ∼5 keV, so we will not discuss them further in the model fits described in the following section.

The iron K line complex contains more flux than any other line in the spectrum. This is consistent with the large normalization for the low ionization component that is required to fit the spectrum. This was discussed by Krolik & Kallman (1987) and Nandra (2006). We have not attempted to fit for a Compton shoulder on the on the red wing of the neutral-like iron line, although such a component may be present in Figure 10 (Matt et al. 2004). We also test for the presence of Ni Kα in our spectrum and find an upper limit of 1.5 × 10⁻⁵ cm⁻² s⁻¹. This is less than the flux for Ni Kα claimed by Matt et al. (2004), which was 5.6$^{+1.8}_{-1.0} \times 10^{-5}$ cm⁻² s⁻¹. One possible explanation could be that their data was extracted from a much larger spatial region, ≃40 arcsec, and is also affected by a nearby point source.

2.5. Iron L Lines

2.6. Fluorescence Lines

A goal of interpreting the X-ray spectrum of NGC 1068 is to understand the distribution of gas: its location relative to the central continuum source, its velocity, density, element abundances, temperature, ionization state. Fitting of models which predict these quantities to the data are the most effective way to derive this information, although there is a range of levels of detail and physical realism for doing this. We have examined this through the use of photoionization table models in xspec, which are described in the Appendix.

Based on what was presented in Section 2 we can summarize some of the criteria for a model for the NGC 1068 HETG spectrum: (1) the model must account for the H- and He-like line strengths from abundant elements from N to Ca (we exclude the iron region from this discussion); (2) the model must account for the G ratios from the elements O, Ne, Mg, Si, which have G ≃ 3 for O and Ne but G ≃ 1 for Mg and Si; (3) the model must account for 2p–3d/2p–3s ratio in Fe xvii which is indicative of recombination rather than radiative excitation; (4) the model must account for the strengths of the RRCs, which have strengths comparable to the corresponding resonance lines for O vii. From criterion (1) it is clear that the emitting gas must have a range of ionization parameter; no single ionization parameter can provide these ions. From criterion (2) it is likely that a range of gas column densities are needed. The hypothesis that the lower ion column densities associated with lower abundance elements (Mg and Si) can account for the diversity in G ratios can be tested in this way. Criterion (4) implies that some of the gas must have high column in order to have a detectable contribution from recombination.

Our fit uses three components with varying ionization parameters; these are spread evenly between log(ξ) = 1 and log(ξ) = 2.6. The HETG spectrum is relatively insensitive to gas outside this range of ionization parameters, except for the iron K lines. Given the strength of both the neutral-like and the H- and He-like iron K lines, it is plausible that there is an approximately continuous distribution of ionization parameters which extends beyond the range considered here. The three ionization parameter components in our fit correspond crudely to the regions dominating the emission for the elements of varying nuclear charge: the low ionization parameter component dominates the emission for N and O, the intermediate component dominates the Fe L lines and Ne and Mg, and the high ionization parameter component dominates for Mg and Si. For each ionization parameter, we include models with two column densities: 3 × 10²² cm⁻² and 3 × 10²³ cm⁻². In this way the models span the likely ionization parameter and column density for which the strongest lines we observe are emitted. We use a constant turbulent line width (sigma) of 600 km s⁻¹ which adequately accounts for the widths the majority of features. Results of our fit are summarized in Table 4.

The red curves in Figures 1–4 show that this model fits the strengths of almost all the strong features in the spectrum. This fit adopts a net redshift for the line emitting gas of 0.0023, corresponding to an outflow velocity of 450 km s⁻¹. This is marginally less than typical speeds from UV and optical lines of ⩽600 km s⁻¹ (both blueshifted and redshifted) (Crenshaw & Kraemer 2000a), but is consistent with the velocities plotted in Figure 5.

Free parameters of the fits include the elemental abundances and the normalizations for the six model components. Experimentation shows that the fit is not substantially improved when all the abundances are allowed to vary independently, so we force the abundances of N, Ne, Mg, Si, S, Ar and Ca to be the same for all the models. We allow the abundances of O and Fe to vary independently for all the models. Our best fit is shown in Figures 1–4 has χ²/ν = 2.34 for 2561 degrees of freedom. More results are shown in Figures 11 and 12, which show the elemental abundances from the various model components, and the masses of the various components (defined below) versus ionization parameter.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to account for the diversity in the importance of radiative excitation between various ions and elements, in the low ionization parameter log(ξ) = 1 component, which is responsible for most of the O vii emission, the oxygen abundance in the low column density component (component 4 in Table 4) is smaller than the oxygen abundance in the high column component (component 1 in Table 4) by a factor 2.5. If not, the recombination emission into O vii would be stronger than observed. Iron abundances are similar between the low and high column density components, and are highest (∼10) at low ionization parameter (green in Figure 12) and lowest (∼1) at high ionization parameter log(ξ) = 2.6. The normalizations and masses of the low column density components (components 4–6) is generally lower than for the high column components (components 1–3). The normalization for component 5 is consistent with zero, though we include it for completeness. Component 6 is needed to provide photoexcitation-dominated lines of Fe xix–xxi.

None of the abundances is significantly subsolar, except for O at log(ξ) = 2.6 and column = 3 × 10²² cm⁻² (component 6). This is necessary to avoid over-producing O viii RRCs. It is worth noting that X-ray observations are not generally capable of determining abundances relative to H or He, but rather only abundances relative to other abundant metal elements such as O or Fe. Thus, an apparent underabundance of O, for example, is equivalent to an overabundance of other elements such as Si, S, and Fe, relative to O.

The highest ionization parameter component is responsible for the Si and Mg emission. Although its mass is dominated by the 3 × 10²³ cm⁻² component (component 3 in Table 4), the abundances of these elements are lower than for lower Z elements, so the optical depths in the He-like resonance lines are small enough that radiative excitation dominates the He-like lines. It is notable that the differences in the importance of radiative excitation between the various He-like ions, and also Fe xvii, can essentially all be accounted for by the effects of varying amounts of resonance line optical depth which is produced solely as a result of ionization and elemental abundance effects. This is consistent with the hypothesis suggested by Ogle et al. (2003). An exception is for the lines from Ne ix, for which a large ion column density is needed in order to suppress the radiative excitation. Models with column less than 3 × 10²² cm⁻² produce G ⩽ 1.6, while the observations give $G\simeq 4^{+6}_{-2}$. Our model produces the Ne ix f line with approximately half the observed strength, while producing r and i lines which are close to the observed strengths. The Ne ix RRC is also too weak in our models. We have found, through experimentation, that an additional component with pure Ne and a column of 3 × 10²⁴ cm⁻³, and this is able to provide the observed Ne recombination features, but we have not included this component in the fits in Figures 1–4 or Table 4.

The table models are calculated assuming that the emitting gas fills a spherical shell defined by an inner radius R_i = (L/nξ)^1/2 and outer radius R_o defined by $\int _{R_i}^{R_o}n(R)dR=N$. Table 4 provides the normalizations κ_i of the table models as derived from the xspec fits. This normalization can be interpreted in terms of the properties of the source as follows: $\kappa =f ({L_{38}}/{D_{{\rm kpc}}^2})$ where L₃₈ is the source ionizing luminosity integrated over the 1–1000 Ry energy range in units of 10³⁸ erg s⁻¹, and D_kpc is the distance to the source in kpc. The normalization also includes a filling factor f which accounts for the possibility that the gas does not fill the solid angle, and f ⩽ 1. We assume the density n scales as n = n_i(R/R_i)⁻². With this assumption the ionization parameter is independent of position for a given component and the total amount of gas in each emission component j, is

$$\begin{equation} M_j=4\pi f_j m_H n_iR_i^2(R_o-R_i). \end{equation} \tag{ 1 }$$

Similarly, the emission measure for each component can be written:

$$\begin{equation} {\rm EM}_j=4\pi f_i n_i^2 R_i^4 \left(\frac{1}{R_i}-\frac{1}{R_o}\right). \end{equation} \tag{ 2 }$$

Since $n_iR_i^2=L/\xi$ both the mass and the emission measure depend on the physical size of the emission region, but not directly on density. They also depend on the ionizing luminosity. Table 4 includes values for f_i, M_i and EM_i derived using these expressions, assuming L = 10⁴⁴ erg s⁻¹. The Chandra HETG spectrum of NGC 1068 constrains the density to be ⩽10¹⁰ cm⁻³, from the strength of the O vii forbidden line. In the results given in Table 4 we have chosen the inner and outer radii of the emission region to reflect simple conclusions from our data: the inner radius is R_i = 1 pc, based on physical arguments about the location of the obscuring torus (Krolik & Begelman 1988), and the outer radius is chosen to be R_o = 200 pc in order to reflect the observed extent of the X-ray emission, and this is roughly consistent with the extraction region width. The density is chosen such that the inner radius of the line emitting gas is at 1 pc for all three emission components. Total masses and emission measures are also given. The total emission measure is log(EM) = 66.3. This is greater than would be required to simply emit the observed lines if the gas had an emissivity optimized for the maximum line emission, owing to the fact that the table models include the gas needed to create the column density responsible for shielding the photoionization dominated gas from the effects of photoexcitation. Furthermore, most of the mass and emission measure comes from hydrogen, while the Chandra HETG spectrum constrains only metals. Values for mass and emission measure are calculated using the elemental abundances given in Table 4.

It is also apparent that the assumptions described so far are not fully self-consistent. That is, if density ∝R⁻² and if R_o ≫ R_i then the column density is ≃  n_iR_i, and an inner radius of 1 pc cannot produce an arbitrary column density and ionization parameter. A self consistent model with density ∝R⁻² requires different inner radii for the different ionization parameter and column density components, and results in a total mass which is greater than that given in Table 4 by approximately a factor of 12. The emission measure is insensitive to these assumptions. Our subsequent discussion is based on the simple scenario given in Table 4: R_i = 1 pc, R_o = 200 pc.

Figures 1–4 show general consistency between most of the lines and a single choice for outflow velocity. This corresponds to a net redshift of 0.0023$^+_-$0.0002, or an outflow velocity of 450$^+_-$50 km s⁻¹ relative to the systemic redshift of 0.00383 from Bland-Hawthorn et al. (1997). There are apparent discrepancies between some model line centroids and observed features when this redshift is adopted. An example is Mg xii Lβ, which has an apparent central wavelength of 7.150 Å, as shown in Table 2. The laboratory wavelength of this line is 7.110 Å (Drake 1971), which would correspond to an outflow velocity of 1688 km s⁻¹. This is greater than other Doppler shifts in the spectrum, and the difference cannot be easily explained from uncertainties in the rest wavelength. Other possible line identifications include Kα lines from Si ions such as Si iii or Si iv. These would be associated with gas at lower ionization parameter than the majority of the other lines in the spectrum and we have not attempted to include them in our model.

One result of global model fits is revealing the presence of features which are not obviously apparent from examination or Gaussian fits to the spectrum. These are features, primarily from Fe L shell ions, which are in the model, but which are blended or are not sufficiently statistically significant in the spectrum to require their inclusion in Gaussian fits. Therefore, these lines are not identified in Table 2. Some are indicative of recombination versus radiative excitation, in a sense similar to the Fe xvii lines discussed above. Examples can be seen from examination of Figures 1–4: RRCs from Fe xvii and xviii at 9.15 and 8.55 Å, respectively, and Fe lines at 10.55 (Fe xvii 2p–5d 10.523 Å in the lab), 10.62, 10.85 (Fe xix 2p–4d 10.6323, 10.8267 Å), 11.05 (Fe xvii 2s–4p 11.043 Å), 11.18, 11.28 (Fe xvii 2p–5d, 11.133, 11.253 Å) 11.45 (Fe xviii 2p–5s, 11.42 Å) 12.15 (Fe xvii 2p–4d 12.12 Å), 12.95 (Fe xix 2s–3p, 12.92 Å) 13.49, 13.53 (Fe xix 2p–3d, 13.46, 13.5249 Å) 14.3 (Fe xviii 2p–3d, 14.208 Å) 15.88 (Fe xviii 2p–3s, 15.83 Å).

It is possible to use the observed spatial extent of the X-ray emission to make further inferences about the X-ray emitting gas. For NGC 1068, the distance is such that 1 arcsec = 72 pc (Bland-Hawthorn et al. 1997). The brightness profile in the zero-order image is such that a significant amount of flux is coming from outside an extraction region which is half the standard size, i.e., 2.4 arcsec. Thus, an approximate, but convenient, size scale for the extent of the observed X-rays is ≃2.4 arcsec ≃ 200 pc. If so, we can compare with the sizes derived from Table 4. These estimates are based on an assumed ionizing luminosity for the nucleus of 10⁴⁴ erg s⁻¹; this is uncertain by a factor of a few, but is bounded by L_bol = 4 × 10⁴⁴ erg s⁻¹ (Bland-Hawthorn et al. 1997).

The mass of material we infer based on the Chandra HETG X-ray spectrum is ≃ 3.7 × 10⁵ M_☉, though this depends on the assumptions about the inner and outer radii of the emission region and on the assumption that density scales with R⁻². This can be compared with typical masses for the entire NLR, which are likely to be ∼10⁶ M_☉, although gas with ionization parameters in the range log(ξ) 1–3 will not radiate efficiently in typical optical/UV narrow lines. If this gas flows uniformly over a distance ≃200 pc at the speed we derive, then the flow timescale is 4.4 × 10⁵ yr, and the mean mass loss rate is 0.3 M_☉ yr⁻¹. This is considerably greater than the mass flux needed to power the nucleus, which is ≃1.7 × 10⁻² M_☉ yr⁻¹ (η/0.1)⁻¹ L/(10⁴⁴ erg s⁻¹) where η is the efficiency of conversion of accreted mass into continuum luminosity.

The filling factors in Table 4 can be interpreted formally as the fraction of the total available volume in the spherical shell bounded by R_i and R_o which is filled with line emitting gas. The angular distribution of this gas is not constrained, so this can be interpreted as a conical region with fractional solid angle given by f or as a spherically symmetric distribution in which the fractional covering of each component is given by f when averaged over solid angle. In either case it is clear that the total covering fraction of the reprocessing gas is at most a few percent. This suggests that, if NGC 1068 were viewed from an arbitrarily chosen angle, and if the line of sight to the central continuum source were not blocked by the Compton thick torus, then the probability of observing a warm absorber similar to those seen in many Seyfert 1 galaxies would be small. If the fractional solid angle of the obscuring torus as seen from the nucleus is ∼50% of the total, then the probability of seeing a warm absorber would be ∼10%. This estimate is inversely proportional to the value we have assumed for the ionizing luminosity of the nucleus in NGC 1068, which is 10⁴⁴ erg s⁻¹. A smaller value for this quantity, which is possible depending on the true spectral shape, would increase the inferred probability of seeing a warm absorber toward a value ∼50%.

Our model fits to the NGC 1068 HETG spectrum allow a simple phenomenological test of the scenario which has been widely used to interpret the spectra of Seyfert 1 objects. That is, Seyfert 1 X-ray spectra, which show prominent blueshifted line absorption, have been analyzed assuming that the absorber lies solely along the line of sight and neglecting any effect of emission filling in the lines (McKernan et al. 2007). This can only be exactly correct if the absorber subtends negligible solid angle as seen from the central source. This would in turn conflict with the apparent presence of absorption in ∼half of known Seyfert 1 objects (Reynolds 1997). On the other hand, examination of the fits shown in Figures 1–4 and the results in Table 2 show a typical emission line flux for a strong line we measure is ≃2 × 10⁻⁴ s⁻¹ cm⁻² for O viii Lα. This can be compared with the amount of energy absorbed in the same line from a Seyfert 1 galaxy; for NGC 3783 this quantity is ≃5 × 10⁻⁴ (Kaspi et al. 2002). Correcting for the difference in distance, this would imply that, if the same emission line gas seen in NGC 1068 is also present in NGC 3783, then the ratio of line flux absorbed from the observed spectrum to emission is ∼8. That is, the apparent flux absorbed in the O viii Lα line, and likely other lines as well, is actually greater than observed by ∼10%–15%. The O viii Lα line in NGC 3783 shows signs of this with an apparent weak P-Cygni emission on the red edge of the trough, though this is not apparent in most other lines in that spectrum. Although O viii Lα is one of the strongest lines in the NGC 1068 spectrum, it is possible that other lines in Seyfert 1 absorption spectra are more affected by filling in from emission. Accurate fitting of these spectra, used to derive mass outflow rates and other quantities, would need to account for this process.

It is interesting to compare our results with those from UV and optical imaging and spectroscopy with the Faint Object Spectrograph and the Space Telescope Imaging Spectrograph (STIS) on the Hubble Space Telescope (HST). These properties have been discussed by Kraemer et al. (1998), Crenshaw & Kraemer (2000a), Kraemer & Crenshaw (2000a), Crenshaw & Kraemer (2000b), and Kraemer & Crenshaw (2000b). The strongest UV emission lines observed by the HST instruments, e.g., C iv λ1550, Si iv λ1398, and the He ii and hydrogen lines, comes from gas with a lower ionization parameter than the inferred from X-ray emitting gas. STIS provides spatially resolved kinematic information, which shows that the UV gas speed increases from the nucleus out to 130 pc, then the speed decreases at larger distances. Also seen are coronal lines from S xii, Fe xiv, Ne v, Ne iv. These ions can exist at the same ionization parameter as the lowest ionization X-ray gas. The intensity as a function of distance from the nucleus shows apparent dilution as $r^{_-2}$, and the emission region appears conical. The continuum spectrum most closely resembles a nuclear power law, implying that it is scattered light from the obscured nucleus. The contribution from stars to the continuum is smaller than the non-thermal power law.

The gas observed by the HETG extends to ∼200 pc from the nucleus, while the obscuration has a much smaller size ∼1 pc. The gas which extends beyond the torus may originate near the torus, resembling a warm absorber flow in a type 1 object, and flowing ballistically to larger distances. If so, the outflow speed would be expected to be constant or decrease with distance owing to gravitational forces, and the density would be subject to purely geometric dilution. Alternatively, the X-ray and UV gas could entrain gas from narrow line clouds, or gas evaporated from the narrow line clouds, or from some other source. Radiation pressure could play a role in energizing the flow. The X-ray gas could play a role in the apparent deceleration of the UV gas seen at ∼100 pc. There is little strong evidence for these latter scenarios from the HETG spectra, since we see no evidence for large changes in the ionization balance or speed of the gas in the spectra extracted from various width regions. Our results appear most nearly consistent with simple geometric dilution of the outflow, leading to constant ionization parameter with distance, and nearly constant outflow speed.

The results presented in this paper can be summarized as follows: (1) no single ionization parameter and excitation mechanism can fit to all the lines in the NGC 1068 X-ray spectrum obtained with the Chandra HETG. The table model fits require three components at ionization parameters ranging from log(ξ) = 1 to log(ξ) = 3. We are able to adequately fit all the strong identified lines in the spectrum with the exception of the lines from He-like Ne, which show a stronger recombination component than are produced by our models. (2) The abundances required by all the models are approximately solar (Grevesse et al. 1996) or slightly greater. A notable exception is a large suppression of the abundance of oxygen at the highest ionization parameter. Otherwise the O viii RRCs would be stronger than observed. (3) The masses and emission measures of gas are greater than expected for optically thin emission, owing primarily to the need for shielding to provide the recombination dominated gas seen in the lines of Mg and Si. (4) The combined constraints on the ionization parameter and column density constrain the location of the emitting gas relative to the continuum source, while the line strengths constrain the amount of emitting gas. Taken together, these require that the emitting gas have a volume filling factor less than unity; the values differ for the various components and range up to ∼0.01. The mass flux through the region included in the HETG extraction region is approximately 0.3 M_☉ yr⁻¹ assuming ordered flow at the speed characterizing the line widths. (5) Limited experimentation with extracting spectra from various positions on the sky does not reveal a clear pattern which separates the various emitting components spatially. It appears that the emitting components coexist in the same physical region.

We make use of calculations using the xstar (Kallman & Bautista 2001) modeling package.⁷ xstar is freely available and distributed as part of the heasoft package. Models can be imported into xspec and other fitting packages as tables, or via the "analytic" model warmabs/photemis/scatemis.

Our models are based on the assumption that the most plausible energy source for the X-ray lines and RRCs observed from NGC 1068 is reprocessing of the continuum from the innermost regions of the AGN. The continuum is presumably associated with the accretion disk and related structures close to the black hole, i.e., within a few R_G ≃ 3 × 10¹¹M₆ cm, where M₆ is the mass of the black hole in units of 10⁶ M_☉. The line and recombination emission we observe is likely formed at distances >10⁶R_G, based on the geometry of the obscuring torus as indicated by high spatial resolution IR imaging (Jaffe et al. 2004). If it is assumed that the line emission is from gas where the heating, excitation and ionization are dominated by the continuum from the black hole, and that there is a local time steady balance between these processes and their inverses (radiative cooling, radiative decay and recombination), then the temperature, atomic level populations and ion fractions can be calculated. xstar carries out such a calculation and also calculates the associated X-ray emission and absorption. These emissivities and opacities are then applied to a simple one-dimensional solution of the equation of radiative transfer to derive the spectrum at a distant observer. Under the simple assumption that the gas is optically thin the important physical quantities depend most sensitively on the ratio of the ionizing X-ray flux to the gas density. We adopt the definition of this ionization parameter which is ξ = 4πF_X/n where F_X is the ionizing (energy) flux between 1 and 1000 Ry and n is the gas number density. An important part of our results concerns the fact that the gas is not optically thin; this is key for explaining the varying amounts of radiative excitation seen in the line ratios. This in turn means that the column density, or the ion column densities, are also free parameters. We implement the models by direct fitting, in which xstar models are used to generate synthetic spectra which are compared to the observations within the fitting program xspec.

xstar explicitly calculates the populations of all atomic levels associated with emission or absorption of radiation; it does not rely on the traditional "nebular approximation" which assumes that every excitation decays only to ground. Therefore it is straightforward to include radiative excitation, and this has been done. Interactive, iterative calculation of full xstar models within xspec and simultaneous fitting to data is not practical owing to computational limitations.

We point out, parenthetically, that the simplest way to use xstar results within xspec is to use the associated xspec "analytic" models. When called from xspec, these read stored tables of ionic level populations as a function of ionization parameter, and then calculate the opacity and emissivity "on the fly." The physical quantity used the xspec model fit is an emitted flux or a transmission coefficient as a function of energy for a photoionized slab of given column density, abundances and ionization parameter. These are calculated from the opacity and emissivity by assuming these quantities are uniform throughout. Advantages of this procedure compared with the use of table models (described below) include: Ability to account for arbitrary element abundances, arbitrary spectral resolution, and arbitrary turbulent broadening. Limitations include the fact that it uses a saved file of level populations calculated for a grid of optically thin models for a fixed choice of ionizing spectrum rather than calculating the ionization balance self-consistently. It implicitly assumes that the absorber has uniform ionization even if the user specifies a large column, and that all emission freely escape. In this sense it is not self-consistent. In fact, analytic models fail when applied to the analysis of the NGC 1068 spectrum because, for column densities large enough to suppress radiative excitation, as demanded by, e.g., the Fe xvii 17 Å line, the associated RRCs for Fe xvii are assumed to escape freely and are greatly over predicted. For this reason we will not discuss analytic models further in this paper.

In order to self-consistently include the effects of absorption of the incident continuum in models for photoionized emission from NGC 1068 it is necessary to use optically thick models. These are calculated using xstar, by performing multiple model runs and using varying ionization parameter and cloud column density. The resulting emission spectra and transmission functions are stored as fits tables, binned in energy, using the table format prescribed by xspec. These can be read by xspec, selected according to parameter values, and interpolation between modeled parameter values is performed. xspec can fit the table models to the observed spectrum and thereby yield values for the cloud column density, ionization parameter, and the normalization, which is described in more detail below. Owing to the fact that the spectra are stored binned in energy the line width cannot be conveniently treated as a free parameter, and therefore it is important that the energy grid spacing used in constructing the table be no greater than the intrinsic instrumental resolution. It is also possible to treat the element abundances as free parameters, by assuming that the escaping line and RRC fluxes depend linearly on the abundances. This ignores the coupling of the cloud ionization structure to abundance, via the temperature and the opacity of the model. In our fits we use this approximation when deriving element abundances.

Within a given optically thick model, xstar calculates the effect of attenuation of the incident continuum using a simple single stream treatment of the radiative transfer. This takes into account the varying ionization and excitation through the model. The inclusion of radiative excitation necessitates a finer spatial grid than for continuum absorption alone, since the transfer must resolve the attenuation length of the photons in the resonance lines. Also, line opacity varies over a large range in a very narrow energy range and it is not feasible to accurately sample the relevant energies for all of the many lines in a photoionized plasma. We adopt a very simple treatment in which the line opacity is binned into continuum bins and used to calculate the attenuated flux. This flux in continuum bins is then used in the calculation of photoexcitation. The continuum bin size is typically Δε/ε = 1.4 × 10⁻⁴ which corresponds to ≃ 41 km s⁻¹ Doppler width. This is comparable to the thermal Doppler width for, e.g., oxygen at a temperature ∼10⁶ K. Thus, we are not fully resolving most lines relevant to this study. The tables used to fit to the Chandra HETG spectrum of NGC 1068 include emission from both the illuminated and unilluminated cloud faces, corresponding to the assumption that the line emitting material is arranged with approximate spherical symmetry around the central continuum source. They are calculated assuming an ionizing spectrum which is a single power law with (photon number) spectral index Γ = 2 and constant density 10⁴ cm⁻³. These results are not expected to depend sensitively on the latter two assumptions; in particular Γ values in the range 1.7–2.3 produce very similar ionization balance distributions.

xstar makes use of atomic data compiled from various sources and described by Bautista & Kallman (2001). This has been extensively updated since that time, as described by Witthoeft et al. (2007, 2009, 2011), Palmeri et al. (2002, 2003a, 2003b, 2008a, 2008b, 2011, 2012), Mendoza et al. (2004), García et al. (2009), and Bautista et al. (2003, 2004). The results presented here which flow from the model calculations are dependent on the assumptions, computational implementation and atomic data used by the models. The atomic data affects the line identifications and the outflow speeds, via the line rest wavelengths. It also affects the ionization balance, via the photoionization and recombination rates, and the line strengths. The quantitative uncertainty associated with the atomic data and the resulting uncertainty in the synthetic spectrum are difficult to estimate; this is a topic of interest for many problems beyond this one. We can point out that the emissivities of most of the lines from the H- and He-like ions in our study depend on two processes: recombination and photoexcitation. These depend in turn on photoionization cross sections and line oscillator strengths, respectively. These are quantities associated with radiative processes in relatively simple ions and therefore are generally more reliable than collisional rate coefficients which are important for ions with partially filled L shells in our models and in coronal plasmas. Foster et al. (2010) have shown that uncertainties of, say, 20% on the collision strengths for O vii can result in a range of a factor of ≃3 in the inferred temperature in a coronal plasma based on the G ratio. We do not anticipate similar sensitivity to atomic data uncertainties in the NGC 1068 models, since the rates for recombination and photoexcitation do not have the strong non-linear dependence on temperature which is intrinsic to collisional rates. In particular, we do not consider the rates embodied in our models to be sufficiently uncertain to change our results qualitatively, or to allow significantly better fits with fewer components or parameters. A more likely shortcoming of our models is the neglect of some important physical mechanism. Examples include: charge transfer; photo-excitation by some unseen source of radiation; non-ionization equilibrium effects; or inhomogeneities in elemental compositions associated with formation or destruction of dust.