COMPTON-THICK AGN IN THE NuSTAR ERA VII. A JOINT NuSTAR, Chandra and XMM-Newton ANALYSIS OF TWO NEARBY, HEAVILY OBSCURED SOURCES

We use Chandra archival observations with an exposure time of $\sim$ 29 ks for NGC 3081 and $\sim$ 10 ks for ESO 565-G019. The images of NGC 3081 and ESO 565-G019 in the 0.3–7.0 keV energy range are shown in Figure 1. The sources in the Chandra images are not point-like, given how the emission is extended beyond the Encircled Energy Fraction (EEF, i.e., the circular region containing a certain fraction of the counts) radius. This is due to the excellent angular resolution of the Chandra telescope, that allows to distinguish the nuclear emission from the extended one. The spectra extraction has been done using the specextract task. This task requires the selection of an extracting region for the sources and for the background (the background region has to be unaffected by the presence of other sources). The source extraction has been chosen on the basis of the EEF for a point like source, at a fixed energy (in the case of the Chandra HRMA a circle with a radius of about 2 arcsec contains $90\%$ of the total energy at 5 keV), so to minimize the contamination from non-nuclear emission we have chosen the energy centroid of the Chandra image in the $E=2-7$ keV to better define the AGN. The source regions have been chosen with a radius of 2” and the background regions have a radius of $\sim 25"$.

Finally we bin the spectra with the grppha task to have at least 15 counts per bin to apply $\chi^{2}$ statistics. Because of the small number of spectral counts, in the case of ESO 565-G019, we have used C-stat instead of $\chi^{2}$ statistics in this analysis (i.e., the spectral analysis has been carried out in the Poisson regime).

The XMM-Newton observations, which are quasi simultaneous to the NuSTAR ones (see Table 1), have an exposure time of $\sim$ 30 ks and $\sim$ 27 ks, for NGC 3081 and ESO 565-G019, respectively.

We create the lightcurve at energies $E>10$ keV to select a threshold to remove part of observation particularly affected by noise (i.e., particle background noise). In the case of pn, we select a value of $0.4$, and $0.35$ counts/s for MOS1 and MOS2 for ESO 565-G019 and, in the case of MOS2, we choose $0.25$ counts/s in order to remove a bright background flare for NGC 3081. In both cases the lightcurve were extracted from the whole field. We use extraction regions corresponding to an aperture which contains the $90\%$ of total energy at 5.0 keV: $40"$ in the case of MOS1, $45"$ for MOS2 and $35"$ in the pn case. For both sources, the background regions have sizes of 70, 90 and 60 arcseconds, for the MOS1, MOS2 and pn, respectively.

Finally, we grouped the spectra to have at least 20 counts per bin, in order to apply the $\chi^{2}$ statistics.

The NuSTAR observations have exposures of $\sim$ 56 ks and $\sim$ 50 ks, for NGC 3081 and ESO 565-G019, respectively. The first step of NuSTAR data reduction is the creation of the calibrated events files which will be cleaned and used to produce an exposure map. This process can be carried out using nupipeline.

The choice of the source extraction regions was made selecting four circles with different radii and inspecting the background counts and signal-to-noise ratio (S/N) variations. In the first case (NGC 3081), the data show a linear increase of the S/N with the region diameter (almost until 80”, with increasing background contribution), so the choice fell on the Half Power Diameter (HPD) which contains $50\%$ of the encircled energy fraction and corresponds to a radius of 60”. In the second case (ESO 565-G019), exceeding 40” leads to an increase of the background contribution. Lastly, we extract the spectra, producing the ARF and RMF matrix, and we bin them to have at least 20 counts per bin.

In table 2 we show the spectral information related to the extracted spectra.

Finally, we find no significant evidences of variability between the Chandra, XMM-Newton and NuSTAR observations.

