© ESO 2021

1. Introduction

The unification model for active galactic nuclei (AGN; Miller & Antonucci 1983; Antonucci & Miller 1985) is able to explain the paradigm that some spectroscopic observations of galaxies display broad emission lines, while others do not. The reason, quite simply, has to do with the orientation and our line-of-sight (LOS) towards the AGN of interest. How we perceive to view the source is the cause of this observational difference, whereby type 1 AGN are viewed from ‘face on’, allowing for both narrow and broad line emission to be seen, whereas type 2 AGN are seen from ‘side on’, so only narrow emission is detected. In type 2 AGN, there is some obscuring material (known as the dusty torus) blocking our LOS to the nucleus, meaning emission from the broad line region (BLR) is invisible; only emission from the narrow line region (NLR) is seen. However, broad Balmer emission lines in polarised light were observed within NGC 1068. They were explained as being due to photons from the BLR scattered off free electrons or possibly dust, and, therefore, polarised, into our LOS (Antonucci & Miller 1985).

In addition, the covering factor of the dusty torus must also be considered in the unification model, as determining between type 1 or type 2 AGN also depends on the probability that a photon can escape the obscuring material (Ramos Almeida & Ricci 2017). This suggests that the larger the covering factor, the more clumps are present within the torus, resulting in a large optical depth such that the photons are less likely to escape. Therefore, the probability of a photon escaping depends on the size of the torus (given by the covering factor), increasing the likelihood of the AGN being obscured (Ramos Almeida & Ricci 2017).

A further consequence of an obscuring torus means the soft X-ray spectra of type 2 AGN do not show absorption features, but they rather display strong emission features above an absorbed continuum. This can be explained by the ionisation cone model (e.g. Kinkhabwala et al. 2002), suggesting that the outflowing wind, which absorbs so much of the X-ray flux seen in Seyfert 1 AGN, re-emits the X-rays through all angles relative to the direction of outflow, producing the emission features that are seen in the spectra of type 2 AGN. In other words, the emission seen in Seyfert 2 AGN is the re-emission of the plasma from the warm absorber (WA) wind in Seyfert 1 AGN (e.g. Kaastra et al. 2002, 2012; Behar et al. 2017; Mehdipour et al. 2018), from the evidence that the ionic column densities of the emission lines were consistent with the values of the WAs in Seyfert 1 AGN (Kinkhabwala et al. 2002).

NGC 1068 is a heavily studied type 2 Seyfert galaxy, which supports perfectly the unification (Antonucci & Miller 1985) and the ionised cone (Kinkhabwala et al. 2002) models, as only narrow emission lines are observed in its soft X-ray spectra. NGC 1068 is at a redshift of z = 0.0038 (Huchra et al. 1999), with a super massive black hole (SMBH) mass of M_BH = 1.6 × 10⁷ M_⊙ (e.g. Panessa et al. 2006).

The very first high resolution X-ray spectra of NGC 1068 were taken by XMM-Newton (Kinkhabwala et al. 2002) and Chandra (Brinkman et al. 2002; Ogle et al. 2003), in the early 2000s, both showing many narrow emission lines originating from photoionised equilibrium (PIE) plasma. Kinkhabwala et al. (2002) found that excess emission in the resonance lines was explained by photoexcitation, rather than emission from collisionally ionised plasma. They suggested that the addition of plasma in collisionally ionised equilibrium (CIE) could explain the emission excess in the lower order resonance lines, but would enhance the other lines in the spectrum, in particular over predicting the intercombination and forbidden lines of the He-like triplets. Kraemer et al. (2015) reanalysed the 2000 RGS spectrum with a two component model, and found the emission lines were blueshifted by ∼160 km s⁻¹, consistent with the optical [O III] λ5007 results (Crenshaw et al. 2010). Both Brinkman et al. (2002) and Kraemer et al. (2015) suggested that the optical and X-ray emission line regions are part of the same outflowing wind, because they have consistent velocities; this is a general property of type 2 Seyfert AGN (Bianchi et al. 2006).

Due to the spatial resolution of Chandra, Brinkman et al. (2002) and Ogle et al. (2003) were able to differentiate between two emission regions within the nucleus of NGC 1068, made up of a central, primary region and an off-centre, secondary region. The secondary region is a bright X-ray source that coincides with bright radio and optical emission, indicating a correlation between the X-ray and optical regions, which may arise from old jet material (Wilson & Ulvestad 1987; Young et al. 2001). The secondary region was found to have an ionic column density three times smaller compared to the primary region, and plasma diagnostics were unable to conclude if the emission was produced in PIE or CIE plasma (Brinkman et al. 2002).

The hard (2−10 keV, and above) X-ray spectrum for NGC 1068 is also very complex (see e.g. Guainazzi et al. 1999; Bianchi et al. 2001) whereby the emission features can be explained by multiple components (e.g. Bianchi et al. 2001; Matt 2002; Matt et al. 2004; Pounds & Vaughan 2006; Bauer et al. 2015). From NuSTAR observations, three reflection components were necessary to model the Fe Kα line and Compton hump with column densities of 10²⁹, 1.5 × 10²⁷, and 5 × 10²⁸ m⁻². The third component bridged the gap between the cold and warm reflectors which dominated the Fe Kα emission line, and the other two components produced the Compton hump (Bauer et al. 2015).

The Atacama Large Millimeter Array (ALMA) observations of NGC 1068 have enabled mapping molecular tracers to study both the obscuring torus, circumnuclear disk (CND) and star burst ring (SBR; e.g. García-Burillo et al. 2014, 2016, 2017; Viti et al. 2014; Imanishi et al. 2018). The outflow rate within the CND was found to be an order of magnitude larger than the star formation rate, suggesting the ionised gas outflow is AGN driven (García-Burillo et al. 2014). The SBR is colder and less dense than the CND (Viti et al. 2014) as it is much further out from the black hole (r_CND <  200 pc, r_SBR ∼ 1.3 kpc; García-Burillo et al. 2017). From ALMA observations, it was found that the torus, with gas mass $M_{\mathrm{gas}}^{\mathrm{torus}}\sim 1\times {10}^5{M}_{\odot }$ and a distance of R_torus ∼ 3.5 pc from the black hole, is inhomogeneous, with a strong non-circular molecular motion, possibly inclined at larger radii (García-Burillo et al. 2016).

With infrared observations, hydrogen emission structures have been observed to infall towards the central SMBH, where one molecular cloud is close enough (∼1 pc) such that a tidally distributed, optically thick stream may possibly makeup the outer part of a clumpy, dusty torus (Müller Sánchez et al. 2009). Evidence for this clumpy structure, in the high energy X-rays, was found during an unveiling event in August 2014 when a flux excess was observed above 20 keV, compared to December 2012 and February 2015 (Marinucci et al. 2016). This was caused by a decrease in column density by about ΔN_H ∼ 2.5 × 10²⁴ cm⁻², which allowed Marinucci et al. (2016) to infer an intrinsic 2−10 keV luminosity of L = 7 × 10⁴³ erg s⁻¹ for the central, obscured AGN.

The goal of this work is to attain a self-consistent best fit photoionisation model for the 2000 and 2014 XMM-Newton observations of NGC 1068 in order to achieve an in depth understanding of the emitting plasma regions within this AGN. We wish to determine if there are any changes regarding the emitting plasma between the epochs, to then obtain other properties which can be derived from our photoionisation modelling. In this paper, we analyse both the reflection grating spectrometer (RGS; den Herder et al. 2001) and European photon imaging camera pn-CCD (EPIC-PN; Strüder et al. 2001) spectra simultaneously from each observation in great detail, with a sophisticated and self consistent photoionisation model.

The layout of this paper is as follows. In Sect. 2 we describe how we reduced the data, while in Sect. 3 we explain how we set up the ionising SED, and how this differs from the observed continuum. Section 4 outlines how we construct the model to analyse the spectra, before presenting our spectral fitting results in Sect. 5. We compare our findings from both epochs in Sect. 6, and then examine the locations and thermal properties of the emission line regions in Sect. 7. Finally, we present our conclusions in Sect. 8.

2. Data reduction and spectral modelling

NGC 1068 was observed twice in 2000 and four times in 2014 with XMM-Newton; the observation log is displayed in Table 1. We reduced the RGS and EPIC-PN data from each observation of NGC 1068 via the SAS (v 17.0.0) pipeline. The RGS data were reduced with RGSPROC and the EPIC-PN data were reduced using the EPPROC command. The 2000 observations were taken during the same XMM-Newton orbit so we stacked both the RGS and EPIC-PN spectra (separately) using RGSCOMBINE and EPICSPECCOMBINE, respectively. For the 2014 observations, separated by different time scales (see Table 1), there was little variability in the RGS spectra, so all four observations were combined. However, due to multiple flares (too close together to filter out the background noise) in the final 2014 observation (February 2015) we were only able to use and combine the first three (July and August 2014) EPIC-PN spectra; only the RGS data for the final 2014 observation were used (see Table 1).

At visual inspection, the Fe Kα and higher energy emission lines in the 2014 EPIC-PN spectrum appear to be shifted with respect to the 2000 observation. The top panel of Fig. 1 compares the 2000 data (red crosses) with the 2014 data (blue crosses) for the Fe Kα line (at around 6.4 keV), where there is an energy difference of ΔE = 50 eV between the observed line peaks. This blueshift (and broadening; σ(E)∼20 eV relative to Chandra observations from 2012) of the Fe Kα line was also found in NGC 5548, where the problem was attributed to uncertainties in the long-term charge transfer inefficiency and instrumental gain correction of the EPIC-PN CCDs, causing modest but significant variations (Cappi et al. 2016, and references within). Similar to Cappi et al. (2016), we compared the 2014 PN spectrum to the 2000 and 2014 MOS spectra, as well as to some 2012 Chandra observations, and found that all Fe Kα lines have energies consistent with 6.4 keV, except for the 2014 PN Fe Kα line. Therefore, this blueshifted iron line is not a new and interesting result, but rather an artefact of the long life of the PN instrument.

Fig. 1. Comparing the Fe Kα line in 2000 (red crosses) and 2014 (blue crosses). Top: there is a clear shift in the Fe Kα line between these two epochs, denoted by the red and blue dashed vertical lines at the peak energies of the observed data. The magenta arrow signifies this energy shift between the observations, equivalent to ΔE = 50 eV. This shift is due to uncertainties in the long-term charge transfer inefficiency and instrumental gain (Cappi et al. 2016). Bottom: we account for this shift manually by multiplying the energies (PI) of the event lists in each 2014 observation by 0.993. Now the 2014 PN data are consistent with the 2000 data and the correction is sufficient for modelling.

To compensate for this energy shift, we manually changed the gain correction of the PN instrument for the 2014 observations. Taking the event list files for each of the three observations in 2014 (the latter was not used due to background flaring), we multiplied the energy (PI) column by 0.993, before restacking the spectra. This procedure corrected the emission line energies in the 2014 PN data to be consistent with the 2000 data (see bottom panel of Fig. 1).

