X-ray analysis of SDSS J165202.60+172852.4, an obscured quasar with outflows at peak galaxy formation epoch

SDSSJ1652 was observed with XMM-Newton (Jansen et al., 2001) on 2019 March 5-6 for an exposure of 130 ks (ObsID: 0830520101). Exposures were taken with the European Photon Imaging Camera (EPIC), which carries a set of three X-ray CCD cameras: MOS-1, MOS-2, and PN (Strüder et al., 2001; Turner et al., 2001). The EPIC CCDs are sensitive from 0.15-15 keV with moderate spectral energy resolution ($E/\Delta E\sim$ 20-50) and point-spread function with 6” full width at half maximum (FWHM).

Each EPIC images was reduced with the XMM-Newton Science Analysis Software (XMMSAS; Gabriel et al. 2004) v18.0.0 and the Current Calibration Files (CCF) v3.12. The calibrated events were produced using the EPPROC task. A rate-filter was applied to remove minor background flares, such that the effective exposure time was cut to approximately 50 ks.

We extracted the source spectra from 15” radius apertures. Background spectra were extracted from the same chip as the source; however the aperture was treated differently depending on the source’s location on the chip. For the EPIC-PN, we took a circular 110” radius aperture near the source aperture. For the EPIC-MOS, we took a 60” to 200” annular aperture centered on the source. Since the background regions also overlapped with neighboring point-source detections, we excluded 15” radii apertures for each contaminant. Images of the source with the 15” extraction apertures for the PN detector at three energy bands are shown in Figure 1. The redistribution matrix file (RMF) and auxiliary response file (AMF) were generated with the RMFGEN and ARFGEN respectively. The final background-subtracted spectra were grouped at 1 photon per bin using grppha. Photometry results are summarized in Table 2.

The NuSTAR (Harrison et al., 2013) observatory is the first focusing satellite with sensitivity over a broad high-energy 3-79 keV energy band. It consists of two focal-plane modules (FPMA and FPMB) with a 12’$\times$12’ field of view. Since we expected low-photon detections due the predicted Compton-thick nature of the SDSSJ1652, a much longer 218 ks observation was performed with NuSTAR on 2019 February 23 (obsID: 60401028002). The NuSTAR data were reduced with the NuSTAR Data Analysis Software (NuSTARDAS) v20160502 and the NuSTAR Calibration Database (CALDB), using the nupipeline tasks. The spectra and light-curves were extracted with the nuproducts task.

Photometry results in Table 2 were extracted from a 30” radius circular region centered on the target coordinates from both detectors. Despite the long exposure, the FPMB detector did not detect any significant signal from the source. Thus we only use data from the FPMA detector for spectral analysis. We calculate the 3$\sigma$ upper limit for FPMB in three bands: 6-10 keV, 10-40 keV, and 40-79 keV. The NuSTAR photometry results are shown in Table 2. The non-detection can be attributed to a few reasons: high attenuation of X-ray signals due to the obscuring material, limited exposure time, and stray light, which complicated background subtraction.

We first fit the SDSSJ1652 spectra with a power-law. The initial fit produces a small photon-index $\Gamma\sim 0.8$. Since the typical value of the power-law slopes measured in active galactic nuclei is $\Gamma\approx 2$ (Nandra & Pounds, 1994; Reeves & Turner, 2000; Page et al., 2005; Piconcelli et al., 2005), the flat spectral slope of $\Gamma\sim 0.8$ is unlikely to reflect the intrinsic properties of the source, but rather indicates the presence of significant absorption (Alexander et al., 2001). The resulting observed spectrum is harder than the intrinsic spectrum of an unobscured quasar (Nandra & Pounds, 1994; Goulding et al., 2018).