In this section we will discuss the main properties and the use of the MYTorus model [35]. The MYTorus model was developed to be used in the XSPEC [5] environment as a combination of additive and multiplicative tables, which represent different components of the nuclear emission. These components are:the zeroth-order emission component (MYTZ), the scattered continuum (MYTS) and the iron line emission component (MYTL). MYTorus models the observed spectrum taking into account the absorbed and the scattered component of the emission, modeling also the presence of fluorescent Fe K$\alpha$ and K$\beta$ emission lines at $6.4$ keV and $7.06$ keV respectively, which are thought to be almost ubiquitous in heavily obscured AGN spectra.

The MYTorus model can be used in two different settings, namely the coupled and the de-coupled configuration. In the first mode, the column density and the inclination angle of the three components (MYTZ, MYTS and MYTL) are tied together. Thus, in this configuration, all the components are produced in the same medium.

The MYTorus model simulates the interaction between input spectrum photons and the obscuring material. The circumnuclear environment is simulated as the classical doughnut-like and azimuthally symmetric structure. The distance from the black hole to the center of the torus section is indicated as $c$, and $a$ is the radius of the section.

$\theta_{obs}$ is the inclination angle: the angle between the torus symmetry axis and the observer line of sight (l.o.s.). It can vary in the range [0^∘ - 90^∘], allowing to reproduce both the face-on ($\theta_{obs}$=0^∘, i.e., the observer looks directly at the nucleus) and the edge-on ($\theta_{obs}$=90^∘,i.e., the observer line of sight intercepts the torus equator).

The torus half-opening angle, which represents the fraction of the sky as seen from the center, is defined as $\alpha=[(\pi-\psi)/2]$=60^∘ (with $\psi$ being the angle subtended by the internal surface of the torus) corresponding to a covering factor C${}_{TOR}=0.5$. The fixed value for the covering factor is linked to the assumptions made on the fraction of obscured AGN with respect to the unobscured ones. Finally, N_H is the equatorial column density (i.e. the column density through the torus diameter); the line of sight column density can be computed as:

The first component in the MYTorus model is the so-called zeroth-order continuum or direct component. This component represents the photons escaping the absorbing medium (i.e., the torus) without being absorbed or scattered.

The second component is the scattered or reprocessed continuum, which represents the photons that escape the medium after being scattered one or more times.

The interaction is via Compton-scattering, thus the energy of the photon after the scattering will be lower with respect to the input photon’s energy. The second component is responsible for the production of the feature observed at $\sim$ 30 keV (i.e. the Compton hump). Moreover, in the MYTorus model the termination energy of the scattered component is variable between 160 keV and 500 keV (in our analysis we use a table with intrinsic continuum extending up to 500 keV). The value of the cut-off energy has been chosen to be consistent with previous similar works [33, 54, 52, see e.g.,]. Moreover, recent works [10, e.g.,] show that it is a reasonable value for the extension of the continuum.

The last component is the fluorescent emission. It takes into account the possibility to have fluorescent emission iron lines, produced in the reprocessing medium. The emission lines that MYTorus models are the Fe K$\alpha$ and K$\beta$ only.

On the one hand, the line photons that escape after being produced by the fluorescence process, constitute the zeroth-order fluorescent emission component. On the other hand, if these photons interact with the reprocessor by scattering processes, they can contribute to form the Compton shoulder.