Figure 2 displays the RGS (left) and PN (right) spectra of NGC 1068, in the observed reference frame, comparing both epochs, where the main emission lines are labelled with their respective ions. The spectra in Fig. 2 show very little to no variability between the two observations (as noted by Marinucci et al. 2016). In order to conclusively determine the emission feature properties and their origins across the 0.35−10 keV energy range, we fitted both RGS and PN spectra simultaneously in each epoch (see Fig. 3), thus allowing us to obtain better continuum constraints compared to using just RGS or PN data individually. Although the EPIC-PN data could be used to search for subtle variations in the soft X-ray band, because of their higher count rate, we found possible pile-up in the 2000 data. There is a flux loss around 10% in the 2000 PN spectra, and a spectral distortion around 2% (estimated as in Jethwa et al. 2015). However, as already noted by Matt et al. (2004), these effects are negligible above 4 keV, while they become important at lower energies. As a result of this, and because there are almost no features in the two PN spectra below 2 keV, we ignored the PN spectra below 1.7 keV. This also reduced the weighted fitting bias as the PN continuum has higher flux counts at lower energies compared to the RGS spectrum.

Fig. 2. Comparing the RGS spectra (left) and EPIC-PN (right) spectra of NGC 1068 in 2000 (red) and 2014 (blue), in the observed reference frame. All spectra show a striking similarity in shape, emission features and flux, suggesting no change between the two epochs.

Fig. 3. Combined RGS and PN spectra for the 2000 (blue and red) and 2014 (green and purple) observations, in the observed reference frame. We simultaneously fit both RGS and PN spectra for each epoch to model all emission features and observed continuum over the 0.35−10 keV energy range.

As the same continuum model goes through both RGS and PN data points, covering the whole X-ray band, it is acceptable to have no or little overlap between RGS and PN data. We wanted the high-resolution spectral features to be fitted well with RGS and the fitting not be affected by the statistically dominant PN data. This has been done in previous papers whereby the PN data were excluded in the RGS range (see e.g. Kaastra et al. 2014; Mehdipour et al. 2017). Therefore, the RGS spectrum covers the energy range between 0.35 and 1.7 keV (∼7−37 Å) and the PN spectrum accounts for the 1.7−10 keV range. We binned the EPIC-PN data using the optimal method as described by Kaastra & Bleeker (2016), whereas the RGS data were binned by a factor of 3. This corresponds to an average bin size of 0.04 Å, which still over-samples the RGS resolution element of ∼0.07 Å FWHM. This is shown in Figs. 2 (right side) and 3.

We modelled the combined spectra of NGC 1068 using SPEX (v 3.04.00; Kaastra et al. 1996). In our model, we set the redshift of this AGN to z = 0.0038 (Huchra et al. 1999) and assumed solar abundances from Lodders et al. (2009). We use the Cash statistic (C-statistic, hereafter; Cash 1979; Kaastra 2017) for statistical significance, with errors at a 1σ confidence.

3. Spectral energy distributions

3.1. Ionising continuum

Due to the nature of NGC 1068, the LOS obscuration from the torus means we are unable to observe the central ionising source directly; this is the origin of the continuum that ionises the outflowing wind which we see either absorbing (Seyfert 1) or emitting (Seyfert 2) in the X-ray energy band (see the cone model in e.g. Kinkhabwala et al. 2002). Therefore, as we cannot directly observe the ionising continuum, we had to assume a broadband continuum model, which we fixed throughout the spectral fitting. It is the shape of the ionising spectral energy distribution (SED) that is important, not the normalisation (Mehdipour et al. 2016), so it is not unreasonable to use a SED from a type 1 Seyfert galaxy. In this case, we adopted the SED of NGC 7469 derived by Mehdipour et al. (2018) as it is a bright and unobscured Seyfert 1 AGN, but we adjusted some of the parameters to match the observed values found for NGC 1068 in the literature (e.g. Bauer et al. 2015). The SED shape of NGC 7469 is a good representation of a typical AGN SED when we do not know the actual ionising SED of NGC 1068, and is a step better than a simple power-law. Mehdipour et al. (2016) compared photoionisation results using the SEDs of NGC 5548 (both unobscured and obscured) and a power-law SED. We compared the SEDs and thermal stability S-curves of NGC 5548 (Mehdipour et al. 2016) and NGC 7469 (Mehdipour et al. 2018), and found that within the X-ray ionisation region where we were probing (log ξ = 1−4), the S-curves did not look significantly different. From this, we are confident that the differences between typical SEDs of the same type (like NGC 7469 and NGC 5548) would not significantly affect our results here. Therefore, NGC 7469 is a reasonable SED to use when modelling the ionising continuum of NGC 1068.

The three main components that contribute to the ionising SED are displayed in Table 2; they are as follows: Firstly, a power-law (POW) was used to model the hard X-ray emission produced when optical/UV disk photons are upscattered in a hot electron corona. To this power-law we applied upper and lower cut-off energies. The lower cut-off energy was fixed at the disk temperature, adopted from NGC 7469 (Mehdipour et al. 2018), to stop the power-law energies becoming lower than the disk photon energies; the high energy cut-off was fixed at 500 keV (Bauer et al. 2015). The next SED component was a reflection component (REFL; Magdziarz & Zdziarski 1995) which modelled the Fe Kα line and Compton hump. The iron abundance (A_Fe = 2.45¹), inclination (θ_inc = 60°) and photon index (Γ = 2.1, coupled to the photon index of POW) were fixed at the values found from previous observations of NGC 1068 (e.g. Matt et al. 2004; Pounds & Vaughan 2006; Bauer et al. 2015). Here we assumed that the X-rays we observed (after being reflected) were the same reflected X-rays that ionise the emitting plasma. This means that we had one REFL component that accounted for both ionising and observed X-ray reflection, so we fitted the normalisation parameter of REFL (N_refl) in the model. This was the only fitted parameter in the SED, all other parameters were fixed to the values in Table 2. Finally, a warm, optically thick medium has been found to best fit the soft X-ray excess and explain the UV/optical emission from the disk in Seyfert 1 AGN (e.g. Petrucci et al. 2013; Middei et al. 2018), which we modelled with the Comptonisation component (COMT; Titarchuk 1994). All the values for COMT were adopted and fixed from the ionising SED in NGC 7469 (Mehdipour et al. 2018).

3.2. Observed continuum

The X-ray spectrum we do see is the result of the ionising continuum reflecting off the circumnucleur material within or surrounding the nucleus of NGC 1068 (e.g. scattering off the dusty torus, or off free electrons; Antonucci & Miller 1985), interacting with and ionising the outflowing plasma wind. This observed continuum was modelled with a reflection component (same REFL as the ionising SED; see above) for the hard X-rays, with the addition of a simple power-law component for the soft X-ray band, for completeness. As we are only interested in the photoionised emission lines within NGC 1068, an ad hoc continuum, in the form of a simple observed power-law plus reflection, is justified.

4. Data analysis

In both 2000 and 2014 spectra, we modelled the observed continuum by fixing both the power-law (POW) and reflection (REFL) parameters to the values in Table 2. However, we did fit the normalisation parameters of both POW (N_pow) and REFL (N_refl; same as the ionising reflection component), and the photon index (Γ) for the power-law, which we initially set to 2.1; Γ in REFL was coupled and fixed at this value (Matt et al. 2004; Pounds & Vaughan 2006; Bauer et al. 2015). We also coupled a VGAU component to REFL to obtain a measurement for the broadness (σ_v) of the Fe Kα line, measured as a velocity.

The Galaxy absorption was modelled with the HOT component, setting the total hydrogen column density to $N_{\mathrm{H}}^{\mathrm{Gal}}=3.34\times {10}^{24}$ m⁻², made up from both H I (Kalberla et al. 2005) and H₂ (Wakker 2006). For a neutral Galactic gas, we fixed the temperature and turbulent velocity to T_Gal = 0.5 eV and ${\mathit{v}}_{\mathrm{turb}}^{\mathrm{Gal}}=5.62\pm 3.19$ km s⁻¹ (for details of obtaining this turbulent velocity result, see the appendix in Grafton-Waters et al. 2020), respectively. All fitted components were redshifted (z = 0.0038; Huchra et al. 1999) and then absorbed by the Galaxy medium in this model.

We modelled the photoionised emission lines in NGC 1068 using the sophisticated and self-consistent photoionisation model PION (Mehdipour et al. 2016) in SPEX. PION uses the full broadband SED continuum (in this case adopted from NGC 7469; Mehdipour et al. 2018) to calculate the emission spectrum and the ionisation/thermal balance of the plasma, simultaneously. To fit the emission lines in both epochs, we fitted one component at a time, fixing the previous component before adding the next, until no further improvement on the overall fit statistic was achieved. In all PION components, we fitted the column density (N_H), ionisation parameter (ξ), covering fraction (C_cov), and outflow and turbulent velocities (v_out and v_turb, respectively).

After fitting the photoionised emission (PION) components, we still found residuals in some of the emission features (mostly in the region between 14 and 18 Å). Although the large consensus regarding CIE plasma is that it is not the dominant source for X-ray emission in the nucleus of NGC 1068 (e.g. Young et al. 2001; Kinkhabwala et al. 2002; Ogle et al. 2003), there is evidence of a star burst region (e.g. García-Burillo et al. 2014). Therefore, we introduced a CIE component to the model to see if this could explain the residuals in the spectra, fitting the electron temperature (T_e), emission measure (EM), and outflow (v_out) and turbulent (v_t) velocities. This was not to account for the lower order resonance lines (expected to be produced via photoexcitation processes; Kinkhabwala et al. 2002), as PION takes these into consideration, but rather to fit the residuals in the emission spectrum between 14 and 18 Å. Figure 4 shows that some of the emission lines in the 2000 RGS spectrum could not be explained by PIE plasma (red line), but were accounted for when we fitted for CIE plasma (orange line). In NGC 7469, Behar et al. (2017) found residuals in the RGS spectrum at 15.3 and 17.4 Å (Fe XVII) which they interpreted as emission in CIE plasma from the star burst region. Therefore, we included a CIE component in our modelling.

Fig. 4. RGS spectrum (grey points) between 8 and 20 Å from the 2000 observation of NGC 1068. The red line shows the best fit model and the blue line displays the photoionisation model made up of four different PION components, fitting the majority of the emission lines in both the RGS and PN spectra, except for the features at 15 and 17 Å. We therefore introduced a collisionally ionised component (CIE; orange line) to account for these emission lines from Fe XVII.

The elemental abundances fitted in this model were for C, N, O, Ne, Mg, Si, S, Ar, Ca and Ni, all with respect to Fe (the abundance of Fe was found to be A_Fe = 2.45², relative to Solar, in NGC 1068; Matt et al. 2004; Bauer et al. 2015) as Fe emission lines were found in both RGS and PN spectra (oxygen is the strongest feature only in the RGS energy band). However, as C, N, O, Ne and Mg were only found in the RGS energy range, and Si, S, Ar, Ca and Ni were only present in the PN spectrum, we split the abundances up into two PION components when fitting. The first PION component fitted the high energy emission lines of the five higher Z elements and the second PION component (RGS band) accounted for the five lower Z elements. The remaining PION and CIE components had their abundance values coupled to these two PION components. This meant that the model fitted the abundances for all the lines from each component, but with each element free only once. In order to do this, we assumed that the chemical enrichment is the same in all emitting plasma regions (photoionised or collisionally ionised) throughout the AGN nucleus and star burst region of NGC 1068.