Next, we model the spectra using a power-law with photoelectric absorption in the rest-frame. We refer to this as the absorbed power-law model, and the corresponding XSPEC model is listed in Table 3. When all model parameters (column density $N_{H}$, normalization $K_{norm}$, and the photon index $\Gamma$) are set as free parameters, the best-fit model also produce a low power-law index $\Gamma\sim 0.8$, which is likely due to the high soft-energy excess, hinting complex attenuation mechanism. Instead, we adopt the approach taken by Goulding et al. (2018), among others, and force a more physically motivated model. We fix the photon index at a typical slope of $\Gamma=1.9$ and freely fit for the column density $N_{H}$ and energy normalization $K_{norm}$. The observed spectra with the best-fit model are shown in Figure 2. The best-fit parameters are $N_{H}=(9^{+4}_{-3})\times 10^{23}\textrm{ cm}^{-2}$ and $K_{norm}=(6\pm 2)\times 10^{-5}\textrm{ phot keV/cm}^{2}\textrm{/s}$ at 1 keV with an improved $C_{stat}=633$ for 693 degree of freedom (dof). The confidence ranges of the calculated $N_{H}$ and $K_{norm}$ are produced with the XSPEC steppar routine. The confidence contours of $N_{H}$ and $K_{norm}$ fits is shown in Figure 2. The power-law model with $\Gamma=1.9$ yields the same results for the intrinsic luminosity whether it is applied to just the NuSTAR data or to the combined XMM-Newton$+$NuSTAR dataset. Spectral fit results are summarized in Table 4.

The best-fitting model spectrum in Figure 2 displays a noticeable inversion around the observed-frame 1.5 keV. Examining the data-model residual, we see that this absorbed power-law model does not capture this observed-frame soft energy excess in the 0.5-1 keV range. Although the Fe K $\alpha$ emission line (1.6 keV in the observed-frame) is present (Piconcelli et al., 2015), due to the low counts observed, we cannot find an appropriate binning that comfortably shows both the line and the continuum. These observed spectral features combined with the initial column density estimate approaching $N_{H}\sim 10^{24}\textrm{ cm}^{-2}$ suggest significant Compton-thick obscuration with a more complex absorption structure.

We adopt the physically motivated MYTorus spectral model (Murphy & Yaqoob, 2009; Yaqoob, 2012), which is a Monte Carlo model based on a toroidal circumnuclear reprocessor that also accounts for geometry. The increased soft X-ray component seen in Figure 2 suggests a Compton scattering continuum, which MYTorus can also model self-consistently. We fit a model that uses pre-calculated MYTorus tables that include: (a) the zeroth-order reprocessed power-law continuum - MYTZ, (b) the Compton-scattered continuum - MYTS, (c) the scattered fluorescent line emission due to Fe K at $\sim 6.4\textrm{ keV}$ rest-frame - MYTL, and (d) the optically-thin scattered power-law continuum, separate from the MYTorus tables, to account for electron scattering in a warm/hot ionized region surrounding the central engine. This added optically-thin scattering region is larger than the obscuring structure. For simplicity, we refer to this model set as the MYTorus model throughout the paper. This spectral shape is motivated by the stacked analysis in Goulding et al. (2018). The XSPEC implementation of this model is shown in Table 3. All models considered assume that the intrinsic power-law continuum terminates at 500 keV, referred to as the terminal energy.

Similar to the absorbed power-law model treatment, we fix $\Gamma=1.9$ and tie it to each MYTorus model components. We also fix the MYTorus normalization constants at unity and link each of their column densities $N_{H}$. Previous UV spectropolarimetry analysis suggested a broad-line region scattering, in which the likely line-of-sight inclination angle $\theta_{i}$ to the central engine is near edge-on; however, $\theta_{i}$ was not constrained at the time (Alexandroff et al., 2018). Rather than treating $\theta_{i}$ as another free parameter, we fix at a particular value, and freely fit for $K_{norm}$, $N_{H}$, and the scattering fraction $f_{sc}$ corresponding to the optically-thin scattered continuum. We repeat this fitting process over four angles: $\theta_{i}=45^{\circ}$, $60^{\circ}$, $75^{\circ}$, and $85^{\circ}$ (edge-on is defined as $90^{\circ}$). In each iteration of the fit, we initially fix $f_{sc}=1\%$ to localize $N_{H}$, then free the $f_{sc}$ parameter for simultaneous $K_{norm}$, $N_{H}$, and $f_{sc}$ fitting.

