The NuSTAR Extragalactic Survey: Average Broadband X-Ray Spectral Properties of the NuSTAR-detected AGNs

The NuSTAR data have been processed using the standard NuSTAR Data Analysis Software (NuSTARDAS, v1.5.1) and calibration files (CALDB version 20131223), distributed within the NASA's HEASARC software (HEAsoft v.6.17²⁶ ). The source spectra have been extracted from each individual pointing using the task nuproduct, from circular extraction regions of 45'' radius (enclosing ∼60% of the NuSTAR point-spread function (PSF)); however, in crowded regions the extraction radius was reduced (to a minimum of 30'', ∼45% of the PSF) to minimize the contamination from nearby sources. The background spectra were extracted from four large regions ($\approx 150^{\prime\prime}$ radius) lying in each of the four quadrants (CCDs) in each pointing, removing areas at the position of known bright Chandra sources (${f}_{2\mbox{--}8\mathrm{keV}}\,\gt 5\times {10}^{-15}$ erg cm⁻² s⁻¹). The source and background spectra, as well as the ancillary files, were then combined using the task addascaspec.

The Chandra source spectra were extracted in a consistent way for all the fields, using the 250 ks and 4 Ms Chandra data in the E-CDFS (Lehmer et al. 2005; Xue et al. 2016) and CDFS (Xue et al. 2011), the 800 ks data in the EGS (Nandra et al. 2015) field, and the 1.8 Ms Chandra COSMOS (C-COSMOS; Elvis et al. 2009) and COSMOS-Legacy survey data (Civano et al. 2016). We used the ACIS Extract (AE) software package²⁷ (Broos et al. 2010; Broos et al. 2012) to extract the source spectra from individual observations using regions enclosing 90% of the PSF; background spectra and relative response matrices and ancillary files were also extracted and then combined by means of the FTOOLS addrmf and addarf.

We first analyzed all the NuSTAR spectra ($E\approx 3\mbox{--}25\,\mathrm{keV}$, observed frame) for each individual source to obtain a rough indication of the spectral slope. This is an essential step to produce the composite spectra (see Section 3.2). We therefore fitted all the spectra with a simple power-law model to obtain the effective photon index (${{\rm{\Gamma }}}_{\mathrm{eff}};$ Figure 2). Since we do not use any redshift information for this initial analysis, we performed the fit for all the NuSTAR sources in our sample. The FPMA and FPMB spectra were fitted simultaneously, with a renormalization factor free to vary, to account for the cross-calibration between the two detectors (Madsen et al. 2015). Due to the poor counting statistics characterizing most of our data (see Figure 2), the spectra have been lightly binned with a minimum of one count per energy bin, and the Cash statistic (C-stat; Cash 1979) was used. The median of the resulting effective photon index is ${{\rm{\Gamma }}}_{\mathrm{eff}}={1.57}_{-1.26}^{+0.79}$ (the errors correspond to the 5th and 95th percentiles). In Figure 2, we also show the mean ${{\rm{\Gamma }}}_{\mathrm{eff}}$ (and 1σ uncertainties) in various net count bins.

Zoom In Zoom Out Reset image size

We then performed spectral fitting including the lower-energy spectra from Chandra. In this case we only include the sources matched to a Chandra counterpart and with redshift identification (173 sources). The Chandra data were fitted between 0.5 and 7 keV, and an absorbed power-law model including Galactic and intrinsic absorption was used. The Galactic column density was fixed to the mean values of ${N}_{{\rm{H}}}^{\mathrm{Gal}}=9.0\times {10}^{19}$ cm⁻², ${N}_{{\rm{H}}}^{\mathrm{Gal}}=1.05\times {10}^{20}$ cm⁻², and ${N}_{{\rm{H}}}^{\mathrm{Gal}}=1.79\times {10}^{20}$ cm⁻² for E-CDFS, EGS, and COSMOS, respectively (Dickey & Lockman 1990), while the intrinsic absorption, the photon index (Γ), and the relative normalization between Chandra, FPMA, and FPMB spectra were left free to vary. From these spectra we calculated the flux at 3–7 keV (observed frame), which is the overlapping energy range between the three telescopes, to test whether there is agreement between the Chandra and NuSTAR data sets (since they are not simultaneous), or whether variability might be an issue. In general, we found good agreement between the Chandra and NuSTAR fluxes within a factor of two (the mean NuSTAR/Chandra flux ratio and 1σ error is 1.1 ± 0.5); however, at the faint end the NuSTAR fluxes are systematically higher than the Chandra ones. This effect has already been shown by Mullaney et al. (2015) and Civano et al. (2015) and is consistent with the Eddington bias, which affects the fluxes close to the NuSTAR sensitivity limit.²⁸ Moreover, this effect could be partly due to the NuSTAR flux being a blend of multiple sources (see Section 2), while only one Chandra spectrum was analyzed together with the NuSTAR data. For the sources with multiple Chandra sources within the matching region, we performed the spectral fit of the NuSTAR data with each of the possible Chandra counterparts and then chose as a unique counterpart the Chandra source with the closest 3–7 keV flux, which is likely to give the highest contribution to the "blended" NuSTAR flux.

Since in several cases the spectral fit could not provide constraints on both ${N}_{{\rm{H}}}$ and Γ simultaneously, we fit the Chandra and NuSTAR data fixing ${\rm{\Gamma }}=1.8$, to obtain some constraints on the intrinsic column density. In some cases, when significant residuals are present at soft energies, we added another power-law component to the model, with the spectral slope fixed to that of the primary component (to limit the number of free parameters in the fits; e.g., Brightman et al. 2013; Lanzuisi et al. 2015; Del Moro et al. 2016), but not affected by intrinsic absorption, to account for any soft excess. From these spectra we also calculated the intrinsic X-ray luminosity at 10–40 keV (${L}_{10\mbox{--}40\mathrm{keV}}$, rest frame) as we aim to construct composite spectra in different bins of ${N}_{{\rm{H}}}$ and ${L}_{10\mbox{--}40\mathrm{keV}}$ (Sections 4.3 and 4.4). In Figure 3 we show the ${N}_{{\rm{H}}}$ versus ${L}_{10\mbox{--}40\mathrm{keV}}$ distribution for all the analyzed sources. For the composite spectra we will divide the sources into three column density bins (unabsorbed, ${N}_{{\rm{H}}}\lt {10}^{22}$ cm⁻², hereafter "NH1"; moderately absorbed, ${N}_{{\rm{H}}}={10}^{22}\mbox{--}{10}^{23}$ cm⁻², herafter "NH2"; and heavily absorbed, ${N}_{{\rm{H}}}\gt {10}^{23}$ cm⁻², hereafter "NH3"; see Section 4.3) and into three luminosity bins (${L}_{10\mbox{--}40\mathrm{keV}}\lt {10}^{44}$ erg s⁻¹, hereafter "L1"; ${L}_{10\mbox{--}40\mathrm{keV}}={10}^{44}\mbox{--}{10}^{45}$ erg s⁻¹, hereafter "L2"; and ${L}_{10\mbox{--}40\mathrm{keV}}\geqslant {10}^{45}$ erg s⁻¹, hereafter "L3"; see Section 4.4).

Zoom In Zoom Out Reset image size