5. Spectral results

In Fig. 3, the combined RGS and PN spectra show many complex features. This is a result of a mixture of many emission lines, some of which are blended, and an obscured continuum that is not easily modelled. This meant that acquiring a statistically significant best fit, represented by the C-statistic, was very challenging. Therefore, in both epochs, we used the change in C-statistic (ΔC) to signify the improvement on the best fit from each fitted component. Despite this however, ΔC could still be very large due to the complexity of the spectrum and its statistical quality.

5.1. 2000 best fit

The initial fitting of the observed continuum (POW and REFL) gave a C-statistic of C = 58 486 (for 1075 degrees of freedom; hereafter d.o.f.). Throughout the emission component fitting (PION and CIE), the continuum parameters were free to vary. The change in C-statistic was with respect to the previous component, which we fixed when introducing the next component, unless stated otherwise.

The first PION emission component (EM1) was required to fit the high energy emission lines in the PN spectrum (E >  6.4 keV). This highly ionised (log ξ = 3.91) component improved the fit with a fairly modest ΔC = 2200 compared to the continuum. The second PION component (EM2; log ξ = 0.63) accounted for many of the narrow emission lines in the RGS spectrum, improving the quality of the fit by ΔC = 35 600. When we fitted both EM1 and EM2 together, the fit improved by ΔC = 270. A third PION component (EM3; log ξ = 1.84) fitted lines in both RGS and PN spectra, improving the fit by ΔC = 4700. Again, we then freed all (three) PION components, decreasing the statistic by ΔC = 1600. A fourth component (EM4; log ξ = 3.02) further improved the fit, albeit not as significantly as the other three components, by ΔC = 220. All four PION components were then fitted together, improving the model significantly by a further ΔC = 240.

The CIE component (CI1) improved the fit by ΔC = 1440 (after fixing all the PION components), accounting for all the emission lines not produced by photoionised plasma, with large enough statistical significance. We then freed and refitted all four photoionised PION components with the CIE component to obtain a change in the C-statistic of ΔC = 840.

Finally, we fitted the abundances of the ten elements in the spectrum, with respect to iron. Firstly, we fixed all the parameters in the PION and CIE components, in order not to fit too many parameters at the same time when adding some new free parameters. Otherwise, this could cause a degenerate result for the parameter values for the same ΔC. We freed and refitted the emission components after we accounted for the abundances. Therefore, we fitted the abundances of the elements in EM2 for the RGS spectrum (C, N, O, Ne, and Mg) followed by the high energy elements (Si, S, Ar, Ca, and Ni) in EM1. The change in C-statistic for these abundances were ΔC = 4900 and ΔC = 330, respectively. From this, we slowly fitted the emission components, firstly EM1 and EM2 (ΔC = 220) followed by EM3 and EM4 (ΔC = 500) and then CI1 (ΔC = 500), with all abundances left free throughout these last few steps. Here, all parameters were fitted to obtain a new global C-statistic minimum of C = 4970 for 1040 d.o.f.

Despite obtaining a relatively good fit so far, there was still an emission feature not accounted for. This line, at around 7.4 keV, belonging to the neutral Ni Kα, comes from the X-ray reflection off neutral material similarly to the Fe Kα line. As REFL in SPEX only accounts for neutral Fe, we fitted this line with a Gaussian (GAUS; for simplicity), freeing the normalisation, energy and line broadening, which improved the fit by ΔC = 130. A similar approach was taken by Semena et al. (2019) for NGC 5643.

All the parameters were then refitted. The best fit for our model, simultaneously applied to both the RGS and PN spectra, was C = 4530 for 1037 d.o.f. The final best fit model, including more GAUS components (see Sect. 5.3), over imposed on the spectrum is shown in Fig. 5 and the best fit parameter values for the continuum, PION, CIE, and abundances are displayed in Tables 3–6, respectively. Table 7 shows the best fit parameter values of the Ni Kα line. All parameter errors are obtained with all fitted parameters free.

Fig. 5. Best fit model to the 2000 spectrum of NGC 1068. We plot the soft X-ray (RGS) band in wavelength units and the hard X-ray (EPIC-PN) band in energy units for display purposes. The red line shows the best fit to the data points (grey crosses) and the other coloured lines (labelled in the legends) represent each component in the model. Bottom panels: residuals between the best fit model and the observed data points.

5.2. 2014 best fit

We fitted first the continuum model (POW and REFL) to the combined spectrum from the 2014 observations, where the initial global fit statistic was C = 125 643 (for 1065 d.o.f.). All continuum parameters were freed when fitting the emission lines. The change in C-statistic was with respect to the previous component, that was fixed, when introducing the next one, unless stated otherwise.

For the 2014 spectrum, we also fitted four PION emission components. The first two components, EMA (log ξ = 3.96) and EMB (log ξ = 0.67), improved the fit by ΔC_EMA = 8600 and ΔC_EMB = 76 100, respectively. Fitting these two components together gave a ΔC = 2700. We then fitted a third component (EMC; log ξ = 1.91) which yielded ΔC = 5740. A further ΔC = 4670 was achieved when we freed the parameters in EMA and EMB together with EMC. The fourth, and final, component (EMD; log ξ = 3.02) improved the fit by ΔC = 1000. When we fitted all four PION components together, we achieved a further improvement in the fit of ΔC = 540.

Similarly to the 2000 spectrum, not all the emission lines were fitted with PION components. We therefore added a CIE component (CIA) to account for the lines produced by a CIE plasma. CIA significantly improved the fit by ΔC = 1960, after we had fixed all of the PION components. Refitting the four PION components with the CIE component further improved the fit by ΔC = 760.

Next we fitted the abundances, subsequently to fixing all PION and CIE components. The five lower Z elements (C, N, O, Ne, Mg) were fitted in EMB (accounting for the RGS lines) obtaining ΔC = 10440, while the five higher Z elements (Si, S, Ar, Ca, Ni) fitted in EMA improved the fit by ΔC = 420. We then fitted the parameters in each component together, with the abundances free, starting with EMA and EMB, followed by EMC and EMD, which gave changes in the C-statistic of ΔC = 790 and ΔC = 2100, respectively. We then freed the CIA parameters too and obtained a new best fit with a decrease in C-statistic of ΔC = 800.

Finally, we accounted for the neutral Ni Kα line at 7.5 keV with a Gaussian component (GAUS), fitting the normalisation, energy and broadening velocity. Adding this line to the model improved the best fit by ΔC = 440. Table 7 shows the best fit parameter values of the Ni Kα line; we were unable to constrain these parameters so we fixed them to their initially fitted values in the subsequent fits and error searches.

All the parameters were fitted together again to obtain a best fit of C = 8510 for 1027 d.o.f. The spectrum and final best fit model, after including more GAUS components (see Sect. 5.3), are shown in Fig. 6. The best fit parameter values for the continuum, PION, CIE, and abundances are displayed in Tables 3–6, respectively. Again, all errors are obtained with all fitted parameters free.

Fig. 6. Best fit model to the 2014 spectrum of NGC 1068. We plot the soft X-ray (RGS) band in wavelength units and the hard X-ray (EPIC-PN) band in energy units for display purposes. The red line shows the best fit to the data points (grey crosses) and the other coloured lines (labelled in the legends) represent each component in the model. Bottom panels: residuals between the best fit model and the observed data points.

5.3. Fitting the RGS and PN spectra separately

Unfortunately, the fitting statistics for the 2000 and 2014 data are not very good: $\frac{C}{\mathrm{d}.\mathrm{o}.\mathrm{f}.}=4.37$ for 2000, and $\frac{C}{\mathrm{d}.\mathrm{o}.\mathrm{f}.}=8.29$ for 2014. The complex spectrum of NGC 1068, including line blending (O VII forbidden line with the N VI RRC between 22 and 23 Å) and unresolved (for example, Mg XI triplet at 9 Å) emission lines, means we are unable to fit all the lines well. For many lines, in particular the He-like triplets, the model underpredicts these features, such that we obtain line ratios that are inconsistent with pure photoionised plasma. This comes down to PION fitting multiple emission lines (with certain ξ and N_H values) simultaneously, rather than individually, making it difficult to achieve a fit that is fully consistent with a photoionised plasma. It is likely that the line-emitting regions in NGC 1068 are more complex than our model can completely match. The geometry and variations in the density and ionisation of the gas, such as clumpy or stratified gas, all affect the overall line spectrum that we get. Nevertheless, this PION model still gives us a more useful insight than an empirical model with Gaussian line fitting. Also, some of the bad fit may be attributed to the instrumental calibration of RGS, which is not perfect. Therefore, we took the best fit models from the simultaneous fits, and folded them to the RGS and PN spectra individually in each epoch, to see if the model is a good fit to individual lines in the separate spectra.

When folding the combined best fit to the PN spectra from each epoch, the statistics were C = 298 for 54 d.o.f. ($\frac{C}{\mathrm{d}.\mathrm{o}.\mathrm{f}.}=5.52$) in 2000, and C = 502 for 54 d.o.f. ($\frac{C}{\mathrm{d}.\mathrm{o}.\mathrm{f}.}=9.30$) in 2014. This clearly shows a poor fit to the PN spectra only, which count only a total of 94 and 83 bins for 2000 and 2014, respectively. There are residuals between 4 and 6 keV in both epochs, and above 9 keV in 2014. There is also a possible feature at around 5.9 keV (just before the Fe Kα), which is also present in the spectrum of Pounds & Vaughan (2006; their Fig. 3).

In the RGS band, the C-statistic for 2000 was C = 3921 for 994 d.o.f. ($\frac{C}{\mathrm{d}.\mathrm{o}.\mathrm{f}.}=3.94$), whereas for the 2014 data C = 7392 for 994 d.o.f. ($\frac{C}{\mathrm{d}.\mathrm{o}.\mathrm{f}.}=7.44$). In the RGS spectra, there were still many residuals where neither PION nor CIE fitted the data well, in particular in the region between 32 and 36 Å. Therefore, we introduced some Gaussian components (GAUS) to the RGS data only, to improve the fit for these lines. Four GAUS components were fitted to the lines at 35.6, 34.9, 34.3 and 32.2 Å, while keeping the other model parameters fixed, improving the best fit in the 2000 model by ΔC = 48, 219, 106, and 66, respectively. In the 2014 data, adding these same four GAUS components, with the other parameters fixed, improved the fit by ΔC = 105, 466, 169, and 80, respectively. We note that there are no likely identifications with ion species for the lines fitted in the 32 to 36 Å wavelength range. There was also a large residual, in both spectra, at 8.4 Å, corresponding to the Mg XII Lyα line, that the PION components only half fitted. We therefore added a GAUS component for this line, improving the fit by ΔC = 56 in 2000 and ΔC = 180 in 2014. Although refitting the two low ionisation PION components and the CIE component to the RGS spectrum (with the GAUS parameters free) did improve the fit in each epoch, there were still many residuals, such as the O VII triplet between 21.5 and 22.5 Å. The best fit statistic for the RGS data only was C = 3411 for 973 d.o.f. ($\frac{C}{\mathrm{d}.\mathrm{o}.\mathrm{f}.}=3.51$) for 2000, and C = 6356 for 973 d.o.f. ($\frac{C}{\mathrm{d}.\mathrm{o}.\mathrm{f}.}=6.53$) for 2014.