The best-fit parameters for the four MYTorus models are listed in Table 4, and the corresponding confidence contour maps of the $N_{H}$ and $f_{sc}$ fits are shown in Figure 3. We find a strong dependence on $\theta_{i}$, in which greater inclination angles are preferred. In fact, the uncertainties for the $\theta_{i}=45^{\circ}$ and $60^{\circ}$ could not be constrained, reaching extreme estimates for $N_{H}$ and $f_{sc}$ as indicated by the $68\%$ contour. From the fit values and the contour maps, we conclude that the $\theta_{i}=75^{\circ}$ and $85^{\circ}$ geometries fit the data better. Since the estimated values for $N_{H}$ and $f_{sc}$ at these two angles are nearly identical with overlapping error bounds, it is difficult to conclusively pick one model over the other. Neither were we able to constrain a lower limit for $f_{sc}$ to high confidence. However, the confidence contour maps in Figure 3 suggest that the model parameters are better constrained at $\theta_{i}=85^{\circ}$ with $N_{H}=$ $($$1.02$${}^{+0.76}_{-0.41}$$)\times 10^{24}$$\textrm{ cm}^{-2}$ and $f_{sc}=(3\pm 2)$ % at $90\%$ confidence with $C_{stat}=610$ for 667 degrees of freedom.

The unfolded energy spectrum with the best-fit $\theta_{i}=85^{\circ}$ MYTorus model is shown in Figure 4. We compare the different best-fitting unfolded MYTorus $\theta_{i}$ models with the observed data. We also plot the models and their respective components in Figure 4. Despite the differences in fits, all of the MYTorus model fit results suggest that SDSSJ1652 is a highly obscured, Compton-thick quasar. The calculated model (redshift and absorption-corrected) rest-frame luminosity for the $\theta_{i}=85^{\circ}$ model is $L_{\textrm{2-10}}$$=$ $($$1.4$${}^{+1}_{-1}$$)\times 10^{45}$ $\textrm{erg s}^{-1}$.

The C-stat values of the best fitting the absorbed power-law and MYTorus models both indicate similar goodness-of-fits. Here we focus on the data-model residuals in Figures 2 and 4. In general, the MYTorus model-fits capture the spectral features across all energies better than absorbed power-law model (Figure 2 left vs. Figure 4 top left). We argue that the MYTorus model fit is mostly driven by the spectral shape in $E>2$ keV observed frame. In fact, according to just the C-stat values, the $\theta_{i}=45^{\circ}$ and $60^{\circ}$ MYTorus models would be preferred over the $\theta_{i}=75^{\circ}$ and $85^{\circ}$ models, which is not supported by the parameter confidence contours. An obscuring structure with a more inclined geometry is preferred to explain the medium/hard energies.

However, there are some uncertainties within the MYTorus model fits that do not fully model the spectral features at $<1.5\textrm{ keV}$ and 3-5 keV, which are more pronounced in the XMM-Newton/PN data (blue data points). The background spectrum in Figure 2 do not show any significant features at 3-5 keV, suggesting that the unmodeled features are not likely to be background driven. Unmodeled components in 3-5 keV may suggest that we are not fully modeling the peak of the X-ray emission, implying underestimated $N_{H}$, $L_{\textrm{2-10}}$, and $f_{sc}$.

To understand the soft energy fits, we also considered MYTorus models without the warm electron scattering, which produce similar best-fit $\theta_{i}$, $N_{H}$, $K_{norm}$, and C-stat values. With this modification, the best-fit spectrum displays deviations in the soft-energies similar to that of the absorbed power-law model, suggesting the need for a soft-energy scattering component. Therefore, we argue that the inclined MYTorus geometry with the scattering region is the preferred model. A possible explanation for the systematic model residuals in our best-fit MYTorus model is a presence of a more complex thermal structure, not captured by MYTorus or the electron scattering components. However, adding a thermal component did not change the parameters nor improve the fit quality, so we cannot preferentially constrain the model. Another possibility is a presence of a more complex obscuring torus that can be described with variable opening angle model (borus; Baloković et al. 2018) as explored in LaMassa et al. (2019), but this would require higher signal data. Finally, the remaining uncertainties may reflect the poor data quality from severely photon limited studies.

There is an ongoing debate about the relationship between outflows and accretion in luminous quasars: whether near-Eddington accretion is associated with disrupting the X-ray emitting corona or whether the corona must be compact and shielded to enable radiative driving (Luo et al., 2013). The X-ray luminosity is one of the key determinants of the wind launching mechanism.