Given the faintness of the sources and the limited counting statistics, in many cases the spectral analysis of the individual sources does not provide constraints on the spectral parameters (∼40% have ${N}_{{\rm{H}}}$ upper limits; see Figure 3). We therefore produce composite spectra in order to investigate the average properties of the AGNs detected by NuSTAR. We use both the Chandra and NuSTAR data together to produce the broadband composite spectra ($\approx 2\mbox{--}40\,\mathrm{keV}$, rest frame), as the Chandra data help improve the counting statistics at low energies $E\lesssim 7\mbox{--}10\,\mathrm{keV}$ and allow us to obtain better constraints on the spectral properties than using the NuSTAR data alone.

3.2.1. Averaging Method

The composite spectra were produced adopting the averaging method described in Corral et al. (2008). Briefly, using the best-fitting parameters obtained from the spectral fits described in the previous Section 3.1, i.e., an absorbed power law (with both ${N}_{{\rm{H}}}$ and Γ free) for the Chandra+NuSTAR spectra, we applied these models to the unbinned, background-subtracted spectra and saved the unfolded spectra in XSPEC (v.12.9.0) in physical units (keV cm⁻² s⁻¹ keV⁻¹) in the energy range 3–25 keV for the NuSTAR data and between $E=\max (0.5,1.0/(1+z))$ and 7 keV for the Chandra data. We limit the Chandra spectra to energies above 1 keV (rest frame) to minimize the contribution from any soft component, which can create distortions in our composite spectra. Each spectrum was then shifted to the rest frame. To combine the Chandra and NuSTAR data, we renormalize the Chandra spectra to the NuSTAR flux at 3–7 keV first, in order to correct for any flux differences, and apply a correction for the Galactic absorption. We then created a common energy grid for all the spectra, with at least 1200 summed counts per energy bin,²⁹ renormalized all the spectra to the flux in the rest-frame 8–15 keV energy range, and redistributed the fluxes in each new energy bin using Equations (1) and (2) from Corral et al. (2008). The renormalization of the spectra is necessary to avoid the brightest sources dominating the resulting composite spectra. We note that the resulting rescaled spectra preserve their spectral slopes and features. Instead of using the arithmetic mean to calculate the average flux in each new energy bin, as done by Corral et al. (2008), we took the median flux in each bin, as the median is less sensitive to outliers of the flux distribution (e.g., Falocco et al. 2012). In the Appendix we analyze in detail the differences between various averaging methods. To estimate the real scatter of the continuum, we performed a resampling analysis and produced 1000 composite spectra drawing random subsamples from the data (e.g., Politis & Romano 1994), excluding some of the spectra (at least one) each time (see Figure 4). For the final composite spectra we took the mean and standard deviation (1σ) of the distribution of median fluxes obtained from the 1000 composite spectral realizations. In Figure 4 we also show the range of composite spectra (1st and 99th percentiles of the distribution) derived from the resampling analysis (shaded areas).

Zoom In Zoom Out Reset image size

We note that the unfolding and the averaging process can distort the shape of the spectrum (e.g., Corral et al. 2008; Falocco et al. 2012). We therefore performed extensive simulations, which are described in the Appendix, to explore the effects of these distortions on the intrinsic average continuum and therefore to derive reliable results from our analyses.

Figure 4 shows the composite spectra from 2 keV to ∼40 keV (rest frame) of the NuSTAR+Chandra data. We produced a composite spectrum for all of the 173 sources with redshift identification (spectroscopic or photometric), and also for BL and NL AGNs separately. Since the number of sources involved in the BL and NL AGN composite spectra is low compared to the whole sample (see Table 1), the scatter in the spectra is larger. The Chandra data significantly improve the counting statistics at $E\lesssim 10\,\mathrm{keV}$ over the NuSTAR-only spectra and allow us to extend the composites to lower rest-frame energies in order to place constraints on the absorbing column densities, as well as on the intrinsic spectral slope. We initially fit these spectra with our baseline model, with both ${N}_{{\rm{H}}}$ and Γ free to vary. The model also includes a Gaussian emission line, as all three composite spectra clearly show an iron Kα emission line at $E\approx 6.4\,\mathrm{keV}$. We fixed the line width to $\sigma =0.1\,\mathrm{keV}$ (which is consistent with the values found for an unresolved Gaussian line in stacked spectra; see, e.g., Iwasawa et al. 2012; Falocco et al. 2013), while the central energy of the line was left free to vary, to account for possible scatter due to the use of photometric redshifts for some of our sources.³⁰ We caution that the hydrogen column densities derived from these spectra do not represent true median values, due to the nonlinear nature of the photoelectric absorption. The results of our spectral fits are reported in Table 2. We find a slightly flatter photon index for the composite of all the sources and for the NL AGNs (${\rm{\Gamma }}={1.65}_{-0.03}^{+0.03}$ and ${\rm{\Gamma }}={1.61}_{-0.07}^{+0.07}$, respectively) compared to the typical ${\rm{\Gamma }}\approx 1.8\mbox{--}2.0$. On the other hand, the BL AGN spectral slope is in good agreement with the typical values (${\rm{\Gamma }}={1.78}_{-0.02}^{+0.03}$).

Table 2.  Spectral Parameters of the Composite Spectra in Different ${N}_{{\rm{H}}}$ and Luminosity Bins

		WA×PO+GAUSS				TORUS			WA×PO+GAUSS+PEXRAV