As a sanity check, if we model the spectrum locally (between 18 and 23 Å), we can fit the O VII triplet and O VIII Lyα well. But when folding this local model to the global spectrum again (7−37 Å), the emission features outside the 18−23 Å range are very poorly fitted, as the parameter values depend on the lines fitted.

5.4. Final best fit

Finally, we added these 5 Gaussian components to our models, fitted to the simultaneous RGS and PN spectra, in each epoch. The 2000 and 2014 models were improved by ΔC = 520 and ΔC = 1000, respectively. As these Gaussians did not add to the photoionisation modelling and have no physical meaning, but were introduced to achieve a statistical improvement in the fitted model, we kept the GAUS parameters fixed. We also note that these Gaussian lines do not affect the overall parameter values of the photoionisation modelling.

We then refitted all the parameters together from Tables 3–7, and obtained final best fit C-statistic values of C = 3790 (for 1037 d.o.f.) and C = 7252 (1027 d.o.f.), for 2000 and 2014, respectively. Unfortunately, in both epochs the models were still statistically poor: $\frac{C}{\mathrm{d}.\mathrm{o}.\mathrm{f}.}=3.65$ for 2000 and $\frac{C}{\mathrm{d}.\mathrm{o}.\mathrm{f}.}=7.06$ for 2014. However, due to the complexity of such a rich spectrum, we were unable to improve the fit any further.

In our modelling, this poor fit (in both epochs) may be a result of the assumed geometry of PION, which takes a slab of material and produces forward emission facing us. This simple geometry works very well with type 1 Seyfert AGN because we are observing the AGN and plasma from face on. However, for type 2 AGN, where we see the plasma from side on, the modelling becomes more complex, especially if the outflowing wind is emitting in an ionising cone (Kinkhabwala et al. 2002). Therefore, a simple slab geometry may break down here as the geometry is fundamental for radiation transfer issues, where resonance lines are sensitive to the column density and can vary depending on the LOS. This being said, PION is an excellent tool at analysing photoionised plasma self consistently, calculating the photoionisation balance on the fly and computing explicitly the photoionisation properties (ξ and N_H) of all the emission lines in the spectrum. However, if the underlying AGN environment is far more complex than assumed by PION, such as an outflowing cone seen in type 2 AGN, then our fitting of the spectra will be simplistic to some extent, explaining why we cannot fit some of the emission lines very well.

On the other hand, this is a common result for NGC 1068, a Compton-thick AGN, as previous models have not been able to fully explain the Chandra or RGS spectra, showing strong residuals. For example, Kinkhabwala et al. (2002) were unable to model many emission lines between 10 and 16 Å or above 34 Å (see their Fig. 11). Brinkman et al. (2002) were able to fit the Chandra spectrum above 20 Å (in their Fig. 7), but many of the emission lines below 16 Å were not fitted well. Furthermore, Kraemer et al. (2015) were unable to fit the emission lines below 17 Å and above 31 Å (see their Figs. 1 and 4), and Kallman et al. (2014) were unable to model some of the strongest lines such as the O VII and Ne IX triplets, again both in the Chandra spectra.

Moreover, none of these modelling attempts accounted for the Fe XVII lines at 15 and 17 Å, which we model with a CIE component. This suggests either previous models did not account for photoexcitation, or CIE is a valid explanation. However, Kinkhabwala et al. (2002) did take these Fe XVII lines into account as due to photoexcitation, although the lines are not fitted in their plot (Fig. 11), whereas Brinkman et al. (2002) say their model did not include the Fe-L transitions in their model. Gu et al. (2019) tested the contributions to the population of levels by different processes within the Fe XVII ions, but do not consider photoexcitation. Further tests and comparisons are required in SPEX to fully understand these lines in order to reduce the limitations of photoionisation codes, and to allow PION to treat photoexcitation properly. Although a CIE component here is able to account for the Fe XVII lines at 15 and 17 Å, photoexcitation is just as physically viable, and should therefore be considered as a possibility in NGC 1068.

As a result, given the complexity of this spectrum, the models fitted (by Kinkhabwala et al. 2002; Brinkman et al. 2002; Kallman et al. 2014; Kraemer et al. 2015, and here) are unable to account for every feature. This means that in all the models the statistical fit is poor and there are many residuals. Only in this paper do we quantify and compare the C-statistic of each component and its contribution to the final model.

As a further comparison, the type 1.5 AGN NGC 4151 has a significantly similar RGS spectrum to that of NGC 1068, with many of the same emission features. The modelling of the soft X-ray spectra of NGC 4151 also showed strong under prediction of the continuum and some emission features (see Figs. 3 and 15 in Schurch et al. 2004; Beuchert et al. 2017, respectively), implying that it is difficult to obtain an acceptable model statistic, whilst fitting all the emission features in these obscured AGN.

5.5. Spectral line features

We display the complex spectra from the 2000 and 2014 observations in Figs. 5 and 6, respectively. The best fit model is shown by the red lines, while the other components are represented by different colours (see the legends in each figure). Four PION components in each epoch are required to fit the majority of the emission lines we observe, and we find that components with similar ξ values from each epoch account for the same emission features. For example, EM1 and EMA (green lines) fit the same high energy features, whereas EM2 and EMB (purple lines) fit the same lower energy RGS features, in their respective epochs. This makes it easier when comparing epochs to determine which lines are fitted by each component. This is the same for the CIE components.

Components EM1 and EMA fit the high energy Fe XXVI Lyα and Lyβ, and Fe XXV lines, present in the PN spectrum. EM2 and EMB fit emission lines in the RGS band, such as the O VII and N VI triplets, O VII, N VI, C VI and C V RRCs, and C V and N VI H- and He-like species lines. EM3 and EMC fit the Mg XI triplet lines, Ne IX RRC, and the Ne IX He-like species lines, in addition to O VIII Lyα, Si XII, and S XIII lines, in both the RGS and PN spectra. EM4 and EMD account for the Fe XVII and Fe XX lines in the RGS spectrum, and the S XVI, Ar XVII, Ar XVIII, Fe XXV lines in the PN data. The CIE component is required to explain the Fe XVII lines at 15.3 and 17.4 Å (see Fig. 4).

We find that some emission features are accounted for by many components at the same time. These include Ne X, O VIIIβ, N VII and C VI Lyα lines, Si XIII and S XV lines, and the Ne IX triplet. In addition, some of the emission lines (in the RGS band) are over predicted due to multiple components fitting them. These include the O VIII and N VI RRCs, the Mg XI forbidden line and Fe XVIII and Fe XIX lines around 14 Å.

On the other hand, we do account reasonably well for the features between 10 and 12 Å which were very under predicted by Kraemer et al. (2015), suggesting that this was because their model had incomplete atomic data. Furthermore, this wavelength region coincides with the failed CCD 7 in RGS1, meaning the effective area has dropped by a factor of two as only data from RGS2 can be analysed between 10 and 15 Å (den Herder et al. 2001). There are still some residuals between 9 and 10 Å, but this may be a result of multiple lines being blended together.

5.6. Luminosities

We take the ionising luminosity (1−1000 Ryd or 13.6 eV−13.6 keV) for both epochs to be L_ion = 1.54 ± 0.04 × 10³⁷ W, calculated from the assumed ionising SED³; we use this value for further calculations (Sect. 7.1). In addition, we calculate the observed 2−10 keV luminosity (using POW and REFL parameters in Table 3) from our best fit model, where we obtain $L_{2-10}^{2000}=1.{07}_{-0.04}^{+0.06}\times {10}^{34}$ W and $L_{2-10}^{2014}=1.{04}_{-0.03}^{+0.04}\times {10}^{34}$ W, for the 2000 and 2014 observations, respectively. The intrinsic 2−10 keV luminosity unveiled during a temporary decrease in absorption was $L_{\mathrm{obs},2-10}=7_{-4}^{+7}\times {10}^{36}$ W (Marinucci et al. 2016). By taking the ratio between our observed luminosities and the intrinsic luminosity, we find that roughly 0.15% of the intrinsic X-ray source is reflected and scattered into our LOS.

6. Comparing epochs

After fitting the simultaneous RGS and EPIC-PN spectra of NGC 1068 from 2000 and 2014, we find almost no difference in the emission features and component properties between each epoch. Here we discuss our overall findings and results from our photoionisation modelling.

6.1. Observed continuum

Starting with the observed continuum, both the power-law and reflection components differ between epochs (see Table 3). The photon index is flatter in 2014, while there is a decrease in normalisation for both POW and REFL components compared to 2000; they are not consistent within the uncertainties. We are able to constrain the line width (σ_v) of the neutral Fe Kα line (measured by coupling a VGAU component to REFL) in 2000, but obtain an upper limit for 2014.

6.2. Photoionised plasma components

For the four PION components, we find very little difference in the ionisation parameters between the two observations (not significant given the very little change in spectra). This suggests that the components with the same ionisation parameter in each epoch are the same plasma region. On the other hand, we do find the total equivalent column density increases from $N_{\mathrm{H}}^{2000}={131}_{-6}^{+5}\times {10}^{25}$ m⁻² in 2000 to $N_{\mathrm{H}}^{2014}={150}_{-6}^{+7}\times {10}^{25}$ m⁻² in 2014. This change of ΔN_H = 19 × 10²⁵ m⁻² is fairly significant, with an average of ΔC ∼ 3900, when substituting the N_H values from one epoch into the model of the other. Below we describe some possible reasons for this increase in total equivalent column density, by comparing each component of similar ionisation state.

A degeneracy was found between N_H and C_cov by Di Gesu et al. (2017) in their analysis of the Seyfert 1 AGN 1H 0419−577, using PION. This may explain the decrease in column density (ΔN_H = −17 × 10²⁵ m⁻²) between EM1 and EMA, where there is an increase in covering fraction (ΔC_cov = 0.13). For EM4 and EMD, the change in column density (ΔN_H = +19 × 10²⁵ m⁻²) between the two epochs may actually be due to a degeneracy between N_H and v_turb. In this case there is a small change in covering fraction (ΔC_cov = 0.01), but the turbulent velocity has increased by Δv_turb = 670 km s⁻¹ (the uncertainties of each parameters are smaller than this difference). This may be due to the large turbulent velocities being driven by line ratios, rather than line widths, causing changes in N_H between components and epochs. In addition, the curve-of-growth illustrates how both N_H and v_turb (line broadness) can affect the intensity of the lines. For EM2 and EMB, and EM3 and EMC, the changes in column density are ΔN_H = +10 × 10²⁵ m⁻² and +7 × 10²⁵ m⁻², respectively. The change in C_cov is negligible, but the v_turb values increase by 40 and 30 km⁻¹, respectively for EM2/B and EM3/C; this difference is larger than the uncertainties on the parameters for EM2/B. The column density changes for each individual component between epochs overall cause a total increase in equivalent column density of ΔN_H = 19 × 10²⁵ m⁻². The varying N_H values between the two epochs are unlikely to be real, however, as there is very little spectral change between the two epochs (Fig. 2). Therefore, the line column densities cannot vary as greatly as suggested here by the fit, and the most probable explanation for the increase in total N_H from 2000 and 2014 is model degeneracy with v_turb and C_cov.