The MYTorus model in “coupled” configuration allows the inclination angle to vary, but does not permit a variation in the column density and in the geometrical properties of the different components. In this way it is not possible to properly characterize a clumpy torus structure. To overcome this problem, as described by [50], it is possible to decouple the MYTorus components by fixing the zeroth-order continuum inclination angle to 90^∘, generating a pure line of sight component. Then, the column density of the scattered component can be untied with respect to the one of the direct continuum. In this way, the direct continuum column density represents the line of sight column density, whereas the scattered component column density represents the “global average” column density. Thus, the ratio between the “global average” column density and the line of sight column density represents a measure of the patchiness (or clumpiness) of the obscuring material: a ratio $\not=1$ will then suggest a scenario in which the column density along the l.o.s. is higher (or lower) than the average column density of the torus, meaning that the structure could likely be clumpy rather than smooth. Following [51], we then fix the inclination angle of the scattered and fluorescent line components to be either $\theta_{S=L}$=90^∘ or $\theta_{S=L}$=90^∘, reproducing an edge-on and face-on geometry. Using the “de-coupled” mode, it is possible to take into account a scenario in which the different MYTorus components are produced with different interactions of the nuclear emission with the reprocessor. Using $\theta_{S=L}$=0^∘ we model a scenario in which the emission is dominated by reflection in the far-side of the torus; if $\theta_{S=L}$=90^∘, the emission is dominated by a near-side Compton-scattering.

Finally, we use the borus02 model (table borus02_v170323a.fits), developed by [9] as an improvement of the BNtorus model [12]. The borus02 model is based on grids of spectral templates obtained using Monte-Carlo simulations of radiative transfer through a neutral spherical torus with polar cut-outs.

The strength of this model lies in the possibility to fit the spectral data having as free parameters the average column density of the torus and its covering factor. The computation of the covering factor it is not possible by using MYTorus  even in its “de-coupled” configuration. Due to the longer variability timescales (years), the average column density represents a more reliable parameter to characterize the thickness of an AGN in respect to the N_H,los, which variability has shown to be of the order of days and weeks, due to the movement of clouds through the line of sight [44, 42, see e.g.,].

Despite this, a proper characterization of the covering factor is not simple. It can be affected by accretion or feedback phenomena taking place nearby the torus [26, 36, e.g.,]; it can depend on the luminosity [7, e.g.,] as well as on the Eddington ratio [41, e.g.,], and these dependencies could vary with redshift [14, e.g.,]. In this perspective, the borus02 model represents an updated tool to compute the covering factor. borus02 model is composed of a single additive table (instead of the three tables of MYTorus ) which takes into account the reprocessed emission component, that is similar to the MYTorus “reprocessed component”, and the fluorescent line emission component, including K$\alpha$ and K$\beta$ lines.

The main parameters of the borus02 model have the following possible values: the covering factor ranges from $0.1$ to $1$, corresponding to a torus opening angle between 84^∘ and 0^∘; the inclination angle is in the range [18^∘-87^∘]. Also, the cut-off energy is a parameter of the model and we fix it to be 500 keV, for consistency with MYTorus. Finally, the iron abundance is also a free parameter, but we fix it to 1, for consistency with the MYTorus analysis.

borus02 does not include the l.o.s. absorption at the redshift of the source, which we model with combination of XSPEC components $zphabs*cabs$ in order to take into account l.o.s. absorption and the losses out of the l.o.s. due to Compton scattering. The primary power law emission is represented by cutoffpl1 that is multiplied by the previous expression; it is characterized by a photon index, cut-off energy and normalization, that must be tied to those of borus02. Also, to properly describe the l.o.s. column density, the $nH$ parameter of cabs and zphabs must be tied together. The soft emission component and the emission lines are included, when needed, as described for the MYTorus modelling.