Bin	No.	Γ	${N}_{{\rm{H}}}$ a	EW (eV)	${\chi }^{2}/$dof	Γ	${N}_{{\rm{H}}}$ a	${\chi }^{2}/$dof	Γ	${N}_{{\rm{H}}}$ a	R (${R}_{{\rm{\Gamma }}=1.8}$)	${\chi }^{2}/$dof
All	173	${1.65}_{-0.03}^{+0.03}$	${1.8}_{-0.5}^{+0.5}$	${92}_{-22}^{+22}$	86.1/90	${1.66}_{-0.05}^{+0.04}$	${1.4}_{-0.4}^{+0.5}$	88.5/90	${1.65}_{-0.02}^{+0.05}$	${1.8}_{-0.5}^{+0.5}$	$\lt 0.17$ (${0.47}_{-0.14}^{+0.15}$)	86.1/89
BL	76	${1.78}_{-0.02}^{+0.03}$	$\lt 0.3$	${75}_{-21}^{+21}$	94.3/90	${1.79}_{-0.01}^{+0.01}$	$\lt 2.7$	93.7/90	${1.84}_{-0.06}^{+0.08}$	$\lt 0.7$	$\lt 0.64$ (${0.12}_{-0.11}^{+0.10}$)	91.8/89
NL	69	${1.61}_{-0.05}^{+0.05}$	${5.8}_{-1.0}^{+1.0}$	${105}_{-41}^{+41}$	89.5/90	${1.60}_{-0.05}^{+0.08}$	${4.8}_{-0.8}^{+0.8}$	93.0/90	${1.62}_{-0.07}^{+0.05}$	${5.8}_{-1.0}^{+1.0}$	$\lt 0.15$ (${0.46}_{-0.30}^{+0.32}$)	89.5/89
HB	64	${1.62}_{-0.05}^{+0.05}$	${1.6}_{-0.7}^{+0.7}$	${76}_{-25}^{+25}$	110.0/90	${1.61}_{-0.03}^{+0.06}$	${1.2}_{-0.5}^{+0.6}$	113.1/90	${1.62}_{-0.03}^{+0.10}$	${1.6}_{-0.7}^{+0.8}$	$\lt 0.39$ (${0.67}_{-0.21}^{+0.22}$)	110.0/89
SB	79	${1.69}_{-0.05}^{+0.05}$	${1.7}_{-0.7}^{+0.7}$	${97}_{-30}^{+30}$	75.2/90	${1.74}_{-0.15}^{+0.13}$	${1.8}_{-1.4}^{+1.4}$	33.7/90	${1.69}_{-0.03}^{+0.07}$	${1.7}_{-0.7}^{+0.7}$	$\lt 0.25$ (${0.30}_{-0.21}^{+0.22}$)	75.2/89
NH1	86b	${1.78}_{-0.02}^{+0.02}$	$\lt 0.3$	${62}_{-21}^{+21}$	84.0/90	${1.79}_{-0.01}^{+0.02}$	$\lt 0.2$	84.1/90	${1.80}_{-0.03}^{+0.08}$	$\lt 0.5$	<0.43 (<0.23)	83.5/89
NH2	44b	1.68${}_{-0.08}^{+0.08}$	${2.4}_{-0.9}^{+0.9}$	${101}_{-43}^{+39}$	59.2/90	${1.69}_{-0.10}^{+0.07}$	${2.1}_{-0.8}^{+0.8}$	58.7/90	${1.72}_{-0.09}^{+0.18}$	${2.5}_{-1.1}^{+1.2}$	<1.08 (${0.52}_{-0.36}^{+0.39}$)	59.0/89
NH3	39b	${1.64}_{-0.11}^{+0.12}$	${28.4}_{-3.3}^{+3.6}$	${100}_{-66}^{+60}$	114.1/77	${1.38}_{-0.12}^{+0.12}$	${16.2}_{-2.6}^{+2.2}$	184.7/79	${1.64}_{-0.11}^{+0.14}$	${28.3}_{-3.3}^{+3.7}$	<0.34 (<0.71)	114.1/76
L1	61b	${1.49}_{-0.10}^{+0.10}$	${1.8}_{-1.3}^{+1.3}$	${115}_{-53}^{+53}$	81.9/90	${1.52}_{-0.12}^{+0.07}$	${1.7}_{-1.1}^{+1.1}$	81.5/90	${1.51}_{-0.06}^{+0.17}$	${1.8}_{-1.3}^{+1.5}$	<0.67 (${1.19}_{-0.50}^{+0.57}$)	81.9/89
L2	84b	${1.71}_{-0.04}^{+0.04}$	${1.5}_{-0.6}^{+0.6}$	${100}_{-17}^{+28}$	98.8/90	${1.72}_{-0.05}^{+0.04}$	${1.2}_{-0.5}^{+0.6}$	102.6/90	${1.71}_{-0.03}^{+0.10}$	${1.5}_{-0.6}^{+0.8}$	<0.34 (${0.29}_{-0.16}^{+0.17}$)	98.8/89
L3	22b	${1.73}_{-0.03}^{+0.05}$	<0.6	${70}_{-37}^{+38}$	77.4/90	${1.75}_{-0.03}^{+0.04}$	$\lt 0.6$	72.0/90	${1.86}_{-0.13}^{+0.19}$	<1.5	<1.74 (${0.34}_{-0.22}^{+0.22}$)	70.2/89

Notes.

^aThe column density ${N}_{{\rm{H}}}$ is expressed in units of 10²² cm⁻². ^bMedian value of the number of sources in each bin, resulting from 1000 spectral realizations with randomized ${N}_{{\rm{H}}}$ or ${L}_{{\rm{X}}}$ values within their error range.

Download table as:  ASCIITypeset image

From extensive spectral simulations (see the Appendix), performed to assess the distortions and variations of the true spectral shape, which might occur at different stages of the stacking process, and thus to validate our spectral analysis results, we find that by simulating unabsorbed NuSTAR and Chandra spectra with a fixed ${\rm{\Gamma }}=1.8$, the resulting composite spectrum has a photon index of ${\rm{\Gamma }}=1.77\mbox{--}1.83$, in good agreement with the slope of the input simulated spectra. This is true also combining unabsorbed spectra with a range of power-law slopes (${\rm{\Gamma }}=1.6\mbox{--}2.0;$ see the Appendix). On the other hand, combining spectra with different levels of X-ray absorption (and therefore different photoelectric cutoff energies) does affect the intrinsic slope of the composite spectra, which becomes slightly flatter (${\rm{\Gamma }}\approx 1.72\mbox{--}1.76;$ Appendix) than the input ${\rm{\Gamma }}=1.8$ of the simulated spectra. This is because the absorption features, which occur at different energies depending on the ${N}_{{\rm{H}}}$ values, are "smoothed" during the stacking process, producing an artificial flattening of the spectral slope (as the true intrinsic ${N}_{{\rm{H}}}$ cannot be recovered). We note, however, that in every test we performed with our simulations we find that this effect is not large enough to explain the relatively flat Γ values observed in the composite spectra from the real data. We can therefore assess that the flattening of spectral slope of the NL AGN composite is real and is significantly different from the slope seen for the BL AGNs. An absorbed power-law model with ${\rm{\Gamma }}\approx 1.61$, in fact, provides a better fit to the data than a slope of ${\rm{\Gamma }}=1.8$ at the >99.9% confidence level, according to the F-test probability. We then fit the data with the TORUS model (Brightman & Nandra 2011), where we fixed the torus opening angle³¹ to ${\theta }_{\mathrm{tor}}=30^\circ$ for all of the spectra, while we fixed the inclination angle to ${\theta }_{\mathrm{inc}}=60^\circ$ for the BL AGNs and ${\theta }_{\mathrm{inc}}=80^\circ$ (nearly edge-on) for the NL AGNs, as these parameters cannot be constrained in the fit. We note that changing the ${\theta }_{\mathrm{inc}}$ value for the BL AGNs to, e.g., 30° makes little difference to the model and to the resulting spectral parameters (the differences are within the parameter uncertainties). We obtain consistent results with those of the baseline model, with a flattened photon index for the composite spectra of all the sources and of the NL AGNs (${\rm{\Gamma }}={1.66}_{-0.05}^{+0.04}$ and ${\rm{\Gamma }}={1.60}_{-0.05}^{+0.08}$, respectively), compared to that of the BL AGNs (${\rm{\Gamma }}={1.79}_{-0.01}^{+0.01};$ see Figure 5). We note, however, that significant residuals in the spectra suggest that the addition of a Gaussian emission line is necessary for all three composite spectra, as the best-fitting TORUS model does not fully account for the iron line emission seen in our spectra. This suggests that possibly also a Compton reflection component is not fully represented by this model, which could explain the resulting slightly flat indices.

Zoom In Zoom Out Reset image size