The C_cov values for EM1 and EMA are significantly larger than the other three components, suggesting that this plasma region is closer to the central black hole than the other components. This is also consistent with the high ionisation parameter. However, these values are significantly smaller compared to the covering factor values of the two RGS components found by Kraemer et al. (2015). Furthermore, the broadening velocities (v_turb) in Table 4 of both EM1 and EMA are consistent with the broadening velocity of the Balmer lines measured at σ_Balmer ∼ 3200 km s⁻¹ (Antonucci & Miller 1985). These measurements indicate that the highly ionised components are close to the central SMBH, possibly consistent with the BLR, although we cannot observe the BLR directly due to obscuration from the torus. So for now, we assume all plasma regions are part of the NLR. On the other hand, if the X-ray BLR emission is largely reflected (like the optical BLR emission; Antonucci & Miller 1985), then the BLR would be a possible origin for the highly ionised emission lines.

In addition, the emission measures (EM) of each PION component are consistent between the two epochs. The definition of the EM is given by ∫n_en_HdV, but as the number densities of electrons (n_e) and ions (n_H), respectively, are unknown, we use the following equation (e.g. Mao et al. 2019; Grafton-Waters et al. 2020)

$$\begin{array}{c}\hfill \mathrm{EM}=\frac{n_{\mathrm{e}}}{n_{\mathrm{H}}}\frac{4\pi N_{\mathrm{H}}{C}_{\mathrm{cov}}{L}_{\mathrm{ion}}}{\xi },\end{array}$$(1)

where n_e/n_H ∼ 1.2 for fully ionised plasma and L_ion is the ionising SED luminosity L_ion = 1.54 × 10³⁷ W. For the PION components, if we assume a constant L_ion and C_cov, then EM depends on N_H and in particular ξ. If ξ is small, then the number density in the plasma is larger (for a fixed distance), which means EM (from the definition) increases, and thus explains why the lower ionised plasmas have the highest EM values. The EM, ξ and N_H values for these four components in each epoch (Table 4) are consistent with the model components from Kallman et al. (2014).

6.3. Collisionally ionised plasma

The electron temperatures (T_e) and line width (v_t) for each CIE component have not changed (within errors) between the two epochs, but the EM and v_out values are not consistent within the errors (see Table 5). However, much like with the PION components, it is the outflow velocities that differ the most. In 2000, the outflow velocity is negligible with ${\mathit{v}}_{\mathrm{out}}=-{40}_{-40}^{+30}$ km s⁻¹, whereas in 2014 the CIE plasma appears to be travelling at −165 ± 30 km s⁻¹.

Although the CIE components in each epoch account for the Fe XVII lines that PION cannot fit⁴, this is not an overall convincing conclusion. One argument against CIE is that the RGS spectra show strong, narrow RRC features, implying the X-ray emission is only from PIE plasma, making CIE unlikely in NGC 1068. However, RRC emission could come from outflowing PIE plasma closer to black hole, while the CIE emission could originate from further out, such as the secondary region or star burst region. In terms of the Fe-L lines themseleves, the Fe XVII line at 17 Å is populated by recombination rather than direct excitation, so PIE plasma should dominate this line. As for the Fe XVII resonance line at 15 Å, the 3d shell is highly populated if the plasma density is high (in a similar way the resonance line in He-like triplets is more dominant over the forbidden line if the density is large), suggesting that CIE plasma can account for this. Figures 1 and 2 from Liedahl et al. (1990) show that the Fe XVII at 17 Å should be accounted for by PIE plasma, whereas the 15 Å line is emitted in CIE plasma, reiterating the arguments above. However, photoexcitation was not considered by Liedahl et al. (1990), and would explain the 15 Å resonance line (Sako et al. 2000; Kinkhabwala et al. 2002); in this case the explanation for the lack of PIE emission at 17 Å is, however, unknown. Therefore, the under-prediction of these lines by PION suggests we are not getting the full picture, and although CIE plasma can produce these Fe XVII lines (Fig. 4), photoexcitation is a more natural process in the outflowing wind (Kinkhabwala et al. 2002). Although we fit a CIE component to each epoch to account for the Fe XVII lines, photoexcitation in PIE plasma cannot be dismissed.

6.4. Abundances

We fitted abundances for ten elements that are present in both RGS and PN spectra of NGC 1068 (from Kinkhabwala et al. 2002; Pounds & Vaughan 2006, respectively). Between the two epochs, there is very little difference in the ratios with respect to Fe as shown in Table 6. Kinkhabwala et al. (2002) and Brinkman et al. (2002) found the N abundance to be 2−3 times Solar, however we do not find any evidence for this (see Table 6).

6.5. The Ni Kα line

Table 7 shows the best fit values for the Ni Kα, modelled with a GAUS component. However, for the 2014 line, we are unable to constrain the parameters of this GAUS, so we fix them to the values initially fitted. As a comparison, Rahin & Behar (2020) constrained the width of the Fe Kα and Ni Kα lines to be of the order of 300 km s⁻¹, lower than what we quote here from EPIC-PN. This could be due to the HETGS instrument having a higher spectral resolution than EPIC-PN.

7. Discussion

In this section we discuss the results from our spectral modelling of NGC 1068 in 2000 and 2014. Here we use our data to infer and determine properties of the emitting plasma that surrounds the nucleus in this AGN. In particular, we investigate the distances of each photoionised component from the SMBH, and determine the thermal stability of the plasma regions for the two epochs, before studying the validity of the CIE emission.

7.1. Photoionised plasma distances

After determining our best fit photoionisation models for the NGC 1068 spectra in both epochs, we now want to establish the locations of the emitting plasma regions, with respect to the central SMBH. Unfortunately, as seen in Fig. 2 and Table 4, there is no variability between the two epochs, either in terms of the spectral features or the ionisation state of the plasma components, suggesting that the emitting plasma has not changed over the 14 year period between observations. This means that we are unable to obtain distance measurements from variability arguments where a change in the ionisation parameter is due to a change in the SED shape (see e.g. Mehdipour et al. 2018, for NGC 7469). In addition, as we are unable to directly observe the ionising X-ray source, attaining an intrinsic SED is practically impossible. Therefore, we have to investigate alternative ways to calculate the distances of the plasma regions, discussing the physicality, as well as the advantages and negatives, of each method. These measurements depend on the geometry and our LOS view of NGC 1068, which adds complexity in comparison to Seyfert 1 AGN.

The uncertainties on the distance estimates here are obtained from the parameters and errors of our best fit modelling using PION (Table 4). The hydrogen number densities for each component, and their respective errors, are shown in Table 8; we multiply these values by 1.2 to get the electron number density that we use in our calculations. Ideally, the plasma density should be constrained through the recombination timescale using the variability between observations, but as this is not the case here, we have used the calculated densities from PION.

7.1.1. Ionisation parameter distances

We start by estimating the locations of these emission components using the definition of the ionisation parameter (ξ)

$$\begin{array}{c}\hfill \xi =\frac{L_{\mathrm{ion}}}{n_{\mathrm{e}}{R}^2},\end{array}$$(2)

where n_e is the electron number density, L_ion is the ionising luminosity (measured between 1 and 1000 Ryd) and R is the distance of the plasma from the black hole. From the fits with PION, we are able to derive the hydrogen number densities (n_H)⁵ for each emission component (displayed in Table 8). Obtaining the densities allows us to estimate the distances of each component; alternatively, Peretz et al. (2019) used the R ratios and plots from Porquet & Dubau (2000) to obtain the number density for the BLR in NGC 4051.

Unfortunately, in both these methods (and the method below in Sect. 7.1.2), the number density values of EM1 and EMA cannot be accurately determined. Therefore the distance estimates cannot be obtained for these two components. For EM4 and EMD, we are unable to constrain the values and instead have density upper limit values. For the other two components in each epoch, the number densities are constrained.

Using the ionising luminosity of L_ion = 1.54 × 10³⁷ W (for both epochs) and the ionisation parameters from Table 4 we can estimate the distances for each component using Eq. (2). Alternatively, we can obtain a similar answer using the EM from Eq. (1) and substituting it into the following equation (Mao et al. 2018)

$$\begin{array}{c}\hfill R^2=\frac{\mathrm{EM}}{4\pi n_{\mathrm{e}}{C}_{\mathrm{cov}}{N}_{\mathrm{H}}},\end{array}$$(3)

where N_H and C_cov are values from Table 4. However, this equation cancels down to Eq. (2) when we substitute EM from Eq. (1), but with a factor of 1.2 difference between n_e and n_H.

These methods provide a first distance indication, but work best when the intrinsic SED is known and when there is large variability between observations. The resulting distances (R_ξ and R_EM, respectively) are listed in Table 8.

7.1.2. Using the size of the component

The second method considers the size of each component estimated from the ratio between the column density (N_H) and electron number density (n_e = 1.2n_H for fully ionised plasma). We then use the relation between the covering fraction (C_cov) and the solid angle (Ω) to find the area A of each component as seen from the black hole:

$$\begin{array}{c}\hfill C_{\mathrm{cov}}=\frac{\mathrm{\Omega }}{4\pi }=\frac{A}{4\pi R^2}\end{array}$$(4)

where R is the distance from the black hole. Here, we assume that each emitting component is spherical and the black hole sees only half of the sphere with area A = 2πr², where r is the radius of the component. However, the radius of each sphere, which is its characteristic thickness, is $\mathrm{\Delta }r\approx \frac{N_{\mathrm{H}}}{n_{\mathrm{e}}}$ (identical to Eq. (6) in Peretz et al. 2019), and therefore A = πΔr².

Substituting A into Eq. (4) and rearranging for R gives the distance from the black hole, given the size of each component

$$\begin{array}{c}\hfill R^2\approx \frac{\mathrm{\Delta }{r}^2}{4C_{\mathrm{cov}}}·\end{array}$$(5)

With this method it is fairly easy to obtain distance measurements based on the geometry between the black hole and plasma region. In addition, Eq. (5) does not depend directly on the SED, although N_H and n_e are measured with PION which does require the correct SED input. However, the N_H/n_e ratio gives an approximation for the thickness of a slab of emitting plasma, whereas we assume each component is spherical. Our assumption that each component is spherical is a crude one, and an over simplification of the plasma geometry that does not relate easily to the thickness. Maybe a more appropriate geometry for the plasma is a cube, with dimensions equal to the thickness N_H/n_e. In addition, the column density (and number density) are measured with respect to our LOS, both of which are dependent on the viewing angle. This model relies heavily on the size and geometry of each plasma component, which is difficult to calculate with certainty, especially when viewed from side on. Distance estimates obtained with this method (R_Δr) are shown in Table 8, except for EM1 and EMA due to inaccurate n_e values.

7.1.3. Escape velocity distances

Next we consider the escape velocity of each component, and assume the outflow velocities from Table 4 are large enough to allow the plasma components to escape the gravitational potential of the black hole. Therefore, setting v_out = v_esc we can obtain distances from