X-ray emission and mid-infrared emission are both tracers of black hole accretion activities. X-rays originate from the corona over the inner accretion flow and the mid-infrared emission is re-radiated by the obscuring dust on larger scales. A common analysis is to compare the intrinsic 2-10 keV X-ray luminosity, $L_{\textrm{2-10}}$, against the mid-infrared rest-frame $6\mu\textrm{m}$ luminosity, $L_{6\mu\textrm{m}}$. Stern (2015) and Chen et al. (2017) showed an empirically derived sub-linear $L_{6\mu\textrm{m}}$-$L_{\textrm{2-10}}$ relationship for luminous Type I quasars with $L_{6\mu\textrm{m}}\gtrsim 10^{45}$ $\textrm{erg s}^{-1}$, in which the measured $L_{\textrm{2-10}}$ at fixed $L_{6\mu\textrm{m}}$ is an order of magnitude below the extrapolations expected from the linear Gandhi et al. (2009) relationship, as shown in Figure 5. The origin of this sub-linear relationship remain unresolved.

One possibility is that luminous quasars with strong X-rays would overly ionize the circumnuclear gas and suppress the production of radiatively-driven winds (Luo et al., 2013). The physical connection would then lead to the observed Stern (2015) and Chen et al. (2017) relations through the coronal quenching phenomenon (Luo et al., 2013), in which X-ray emission is suppressed by a massive accretion flow such as a failed disk wind that falls into the accretion disc (Proga, 2005; Leighly et al., 2007; Lusso et al., 2010; Jin et al., 2012). This would suggest that high accretion efficiency leads to a decrease in radiative efficiency.

Interestingly, there is evidence for other luminous quasars with even weaker observed X-ray luminosities, falling below the Stern (2015) and Chen et al. (2017) relations even after accounting for absorption. Examples of quasars with a very weak X-ray emission include a Type I quasar with emission line outflows (Leighly et al., 2007), broad absorption line quasars (BAL; Luo et al. 2013), ultraluminous infrared galaxies (ULIRGs; Teng et al. 2014, 2015), and hot dust-obscured galaxies (HotDOGs; Ricci et al. 2017). Either these quasars are intrinsically X-ray weak or they have a complex inner geometrical structure with multiple absorbers which cannot be identified in a single low-quality spectrum.

Combining the WISE photometry data and the X-ray analysis from this paper, we can perform a similar analysis to probe the outflow mechanisms of SDSSJ1652. In Figure 5, we compare the results for SDSSJ1652 against the Stern (2015) and Chen et al. (2017) relationships and the stacked Chandra measurements of ERQs, which includes SDSSJ1652, by Goulding et al. (2018). Given the uncertainties of individual sources using the stacking technique, the stacking estimate and our estimate of the intrinsic luminosity of SDSSJ1652 are consistent, while deeper observations of SDSSJ1652 places a better constraint on the absorption-corrected luminosity. This result supports the Goulding et al. (2018) observation that there is no strong evidence that ERQs are intrinsically X-ray weak: they appear to be consistent with Type I quasars once the intervening absorption is accounted for. These results would suggest that despite the presence of extreme outflows, the heavily obscured SDSSJ1652 and other ERQs are not dominated by extreme coronal quenching that would otherwise result in extreme X-ray suppression, hence they are not intrinsically X-ray weak.

It is also possible for the measured $L_{6\mu\textrm{m}}$ to be underestimated due to obscuration or due to anisotopic emission (Zappacosta et al., 2018), which would imply $L_{\textrm{2-10}}$ deficit compared to the standard relationships and strong X-ray suppression. This would require a correction of nearly an order of magnitude to impose a significant departure from the observed $L_{6\mu\textrm{m}}$-$L_{\textrm{2-10}}$ relationship, which is unlikely. Goulding et al. (2018) notes that the observed high luminosities and outflows in ERQs may be due to a combination of color selection and orientation effects. Current observations of SDSSJ1652 and other ERQs suggest that they are dissimilar from X-ray weak BAL quasars in Luo et al. (2013) in terms of their relationship between extreme accretion, coronal quenching, and outflows. Therefore, it is possible that the ERQs have an outflow mechanism that is different from those of BAL quasars.