Indeed, there are two main effects that can produce a flattening in the spectral slope: (1) photoelectric absorption at soft X-ray energies, which can be underestimated and therefore compensated in the fit by a flatter Γ (e.g., Mateos et al. 2005); and (2) Compton reflection contributing to the flux at $E\gt 10\,\mathrm{keV}$ (e.g., Bianchi et al. 2009; Ballantyne et al. 2011). The reflection component may arise from either the accretion disk or cold, dense gas at large distances from the nucleus, such as from the inner part of the putative obscuring torus (e.g., Ross & Fabian 2005; Murphy et al. 2009). To verify how much the reflection can be contributing to our composite spectra, we then used the pexrav model, as described above. In this model we fixed the photon index of the pexrav component to be the same as that of the primary power law ${\rm{\Gamma }}={{\rm{\Gamma }}}_{\mathrm{ref}}=1.8$ (assuming that the NL AGNs have intrinsically the same spectral slope as the BL AGNs), the inclination angle to a mean value of $\cos \,\theta =0.45$, and the power-law cutoff energy to ${E}_{c}=200\,\mathrm{keV}$ (e.g., Ballantyne 2014). We constrain the reflection parameter³² to be $R={0.46}_{-0.30}^{+0.32}$ for the NL AGNs, while for the BL AGNs we obtained $R={0.12}_{-0.11}^{+0.10}$ (see Figure 6). From the composite of all the sources we obtained $R={0.47}_{-0.14}^{+0.15}$ (${\rm{\Gamma }}=1.8$ fixed). Although the scatter on the constraints for the NL AGNs is relatively large, this shows that a larger contribution from a reflection component in the NL AGNs compared to the BL AGNs could indeed be responsible for the flattening of the spectral slope.

Zoom In Zoom Out Reset image size

From the composite spectra we also derived the EW of the iron Kα emission line at $E\approx 6.4\,\mathrm{keV}$. As stated above, the central energy was left free to vary in the fits to allow for the uncertainties that might arise from using photometric redshifts, while the line width was fixed to $\sigma =0.1\,\mathrm{keV}$. The resulting centroid of the emission line is always at $E\approx 6.4\,\mathrm{keV}$ with typical scatter of ${\rm{\Delta }}E\lesssim 0.05\,\mathrm{keV}$. The obtained EWs (reported in Table 2) from our composite spectra are broadly in agreement with the results from previous works that have investigated the composite X-ray spectra of AGNs (e.g., Corral et al. 2008; Falocco et al. 2012; Liu et al. 2016). For the NL AGNs we find slightly higher EW values than for the BL AGNs, which is consistent with the trend seen for the strength of the reflection component; however, given the uncertainties, the trend for the iron line EW is only tentative. Moreover, although some broadening and complexities of the line are visible in the spectra, we only fit a single narrow emission line, since it has been shown that the averaging process and the rest-frame shifting of all the spectra can cause broadening and distortions of the emission line and of the underlying continuum (e.g., Yaqoob 2007; Corral et al. 2008; Falocco et al. 2012). Spectral simulations of the emission line would be required to assess the true nature of these complexities. However, since the detailed properties of the emission line are not the main focus of our paper, we do not attempt more complex analyses.

With our analyses we also aim to test whether there are significant differences between the intrinsic spectral properties of the sources detected in the NuSTAR HB (8–24 keV) and those that formally are not (i.e., they are undetected according to the threshold used by Mullaney et al. 2015; Civano et al. 2015). We then produced a composite spectrum for all the 64 HB-detected sources (see Table 1) and for the soft-band (SB; 3–8 keV) detected sources that are HB undetected (79 sources) and fit the spectra with the models described in the previous section. The spectra are shown in Figure 7. Using our baseline model, the resulting spectral parameters for the HB composite spectrum are ${\rm{\Gamma }}={1.62}_{-0.05}^{+0.05}$ and ${N}_{{\rm{H}}}=(1.6\pm 0.7)\,\times {10}^{22}$ cm⁻², with a fairly weak emission line ($\mathrm{EW}=76\,\pm 25$ eV). For the SB composite spectrum the best-fitting parameters are ${\rm{\Gamma }}={1.69}_{-0.05}^{+0.05}$ and ${N}_{{\rm{H}}}=(1.7\pm 0.7)\,\times {10}^{22}$ cm⁻², with an EW of the iron Kα emission line of $\mathrm{EW}=97\pm 30\,\mathrm{eV}$. Adopting the TORUS model yields consistent results for the spectral slope and intrinsic ${N}_{{\rm{H}}}$ with our baseline model (see Figure 5 and Table 2).

Zoom In Zoom Out Reset image size

Using the pexrav model, as described in the previous section, we constrain ${\rm{\Gamma }}={1.62}_{-0.03}^{+0.10}$ and $R\lt 0.39$ ($R\,={0.67}_{-0.21}^{+0.22}$ for ${\rm{\Gamma }}=1.8$ fixed) for the HB-detected sources; the spectral parameters are consistent within the uncertainties with those obtained for the SB-detected sources, i.e., ${\rm{\Gamma }}={1.69}_{-0.03}^{+0.07}$ and $R\lt 0.25$ ($R={0.30}_{-0.21}^{+0.22}$ for ${\rm{\Gamma }}=1.8$ fixed; Figure 6). If we attempt to account for the "artificial" flattening due to the stacking process, e.g., by fixing ${\rm{\Gamma }}=1.76$ (instead of ${\rm{\Gamma }}=1.8;$ see the Appendix), the constraints on the reflection parameters become $R={0.50}_{-0.19}^{+0.20}$ for the HB composite spectrum and $R\lt 0.37$ for the SB composite spectrum. Although there is a hint of an increase of the reflection strength in the composite spectrum of the HB-detected sources, the spectral parameters of the HB- and SB-detected sources are consistent within the uncertainties. These results therefore suggest that the two subsamples have similar characteristics, with no significant biases toward more obscured or reflection-dominated objects in the HB-detected samples compared to the SB-detected samples (the fraction of X-ray-obscured sources in the two subsamples is ∼52% and ∼46%, respectively).

To place better constraints on the contribution from Compton reflection to the average spectra of the NuSTAR sources and disentangle its effect from that of absorption in flattening the spectral slope (see Section 4.1), we produce composite spectra in different ${N}_{{\rm{H}}}$ bins. In this way, by knowing a priori the median ${N}_{{\rm{H}}}$ of the spectra, we can attribute any hardening of the spectral slope just to the strength of the Compton hump. As an estimate of the intrinsic ${N}_{{\rm{H}}}$ of the sources we used the results from the fitting of the Chandra and NuSTAR spectra for individual sources with ${\rm{\Gamma }}=1.8$ fixed (see Figure 3 and Section 3.1). We caution that this is a very simple model and provides a crude estimate of the column densities, as more complex models might be needed to fully characterize the individual source spectra (e.g., Del Moro et al. 2014; L. Zappacosta et al. 2017, in preparation). However, such analysis is not feasible for many of the sources given the limited counting statistics.