$$\begin{array}{c}\hfill R=\frac{2G{M}_{\mathrm{BH}}}{v_{\mathrm{out}}^2},\end{array}$$(6)

where G is the gravitational constant and M_BH = 1.6 × 10⁷ M_⊙ is the mass of the black hole (Panessa et al. 2006).

This method is independent on the SED. However, the outflow velocity, v_out, is along our LOS and therefore depends on the viewing angle with respect to a face on view (type 1) AGN. Although the outflowing wind makes a cone-like shape, the v_out we measure is not the net outflow velocity and depends strongly on wavelength accuracy of the RGS spectrum (δλ = 8 mÅ; den Herder et al. 2001). In addition, while Eq. (6) was used by Blustin et al. (2005) to obtain minimum distance limits on the WA components in AGN, the R_vout values we derive are not the lowest distance measurements for each component in Table 8.

7.1.4. Estimated volume filling factor distances

In NGC 7469, both the WA (Blustin et al. 2007) and emission line region (Grafton-Waters et al. 2020) distances were estimated using the volume filling factor f_v. The f_v values were calculated for the WA components, but the values for the emission line regions had to be fixed at some arbitrary values for the NLR and BLR (see Grafton-Waters et al. 2020, and references within). Ideally, we would want to calculate the f_v values for each individual emission component within NGC 1068.

The volume filling factor equation was derived from the mass outflow rate by Blustin et al. (2005). They calculated the distances of the WA in many type 1 AGN by assuming that the momentum of the outflowing wind is related to the momentum of radiation being absorbed (P_abs) plus the momentum of the scattered radiation (P_scat), which depends on the size of each plasma region (Blustin et al. 2005). The equation is given by

$$\begin{array}{c}\hfill f_v=\frac{({\dot{P}}_{\mathrm{abs}}+{\dot{P}}_{\mathrm{scat}})\xi }{1.23m_{\mathrm{p}}{L}_{\mathrm{ion}}{v}_{\mathrm{out}}^{2}4\pi C_{\mathrm{cov}}},\end{array}$$(7)

where m_p is the proton mass, and the remaining parameters are calculated in the best fit models for each epoch.

This equation works well for the WA components in type 1 AGN, such as NGC 7469, because our LOS is approximately aligned with the wind’s direction of motion, and the radiation, which drives the wind (giving it momentum), can be measured directly. However, in type 2 AGN, the direction of the plasma motion is almost perpendicular to our LOS, and we are unable to measure the intrinsic luminosity in order to determine, for example, the two momenta in Eq. (7). This, therefore, means we are unable to obtain accurate f_v values for the components within NGC 1068.

Instead we set values for the volume filling factor. The volume filling factor has been quoted to be between 0.001 and 0.01 (e.g. Osterbrock 1991; Snedden & Gaskell 1999) for the BLR, but there is no information regarding the NLR. In NGC 7469, Grafton-Waters et al. (2020) set f_v = 0.1 for the two NLR components, such that the emission line regions were further from the black hole than the WA components, and f_v = 0.001 for the BLR component. In NGC 1068, we do not measure any absorption due to our LOS view (e.g. the cone model in Kinkhabwala et al. 2002), but we initially set f_v = 0.1 for all PION components, although Kallman et al. (2014) suggested a volume filling factor of 0.01.

Components EM1 in 2000 and EMA in 2014 have large C_cov values relative to the other components, suggesting they are more extended as seen by the SMBH, and are therefore closer (see Eq. (5)). The broadening velocities (v_turb) in Table 4 of these components are consistent with those of the Balmer lines in NGC 1068 (σ_Balmer ∼ 3200 km s⁻¹; Antonucci & Miller 1985), and with the line broadening velocities for BLRs found in type 1 AGN (e.g. Kollatschny & Zetzl 2013). Furthermore, the ionisation parameters of EM1 and EMA are very large, and these components only fit the high energy lines in the PN spectra for both epochs, suggesting that these lines are produced from highly ionised plasma, closer to the central black hole than the other components.

The parameter values point us to the conclusion that EM1 and EMA could be part of the BLR. However, we cannot see the emission from the BLR because the nucleus in NGC 1068 is obscured by the torus. In optical observations, the BLR is reflected off the far side of the dusty torus, or scattered off free electrons, where the light we observe is polarised, relative to the unreflected emission (Antonucci & Miller 1985). The amount of polarisation of the BLR emission is consistent to that of the underlying continuum (in our case modelled with REFL), whereas the NLR shows negligible polarisation (Antonucci & Miller 1985). This could also be true for the X-ray BLR whereby the emission lines are reflected into our LOS. However, there is no significant evidence in our modelling that suggests the emission is being reflected into our LOS, so we cannot say conclusively that we are seeing emission from the BLR. Therefore, we assume both EM1 and EMA are observed directly and are part of the NLR.

There are two methods that can be used to obtain distance estimates using the volume filling factor. Both methods start with N_H = n_ef_vΔr and ξ = L_ion/n_eR², where n_e is the electron number density, Δr is the thickness of the plasma region and R is the plasma distance from black hole; L_ion, ξ, and N_H are from our modelling. The first derivation, from Blustin et al. (2005) (and used on NGC 7469 by Blustin et al. 2007), assumed a thin layer (of thickness Δr) for each plasma component containing most of the mass, with an ionisation ξ. This method led to a maximum distance limit because Blustin et al. (2005) assumed Δr ≤ R. Alternatively, Grafton-Waters et al. (2020; see also Behar et al. 2003; Whewell et al. 2015) derived the distance estimate by integrating n_e (substituting for $n_{\mathrm{e}}=\frac{L_{\mathrm{ion}}}{\xi R^2}$; Eq. (2)) over the thickness of the plasma region such that $N_{\mathrm{H}}={\displaystyle \int \frac{L_{\mathrm{ion}}{f}_{\mathit{v}}}{\xi R^2}}\mathrm{d}R$, where they assumed the plasma outer distance to be further from the black hole than the inner distance. Therefore, this second method measured a minimum distance value. The values calculated with these two methods are the same, it just depends on the method used as to whether an upper distance limit or a minimum distance is achieved. Therefore, for ease, we shall assume this equation gives the locations of the emission components,

$$\begin{array}{c}\hfill R=\frac{L_{\mathrm{ion}}{f}_v}{\xi N_{\mathrm{H}}}·\end{array}$$(8)

For NGC 1068, f_v cannot be measured for the plasma regions, therefore we set f_v to the same arbitrary values for each component in 2000 and 2014. Consequently, there are very large uncertainties associated with f_v, which means the errors on the distances obtained from this method are significantly greater than what is derived from our fitting and what is quoted in Table 8. The evident discrepancies between the result of this distance method and the others are likely to come from these large f_v uncertainties, which may be up to orders of magnitude in difference. Furthermore, Eq. (8) requires two PION parameters (N_H and ξ) as well as the ionising luminosity, none of which are secure because we do not know the true shape of the SED, adding to the uncertainties on the distance values. Distance estimates obtained from this method (R_fv) are shown in Table 8.

7.1.5. Comparing distance estimates

Figure 7 illustrates the estimated distances from the central SMBH, as a function of ionisation parameter of each component, while Table 8 displays them quantitatively. The distances vary because the parameters and assumptions used in each method are different, in addition to using a non-intrinsic SED in our photoionisation modelling. Although we are unable to constrain the distances shown in Fig. 7, and a range of estimated values for each component is found, these findings are very useful, with some components being within the same distance ranges as each other, suggesting they are part of the same emitting region (as found by Kallman et al. 2014). Overall, the distances for each component are consistent with the expectations based on the values in Table 4.

Fig. 7. Distances of each component are plotted for the different methods in Sect. 7.1, at the respective ionisation parameters for 2000 (top) and 2014 (bottom). The data point shapes are shown in the legend, as are the colours for each component. The dashed line in each plot represents the location of the torus, at rtorus = 3.5 pc (García-Burillo et al. 2016), while the dotted line is at 1000 gravitational radii ($R_{\mathrm{g}}=\frac{2G{M}_{\mathrm{BH}}}{c^2}$) from the SMBH, corresponding to a typical outer disk radius (Netzer 2013).

Figure 7 shows that the most ionised components, EM1 and EMA (log ξ ∼ 4), and EM4 and EMD (log ξ ∼ 3) are closest to the black hole, compared to the low ionised components. For EM1 and EMA, due to unconstrained number densities, we only have two distance estimates. The larger column density and outflow velocity of EM1 puts it closer to the black hole compared to EMA. In both epochs, the distances are within the torus radius, suggesting they could be part of the BLR. EM4 and EMD have the largest range in distance estimates, mostly due to the upper limits of the number density and therefore a lower distance limit. The method using the plasma size (R_Δr) places EM4 and EMD close to the black hole, with the number density upper limit consistent with the BLR (n_BLR ∼ 10¹⁵ m⁻³). However, this location is probably unphysical as the high density limit is a result of EM4 and EMD not fitting many emission lines (Figs. 5 and 6) and statistically improving the C-statistic the least. The consistencies in distances between high ionisation components, EM1 and EM4, and EMA and EMD, suggest that they are at similar locations; this is also shown by the stability plots in Fig. 8 (see Sect. 7.2).

Fig. 8. Top panel: electron temperature (Te) of the plasma as a function of ionisation parameter (ξ). The labelled circles correspond to the four PION components for each epoch (red for 2000 and blue for 2014) at their respective ionisation parameter values from Table 4. Bottom panel: thermal stability curve of electron temperature (Te) as a function of the pressure form of the ionisation parameter (Ξ). On the curve, the rate of cooling equals the rate of heating. Left of the curve, cooling dominates over heating and to the right, heating dominates over cooling. The labelled circles correspond to the Ξ values of the four PION components for each epoch: red for 2000 and blue for 2014.

For the low ionised components in each epoch, EM2 and EMB are furthest from the black hole, and have a large range of distance estimates, placing them either side of the torus, depending on the method used. For example, the large outflow velocity (compared to the other components) places them close to the black hole, while the low ionisation state indicates they are further away. EM3 and EMC, on the other hand, have the most constrained distance estimates, placing the plasma region consistent with the torus distance (R_torus = 3.5 pc; García-Burillo et al. 2016). These components have the lowest outflow velocities in each epoch, along with a low column density for a moderate ionisation state of log ξ ∼ 2. In addition, we can constrain the number density of EM3 and EMC.

None of the components are found at the location of the secondary region of X-ray emission at 290 pc (which extends to about 400−500 pc from the central black hole; Brinkman et al. 2002; Ogle et al. 2003; García-Burillo et al. 2017). In addition, the secondary region has a column density three times smaller compared to the primary region (Brinkman et al. 2002; Ogle et al. 2003), but none of the equivalent column densities of the four components (Table 4) are a third of the size of any other, again suggesting none of these components, in either epoch, are situated outside the circumnuclear disk.

7.2. Thermal stability of the plasma