The obscuration geometry of the best-fit MYTorus model is broadly comprised of two components: the obscuring torus and optically-thin scattering region. The $\theta_{i}=85^{\circ}$ MYTorus result is consistent with the near-edge-on line-of-sight geometry described by Alexandroff et al. (2018), which suggest the observed X-ray absorption must occur in same wind as the observed UV scattering (Goulding et al., 2018). In the Alexandroff et al. (2018) model, the polarized UV wind in SDSSJ1652 is likely due to dust scattering somewhere in the broad emission line region (BELR) and broad absorption line region (BALR) at scales greater than the obscuration. The inclined viewing angle inferred in the UV is based on the observed shape of the infrared SED. With the inclusion of the X-ray results, our MYTorus model-fits adds the absorbing torus and the warm/hot ionized electron scattering region surrounding the central engine. A cartoon picture of the combined SDSSJ1652 model is shown in Figure 6.

The optically-thin scattered power-law continuum is commonly observed as a rise in the X-ray spectrum towards the low energies. The scattered spectrum is assumed to have the same shape (i.e. $\Gamma$) as the incident power-law spectrum and is assumed to be due to electron scattering. This means that the scattering fraction $f_{sc}$ is the ratio of the scattered flux to the incident flux. This extended scattering region needs to be on scales larger than the central engine and the obscuring torus, which in turn must be larger than the derived dust sublimation region. Using the dust sublimation estimate (Barvainis, 1987) and the extrapolated UV luminosity, we find that the dust sublimation distance is $\sim 1\textrm{ pc}$ (Alexandroff et al., 2018).

We can describe $f_{sc}$ modeled in MYTorus in terms of the size of the scattering region and the density of the scattering particles (Zakamska et al., 2005). We assume scattering is dominated by Thomson electron scattering with scattering cross-section $\sigma_{T}$ in a region of uniform electron density $n_{e}$ with differential cross section $d\sigma/d\Omega\sim\sigma_{T}/4\pi$ through a cone subtending $\Delta\Omega=2\pi$, defined by MYTorus in Figure 6:

Interestingly, $f_{sc}\sim 3\%$ is similar to the UV scattering estimated from Alexandroff et al. (2018), which was obtained by extrapolating the IR emission of ERQs using a Type I quasar SED to the observed UV flux. The similar geometry suggested by both X-ray modeling and UV spectropolarimetry, which is dominated by dust scattering, also suggests that the X-ray absorption occurs in the same wind as seen in spectropolarimetric observations. This is consistent with the Goulding et al. (2018) model. Thus, we consider the dust scatter scale height limits $h\sim 1-10\textrm{ pc}$ from Alexandroff et al. (2018) in Equation 1 to calculate the electron scattering density ranging between $n_{e}\sim 7.5\times(10^{2}-10^{3})\textrm{ cm}^{-3}$.

We also use the scale-height to approximate the line-of-sight mass of the obscuring material. Gas dominates the absorption at X-ray energies (Hickox & Alexander, 2018), so we consider a hydrogen cloud at fixed metallicity with a constant density $n_{H}$ distributed over a spherical volume with a scale height $h$, based on the assumption that the scattering zone would be of scales similar to those of obscuration. To match a torus-like distribution, we take the volume $V=\frac{1}{3}h^{3}\Delta\Omega$, where $\Delta\Omega=2\pi$ is the solid angle subtended by the torus, based on the MYTorus model in Figure 6. The resulting gas mass is therefore $M_{obsc}=m_{H}n_{H}V$. In terms of the measured column density $N_{H}=n_{H}h$, we have

Given the uncertainty in $h$ and in the torus geometry, we calculate the mass range of the line-of-sight obscuration to $M_{obsc}\sim 1.7\times(10^{4}-10^{6})M_{\astrosun}$. Although it is more appropriate to properly model the the torus structure, we cannot set any limits on its size from MYTorus. At most, we expect a small correction by a factor of a few.