We divided the sources into three ${N}_{{\rm{H}}}$ bins (see Section 3.1): ${N}_{{\rm{H}}}\lt {10}^{22}$ cm⁻² (NH1), ${N}_{{\rm{H}}}={10}^{22}\mbox{--}{10}^{23}$ cm⁻² (NH2), and ${N}_{{\rm{H}}}\gt {10}^{23}$ cm⁻² (NH3), including upper limits. To distribute the sources in each bin, we approximated the ${N}_{{\rm{H}}}$ and errors of each source to a Gaussian distribution and performed 1000 spectral realizations, randomly picking an ${N}_{{\rm{H}}}$ value from the distribution and assigning the source to one of the three ${N}_{{\rm{H}}}$ bins accordingly. For the sources with an ${N}_{{\rm{H}}}$ upper limit (∼40% of the sample; see Figure 3) we assumed a constant probability distribution for the column density values ranging from log ${N}_{{\rm{H}}}$ = 20.5 cm⁻² and the ${N}_{{\rm{H}}}$ upper limit of the source. We excluded two (absorbed) sources owing to their particularly strong soft component and/or spectral complexity, which would further increase the scatter of the composite spectra at $E\lesssim 4\,\mathrm{keV}$ and around the iron line. The median number of sources in each bin, resulting from the 1000 spectral realizations, is 86 in NH1, 44 in NH2, and 39 in NH3 (see Table 2). As described in Section 3.2.1, the final average spectra are obtained taking the mean and the 1σ standard deviation of the distribution of median fluxes in each energy bin, resulting from the 1000 spectral realizations described above, and resampling analysis. We note that the NH3 bin has typically a smaller number of sources compared to the other two, and therefore the scatter in the composite spectrum is larger. The composite spectra in the three different bins are shown in Figure 8 (left panel).

Zoom In Zoom Out Reset image size

We again fit the data with our baseline model and with the TORUS model, including a Gaussian line. The spectrum in the NH3 bin was fitted between 3.5 and 30 keV, as the spectrum at $E\lt 3.5\,\mathrm{keV}$ has large scatter (see Figure 8, left panel), likely due to the presence of a soft-scattered component in some of the individual source spectra. The two models yield fairly consistent results, which are summarized in Table 2. From the fit we obtained values of ${N}_{{\rm{H}}}$ consistent with the median values in each bin for all three composite spectra, while we find a significant flattening of the spectral slope for the NH3 sources (with the highest ${N}_{{\rm{H}}}$), especially with the TORUS model, compared to the unabsorbed and moderately absorbed sources (NH1 and NH2), for which the spectral slopes are consistent with the typical values of ${\rm{\Gamma }}=1.8\pm 0.2$ (see Figure 5). We note, however, that the fitting statistic for the NH3 composite spectrum is poor (see Table 2), and therefore the constraints on the spectral parameters for the sources in this ${N}_{{\rm{H}}}$ bin are less reliable than for the sources with lower X-ray absorption.

To constrain the contribution from the Compton reflection component to the composite spectra, we used the pexrav model, as in the previous sections (see Sections 4.1 and 4.2). For the NH1 spectrum the best-fit parameters are ${\rm{\Gamma }}={1.80}_{-0.03}^{+0.08}$ and ${N}_{{\rm{H}}}\lt 0.5\times {10}^{22}$ cm⁻², with a relative reflection fraction upper limit of $R\lt 0.43$ (${\chi }^{2}$/dof = 83.5/89; see Table 2). For the NH2 spectrum the best-fit parameters are ${\rm{\Gamma }}={1.72}_{-0.09}^{+0.18}$ and ${N}_{{\rm{H}}}=({2.5}_{-1.1}^{+1.2})\times {10}^{22}$ cm⁻², with a relative reflection fraction upper limit of $R\lt 1.08$ (${\chi }^{2}$/dof = 59.0/89). For the NH3 spectrum the best-fitting solution still favors a slightly flatter photon index of ${\rm{\Gamma }}={1.64}_{-0.11}^{+0.14}$ (although still consistent with typical values within the errors) and ${N}_{{\rm{H}}}=({28.3}_{-3.3}^{+3.7})\,\times {10}^{22}$ cm⁻², with a low reflection fraction upper limit of $R\lt 0.34$ (${\chi }^{2}$/dof = 114.1/76).

Since the spectral parameters Γ and R are somewhat degenerate when the scatter in the spectra is large (e.g., Del Moro et al. 2014), we then fixed the photon index to ${\rm{\Gamma }}=1.8$, i.e., the intrinsic value found for the unobscured and BL AGNs, to obtain better constraints on the reflection fraction. From these spectral fits we obtained $R\lt 0.23$ for the NH1 spectrum, $R={0.52}_{-0.36}^{+0.39}$ for the NH2 spectrum, and $R\lt 0.71$ for the NH3 spectrum (see Figure 6; $R\lt 0.51$ if we fix ${\rm{\Gamma }}=1.75$ to account for the "artificial" flattening of the composite spectra with high ${N}_{{\rm{H}}}$). The contribution from reflection is typically low in unobscured sources (NH1), while it seems to increase in obscured sources (e.g., Ricci et al. 2011; Vasudevan et al. 2013). The relatively poor fit obtained for the NH3 composite spectrum, however, does not allow us to place tight constraints on R and therefore to securely assess whether there is a clear dependence between the strength of the Compton reflection and ${N}_{{\rm{H}}}$. This is true also for the EW of the iron Kα line (see Table 2 and Figure 9, left panel). Although we would expect an increase of EW with absorption (as the EW is measured against the absorbed continuum), we cannot find a significant dependence, due to the large errors on our EW measurements. We note, however, that the EW of the iron line in the three ${N}_{{\rm{H}}}$ bins has a similar trend to that of the reflection strength (see Figure 6), despite these two components being fitted independently in our models. This is expected, since the iron line and the Compton hump are two features of the same reflection spectrum.

Zoom In Zoom Out Reset image size

To test whether the average spectral parameters of the NuSTAR sources change as a function of luminosity, we constructed composite spectra in three different luminosity bins: L1, i.e., sources with ${L}_{10\mbox{--}40\mathrm{keV}}\lt {10}^{44}$ erg s⁻¹; L2, with ${L}_{10\mbox{--}40\mathrm{keV}}\,={10}^{44}\mbox{--}{10}^{45}$ erg s⁻¹; and L3, with ${L}_{10\mbox{--}40\mathrm{keV}}\geqslant {10}^{45}$ erg s⁻¹. Similarly to the method adopted in the previous section, we have approximated the luminosity ${L}_{10\mbox{--}40\mathrm{keV}}$ and its uncertainties for each source with a Gaussian probability function and performed 1000 realizations randomly picking a ${L}_{10\mbox{--}40\mathrm{keV}}$ value from the distribution and assigning the source to each luminosity bin. For the luminosity upper limits we assumed again a constant probability function ranging from ${{logL}}_{10\mbox{--}40\mathrm{keV}}=42.0$ to the corresponding upper-limit value of the source. From the 1000 realizations the median number of sources in each bin is 61 in L1, 84 in L2, and 22 in L3 (see Table 2). In this case, although the L3 bin has a small number of sources, the composite spectrum has a fairly high signal-to-noise ratio (S/N), as the sources in this bin are brighter and have typically good counting statistics in each individual spectrum. On the other hand, the L1 composite spectrum has relatively high scatter since it comprises the least luminous sources in the sample and therefore the S/N of the individual spectra is typically lower than those in the other luminosity bins. Figure 8 (right panel) shows the composite spectra in the three luminosity bins.