From PION, we are able to extract the temperature of the plasma components using ionisation and thermal balance equations, by iterating over a range of ionisation parameters, assuming photoionisation and thermal equilibrium, for each ionisation state (Mehdipour et al. 2016). Figure 8 displays the thermal stability of the emitting plasma, shown as a plot of electron temperature (T_e) as a function of either ionisation parameter (ξ; top panel) or the pressure form of the ionisation parameter (Ξ; bottom panel), within NGC 1068. The labelled circles show the locations of the four PION components for 2000 (red) and 2014 (blue), respectively. In the top panel of Fig. 8, T_e increases with an increase in ξ (same for both epochs), which is to be expected as the higher the plasma ionisation, the more energy the electrons within it have, and therefore the hotter the gas is.

The bottom panel of Fig. 8 provides an effective tool to determine which plasma regions are thermally stable or unstable, depending on the gas being in thermal equilibrium. This plot, known as a thermal stability curve, shows the electron temperature (T_e) as a function of the pressure form of the ionisation parameter (Ξ; Krolik et al. 1981) which is defined as

$$\begin{array}{c}\hfill \mathrm{\Xi }=\frac{F}{n_{\mathrm{H}}ck{T}_{\mathrm{e}}}=\frac{\xi }{4\pi ck{T}_{\mathrm{e}}},\end{array}$$(9)

where F = L_ion/4πr² is the flux of the ionising source (between 1 and 1000 Ryd), n_H is the hydrogen number density, c is the speed of light, k is the Boltzmann constant, T_e is the electron temperature, $\xi =\frac{L_{\mathrm{ion}}}{n_{\mathrm{H}}{r}^2}$ is the ionisation parameter, and r is the distance of the plasma from the central SMBH. On the thermal stability curve (black line in Fig. 8) the rate of heating equals the rate of cooling; to the left, cooling dominates over heating, and to the right heating dictates over cooling. Over-plotted onto the curve are the temperatures of the four PION components from each epoch (red for 2000 and blue for 2014), at their respective Ξ values (from Eq. (9)). EM2, EM3, and EM4 (in 2000) and EMB, EMC and EMD (in 2014) lie on the curve where the gradient is positive. This means that these components are thermally stable and can reach thermal equilibrium: should a small change in temperature occur, then the gas will move towards a heating or cooling process to balance out this change. EM1 (2000) and EMA (2014), however, lie on the curve where there is a negative gradient, although at the point where the curve starts to turn towards a positive gradient, which could be the result of an inappropriate SED we have used in the modelling. This means that the two components are thermally unstable and a little perturbation is enough to move the gas away from equilibrium. Instead, the gas may exist as two plasma phases, one cold with high density, the other hot with a low density, which together are in pressure equilibrium (Krolik et al. 1981).

Components EM1 and EMA and EM4 and EMD have similar Ξ values, suggesting they may be part of the same region within the outflowing wind, and are far from the locations of the lower ionised components, implying they are not associated with these. This is shown in Fig. 7, where the distances of EM1, EM4, EMA, and EMD are similar to each other. However, three of the distance methods for EM1 and EMA are not available, whereas three of the distance estimates for EM4 and EMD are just lower limits.

The other two components in each epoch are possibly discrete phases, rather than possessing a continuous distribution of ξ, as they lie on different stable parts of the curve (Ebrero et al. 2011). Furthermore, wind instabilities will cause clumping, but the lifespan of a clump is dependent on pressure equilibrium or stability of the plasma, as shown by its location on the curve (Fig. 8; Steenbrugge et al. 2005). Hence, EM1 and EMA cannot be in equilibrium with EM2 and EMB as their values of ξ (and therefore Ξ) are very different. This reiterates our findings from the distances for these two components, where we can conclude that EM1 and EMA, and EM2 and EMB are not part of the same outflowing region and the ionised clouds have different properties, thus placing them at the different locations in the wind.

Bianchi et al. (2019) proposed an alternative scenario for obscured AGN where radiation pressure compression (RPC) in the outflowing gas leads to a well defined ionisation distribution and differential emission measure (DEM). However, a constant gas pressure multi-phase scenario, as adopted here, was not excluded by their results. Their alternative scenario comes from the evidence that gas surrounding the AGN is photoionised by the continuum, producing soft X-ray and optical emission from the same plasma region in AGN, suggesting high and low density gas can coexist together (Bianchi et al. 2006). The AGN continuum ionises and heats up the plasma, compressing the gas and causing a gradient gas pressure.

Therefore, a DEM distribution is required to explain and fit the slope of a RPC gas, meaning we see X-ray emission from high-ξ and low-density gas, and optical emission from low-ξ and high-density gas, further away from the ionising continuum, in the same cloud (Bianchi et al. 2019). As Bianchi et al. (2019) point out, density diagnostics with forthcoming high-throughput and high spectral resolution instrumentation will be able to distinguish between the two interpretations.

In Fig. 9, we plot the emission measure (EM) of each PION component against log ξ. We then plot the 12 ions (purple squares) from Fig. 3 in Bianchi et al. (2019), adapted for the EM axis by multiplying the DEM values by the width of each log ξ bin. Finally, the DEM distribution (brown line) predicted for RPC gas is overlaid on top to test whether PIE plasma within a RPC gas is a solution. The overall trend shows a decreasing EM for an increase in log ξ, consistent with the relation for a DEM ($\frac{\mathrm{d}(\mathrm{EM})}{\mathrm{d}(log\xi )}$) in RPC (Bianchi et al. 2019). Both the distribution and overall trend suggest that the plasma is in RPC, which represents a natural solution. Here, however, the highest ξPION component does not fit the RPC trend, because it comes from the region that is almost transparent (almost fully ionised), so adding more high ξ gas makes very little difference for the RPC solution. This can also be seen in Fig. 11 in Kallman et al. (2014), where a plot of mass distribution as a function of ξ shows an inverse relation (increase in ξ for a decrease in mass). This suggests that a high ξ gas has less mass and a low optical depth, and is therefore located closer to the black hole, compared to a high-density, low-ξ gas, which we see in the RPC scenario.

Fig. 9. Emission measure (EM) of the four PION components in each epoch plotted versus their respective ionisation parameter (log ξ). The purple squares are the ions from Fig. 3 in Bianchi et al. (2019), whereby the EM values decrease for increasing log ξ, consistent with the PION components (same colours and data point styles from Fig. 7). We also plot the expected radiative pressure compression (RPC) curve (brown line) that follows the trend of both the ions (purple squares) and PION components suggesting that the photoionised plasma within NGC 1068 is consistent with RPC gas.

In a similar way to calculating the temperatures of each component, PION can also calculate the rates of the various heating and cooling mechanisms, for a given ionisation parameter. Figure 10 displays the total heating rate (top panel) and total cooling rate (bottom panel) within the PIE plasmas in NGC 1068, with the 2000 and 2014 PION components overplotted as red and blue circles, respectively. We also display the individual processes that make up these heating and cooling rates. The main heating processes are: Compton scattering (CS), Auger electrons (AE) and photoelectrons (PE). In Fig. 10, the largest contribution to heating the plasma comes from photoelectrons up to about log ξ = 3.5, where Compton scattering starts dominating. This is due to X-ray ionisation increasing the number of photoelectrons in the plasma with log ξ <  3.5, while above this value the free electrons interact and gain energy from the X-rays. The main cooling processes are: inverse-Compton scattering (ICS), bremsstrahlung (FF; free-free), collisional excitation (CE) and recombination (RR). The dominant cooling processes are collisional excitation for log ξ <  2.5 and bremsstrahlung for log ξ >  2.5. The latter comes about because the ions in the highly ionised plasma have lost all their outer shell electrons, such that when a free electron passes by, the large charge difference causes the electron to be deflected greatly, losing its kinetic energy to radiation and therefore cooling the gas. These results are consistent with the heating and cooling processes for various SEDs (Mehdipour et al. 2016).

Fig. 10. Top panel: total heating rate as a function of ionisation parameter (ξ) for NGC 1068 (black solid line). The contributions from the different heating processes are also shown: Compton scattering (CS; green dot-dashed line), photoelectrons (PE; purple dotted line), and Auger electrons (AE; orange dashed line). Bottom panel: total cooling rate as a function of ξ for NGC 1068 (black solid line). The contributions from the individual cooling processes are also displayed: inverse-Compton scattering (ICS; orange dot-dashed line), Bremsstrahlung (FF; green dashed line), collisional excitation (CE; purple dot-dashed-dotted line), and recombination (RE; magenta dotted line). Both panels: the labelled circles correspond to the four PION components in each epoch, shown as red for 2000 and blue for 2014, at their respective ionisation parameters from Table 4.

7.3. Origin of collisionally ionised plasma

Finally, if CIE plasma is the origin of the 15 and 17 Å lines (not considering photoexcitation), then we need to know the location of this emission. The most likely origin is from the star burst region, situated up to 1.3 kpc (García-Burillo et al. 2017) from the central black hole, which appears to be made up of two arcs, that together form a misshapen ring (Fig. 1a in García-Burillo et al. 2014). The SBR (and CND) appears to be a dynamic structure, as shown in Fig. 10 of García-Burillo et al. (2014), where the west of the SBR is blueshifted (down to −180 km s⁻¹), while the east is redshifted (up to +180 km s⁻¹). The RGS instrument has a spatial resolution (∼5′) large enough to cover the full SBR which is 0.8′ (∼50″) in diameter, meaning we can observe the full SBR and therefore emission from the blue and redshifted sides.

Chandra LETGS images of NGC 1068 show bright X-ray emission at the nucleus (including the secondary region), but only diffuse emission further out (Brinkman et al. 2002). On the other hand, Chandra HETGS images show X-ray emission that corresponds to the star burst ring of the host galaxy, south-east of the nucleus, in the energy range between 0.8−1.3 keV (Ogle et al. 2003). In addition, ACIS images show X-ray emission between 0.25 and 2 keV as far out as the SBR, in both north-west and south-east directions (Young et al. 2001). The X-ray energy ranges in which emission from the SBR is observed are consistent with the RGS band where we model the CIE plasma, suggesting that the SBR is a good candidate for the origin of collisionally ionised emission, although direct detection is yet to be confirmed (e.g. Ogle et al. 2003).

To confirm the legitimacy of this CIE component, more so than just by statistical significance in the model, we compare the calculated CIE luminosity (L_CIE) with the X-ray luminosity inferred from the far infrared (FIR) luminosity from the star burst region in NGC 1068. The star formation rate (SFR) in NGC 1068 was estimated between 0.4 and 0.7 M_⊙ yr⁻¹ (García-Burillo et al. 2014). From here, we can estimate the FIR luminosity of the SFR using $L_{\mathrm{FIR}}=\frac{\mathrm{SFR}}{4.5\times {10}^{-44}}$ erg s⁻¹ (Kennicutt 1998): we find the FIR luminosity to be between L_FIR[0.4]=8.9 × 10³⁵ and L_FIR[0.7]=1.6 × 10³⁶ W. We then use the relation $log\left[\frac{L_{\mathrm{X}}}{L_{\mathrm{IR}}}\right]>-4.3$ (Symeonidis et al. 2014), where L_IR is the infrared luminosity and L_X is the soft X-ray (0.5−2 keV) luminosity (as CIE emission is only seen in the RGS band), to infer an X-ray luminosity based on the IR luminosity. We therefore estimate a lower limit for the X-ray luminosity to be within the range L_X >  4.46 × 10³¹ − 7.82 × 10³¹ W. The measured (0.5−2 keV) luminosities of the CIE components in 2000 and 2014 are $L_{\mathrm{CIE}}^{2000}=4.{89}_{-0.11}^{+0.14}\times {10}^{33}$ W and $L_{\mathrm{CIE}}^{2014}=4.{04}_{-0.07}^{+0.11}\times {10}^{33}$ W, respectively, which are significantly larger than the X-ray luminosities inferred from the FIR. To double check this, we use the relation L_X = L_FIR10^−3.70 (Eq. (8) from Ranalli et al. 2003) and obtain soft X-ray luminosities of L_X = 1.78 × 10³² − 3.19 × 10³² W, but these are still an order of magnitude less than L_CIE.