We calculate the order of magnitude approximations for the central black hole accretion properties. The directly observed $L_{6\mu\textrm{m}}$ of $10^{47.19}$ $\textrm{erg s}^{-1}$ suggests that the bolometric luminosity may be on the order $L_{Bol}\sim 10^{47-48}$ $\textrm{erg s}^{-1}$, which is in agreement with the Perrotta et al. (2019) measurement of $L_{Bol}\sim 10^{47.73}$ $\textrm{erg s}^{-1}$. This would correspond to a X-ray bolometric correction $\kappa_{X}\approx 380$ for $L_{Bol}=\kappa_{X}L_{\textrm{2-10},obs}$. This is roughly in agreement with values found for $z=2-4$ hyperluminous Type I quasars (Martocchia et al., 2017).

While there are many scaling relationships for black hole mass $M_{BH}$ measurements (e.g. host stellar mass, $H\beta$ kinematics) there are many caveats that complicate direct comparison to ERQs, as noted by Zakamska et al. (2019). Without reliable X-ray to $M_{BH}$ scaling relations, it is also difficult to independently measure $M_{BH}$ from this study. With this in mind, we consider the $H\beta$ virial mass estimate of the SDSSJ1652 supermassive black hole by Perrotta et al. (2019): $M_{BH}=3\times 10^{9}M_{\astrosun}$, which corresponds to an Eddington luminosity of $L_{Edd}=4\times 10^{47}\textrm{erg s}^{-1}$. Combined with the Perrotta et al. (2019) bolometric luminosity, the Eddington ratio is estimated to be $L_{Bol}/L_{Edd}\sim 1.2$. Zakamska et al. (2019) explained that a possible argument against high Eddington ratios is an overestimated $L_{Bol}$. However, for obscured populations, underestimated $L_{Bol}$ is more likely especially when the observed $L_{6\mu\textrm{m}}$-$L_{\textrm{2-10}}$ relationship is consistent with Stern (2015); Chen et al. (2017). Therefore, SDSSJ1652 is consistent with being a near- or super-Eddington, Compton-thick source, exhibiting high $L_{2-10}$ and powerful outflows.

Assuming the material scattering X-rays has an outflowing component, we can also calculate its maximum kinetic energy power using the $n_{e}$ estimates of the warm scattering region. We assume the same velocities $v\approx 3000$ km s^-1 seen in Alexandroff et al. (2018) launched at scale-height $h=10\textrm{ pc}$:

which makes up a small fraction of the bolometric output. Thus, the phase of gas responsible for producing the scattered X-rays is not an energetically important component of the quasar ouput.

While there are X-ray studies of highly obscured quasars, most have been limited to low-redshift samples. Another class of reddened quasars are the HotDOGs, which have properties similar to those of ERQs in terms of spectral energy distribution, luminosity, and redshift (Eisenhardt et al., 2012; Assef et al., 2015). Despite these similarities, there are some interesting differences stemming from somewhat different selection methods. One such difference is that HotDOGs are well below the Stern (2015); Chen et al. (2017) relationships (Ricci et al., 2017; Vito et al., 2018), while also showing Compton thick obscurations (Piconcelli et al., 2015), in contrast to ERQs explored in Goulding et al. (2018) and in this paper. We compare their $L_{6\mu\textrm{m}}$-$L_{\textrm{2-10}}$ properties in Figure 5.

A recent NuSTAR observation of a $z=2.298$ HotDOG (W1835+4355) showed Compton-thick obscuration, $N_{H}\approx 0.9\times 10^{24}\textrm{ cm}^{-2}$, and absorption-corrected luminosity of $L_{2-10}=1.6\times 10^{45}\textrm{ erg s}^{-1}$ (Zappacosta et al., 2018) similar to those of SDSSJ1652 presented here. Compared to other HotDOGs, this target also exhibited higher X-ray luminosity, nearly consistent with the Stern (2015); Chen et al. (2017) relations, shown in Figure 5. In contrast to SDSSJ1652, their target was modeled with $\theta_{i}=70^{\circ}$ and a larger model-dependent $f_{sc}\approx 5-15\%$. This suggest remarkable similarities between ERQs and HotDOGs with the only differences in the obscuration properties: $\theta_{i}$ and $f_{sc}$. This indicates a variety in HotDOG properties with some that have overlaping $L_{6\mu\textrm{m}}$-$L_{\textrm{2-10}}$ properties to ERQs.