We have analyzed the “soft” Chandra and XMM-Newton spectra alone in order to quantify the effect of adding the information from the NuSTAR data. From the 0.6–10 keV analysis of NGC 3081, we obtain a best fit model with $\Gamma=1.91_{-0.16}^{+0.15}$ and $N_{H,z}=0.74_{-0.07}^{+0.06}\times 10^{24}$ cm^-2. The spectral slope is in agreement with typical values observed in AGN. However, in order to increase the statistics and properly constrain the spectral parameters (in particular, the photon index and the column density), we have combined the NuSTAR data with the “soft” spectrum, obtaining the 0.6–70 keV NGC 3081 spectrum. As expected, the uncertainties on the main spectral parameters significantly decrease: the errors associated to the photon index decrease from $\sim 8\%$ to $\sim 4\%$ and those on the l.o.s. column density decrease from $\sim 8\%$ to $\sim 3\%$. Also, in accordance with [30], we find a shift in both the photon index and the l.o.s. column density values. The first is reduced by $\sim 5\%$ [30, the average decrease in $\Gamma$ value found by] is $\sim 13\%$ and the $N_{H,z}$ one decreases by $\sim 19\%$ [30, the average value in ] is $\sim 32\%$; however, several sources in that sample only had low-count statistic Swift-XRT coverage in the 0.5–10 keV band). The smaller errors allow us to break the degeneracy between $\Gamma$ and N_H,z. In Figure 6 we show the comparison between the 0.6–10 keV and 0.6–70 keV $\Gamma$-N_H,z confidence regions. It is clear that, when adding NuSTAR data to the 0.6–10 keV spectrum, there is a shift in the spectral parameters to lower values and, also, a significant decrease of their uncertainties; this result highlights the strength of the X-ray broad-band approach to characterize candidate CT-AGN, and the key role played by NuSTAR to achieve this goal.

We investigate possible variability between the Chandra data and the NuSTAR + XMM-Newton data. While NuSTAR and XMM-Newton observations are simultaneous, Chandra targeted NGC 3081 about one year before. In Figure 6 we show the contour plot between the normalization (indicator of the flux) and the column density (which accounts for the absorption properties). Variability can either be due to intrinsic variation of the emission from the central engine or of the absorbing structure. If there were a difference in both the normalization and in the column density (between the two sets of spectra), we may affirm that it can be due to variation in the geometrical properties of the torus through time (or variation of the accretion efficiency in the case of the normalization). In Figure 6 the superposed Chandra and the XMM-Newton + NuSTAR Normalization-N_H contour plots are shown. The XMM-Newton + NuSTAR contours plot (solid lines) is much smaller than the Chandra one (dashed lines); however, they are consistent, and no variability effects can be attested.

We also search for variability effects, for ESO 565-G019, between the Chandra and XMM-Newton+NuSTAR observations through the Normalization-N_H,z contour plots. We do not find indications of variability, although we note that the errors on the parameters, due to the limited photon statistics, are large.

Given the presence of Chandra observations, it is possible to study the properties of the extended emission in both NGC 3081 and ESO 565-G019, thus to establish a better portrait of the soft band spectrum.

Following the approach used in [21, 27, 29], we obtain the soft (i.e. 0.3–3.0 keV) and hard (i.e. 3.0–7.0 keV) Chandra images of both sources (see Figures 7 and 8). We can notice that the diffuse emission extends in the NW-SE and N-S direction, showing an elongated structure on projected scales of about 2 kpc and 3.5 kpc for NGC 3081 and ESO 565-G019, respectively.

To quantitatively assess the presence (or lack) of extended emission in our sources, we compare the radial distribution of their surface brightness with the one obtained from a simulated PSF in the two energy ranges using the ChaRT and MARX 5.5.1 tasks [see e.g., 21, 27, 29, for a detailed description of this technique]. Figure 9 shows the radial profiles of the emission vs. PSF expectations in the 0.3–3.0 keV and 3.0–7.0 keV energy bands. These profiles have been obtained from an annular region comprising 8 annuli from $\sim 0.7$ to $\sim 8$ arcsec. The 0.3–3.0 keV profiles of NGC 3081 and ESO 565-G019 show a significant emission above the PSF values up to $\sim 7$ arcsec, whose origin could be related to non-nuclear processes like star formation or diffuse emission on the host-galaxy scales. In Figure 10 the Chandra contours plotted over the optical DSS images are shown. As mentioned in Section II, ESO 565-G019 has a SFR on the order of $\sim$3–4 $M_{\odot}/yr$ [23], thus, its X-ray diffuse emission could be ascribed to a thermal emission on the scales of the host with a possible contribution of star formation processes. Given the 0.5–2 keV spectrum, it is possible to compute the X-ray SFR for these sources following the relation between the 0.5–2 keV luminosity and the SFR [40, e.g.,]. We find that ESO 565-G019 has a SFR$=4.4_{-0.7}^{+0.8}$ M${}_{\odot}/yr$ (the errors correspond to the dispersion of the relation we adopted), consistent with literature. For NGC 3081 there is little [20, nuclear SFR $<$0.05 M_⊙/yr,] to no evidences for nuclear star formation activity [19, 22, see, e.g.,], therefore the diffuse emission could be produced by galaxy-scale processes (e.g., hot gas in the nuclear region of the galaxy). However, we cannot rule out the possibility that star formation also contributes to some extent. Indeed, from the 0.5–2 keV spectrum, we find a SFR$=1.2_{-0.2}^{+0.2}$ M${}_{\odot}/yr$.