Fitting the spectra with our baseline model, we find that for low-luminosity sources the photon index is flatter (${{\rm{\Gamma }}}_{{\rm{L}}1}={1.49}_{-0.10}^{+0.10}$) than those of higher-luminosity sources (${{\rm{\Gamma }}}_{{\rm{L}}2}={1.71}_{-0.04}^{+0.04}$ and ${{\rm{\Gamma }}}_{{\rm{L}}3}={1.73}_{-0.03}^{+0.05}$). All the best-fitting spectral parameters are reported in Table 2. Similar results are obtained when using the TORUS model to fit the spectra (Figure 5). The flattening of the spectrum of the L1 sources might be due to a higher incidence of absorbed sources in this luminosity bin (∼57%, compared to ∼47% and ∼36% in the L2 and L3 bins, respectively). Indeed, the median value of the column density at ${L}_{10\mbox{--}40\mathrm{keV}}\lt {10}^{44}$ erg s⁻¹ is ${N}_{{\rm{H}}}({\rm{L}}1)\approx 2.5\times {10}^{22}$ cm⁻², while for the other two luminosity bins the values are lower: ${N}_{{\rm{H}}}({\rm{L}}2)\approx 8.1\times {10}^{21}$ cm⁻² and ${N}_{{\rm{H}}}({\rm{L}}3)\approx 6.4\times {10}^{21}$ cm⁻². However, a K-S test on the ${N}_{{\rm{H}}}$ distribution in the three luminosity bins (see Figure 3) suggests that the distributions are not significantly different ($D({\rm{L}}1\mbox{--}{\rm{L}}2)\,=0.178$ and $\mathrm{Prob}({\rm{L}}1\mbox{--}{\rm{L}}2)=0.227$, and $D({\rm{L}}1\mbox{--}{\rm{L}}3)=0.240$ and $\mathrm{Prob}({\rm{L}}1\mbox{--}{\rm{L}}3)=0.278$).

Some studies have found a dependence of the photon index on luminosity, where high-luminosity sources show steeper Γ than the low-luminosity ones (e.g., Dai et al. 2004; Saez et al. 2008). Conversely, other studies have found an anticorrelation between Γ and X-ray luminosity (e.g., Corral et al. 2011; Scott et al. 2011). Some of these trends are probably a consequence of the stronger dependence that has been found for the photon index with Eddington ratio (e.g., Brightman et al. 2013; Ricci et al. 2013b). The flat photon index we found for the L1 spectrum might be partly due to these correlations. However, the spectral slopes of the sources in the L2 and L3 luminosity bins are pretty much the same (see Figure 5 and Table 2), and therefore there is no clear dependence of Γ on luminosity in our sample (e.g., Winter et al. 2009; Brightman et al. 2013).

We therefore investigate whether there is any difference in the amount of Compton reflection contributing to the spectra as a function of luminosity. We fit the pexrav model to our spectra, as in the previous sections (Sections 4.1 and 4.3), and we constrain the spectral parameters to be ${\rm{\Gamma }}={1.51}_{-0.06}^{+0.17}$, ${N}_{{\rm{H}}}=({1.8}_{-1.3}^{+1.5})\times {10}^{22}$ cm⁻², and $R\lt 0.67$ for the L1 sources; ${\rm{\Gamma }}={1.71}_{-0.03}^{+0.10}$, ${N}_{{\rm{H}}}=({1.5}_{-0.6}^{+0.8})\times {10}^{22}$ cm⁻², and $R\lt 0.34$ for the L2 sources; and ${\rm{\Gamma }}={1.86}_{-0.13}^{+0.19}$, ${N}_{{\rm{H}}}\lt 1.5\times {10}^{22}$ cm⁻², and $R\lt 1.74$ for the L3 sources, which is basically unconstrained. Fixing the photon index ${\rm{\Gamma }}=1.8$, we obtain tighter constraints on the reflection fraction in the three luminosity bins, i.e., $R={1.19}_{-0.50}^{+0.57}$, $R={0.29}_{-0.16}^{+0.17}$, and $R={0.34}_{-0.22}^{+0.22}$ for the L1, L2, and L3 sources, respectively (see Figure 6). While for the low-luminosity sources (${L}_{10\mbox{--}40\mathrm{keV}}\lt {10}^{44}$ erg s⁻¹) the reflection strength is relatively high, we find that it decreases significantly at luminosities ${L}_{10\mbox{--}40\mathrm{keV}}\geqslant {10}^{44}$ erg s⁻¹ (e.g., Bianchi et al. 2007; Shu et al. 2010; Liu et al. 2016; however, see also Vasudevan et al. 2013). Although this trend is only seen fixing the photon index, due to the degeneracies in the spectral parameters, this decrease is significant also if we fit the L1 spectrum with a flatter photon index of ${\rm{\Gamma }}=1.75$ to correct for the "artificial" flattening discussed in the previous sections, which yields $R={0.93}_{-0.47}^{+0.52}$. We find hints of a similar, decreasing trend also for the EW of the iron Kα line measured from the composite L1, L2, and L3 spectra (see Figure 9, right panel), although, given the large errors, this trend is not statistically significant. An anticorrelation between the strength of the Fe Kα line and the X-ray luminosity has been observed and investigated in several previous studies (e.g., Iwasawa & Taniguchi 1993; Nandra et al. 1997; Page et al. 2004; Bianchi et al. 2007; Ricci et al. 2013a), and it is often referred to as the "X-ray Baldwin effect" or the "Iwasawa-Taniguchi effect" (IT effect). We compared our results with the anticorrelation found by Bianchi et al. (2007) for a sample of radio-quiet, type 1 AGNs (dashed and dotted lines in Figure 9, right panel). We converted their 2–10 keV luminosity into the 10–40 keV luminosity assuming a power-law model with ${\rm{\Gamma }}=1.8$. The slope of the anticorrelation is consistent with our results; however, the EWs we find in our spectra are typically larger compared to those found by Bianchi et al. (2007). This is likely due to the different types of analyses and different types of sources (as we also include absorbed AGNs) used in the two papers, as the measurements from our stacking analyses might be biased toward higher values of EW compared to the distribution found from individual sources. Moreover, the emission line in our composite spectra is broadened owing to the rest-frame shifting of the spectra (see Section 3.2.1), thus likely increasing the measured EW.

If the IT effect is due to a decrease of the covering factor with increasing AGN luminosity (e.g., Bianchi et al. 2007), this could also explain the drop of the reflection strength, as in high-luminosity sources there is less material close to the nucleus obscuring/reflecting the intrinsic X-ray emission, thus producing a weaker Compton reflection spectrum. We caution, however, that we cannot exclude that the drop of R at high luminosities could be partly due to an evolution of the spectral properties with redshift, as the mean redshifts of the sources in the L1, L2, and L3 luminosity bins are $\langle z\rangle =0.55$, $\langle z\rangle =1.17$, and $\langle z\rangle =1.86$, respectively.

