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 ∼29 ks for NGC 3081 and ∼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, which allows us to distinguish the nuclear emission from the extended emission. 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'' 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 ∼25''.

Zoom In Zoom Out Reset image size

Finally we bin the spectra with the grppha task to have at least 15 counts per bin to apply χ² statistics. Because of the small number of spectral counts, in the case of ESO 565-G019, we have used C-stat instead of χ² 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 ∼30 and ∼27 ks, for NGC 3081 and ESO 565-G019, respectively.

We create the light curve 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 light curves were extracted from the whole field. We use extraction regions corresponding to an aperture that 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'', 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 χ² statistics.

The NuSTAR observations have exposures of ∼56 and ∼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 that 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, 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. Last, 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 (Murphy 2009). The MYTorus model was developed to be used in the XSPEC (Arnaud 1996) 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α and Kβ 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 decoupled 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.

The inclination angle is θ_obs, 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 the model to reproduce both the face-on (θ_obs = 0°, i.e., the observer looks directly at the nucleus) and the edge-on (θ_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 α = [(π − ψ)/2] = 60° (with ψ 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 AGN. 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:

$$\begin{eqnarray*}&&{N}_{{\rm{H}},z}={N}_{{\rm{H}}}\left[1-{\left(\frac{c}{a}\right)}^{2}{\cos }^{2}{\theta }_{\mathrm{obs}}\right].\end{eqnarray*}$$

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 ∼30 keV (i.e., the Compton hump). Moreover, in the MYTorus model the termination energy of the scattered component is variable between 160 and 500 keV (in our analysis, we use a table with intrinsic continuum extending up to 500 keV). The value of the cutoff energy has been chosen to be consistent with previous similar works (see, e.g., Marchesi et al. 2017b; Zhao et al. 2019b, 2019a). Moreover, recent works (e.g., Baloković et al. 2021) 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 of having fluorescent emission iron lines, produced in the reprocessing medium. The emission lines that MYTorus models are the Fe Kα and Kβ only.

The line photons that escape after being produced by the fluorescence process constitute the zeroth-order fluorescent emission component. 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 Yaqoob (2012), 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 ≠ 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 Yaqoob et al. (2015), we then fix the inclination angle of the scattered and fluorescent line components to be either θ_S=L = 90° or θ_S=L = 0°, reproducing an edge-on and face-on geometry. Using the decoupled 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 θ_S=L = 0° we model a scenario in which the emission is dominated by reflection in the far-side of the torus; if θ_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 Baloković et al. (2018) as an improvement of the BNtorus model (Brightman 2011). 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 cutouts.

The strength of this model lies in the possibility of fitting the spectral data having as free parameters the average column density of the torus and its covering factor. The computation of the covering factor is not possible by using MYTorus even in its decoupled configuration. Due to the longer variability time scales (years), the average column density represents a more reliable parameter to characterize the thickness of an AGN in respect to the N_H,los, whose variability has shown to be of the order of days and weeks, due to the movement of clouds through the line of sight (see, e.g., Risaliti et al. 2002; Ricci et al. 2016).

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 (e.g., Heckman & Best 2014; Netzer 2015); it can depend on the luminosity (e.g., Assef et al. 2013) as well as on the Eddington ratio (e.g., Ricci et al. 2017), and these dependencies could vary with redshift (e.g., Buchner et al. 2015). In this perspective, the borus02 model represents an updated tool to compute the covering factor. The borus02 model is composed of a single additive table (instead of the three tables of MYTorus) that takes into account the reprocessed emission component, which is similar to the MYTorus "reprocessed component," and the fluorescent line emission component, including Kα and Kβ 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 cutoff 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.

The borus02 model 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, which is multiplied by the previous expression; it is characterized by a photon index, cutoff 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 modeling.

4.1.1.  MYTorus Model

We use MYTorus in both the coupled and decoupled configurations (the last one in the edge-on and face-on mode).

The best fit model consists of the three MYTorus components, the second power-law component, the mekal component (to model the soft thermal emission), and the emission lines. Moreover, we included other two constants to the models, A_S and A_L , to take into account the possible different normalizations of the other two components with respect to the zeroth-order continuum:

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{Model}\,\mathrm{NGC}\_A=\mathrm{pha}\ast (z\mathrm{po}1\ast \mathrm{MYTZ}+{A}_{S}\ast \mathrm{MYTS}\\ \quad +\ {A}_{L}\ast \mathrm{MYTL}+{f}_{s}\ast z\mathrm{po}2+\mathrm{mekal}+3\ast z\mathrm{gauss})\end{array}.\end{eqnarray} \tag{ 1 }$$

The photon index is ${{\rm{\Gamma }}}_{\mathrm{MYT},c,B}={1.59}_{-0.03}^{+0.03}$. The column density, which is ${N}_{{\rm{H}},\mathrm{eq}}={0.62}_{-0.02}^{+0.02}\times {10}^{24}$ cm⁻², is below the CT threshold. We also fit our data leaving the inclination angle free to vary and we find no significant improvement in the fit statistic (χ²/d. o. f. = 1773/1403). Therefore, we fix the inclination angle to 90°. The best fit model and the spectrum are reported in Figure 2. Also, we show the residuals for the iron Kα line in the best fit model without the line component. As can be seen, the line component is required to improve the fit statistic.

Zoom In Zoom Out Reset image size

In order to reach a more complete understanding of the geometrical properties of the NGC 3081 torus, we fit the 0.6–70 keV spectrum with MYTorus in the decoupled mode. The photon indices we obtain are ${{\rm{\Gamma }}}_{\theta ,S=90}={1.81}_{-0.06}^{+0.07}$ and ${{\rm{\Gamma }}}_{\theta ,S=0}={1.75}_{-0.05}^{+0.03}$, for the edge-on and face-on modes, respectively. Also the l.o.s. column densities are: ${N}_{{\rm{H}},z}\,={0.60}_{-0.02}^{+0.02}\times {10}^{24}$ cm⁻² and ${N}_{{\rm{H}},S}={1.59}_{-0.17}^{+0.19}\times {10}^{24}$ cm⁻² for the edge-on mode and ${N}_{{\rm{H}},z}={0.66}_{-0.02}^{+0.02}\times {10}^{24}$ cm⁻² and ${N}_{{\rm{H}},S}={3.00}_{-0.68}^{+0.69}\times {10}^{24}$ cm⁻² for the face-on mode. In both modes, the photon index is steeper than the one found in the analysis with the coupled configuration (see Table 3); also the l.o.s. column densities are lower than the global average column densities. The χ² statistics favors the edge-on configuration, with a χ²/d. o. f. = 1723/1403. We show in Figure 3 (left panel) the unfolded spectrum and the MYTorus model in the edge-on mode.

Zoom In Zoom Out Reset image size

Table 3. Summary Table of the Spectral Results Obtained with MYTorus (Coupled and Decoupled) and borus02 Applied to NGC 3081 Data

	NGC 3081
	MYTorus	MYTorus		borus02
	coupled	decoupled
		edge-on	face-on
χ2/d. o. f.	1774/1404	1723/1403	1747/1403	1753/1403
Γ	${1.59}_{-0.03}^{+0.03}$	${1.81}_{-0.06}^{+0.07}$	${1.75}_{-0.05}^{+0.03}$	${1.80}_{-0.04}^{+0.06}$
NH,eq a	${0.62}_{-0.02}^{+0.02}$	⋯	⋯	⋯
NH,z b	⋯	${0.60}_{-0.02}^{+0.02}$	${0.66}_{-0.02}^{+0.02}$	${0.61}_{-0.02}^{+0.02}$
NH,S c	⋯	${1.59}_{-0.17}^{+0.19}$	${3.00}_{-0.68}^{+0.69}$	${1.51}_{-0.10}^{+0.11}$
AS = AL d	${0.86}_{-0.12}^{+0.13}$	${2.23}_{-0.37}^{+0.46}$	${0.42}_{-0.05}^{+0.05}$	⋯
θobs e	90f f	90f	0f	⋯
fs g	${5.63}_{-0.56}^{+0.62}$	${4.36}_{-0.62}^{+0.57}$	${3.17}_{-0.34}^{+0.53}$	${4.51}_{-0.55}^{+0.44}$
EW h	${0.181}_{-0.003}^{+0.010}$	${0.185}_{-0.007}^{+0.007}$	${0.185}_{-0.007}^{+0.007}$	⋯
CTOR i	⋯	⋯	⋯	${0.73}_{-0.10}^{+0.09}$
kT j	${0.26}_{-0.01}^{+0.02}$	${0.25}_{-0.02}^{+0.02}$	${0.26}_{-0.02}^{+0.02}$	${0.26}_{-0.02}^{+0.02}$
F2–10 keV k	${4.50}_{-0.13}^{+0.11}$	${4.40}_{-0.19}^{+0.10}$	${4.45}_{-0.17}^{+0.11}$	${4.38}_{-0.14}^{+0.10}$
F10–40 keV l	${4.12}_{-0.08}^{+0.07}$	${4.19}_{-0.13}^{+0.05}$	${4.16}_{-0.28}^{+0.03}$	${4.18}_{-0.08}^{+0.10}$
$\mathrm{log}({L}_{2\mbox{--}10\,\mathrm{keV}})$ m	${41.73}_{-0.01}^{+0.01}$	${41.74}_{-0.05}^{+0.05}$	${41.75}_{-0.06}^{+0.06}$	${42.78}_{-0.01}^{+0.01}$
$\mathrm{log}({L}_{10\mbox{--}40\,\mathrm{keV}})$ n	${42.56}_{-0.02}^{+0.02}$	${42.57}_{-0.01}^{+0.01}$	${42.66}_{-0.07}^{+0.07}$	${42.83}_{-0.01}^{+0.01}$

Notes.

^a Equatorial column density in the MYTorus model, in the coupled configuration, in units of 10²⁴cm⁻². ^b Column density along the l.o.s. in units of 10²⁴cm⁻². ^c Global average column density in units of 10²⁴cm⁻². ^d Normalization between the reprocessed MYTorus component and the zeroth-order continuum. ^e Torus inclination angle in degrees. ^f The f indicates that a parameter is fixed. ^g Fraction of the scattered component. ^h Equivalent width of the Kα Iron line in units of keV. ⁱ Covering factor of the torus, in the borus02 model. ^j Temperature of the thermal component in keV. ^k 2–10 keV flux in units of 10⁻¹² erg s⁻¹ cm⁻². ^l 10–40 keV flux in units of 10⁻¹² erg s⁻¹ cm⁻². ^m 2–10 keV intrinsic luminosity in units of erg s⁻¹. ⁿ 10–40 keV intrinsic luminosity in units of erg s⁻¹.

Download table as:  ASCIITypeset image

4.1.2.  borus02 Model

Finally, we model the 0.6–70 keV combined spectrum with the borus02 model. In addition to the main emission component and the reprocessed component, we also included the second power-law component and the thermal component to model the contribution in the soft part of the spectrum. Also, we add the three emission lines previously described.

$$\begin{eqnarray}&&\begin{array}{c}\mathrm{Model}\,\mathrm{NGC}\_B=\mathrm{pha}\ast (\mathrm{borus}02+z\mathrm{pha}\ast \mathrm{cabs}\\ \quad \ast \mathrm{cutoff}{pl}1+{f}_{s}\ast \mathrm{cutoff}{pl}2+\mathrm{mekal}+3\ast z\mathrm{gauss})\end{array}.\end{eqnarray} \tag{ 2 }$$

The best fit model has a χ²/d. o. f. = 1753/1403; the photon index is ${\rm{\Gamma }}={1.80}_{-0.04}^{+0.06};$ the l.o.s. column density is ${N}_{{\rm{H}},z}={0.61}_{-0.02}^{+0.02}\times {10}^{24}$ cm⁻² and the average column density ${N}_{{\rm{H}},S}={1.51}_{-0.10}^{+0.11}\times {10}^{24}$ cm⁻² is consistent with the torus being CT, as we also found using MYTorus in the edge-on configuration. In Figure 3 (right panel), we show the 0.6–70 keV spectrum. The borus02 model allows us to leave the covering factor as a free parameter; in the case of NGC 3081 we obtain ${C}_{\mathrm{TOR}}={0.73}_{-0.10}^{+0.09}$.

4.1.3. Summary of the NGC 3081 Spectral Analysis Results

4.2.1.  MYTorus Model

We first analyze the ESO 565-G019 0.6–40 keV spectrum with MYTorus in its coupled configuration, using the three MYTorus components plus the second power-law component and the thermal component to describe the soft part of the spectrum:

$$\begin{eqnarray}&&\begin{array}{c}\mathrm{model}\,\mathrm{ESO}\_A=\mathrm{pha}(z\mathrm{po}1\ast \mathrm{MYTZ}+{A}_{S}\ast \mathrm{MYTS}\\ +\ {A}_{L}\ast \mathrm{MYTL}+\mathrm{mekal}+{f}_{s}\ast z\mathrm{po}2)\end{array}.\end{eqnarray} \tag{ 3 }$$

In this case, leaving the inclination angle free to vary, it is not possible to obtain a statistically acceptable solution for the fit (χ² > 2). To this purpose, we tried two different configurations: one with the inclination angle fixed to 90°, and the other with θ_obs = 65° (i.e., we are observing through the brink of the torus). We report the results of the spectral fitting in Table 4: both the photon index and the column density are very different in the two models, in particular we obtain ${N}_{{\rm{H}},\mathrm{MYT},c,65}={10}_{-1.29}^{+0.00* }\times {10}^{24}$ cm⁻² that is the MYTorus upper limit. We also try to fit the data fixing the angle to an intermediate value (θ_obs = 77°), but the statistics do not show significant improvement (${\chi }_{\nu }^{2}=297/212$). We report in Figure 4 the 0.6–40 keV spectrum of ESO 565-G019 fitted with MyTorus coupled with θ=65°.

Zoom In Zoom Out Reset image size

Table 4. Summary Table of the Spectral Results Obtained with MYTorus (Coupled and Decoupled) and borus02 Applied to ESO 565-G019 Data. The parameters are reported as in Table 3

	ESO 565-G019
	MYTorus	MYTorus		borus02
	coupled	decoupled
		edge-on	face-on
χ2/d. o. f.	286/212	292/211	256/211	248/209
Γ	${2.08}_{-0.08}^{+0.07}$	${1.56}_{-* }^{+0.16}$	${2.22}_{-0.08}^{+0.12}$	${1.75}_{-0.22}^{+0.04}$
NH,eq	${5.30}_{-0.68}^{{+}^{* }}$	⋯	⋯	⋯
NH,z	⋯	${2.96}_{-0.38}^{+0.53}$	${5.80}_{-2.53}^{+4.20}$	${3.00}_{-0.69}^{+* }$
NH,S	⋯	${0.35}_{-0.05}^{+0.07}$	${3.30}_{-1.80}^{+1.19}$	${3.72}_{-1.15}^{+* }$
AS = AL	1f	1f	1f	⋯
θobs	65f	90f	0f	⋯
fs	${1.24}_{-0.42}^{+0.35}$	${6.61}_{-1.53}^{+2.03}$	${1.73}_{-0.60}^{+1.66}$	${9.14}_{-4.58}^{+4.09}$
EW	${1.90}_{-0.27}^{+* }$	${1.42}_{-0.77}^{+0.18}$	${1.60}_{-\mathrm{0.0.42}}^{+* }$	⋯
CTOR	⋯	⋯	⋯	${0.47}_{-0.08}^{+0.18}$
kT	${0.59}_{-0.03}^{+0.03}$	${0.57}_{-0.05}^{+0.03}$	${0.59}_{-0.04}^{+0.03}$	${0.59}_{-0.03}^{+0.03}$
F2–10 keV	${0.48}_{-0.10}^{+24.25}$	${0.52}_{-0.10}^{+0.03}$	${0.48}_{-0.09}^{+16.42}$	${0.50}_{-0.39}^{+5.07}$
F10–40 keV	${3.89}_{-0.07}^{+13.42}$	${3.61}_{-0.86}^{+0.30}$	${3.44}_{-0.18}^{+9.51}$	${3.68}_{-2.79}^{+2.53}$
$\mathrm{log}({L}_{2\mbox{--}10\,\mathrm{keV}})$	${43.17}_{-0.06}^{+0.05}$	${42.87}_{-0.03}^{+0.03}$	${43.03}_{-0.05}^{+0.04}$	${42.48}_{-0.63}^{+0.25}$
$\mathrm{log}({L}_{10\mbox{--}40\,\mathrm{keV}})$	${42.22}_{-0.07}^{+0.06}$	${42.12}_{-0.04}^{+0.03}$	${42.04}_{-0.06}^{+0.05}$	${42.56}_{-0.56}^{+0.24}$

Download table as:  ASCIITypeset image

These fitting issues may be due to the limitations of the coupled mode, thus, using MYTorus in its decoupled configuration, we expect to achieve a more physical description of the circumnuclear region for ESO 565-G019.

We then model the spectrum using the MYTorus model in the edge-on and face-on modes of the decoupled configuration and adding both the second power-law and the thermal components to reproduce the soft emission. In the edge-on configuration the photon index is ${{\rm{\Gamma }}}_{\theta ,S=90}={1.56}_{-0.16* }^{+0.16};$ in the face-on mode it is instead steeper ${{\rm{\Gamma }}}_{\theta ,S=90}={2.22}_{-0.08}^{+0.12}$. The l.o.s. column densities are ${N}_{{\rm{H}},z}={2.96}_{-0.38}^{+0.53}\times {10}^{24}$ cm⁻² for the edge-on mode and ${N}_{{\rm{H}},z}={5.80}_{-2.53}^{+4.20* }\times {10}^{24}$ cm⁻² in the face-on configuration. The "global average" column density is ${N}_{{\rm{H}},S}={0.35}_{-0.05}^{+0.07}\times {10}^{24}$ cm⁻² in the edge-on configuration and ${N}_{{\rm{H}},S}={3.30}_{-1.80}^{+1.19}\times {10}^{24}$ cm⁻² in the other configuration. For both modes, the preferred scenario is the one in which we are observing through a particularly dense region of the torus, which has a lower global average column density. In Table 4, we report the spectral parameter for the 0.6–40 keV spectra: as it can be seen, the statistically favored scenario is the face-on scenario, whose best fit model with the spectrum is reported in Figure 5.

Zoom In Zoom Out Reset image size

4.2.2.  borus02 Model

We finally analyze the ESO 565-G019 spectrum using the borus02 spectral model. The model consists of the borus02 table, the two power-law components with a cutoff, and the mekal component to take into account the soft emission:

$$\begin{eqnarray}&&\begin{array}{c}\mathrm{model}\,\mathrm{ESO}\_B=\mathrm{pha}\ast \left(\mathrm{borus}02+z\mathrm{pha}\ast \mathrm{cabs}\right.\\ \quad \ast \left.\mathrm{cutoffpl}1+C2\ast \mathrm{cutoffpl}2+\mathrm{mekal}\right)\end{array}.\end{eqnarray} \tag{ 4 }$$

The best fit model (χ²/d. o. f. = 248/209) is characterized by a photon index ${\rm{\Gamma }}={1.75}_{-0.22}^{+0.04}$ and a l.o.s. absorption that is consistent with a CT scenario ${N}_{{\rm{H}},{\rm{l}}{\rm{.}}{\rm{o}}{\rm{.}}{\rm{s}}.}={3.72}_{-1.15}^{+* }\times {10}^{24}$ cm⁻². The average column density is also CT, but unconstrained in its upper bound and the covering factor is ${C}_{\mathrm{TOR}}={0.47}_{-0.08}^{+0.18}$, We show the results in Table 4 and the combined 0.6–40 keV spectrum in Figure 5.

4.2.3. Summary of the ESO 565-G019 Spectral Analysis Results

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 ${\rm{\Gamma }}={1.91}_{-0.16}^{+0.15}$ and ${N}_{{\rm{H}},z}={0.74}_{-0.07}^{+0.06}\times {10}^{24}$ cm⁻². 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 ∼8% to ∼4% and those on the l.o.s. column density decrease from ∼8% to ∼3%. Also, in accordance with Marchesi et al. (2018), we find a shift in both the photon index and the l.o.s. column density values. The photon index is reduced by ∼5% (the average decrease in Γ value found by Marchesi et al. 2018 is ∼13%) and the N_H,z one decreases by ∼19% (the average value in Marchesi et al. 2018 is ∼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 Γ and N_H,z . In Figure 6, we show the comparison between the 0.6–10 keV and 0.6–70 keV Γ–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.

Zoom In Zoom Out Reset image size

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 could 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 plot (dashed lines); however, they are consistent, and no variability effects can be seen.

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 Fabbiano et al. (2017), Jones et al. (2020), and Ma et al. (2020), 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 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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., Fabbiano et al. 2017; Jones et al. 2020; Ma et al. 2020, for a detailed description of this technique). Figure 9 shows the radial profiles of the emission versus 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 eight annuli from ∼07 to ∼8''. The 0.3–3.0 keV profiles of NGC 3081 and ESO 565-G019 show a significant emission above the PSF values up to ∼7'', 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 2, ESO 565-G019 has a SFR on the order of ∼3–4 M_⊙ yr⁻¹ (Gandhi et al. 2013), 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 (e.g., Ranalli et al. 2003). We find that ESO 565-G019 has a $\mathrm{SFR}={4.4}_{-0.7}^{+0.8}$ M_⊙ yr⁻¹ (the errors correspond to the dispersion of the relation we adopted), consistent with literature. For NGC 3081, there is little (nuclear SFR <0.05 M_⊙ yr⁻¹, Esquej et al. 2014) to no evidences for nuclear star-formation activity (see, e.g., Esparza-Arredondo et al. 2018; Fuller et al. 2019), 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 $\mathrm{SFR}={1.2}_{-0.2}^{+0.2}$ M_⊙ yr⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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% ± 1.8% for the 0.3–3.0 keV extended emission of NGC 3081 and 18.9% ± 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 (Fabbiano et al. 2017; Jones et al. 2020; Ma et al. 2020). 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 (see e.g., Ma et al. 2020).

Based on the results of this work along with those reported by Torres-Albà et al. (submitted), 8 out of 10 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 (∼91%) of sources show this discrepancy (see Figure 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²⁴ cm⁻². With the addition of this source, the percentage of NuSTAR-confirmed CT-AGN in the BAT sample at z ≤ 0.05 is ∼8% (32/417), ⁸ still much lower than the predictions. However, with the analysis of the full ∼60 source sample, Torres-Albà et al. (submitted) observed a significant decrease of the CT fraction with redshift. In fact, at z < 0.01, the fraction is 20.0% ± 5.7%, which is much closer to the theoretical predictions.

Despite the many studies carried out over the last 20 yr, 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 (e.g., Alexander et al. 2003; Gilli et al. 2007; Treister et al. 2009). 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 (see, e.g., Chen et al. 2013; Azadi et al. 2015).

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 CT 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 shown. 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 a 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 CT column density both on the l.o.s. and in average, and this is represented by over-dense clumps in the whole structure.

Zoom In Zoom Out Reset image size

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 that characterize obscured AGN at low redshift (z ≤ 0.1). Its contribution mainly consists in decreasing the errors of the parameters of interest and in breaking the 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 that 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⁶ and 10⁷ 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 (e.g., Brinkman et al. 2002; Piconcelli et al. 2011).
• 3.  
  We found that NGC 3081 is best fitted by the MYTorus model in the decoupled 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, CT. ESO 565-G019 is classified as CT 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}_{\mathrm{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}_{\mathrm{TOR}}={0.47}_{-0.08}^{+0.18}$, along with the ΔN_H, suggests a scenario in which the obscuring structure, which is CT, 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 Marchesi et al. (2019), C_TOR ≃ 0.6, for a ∼30 CT-AGN candidates sample.
• 5.  
  The main nuclear emission can be divided into the reprocessed component, which is heavily suppressed in CT AGN, and the component that 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 ∼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.