In the hard-band profiles, where the AGN contribution is dominant and negligible contamination from non-AGN processes is expected, this extension is much more reduced, especially for ESO 565-G019. We also measure the excess fraction, defined as the ratio between the counts above the PSF and the total counts in the analyzed area. We find an excess fraction of $20\pm 1.8\%$ for the 0.3–3.0 keV extended emission of NGC 3081 and $18.9\pm 8.7\%$ for ESO 565-G019. In the 3.0–7.0 keV range, the excess fraction is negative for both NGC 3081 ($<-0.01$) and ESO 565-G019 ($<-0.07$), meaning that the emission is consistent with the PSF in the 3.0–7.0 keV energy range.

The detection of a diffuse axial emission in the 0.3–7.0 keV interval is in agreement with the aforementioned works [21, 27, 29]. However, although we detect an excess in the soft data, we find no significant excess in the hard extended emission, which has been found to be 12–22% and likely due to the existence of reprocessed emission on scales beyond the torus [29, see e.g.,].

Based on the results of this work along with those reported by Torres-Albà et al. (submitted), eight out of ten sources are incompatible, at 90% significance, with having the same line-of-sight and average torus column densities. This follows the overall trend observed for the full sample of 57 obscured AGN, in which the large majority ($\sim$91%) of sources show this discrepancy (see Fig. 4 in Torres-Albà et al., submitted). We link this observational evidence to the presence of a clumpy torus. Moreover, including the results from the 10 sample analysis, 13 out of 57 candidate CT-AGN (one of which is ESO 565-G019) have both l.o.s and average column density larger than $10^{24}$ cm^-2. With the addition of this source, the percentage of NuSTAR-confirmed CT-AGN in the BAT sample at $z\leq 0.05$ is $\sim$8% (32/417) ¹⁵¹⁵ 15 https://science.clemson.edu/ctagn/ctagn/, still much lower than the predictions. However, with the analysis of the full $\sim$60 source sample, Torres-Albà et al. (submitted) observed a significant decrease of the Compton-thick fraction with redshift. In fact, at $z<0.01$, the fraction is $20.0\pm 5.7\%$, which is much closer to the theoretical predictions.

Despite the many studies carried out over the last 20 years, the obscured AGN population is not entirely well characterized and many questions are still unresolved. The study of obscured AGN is relevant for several astrophysical issues concerning galaxy evolution, as well as the CXB content and determination of the accretion history of the Universe [2, 24, 45, e.g.,]. In fact, the quest for and characterization of Type 2 AGN provides a census of the population of galaxies which are thought to be in the phase of the building up of their mass (probably after a merger) or in the first phases of the nuclear activity [16, 8, see, e.g.,].

Combining the capabilities of Chandra and XMM-Newton at E$<10$ keV with those of NuSTAR (3–79 keV), it is possible to properly characterize the spectral parameters (e.g., intrinsic column density and opening angle of the obscuring material) of heavily obscured AGN, allowing to distinguish between the Compton-thin and Compton-thick regime. These combined observations can also break the degeneracy between spectral parameters (e.g., between the photon index and the obscuring material column density) through the use of advanced, physically motivated spectral models (e.g., MYTorus and borus02 ).