Even dividing the sources into bins of column density or X-ray luminosity, it is not possible to fully understand how the average spectral properties, such as the strength of the iron Kα line and the Compton reflection, depend on these two quantities, as ${N}_{{\rm{H}}}$ and ${L}_{{\rm{X}}}$ are somewhat linked. For instance, the NH1 bin contains a larger fraction of high-luminosity sources (∼74%) compared to the NH2 and NH3 bins (∼67% and ∼59%, respectively). Similarly, the L1 bin has a larger fraction of absorbed sources than the high-luminosity bins (see Section 4.4). To further investigate how the Compton reflection might depend on luminosity and/or ${N}_{{\rm{H}}}$ independently, we divided the sources into two luminosity bins, ${L}_{10\mbox{--}40\mathrm{keV}}\,\lt 2\times {10}^{44}$ erg s⁻¹ and ${L}_{10\mbox{--}40\mathrm{keV}}\geqslant 2\times {10}^{44}$ erg s⁻¹; the separation was chosen to have a comparable number of sources within the low- and high-luminosity bins (see Table 3). Within these bins we produce separate composite spectra for the unabsorbed (${N}_{{\rm{H}}}\lt {10}^{22}$ cm⁻²) and absorbed (${N}_{{\rm{H}}}$ ≥ 10²² cm⁻²) sources. We adopted the same randomization method described in Sections 4.3 and 4.4 to assign sources to each bin. Since our aim here is to investigate the reflection component, we only fit the composite spectra with the pexrav model. The resulting parameters are reported in Table 3. In Figure 10 (left panel) we show the reflection parameter (R) derived from this analysis as a function of the rest-frame 10–40 keV luminosity; for each composite spectrum, we plot the median ${L}_{10\mbox{--}40\mathrm{keV}}$ of the sources. The uncertainties on the derived spectral parameters are quite large, especially for the obscured sources, and therefore it is not possible to derive a statistically significant dependence of R on luminosity or ${N}_{{\rm{H}}}$; however, there is an indication of a decreasing trend with luminosity. Moreover, the unobscured sources seem to have typically a smaller reflection fraction than the obscured ones. These trends seem to be confirmed by the strength of the iron line measured from these spectra, as shown in Figure 10 (right panel), where the iron line EW is plotted against R. The strength of the emission line for obscured and unobscured sources seems to follow a similar behavior to the Compton reflection strength. Indeed, a Spearman rank correlation test, performed on 10,000 simulated data sets accounting for the uncertainties of the data points, indicates that there is a strong correlation, as the resulting correlation coefficient is ${\rho }_{{\rm{s}}}=0.80;$ however, given the small number of data points and the large uncertainties, the null hypothesis probability is $P=0.38$.

Zoom In Zoom Out Reset image size

Table 3.  Spectral Parameters of the Composite Chandra+NuSTAR Spectra for the Low- and High-luminosity Sources, Further Divided into Unobscured (unob) and Obscured (obs), Fitted with Our pexrav Model, wabs×pow+Gauss+pexrav

Bin	No.a	log ${L}_{10-40\,\mathrm{keV}}^{b}$	zc	Γ	${N}_{{\rm{H}}}$ d	EW (eV)	R (${R}_{{\rm{\Gamma }}=1.8}$)	${\chi }^{2}/\mathrm{dof}$
Low-L unob	39	43.95	0.698	${1.80}_{-0.07}^{+0.10}$	$\lt 1.4$	${103}_{-38}^{+45}$	$\lt 0.36$ ($\lt 0.24$)	71.2/89
Low-L obs	48	43.84	0.660	${1.42}_{-0.10}^{+0.12}$	${6.5}_{-1.9}^{+1.2}$	${118}_{-69}^{+70}$	$\lt 0.30$ (${1.35}_{-0.60}^{+0.72}$)	91.8/89
High-L unob	46	44.81	1.342	${1.85}_{-0.05}^{+0.06}$	$\lt 0.5$	${68}_{-25}^{+23}$	${0.39}_{-0.31}^{+0.40}$ (${0.12}_{-0.11}^{+0.12}$)	83.4/89
High-L obs	33	44.63	1.375	${1.73}_{-0.17}^{+0.22}$	${8.4}_{-2.2}^{+2.3}$	${106}_{-50}^{+62}$	$\lt 1.65$ (${0.71}_{-0.38}^{+0.43}$)	42.5/89

Notes.

^aNumber of spectra used for the composite; this is the median value obtained from the 1000 spectral realizations. ^bLogarithm of the median X-ray luminosity (10–40 keV) of the sources in each bin. ^cMedian redshift of the sources in each bin. ^dThe hydrogen column density is expressed in units of 10²² cm⁻².

Download table as:  ASCIITypeset image

We initially simulated NuSTAR and Chandra spectra as a simple power law with ${\rm{\Gamma }}=1.8$ fixed, using the response and ancillary files extracted from our survey data and a flux distribution similar to that of the real sources. For practicality, since the simulations are very time-consuming, we only used the files of the sources detected in the E-CDFS field to simulate the spectra (i.e., 45 sources). When more spectra are needed to increase the counting statistics in the composite spectra (see below), we generated multiple simulated spectra for each real source (using photon randomization). We simulated both background and source spectra to test more closely the results from the real data. We applied to the simulated spectra the same procedure used for the real data (see Section 3.2), i.e., we fitted the individual background-subtracted spectra with a power-law model (with Γ free to vary) and saved them in physical units (unfolded spectrum in XSPEC); shifted the spectra to the rest frame, using the redshifts of our real sources; created a new energy grid for all the spectra; normalized them to the flux at 8–15 keV; and redistributed the fluxes in the new energy bins. Our average spectra are obtained by taking the median flux in each new energy bin, performing a resampling analysis to estimate the 1σ errors. We then fitted the obtained spectrum with a power-law model to verify whether we can recover the input photon index of the simulated spectra. We obtained ${\rm{\Gamma }}=1.77\pm 0.03$, which is slightly lower than the input spectral slope, but in good agreement, within the errors, with the input parameters.

We performed the same test simulating spectra with different photon indices (45 spectra for each Γ value), in the range of ${\rm{\Gamma }}=1.6\mbox{--}2.0$, with a step size of ${\rm{\Delta }}{\rm{\Gamma }}=0.5$ and obtaining the average spectrum for each input spectral slope. Subsequently, we also combined spectra with different slopes, to test whether we can recover the median input Γ from our composite spectra or whether significant distortions are affecting the spectra. In Figure 12 we show the results of these tests by plotting the resulting spectral slopes, obtained by fitting the composite spectra with a power-law model, versus the input Γ used to simulate the individual spectra. When spectra with different slopes are combined, we plot their median Γ on the x-axis. In general, the spectral slopes obtained from the composite spectra are in good agreement with the input Γ of the simulated spectra used to produce the composites. However, for steep slopes (${\rm{\Gamma }}\gt 1.9$), the resulting Γ tends to be flatter than the input values. This is likely because when the spectra are steep, there are fewer counts contributing to the spectra at high energies (e.g., in the NuSTAR band), increasing the scatter in the composite spectra at $E\gtrsim 10\,\mathrm{keV}$, thus impairing the constraints on the intrinsic spectral slope. This results in a slightly flatter value of Γ. We note that this issue is affecting the results only when a relatively small number of spectra are used to produce the composites (for instance, 45 spectra in our tests with single Γ values), while when we increase the number of spectra (hence, the counting statistics), as in the case of the composites obtained from spectra with different Γ values, the issue is no longer present.

Zoom In Zoom Out Reset image size

We can therefore conclude that our averaging method does not introduce significant distortions to the final composite spectra when combining spectra with different slopes, in the case where a simple power-law model is assumed. From our tests on these simulated composite spectra we can reliably recover the input Γ of the individual spectra, and/or their median value, when a distribution of Γ is assumed.