With the CIE luminosity in both epochs being much larger than the soft X-ray luminosity limit inferred from the FIR luminosity, we conclude that either something else is contributing to the CIE X-ray luminosity, or some other process is producing this emission. The SFR values estimated from L_CIE in both epochs, using the relation SFR = 2.2 × 10⁻³³ L_CIE M_⊙ yr⁻¹ (modified Eq. (14) from Ranalli et al. 2003), are SFR₂₀₀₀ = 11 M_⊙ yr⁻¹ and SFR₂₀₁₄ = 9 M_⊙ yr⁻¹. These values are significantly larger compared to the SFR from the FIR (García-Burillo et al. 2014), again suggesting that the CIE plasma is unlikely to come from the SBR.

Alternatively, the CIE emission could originate from the secondary region, north-east of the central black hole (Brinkman et al. 2002), assuming CIE plasma is the cause of the Fe XVII lines. The secondary region is bright in both radio and X-rays, possibly due to old jet material (Wilson & Ulvestad 1987). Due to the field of view and spatial resolution of RGS, we cannot resolve the SBR and the primary and secondary regions with XMM-Newton, hence we see all the emission from each of these three regions combined. Therefore, the ionisation and outflow velocities we obtain from our PION modelling may also apply to emission from the secondary region, harbouring similar plasma to the primary region; however, we do not expect the location of this plasma to be as far out as the secondary region (from our distance estimates in Sect. 7.1.5). It is possible, therefore, that CIE plasma may be located in the secondary region if there are shocks generated in the gas at times when the nuclear jet is particularly strong.

8. Conclusions

In this paper, we investigate the type 2 AGN NGC 1068, comparing XMM-Newton observations in 2000 and 2014, modelling the RGS and EPIC-PN spectra simultaneously. We find that the majority of the emission lines are explained by photoionised plasma, while at least two spectral features require a collisionally ionised component. The main conclusions from our analysis are detailed below.

1. Four photoionised components (PION in SPEX) are required to fit the majority of the emission lines in both epochs, across the full X-ray band. We find that there is very little change in the ionisation states (log ξ = 1−4) of these components between 2000 and 2014, implying the four plasma regions are the same.

2. We compare distance methods to get an insight into the possible locations of the emitting plasma regions within NGC 1068. Although we cannot constrain the estimates, these methods are useful in obtaining rough locations for the outflowing wind, even if there are no spectral changes between observations.

3. We investigate the thermal properties of emitting plasma regions, studying their thermal stability and heating and cooling processes. Using the thermal stability curve (Fig. 8) we find that the highest ionisation state plasmas in each epoch are thermally unstable, while the other three PION components are stable and in thermal equilibrium. We also conclude that the highest and lowest ionisation phases cannot be part of the same regions in the outflowing wind.

4. A possible scenario for obscured AGN where radiation pressure compression (RPC) in the outflowing gas leads to a well defined ionisation distribution and differential emission measure was proposed by Bianchi et al. (2019). The multi-phased PIE plasma modelled here follows the trend of decreasing emission measure for increasing ionisation parameter, which traces the expected slope of a RPC gas. This result suggests that the plasma surrounding the central black hole has properties consistent with a RPC gas, with X-ray and optical emissions coming from the same plasma region, depending on the ionisation and density. However, forthcoming high-throughput and high spectral resolution instrumentation will be able to distinguish between a constant gas pressure multiphased plasma, which we consider here, or one within a RPC scenario, using density diagnostics (Bianchi et al. 2019).

5. From our photoionisation modelling, the ionisation parameter (ξ) and velocity broadening (v_turb; Table 4) values, the distances (Fig. 7), and the thermal stability properties (Fig. 8) of EM1 and EMA suggest that these highest ionisation components are very close to the central SMBH, compared to the other components. This indicates that EM1 and EMA could be part of the BLR. We cannot conclude this with certainty as the BLR is obscured by the dusty torus, although we do not fully rule out this idea as it is possible that the X-ray BLR emission is scattered into our LOS.

6. We find evidence for emission from collisionally ionised plasma in NGC 1068, as photoionised plasma (PION) is unable to account for the Fe XVII lines at 15 and 17 Å. However, photoexcitation would be a more natural way to produce these lines in the outflowing gas surrounding the central black hole (Kinkhabwala et al. 2002), and the lack of PION emission in these lines may be due to incomplete information regarding photoexcitation in SPEX. However, if CIE emission was the cause of the Fe-L lines, we find no conclusive evidence for the CIE emission originating from the SBR or the secondary region. This suggests that CIE may not be the underlying emission component to produce these lines, but instead only statistically accounts for the lines PION cannot produce, where photoexcitation is not accounted for.

NGC 1068 has changed very little, if at all, in the 14 year period between XMM-Newton observations. Future monitoring of NGC 1068 would be important in assessing if this ionisation state of the plasma is maintained, or if this AGN is dynamic in a similar way to type 1 AGN. Multiwavelength observations are paramount in establishing as many intrinsic properties as possible of this obscured AGN in order to constrain and develop the underlying SED for accurate photoionisation modelling.

Acknowledgments

We would like to thank Jelle Kaastra for his helpful conversions and the anonymous referee for their useful and insightful comments. This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. S.G.W. acknowledges the support of a PhD studentship awarded by the UK Science & Technology Facilities Council (STFC) under the grants ST/R505171/1 and ST/M503861/1. M.M. is supported by the Netherlands Organisation for Scientific Research (NWO) through the Innovational Research Incentives Scheme Vidi grant 639.042.525. S.B. acknowledges financial support from the Italian Space Agency under grant ASI-INAF 2017-14-H.O.

References

All Tables

All Figures

Fig. 1. Comparing the Fe Kα line in 2000 (red crosses) and 2014 (blue crosses). Top: there is a clear shift in the Fe Kα line between these two epochs, denoted by the red and blue dashed vertical lines at the peak energies of the observed data. The magenta arrow signifies this energy shift between the observations, equivalent to ΔE = 50 eV. This shift is due to uncertainties in the long-term charge transfer inefficiency and instrumental gain (Cappi et al. 2016). Bottom: we account for this shift manually by multiplying the energies (PI) of the event lists in each 2014 observation by 0.993. Now the 2014 PN data are consistent with the 2000 data and the correction is sufficient for modelling.

In the text

	Fig. 2. Comparing the RGS spectra (left) and EPIC-PN (right) spectra of NGC 1068 in 2000 (red) and 2014 (blue), in the observed reference frame. All spectra show a striking similarity in shape, emission features and flux, suggesting no change between the two epochs.
In the text

	Fig. 3. Combined RGS and PN spectra for the 2000 (blue and red) and 2014 (green and purple) observations, in the observed reference frame. We simultaneously fit both RGS and PN spectra for each epoch to model all emission features and observed continuum over the 0.35−10 keV energy range.
In the text

Fig. 4. RGS spectrum (grey points) between 8 and 20 Å from the 2000 observation of NGC 1068. The red line shows the best fit model and the blue line displays the photoionisation model made up of four different PION components, fitting the majority of the emission lines in both the RGS and PN spectra, except for the features at 15 and 17 Å. We therefore introduced a collisionally ionised component (CIE; orange line) to account for these emission lines from Fe XVII.

In the text

Fig. 5. Best fit model to the 2000 spectrum of NGC 1068. We plot the soft X-ray (RGS) band in wavelength units and the hard X-ray (EPIC-PN) band in energy units for display purposes. The red line shows the best fit to the data points (grey crosses) and the other coloured lines (labelled in the legends) represent each component in the model. Bottom panels: residuals between the best fit model and the observed data points.

In the text

Fig. 6. Best fit model to the 2014 spectrum of NGC 1068. We plot the soft X-ray (RGS) band in wavelength units and the hard X-ray (EPIC-PN) band in energy units for display purposes. The red line shows the best fit to the data points (grey crosses) and the other coloured lines (labelled in the legends) represent each component in the model. Bottom panels: residuals between the best fit model and the observed data points.

In the text

Fig. 7. Distances of each component are plotted for the different methods in Sect. 7.1, at the respective ionisation parameters for 2000 (top) and 2014 (bottom). The data point shapes are shown in the legend, as are the colours for each component. The dashed line in each plot represents the location of the torus, at rtorus = 3.5 pc (García-Burillo et al. 2016), while the dotted line is at 1000 gravitational radii ($R_{\mathrm{g}}=\frac{2G{M}_{\mathrm{BH}}}{c^2}$) from the SMBH, corresponding to a typical outer disk radius (Netzer 2013).

In the text

Fig. 8. Top panel: electron temperature (Te) of the plasma as a function of ionisation parameter (ξ). The labelled circles correspond to the four PION components for each epoch (red for 2000 and blue for 2014) at their respective ionisation parameter values from Table 4. Bottom panel: thermal stability curve of electron temperature (Te) as a function of the pressure form of the ionisation parameter (Ξ). On the curve, the rate of cooling equals the rate of heating. Left of the curve, cooling dominates over heating and to the right, heating dominates over cooling. The labelled circles correspond to the Ξ values of the four PION components for each epoch: red for 2000 and blue for 2014.

In the text

Fig. 9. Emission measure (EM) of the four PION components in each epoch plotted versus their respective ionisation parameter (log ξ). The purple squares are the ions from Fig. 3 in Bianchi et al. (2019), whereby the EM values decrease for increasing log ξ, consistent with the PION components (same colours and data point styles from Fig. 7). We also plot the expected radiative pressure compression (RPC) curve (brown line) that follows the trend of both the ions (purple squares) and PION components suggesting that the photoionised plasma within NGC 1068 is consistent with RPC gas.

In the text

Fig. 10. Top panel: total heating rate as a function of ionisation parameter (ξ) for NGC 1068 (black solid line). The contributions from the different heating processes are also shown: Compton scattering (CS; green dot-dashed line), photoelectrons (PE; purple dotted line), and Auger electrons (AE; orange dashed line). Bottom panel: total cooling rate as a function of ξ for NGC 1068 (black solid line). The contributions from the individual cooling processes are also displayed: inverse-Compton scattering (ICS; orange dot-dashed line), Bremsstrahlung (FF; green dashed line), collisional excitation (CE; purple dot-dashed-dotted line), and recombination (RE; magenta dotted line). Both panels: the labelled circles correspond to the four PION components in each epoch, shown as red for 2000 and blue for 2014, at their respective ionisation parameters from Table 4.

In the text