In Figure 11 a schematic sketch of the best-fit configuration found for the two sources is reported. The torus is represented by several clumps distributed following the toroidal structure. In the figure the following properties are present: the covering factor is represented as the angle of the sky free from the torus by the central source point of view; the differences in the column density of the torus is represented by differences in the colors of the clumps, darker colors mean higher column density. Thus, NGC 3081 has an higher covering factor (lower angle) than ESO 565-G019. The l.o.s. column density of NGC 3081 is lower than the “global average” column density: the observer is looking at the central source through an under-dense region with respect to the torus average column density. ESO 565-G019 is characterized by a Compton-thick column density both on the l.o.s. and in average, and this is represented by over-dense clumps in the whole structure.

The conclusions of this work can be summarized as follows:

1. We have verified that the NuSTAR data significantly contribute to the determination of the main spectral parameters which characterize obscured AGN at low redshift ($z\leq 0.1$). Its contribution mainly consists in the decreasing of the errors of the parameters of interest and in breaking of degeneracy between them, thus allowing a better characterization of the spectral emission properties. Moreover, the use of XMM-Newton and NuSTAR simultaneous observations allows one to avoid variability effects.

2. The spectra of both sources present a significant emission at energies $<2-3$ keV which cannot be ascribed solely to the main emission from the nucleus. We have modeled this emission with a thermal component and found that it can be produced by a medium with temperature between $10^{6}-10^{7}$ K. However, we do not place constraints on the origin of this emission. This thermal component can be due either to a thermally emitting gas in the nuclear region or to a population of X-ray emitting unresolved sources (e.g., X-ray binaries) or to a combination of the two phenomena. In the soft part of the spectrum of NGC 3081, we find evidences of emission lines, which are typical of obscured sources [13, 39, e.g.,].

3. We found that NGC 3081 is best fitted by the MYTorus model in the de-coupled mode (edge-on configuration), while ESO 565-G019 is best fitted by the borus02 model. NGC 3081 is Compton-thin along the l.o.s., but with the obscuring material being, on average, Compton-thick. ESO 565-G019 is classified as Compton-thick in both the l.o.s. and average components of the column density.

4. For both sources we were able to compute the torus covering factor through the borus02 modeling. NGC 3081 has C${}_{TOR}=0.73_{-0.10}^{+0.09}$, suggesting that the torus contributes in a significant way to cover the central emission. Moreover, the ratio between l.o.s. and average column densities is typical of a clumpy scenario. For ESO 565-G019 the lower covering factor C${}_{TOR}=0.47_{-0.08}^{+0.18}$, along with the $\Delta$N_H, suggests a scenario in which the obscuring structure, which is Compton-thick, may be distributed in several individual clouds, responsible for the obscuration. The values are consistent, between the uncertainties, with the average covering factor found by [31], C${}_{TOR}\simeq 0.6$, for a $\sim$ 30 CT-AGN candidates sample.

5. The main nuclear emission can be divided into the reprocessed component, which is heavily suppressed in Compton-thick AGN, and the component which is scattered (and unabsorbed) by Compton-thin material and can reach the observer at lower energies (i.e. $<5$ keV). For both AGN presented in this work we find that this component represents a low fraction of the main emission, being $<1\%$ and $\sim 1\%$ for NGC 3081 and ESO 565-G019, respectively. These results suggest that the Compton-thin material is a small fraction of the circumnuclear environment.

6. Thanks to the presence of Chandra data, we were able to investigate the extended emission of NGC 3081 and ESO 565-G019. We found significant diffuse emission in the 0.3–3.0 keV band extending for about 2 (NGC 3081) and 3.5 kpc (ESO 565-G019). However, we were not able to detect any diffuse emission in the 3–7 keV energy range.