An important effect that can create distortions in the composite spectra and modify the results from our analyses is the X-ray column density. To test how the X-ray absorption modifies the composite spectra, we simulated the Chandra and NuSTAR source spectra (and the relative background) with different values of ${N}_{{\rm{H}}}$ (${N}_{{\rm{H}}}={10}^{21},{10}^{22},{10}^{23},5\times {10}^{23}$ cm⁻²) and a fixed ${\rm{\Gamma }}=1.8$ and produced the composite rest-frame spectra for individual values of ${N}_{{\rm{H}}}$. An example is shown in Figure 13. Analyzing the spectra (at $E\approx 3\mbox{--}30\,\mathrm{keV}$, rest frame), in all cases we obtained parameters in good agreement with the input values; the results are summarized in Table 4. We note that for the ${N}_{{\rm{H}}}={10}^{21}$ cm⁻² composite, where we obtained a much lower ${N}_{{\rm{H}}}$ value than the input simulated spectra, the discrepancy is due to the fact that we cannot constrain such a low value of ${N}_{{\rm{H}}}$ from our composite spectrum, as it spans the energy range $E\approx 3\mbox{--}40\,\mathrm{keV}$, rest frame. For high levels of ${N}_{{\rm{H}}}$ the resulting photon index tends to be slightly flatter than the input ${\rm{\Gamma }}=1.8$, but it is always consistent within the errors. In these cases, the flattening is likely due to a decrease of the number of counts in the composite spectra at low energies, due to the absorption, thus yielding poorer constraints on the intrinsic spectral slope.

Zoom In Zoom Out Reset image size

Table 4.  Spectral Parameters of the Simulated Composite Spectra, Obtained Combining Simulated Chandra+NuSTAR Spectra with ${\rm{\Gamma }}=1.8$ and Different Levels of X-ray Absorption

Input ${N}_{{\rm{H}}}$ a	No. of Sources	Γ	${N}_{{\rm{H}}}$ a
0.1	45	1.83 ± 0.01	<0.05
1.0	45	1.82 ± 0.02	0.9 ± 0.1
10.0	45	1.78 ± 0.02	10.1 ± 0.2
50.0	45	1.75 ± 0.07	51.4 ± 2.7
0.1, 1.0 (a)	90	1.79 ± 0.03	0.6 ± 0.1
1.0, 10.0 (b)	90 (180)	1.72 ± 0.05 (1.78 ± 0.03)	6.3 ± 0.7 (6.3 ± 0.4)
0.1, 1.0, 10.0 (c)	135	1.76 ± 0.04	4.0 ± 0.7

Note.

^aThe hydrogen column density is expressed in units of 10²² cm⁻².

Download table as:  ASCIITypeset image

We then produced composite spectra combining simulated source spectra with different levels of absorption, to see whether, also in these cases, the resulting spectral parameters are in agreement with the input distributions. We note that we do not expect to obtain the true median value of ${N}_{{\rm{H}}}$ from the composite spectra, as the absorption is not a linear parameter. However, the resulting ${N}_{{\rm{H}}}$ should be somewhat consistent with (or within the range of) the input values. We therefore combined spectra with (a) ${N}_{{\rm{H}}}={10}^{21},{10}^{22}$ cm⁻²; (b) ${N}_{{\rm{H}}}={10}^{22},{10}^{23}$ cm⁻²; and (c) ${N}_{{\rm{H}}}={10}^{21},{10}^{22},{10}^{23}$ cm⁻². The spectral parameters obtained from these composite spectra are reported in Table 4. For the unabsorbed and lightly absorbed sources (composite spectrum "a") the resulting parameters are in very good agreement with those of the input simulated spectra, while for spectra "b" and "c," which also include more absorbed sources, the resulting spectral slope tends to be flatter than those of the input sources. The flattening is mainly due, as described above, to a decrease of the number of counts in the composite spectra for the absorbed sources. Indeed, we tested that increasing the number of input spectra to produce the composite "b" (e.g., from 90 to 180; see Table 4), thus, increasing the counting statistics in the composite spectrum, we can recover a photon index in good agreement with the ${\rm{\Gamma }}=1.8$ of the input spectra. However, another effect plays a non-negligible role in flattening the photon index of the composite spectra: combining spectra with different values of column densities, and therefore different photoelectric cutoff energies, can lead to a slight flattening in the spectral slope, as the typical absorption features get "smoothed" in the final composites, yielding a flatter Γ and possibly an underestimation of ${N}_{{\rm{H}}}$ (since it is not possible to recover the true median value of ${N}_{{\rm{H}}}$ from the composite spectra). This effect is also seen in the real data (Sections 4.1 and 4.3 and Figure 5). However, our simulations indicate that this "artificial" flattening of the spectra due to the averaging process is relatively small (${\rm{\Delta }}{\rm{\Gamma }}\approx 0.2\mbox{--}0.8$) and cannot be the sole cause of the flat spectral slopes seen in the composites of the NL AGNs and the heavily absorbed sources (NH3), for which the resulting photon indices are ${\rm{\Gamma }}\approx 1.6$.

There are different methods one can use to determine the average spectrum of a population. The arithmetic mean (or median) tends to preserve the relative fluxes of the emission features, while the geometric mean tends to preserve the global continuum shape, when it can be approximated by a power law (Vanden Berk et al. 2001). Since the aim of our analyses is to study the intrinsic AGN continuum, the geometric mean would be the most appropriate choice to derive our composite spectra. However, since the NuSTAR spectra have typically high background levels, several energy bins in the individual spectra result in a null (or negative) flux value after background subtraction, due to the background fluctuations. These energy bins have to be excluded from the geometric mean, therefore biasing the resulting average flux.

We performed several tests to assess the differences between these averaging methods and evaluate the best approach to use for our data. We initially used the simulated spectra (see details in the previous section) with known input spectral slopes and analyzed the resulting composite spectra. When we produce composite spectra for unabsorbed ${\rm{\Gamma }}=1.8$ simulated spectra, the three averaging methods yield comparable results, recovering the input ${\rm{\Gamma }}=1.8$. When we introduce various levels of absorption in our simulated spectra, however, the resulting composite spectra are affected by some distortions (see details below), and the averaging methods used to produce the composite spectra yield different results. In Figure 14 we show the confidence contours for the spectral parameters Γ and ${N}_{{\rm{H}}}$ derived from the composite spectra. On the left the input spectra used for the composite are all simulated with the same ${\rm{\Gamma }}=1.8$ and ${N}_{{\rm{H}}}={10}^{23}$ cm⁻² (which are the values we want to recover from the resulting composite spectra); on the right the input spectra used for the composite have a photon index ${\rm{\Gamma }}=1.8$ and ${N}_{{\rm{H}}}=0,{10}^{22},{10}^{23}$ cm⁻² in equal numbers. The median and the geometric mean provide comparable results in the first case, while the arithmetic mean yields a steeper photon index than the input value. When spectra with various column densities are stacked together to produce the composite, the parameters derived from the median spectrum are consistent with the input ${N}_{{\rm{H}}}$ values of the individual spectra (although we do not expect to recover the true median ${N}_{{\rm{H}}}$, as described in the previous section), while slightly underestimating the photon index, which is flatter than the input value ${\rm{\Gamma }}=1.8$, but still consistent within the errors. On the other hand, Γ tends to be overestimated in the composite generated adopting the geometric mean, as well as the arithmetic mean. We also note that, since the median is less sensitive to outliers or extreme values, the composite spectrum produced with this method is less noisy than those produced using the arithmetic or geometric means. We therefore favor the median as the method to produce our composite spectra (see also Falocco et al. 2012), keeping in mind that the flattening of the photon index seen in the composite spectra of our real data is partially due to the stacking method (see previous section).

Zoom In Zoom Out Reset image size