Relativistic Components of the Ultra-fast Outflow in the Quasar PDS 456 from Chandra/HETGS, NuSTAR, and XMM-Newton Observations

The aim of this work is to look for signatures of UFOs in PDS 456. For this purpose, we chose to focus on three different data sets, as presented in Table 1. The strongest signatures of UFOs are found at high energy, above 7 keV, as found in particular in PDS 456, according to previous studies. We thus selected observations including NuSTAR data, namely XMM-Newton and NuSTAR data from 2013 to 2014 (observations "XN"), and Chandra and NuSTAR data from 2015 (observation "CN"). Other signatures of the UFOs are expected at lower energy, so Chandra/HETGS is an excellent instrument to look for such features. We also considered Chandra/HETGS data from 2003 (observation "C"), data that together with CN data have never been analyzed in detail. Despite the lack of high-energy data for this C observation, our analysis, which is performed consistently along the different data sets, allowed us to constrain the wind parameters in this observation as well. For observations XN, we focused on EPIC/pn data, mostly to look for high-energy signatures of the UFOs. We did not reanalyze RGS data, which have already been studied in detail by Reeves et al. (2016), because we focus on high-velocity and high-ionization outflows whose signatures are expected to be beyond the energy range of RGS. The scope of the present paper is to look for signatures of UFOs at both high and low energies, so we did not examine other archival data.

Table 1.  Observation Log

Label	Instrument	Obs. ID	Obs. Start	Obs. End	${T}_{\mathrm{tot}}$	${T}_{\mathrm{net}}$	Channel	Count Rate	Flux
CN	Chandra/HETGS	17452	2015 Jul 21, 12:52	2015 Jul 23, 03:48	138	136	HEG	0.098	5.16 ± 0.07
							MEG	0.100	6.00 ± 0.12
	NuSTAR	90101008002	2015 Jul 21, 11:01	2015 Jul 23, 07:46	160	74	FPMA	0.048	4.63 ± 0.14
							FPMB	0.044	4.51 ± 0.16
		90101008004	2015 Jul 24, 11:36	2015 Jul 25, 10:51	83	38	FPMA	0.044	4.75 ± 0.81
							FPMB	0.042	4.47 ± 0.20
C	Chandra/HETGS	4063	2003 May 7, 03:29	2003 May 8, 20:08	145	143	HEG	0.070	3.85 ± 0.06
							MEG	0.073	4.23 ± 0.07
XN1	XMM-Newton	0721010201	2013 Aug 27, 04:41	2013 Aug 28, 11:13	110	85.5	EPIC/pn	3.175	10.87 ± 0.04
	NuSTAR	60002032002	2013 Aug 27, 03:41	2013 Aug 28, 11:41	114	44	FPMA	0.134	9.70 ± 0.22
							FPMB	0.128	9.50 ± 0.23
XN2	XMM-Newton	0721010301	2013 Sep 6, 03:24	2013 Sep 7, 10:36	112	92.1	EPIC/pn	2.098	5.63 ± 1.39
	NuSTAR	60002032004	2013 Sep 6, 02:56	2013 Sep 7, 10:51	114	43	FPMA	0.051	3.56 ± 0.38
							FPMB	0.049	3.61 ± 0.17
XN3	XMM-Newton	0721010401	2013 Sep 15, 18:47	2013 Sep 17, 03:57	119	102.0	EPIC/pn	1.974	6.08 ± 0.02
	NuSTAR	60002032006	2013 Sep 15, 17:56	2013 Sep 17, 04:01	119	44	FPMA	0.072	5.20 ± 0.18
							FPMB	0.068	4.59 ± 0.18
XN4	XMM-Newton	0721010501	2013 Sep 20, 02:47	2013 Sep 21, 09:37	111	92.9	EPIC/pn	1.925	6.35 ± 0.03
	NuSTAR	60002032008	2013 Sep 20, 03:06	2013 Sep 21, 11:11	119	44	FPMA	0.075	5.18 ± 0.16
							FPMB	0.073	5.38 ± 0.17
XN5	XMM-Newton	0721010601	2014 Feb 26, 08:03	2014 Feb 27, 22:51	140	103.9	EPIC/pn	1.383	4.47 ± 1.38
	NuSTAR	60002032010	2014 Feb 26, 08:16	2014 Feb 28, 22:56	224	110	FPMA	0.045	3.01 ± 0.53
							FPMB	0.043	3.00 ± 0.58

Note. Summary of PDS 456 Chandra, XMM-Newton, and NuSTAR observations used in this study. Observation CN considers the contemporaneous observations with Chandra/HETGS and NuSTAR from 2015, while observation C refers to the archival Chandra/HETGS observation from 2003. Observations XN1 to XN5 refer to XMM-Newton and NuSTAR observations performed in 2013–2014 and previously analyzed by Nardini et al. (2015). We refer to the combination of observations XN1 to XN5 as observation XN. Observation start and end times are in UT. ${T}_{\mathrm{tot}}$ (in ks) is the total elapsed time, while ${T}_{\mathrm{net}}$ (in ks) is the net exposure after screening and deadtime correction. The count rates (in ${s}^{-1}$) are the background-subtracted count rates. The fluxes (in 10⁻¹² erg cm⁻² s⁻¹) refer to the 0.8–7 keV (HEG), 0.4–7 keV (MEG), 0.4–10 keV (EPIC/pn) and 3–30 keV (FPMA/FPMB) bands, obtained using the ISIS function data_flux(), which calculates the absorbed X-ray flux solely from the spectral data, and estimates errors from the data and information from the instrumental responses.

Download table as:  ASCIITypeset image

PDS 456 was observed with the HETGS (Markert et al. 1994; Canizares et al. 2005) using the Chandra Advanced CCD Imaging Spectrometer (ACIS; Garmire et al. 2003), from 2015 July 21 to July 23 (observation CN, 138 ks of exposure time, see Table 1) and from 2003 May 07 to 08 (observation C, 145 ks of exposure time, see Table 1). The HETGS spectrometer is composed of two grating types: the medium-energy gratings (MEGs), covering the 0.4–7 keV energy band, with an FWHM resolution of 0.023 Å, and the high-energy gratings (HEGs), having an FWHM resolution of 0.012 Å in the 0.8–10 keV band. Data processing was performed with the TGCat software (Huenemoerder et al. 2011), which employs Chandra Interactive Analysis of Observations tools (CIAO v.4.8; Fruscione et al. 2006) and Calibration Data base (CALDB v.4.8.0). Chandra/HETGS data for both periods were reduced in a standard way, using a narrow mask to avoid mask confusion above 6 keV. Plus and minus first-order (m = ±1) MEG and HEG data were extracted from the −1 and the +1 arms of the MEG and HEG gratings for the source and the background for both observations, using the CIAO tool tgextract. Spectral redistribution matrix files and effective area files were generated with mkgrmf and mkgarf.

PDS 456 was simultaneously observed by NuSTAR during observation CN of Chandra, from 2015 July 21 to 23 (74 ks exposure) and from July 24 to 25 (38 ks exposure). NuSTAR data were processed using the NuSTARDAS v1.7.1. The source spectra were extracted using a 45'' circular region centered on the source, and the background from a 45'' circular region clear of straylight in the same detector for both focal plane modules A and B (FPMA and FPMB). Observation reports depict a notable solar activity, which may impact the background event rate during both NuSTAR observations. Passage of solar coronal mass ejections over the Earth induces a temporary increase in the low Earth orbit radiation environment (which can persist for many orbits) and can significantly increase the background event level in the detectors of NuSTAR when the observatory is close to the South Atlantic Anomaly (SAA). We thus optimized the screening out of high-background periods near SAA passages using the module nucalcsaa.

PDS 456 was observed in 2013–2014 by NuSTAR and by XMM-Newton in five observations labeled XN1 to XN5 in Table 1 (see Nardini et al. 2015). We reduced the NuSTAR observations as described above, filtering the SAA passages when required, in order to compare our results to previous observations. We also reduced the XMM-Newton original data files using the XMM-Newton Standard Analysis Software (SAS v16.0.0 - Gabriel et al. 2004) considering the EPIC/pn (Strüder et al. 2001) spectrum for each observation of PDS 456. Events corresponding to flaring particle background were filtered using the SAS standard procedure. Single and double events were selected for extracting spectra. The data were screened for any increased flux of background particles. Spectra were extracted from a circular region of 30'' centered on the source. We checked for pile-up in all observations. The background was extracted from a nearby source-free region of 40'' in the same CCD as the source. Response matrices were generated for each source spectrum using the SAS arfgen and rmfgen tasks.

In order to increase the signal-to-noise ratio and because of the complex shape of the broadband spectrum around 9 keV (as described in Section 1), we combined the data to perform the spectral analysis, using the combine_data sets function in ISIS. For observation CN (see Table 1), we placed the HEG data on the same grid as MEG and combined all HETGS spectra (HEG-1, HEG+1, MEG-1, and MEG+1), and we also combined NuSTAR FPMA and FPMB data of the two observations 90101008002 and 90101008004. We used a combination of binning by signal-to-noise ratio and channel binning in order to obtain an adequate signal-to-noise ratio even below 1 keV (as discussed by Danehkar et al. 2018). We have thus binned Chandra and NuSTAR data with a minimum of one channel per bin (binned to the half-width at half-maximum of the MEG resolution) and a minimum signal-to-noise ratio of 2. We took into account the background. The same binning as for Chandra data has been used for observation C. As our C and CN observations have a Poisson distribution, we used Cash statistics (Cash 1979) to analyze these data. To estimate the error bars and evaluate the significance of our results, we used a ΔC value, which is appropriate for a likelihood estimate with Poisson statistics (corresponding to a given confidence level and a given number of degrees of freedom) that is analogous to Δχ² for the case of Gaussian statistics (see, e.g., the statistical textbook from Breiman 1973, works from Cash 1979 and Wilks 1938, and the pedagogical discussions in Arnaud et al. 2011).

For the simultaneous XMM-Newton and NuSTAR data set, we combined the pn spectra of the five observations XN1 to XN5, as well as the FPMA and FPMB spectra. We refer to this combined data set as XN. Despite the significant variability in the shape of the spectra that was described by Nardini et al. (2015), this combination of data is possible because of compatible fluxes for each observation, and because the absorption through around 9 keV is detected in all observations with similar parameters. We calculated the hardness ratios of the five observations (see Figure 2) using count rates CR in different energy bands:

$$\begin{eqnarray}&&{\mathrm{HR}}_{1}=\displaystyle \frac{{{CR}}_{5-10\ \mathrm{keV}}^{{XMM}}-{{CR}}_{0.5-4\ \mathrm{keV}}^{{XMM}}}{{{CR}}_{5-10\ \mathrm{keV}}^{{XMM}}+{{CR}}_{0.5-4\ \mathrm{keV}}^{{XMM}}},\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&{\mathrm{HR}}_{2}=\displaystyle \frac{{{CR}}_{15-30\ \mathrm{keV}}^{{NuSTAR}}-{{CR}}_{3-8\ \mathrm{keV}}^{{XMM}}}{{{CR}}_{15-30\ \mathrm{keV}}^{{NuSTAR}}+{{CR}}_{3-8\ \mathrm{keV}}^{{XMM}}},\end{eqnarray} \tag{ 2 }$$

in order to study the possible impact of spectral changes among the observations. By their definitions, HR₁ gives information about spectral variability regarding the whole EPIC/pn energy band, while HR₂ represents the variability of the continuum shape around the absorption features at 9 and 11 keV. In the top panel of Figure 2, the hardness ratio of the second XN observation is slightly lower, which is due to a small change in the spectral shape of the EPIC/pn spectrum, as shown in Figure 1 from Nardini et al. (2015). However, by comparing the spectra resulting from the sum of XN observations with and without this second observation, we did not notice any significant change in the spectral shape. In the bottom panel of Figure 2, we see that the hardness ratio that is representative of the variability around the absorption features is consistent among the five XN observations. The hardness ratios presented in Figure 2 show that all XN observations are consistent spectrally and can be combined for our analysis. The stability of HR₂ shows that the absorption features around 9 and 11 keV are not due to changes of the continuum. We thus combined the data in order to increase the signal-to-noise ratio and to obtain an average shape of the P Cygni-like profile that is detected in individual observations, which is to be compared with the CN observation. We binned NuSTAR data to a minimum of two channels per bin and a minimum signal-to-noise ratio of 2, and XMM-Newton data to a minimum of one channel per bin and a minimum signal-to-noise ratio of 5 (during the analysis described in the following sections, we checked that rebinning the data to a minimum of five channels per bin does not change the results on the absorption lines). The binning of these data is sufficient to obtain enough counts per channel to allow us to use χ² statistics for the analysis.

Zoom In Zoom Out Reset image size

To perform a precise analysis of the combined data sets described above, in particular the absorption features, we needed to determine the continuum carefully. Following the modeling of the continuum performed in previous studies, and in particular by Nardini et al. (2015), on the same XMM-Newton and NuSTAR data, we first determined the continuum for observation CN. We thus first fit the hard-energy band above 3 keV with a simple power law, taking into account Galactic absorption (with ${N}_{{\rm{H}}}^{\mathrm{gal}}$ fixed to 2.4 × 10²¹cm⁻² according to 21 cm measurements, Dickey & Lockman 1990; Kalberla et al. 2005), modeled with the tbabs absorption model (Wilms et al. 2000). We used a cross-calibration factor to allow us to simultaneously fit data from different satellites, but the photon index of the power law was the same. The fit resulted in a C statistic of 472.8 for 364 dof. In order to account for the spectral curvature, which has been found in previous studies to range between 2 and 5 keV (e.g., Reeves et al. 2009, 2014, 2018a; Turner & Miller 2009; Behar et al. 2010), we added a partially covering absorption to the power law using the zpcfabs model, which improved the fit significantly, with a ΔC of 15.4 (for two parameters of interest) and an ftest F-value of 4.05 and a p-value of 0.007. This partial covering by a moderately ionized absorber could well describe the continuum spectral variability, due to patchy obscuration, e.g., in NGC 5548 (Kaastra et al. 2014) or NGC 3516 (Turner et al. 2011). We also checked for the presence of neutral distant reflection by replacing the simple power law with a pexmon model. However, the improvement of the fits is less significant, with a ΔC of 11.1 (for three free parameters), an F-value of 2.97, but a p-value of 0.03. Furthermore, the spectra of our data sets do not show an excess around 30 keV or a fluorescence of the Fe Kα line at 6.4 keV. The spectral complexity of our data sets above 3 keV is well illustrated by the data-to-model ratios of observations CN and XN in Figure 3.

Zoom In Zoom Out Reset image size

When the soft part of the spectra is added and the continuum model determined in hard X-rays is extrapolated in the soft band, we clearly detect a soft excess below 1 keV, as well as strong emission and absorption features, mostly broadened, especially around 1 keV and above 7 keV (QSO frame energy). We modeled the soft excess with a broad Gaussian emission line (zgauss), that is, with a phenomenological model. We chose to fit the soft excess with a broad line to be consistent with the analysis from Nardini et al. (2015; because we reanalyzed the same data) and because using a more common blackbody component instead did not improve the fit. The resulting continuum model is

• cst ∗ tbabs∗zpcfabs ∗(power law+zgauss).

As the continuum model of observation CN is totally consistent with the analysis of XMM-Newton and NuSTAR data from Nardini et al. (2015), we applied it to our combined XN data. In order to analyze all our observations in a similar way, we also applied this model to our C observation, for which the continuum is determined on a smaller energy band. The parameters of the broadband fitting with this continuum model are given in Table 2.

Table 2.  Continuum Parameters

Parameters	Obs. CN	Obs. XN	Obs. C
Power law:
Photon index	${2.30}_{-0.012}^{+0.0048}$	${2.33}_{-0.008}^{+0.009}$	${2.47}_{-0.01}^{+0.01}$
${F}_{7-30\mathrm{keV}}$	${1.28}_{-0.011}^{+0.011}$	${1.24}_{-0.001}^{+0.001}$	${1.27}_{-0.013}^{+0.013}$
${F}_{0.4-30\mathrm{keV}}$	${12.27}_{-0.10}^{+0.10}$	${13.31}_{-0.0015}^{+0.0015}$	${16.49}_{-0.16}^{+0.16}$
Partial covering:
NH	${3.03}_{-0.53}^{+0.46}$	${8.46}_{-0.42}^{+0.49}$	${3.42}_{-0.11}^{+0.14}$
Covering factor	${0.32}_{-0.006}^{+0.009}$	${0.33}_{-0.009}^{+0.009}$	${0.75}_{-0.004}^{+0.004}$
Soft excess:
E	${0.060}_{-0.01}^{+0.008}$	${0.69}_{-0.007}^{+0.007}$	${0.14}_{-0.03}^{+0.05}$
σ	${0.10}_{-0.0017}^{+0.0017}$	${0.19}_{-0.004}^{+0.004}$	${0.06}_{-0.006}^{+0.003}$
${F}_{0.4-30\mathrm{keV}}$	${65.59}_{-15.45}^{+13.67}$	${0.97}_{0.007}^{+0.007}$	${9.43}_{-2.87}^{+3.44}$
Cross-calibration	${0.90}_{-0}^{+0.005}$	${0.97}_{-0.007}^{+0.007}$	⋯
C or ${\chi }^{2}$/dof	1693.48/1609	2240.58/1704	1409.60/1331

Note. Parameters of the continuum fit, performed above 0.4 keV, for the three observations, as described in Section 3.2. Fluxes are in units of 10⁻¹² erg s⁻¹ cm⁻².

Download table as:  ASCIITypeset image

After establishing the continuum model, we looked for absorption features above 7 keV. We fit the spectra above 3 keV, adding Gaussian and edge models to the base continuum model in order to constrain the emission and absorption profiles in the Fe–K band, similarly to Nardini et al. (2015). We first added gabs lines, one in emission (using a negative normalization) and two in absorption, and an edge, in order to account for the Fe xxvi Lyα emission and absorption lines (at E_lab = 6.97 keV), the Fe xxvi Lyβ absorption line (at E_lab = 8.25 keV) and the Fe xxvi K edge (at E_lab = 9.28 keV). We tied the widths and shifts of the absorption lines together. The width and shift of the emission line were also free to vary, independently of those of the absorption features. In this way, we were able to measure the outflow velocities of the UFO in emission and in absorption. We performed this fit for both CN and XN observations because the NuSTAR data helped so well to define the continuum. The significance of each line is given by ΔC for observation CN and Δχ² for observation XN, first adding the emission line to the continuum (with three parameters of interest), then the Lyα absorption line (three free parameters), then the Lyβ absorption line (linked in shift and width, so only one free parameter), and finally the edge (linked in shift, leaving one parameter of interest). The parameters of this fit are shown in Table 3. Our results show that the Fe xxvi Lyα emission and absorption line are significantly detected, in addition to the Fe xxvi Lyβ absorption line and the Fe xxvi K edge (at lower significance), with a slight emission blueshifted by −0.01 to −0.02 c and a large absorption blueshift of −0.27c. Note that we assessed the statistical significance of the lines using Δχ² and ΔC, but this method only gives an approximated significance according to Protassov et al. (2002).

Table 3.  Fits Performed above 3 keV with One Set of Fe Lines

Parameters		Obs. CN			Obs. XN
Line		Erest (keV)	EW (eV) or τ	ΔC/Δdof	Erest (keV)	EW (eV) or τ	Δχ2/Δdof
Fe xxvi Lyα	em.	7.04	${233}_{-46}^{+46}$	28.28/3	7.13	${91}_{-16}^{+16}$	96.50/3
	abs.	9.16	${259}_{-71}^{+71}$	15.55/3	9.15	${212}_{-23}^{+23}$	134.18/3
Fe xxvi Lyβ abs.		10.84	${142}_{-88}^{+88}$	3.09/1	10.83	${132}_{-29}^{+29}$	23.58/1
Fe xxvi K edge		12.20	${0.10}_{-0.05}^{+0.05}$	4.31/1	12.18	${0.06}_{-0.03}^{+0.03}$	5.06/1
${v}_{\mathrm{out}}/c$	em.	$-{0.010}_{-0}^{+0.02}$			$-{0.022}_{-0.012}^{+0.012}$
	abs.	$-{0.267}_{-0.023}^{+0.0312}$			$-{0.266}_{-0.0052}^{+0.0051}$
σ(keV)	em.	${0.600}_{-0.0791}^{+0}$			${0.442}_{-0.0769}^{+0.0874}$
	abs.	${0.457}_{-0.1568}^{+0.1432}$			${0.319}_{-0.0390}^{+0.0423}$
C or ${\chi }^{2}$/dof		379.09/358 = 1.059			1150.37/1179 = 0.9757

Note. Energies are in keV, velocities in units of c, and line widths σ are in keV. The equivalent width (EW, in eV) and the depth τ give information about the strength of the Gaussian line (in the first case) and of the edge (in the second case). The significance of each line is given by ΔC or Δχ².

Download table as:  ASCIITypeset image

In order to assess the presence of a second higher-velocity UFO as claimed by Reeves et al. (2018b), we added a second set of absorption lines (using two additional absorption Gaussian lines and another edge), with linked widths and shifts, these parameters being allowed to be different from those of the first UFO and of the emission component. The parameters of this fit are shown in Table 4. The results show that the Fe xxvi Lyα emission and absorption lines for both UFOs (v_out1 = −0.26 c and v_out2 = −0.47 c) are detected significantly, the Fe xxvi Lyβ absorption line from the second UFO being detected at a lower confidence level. However, the Fe xxvi Lyβ absorption line of the first UFO and the Fe xxvi K edges of both winds are not detected significantly. Furthermore, the width of the lines from the second UFO is hardly constrained in observation CN because it can only be established from NuSTAR data that have a lower spectral resolution than data from XMM-Newton. Considering the improvement of the fit when taking into account the second UFO, looking at the statistics (last line of Table 4), and at the residuals (in Figure 4), we assess that a second faster UFO is required by the data, significantly for observation XN (at about 97% confidence level), but only marginally for observation CN (slightly above 1σ confidence level).

Zoom In Zoom Out Reset image size

Table 4.  Fits Performed above 3 keV with Two Sets of Fe Lines

Parameters		Obs. CN			Obs. XN
Line		Erest (keV)	EW (eV) or τ	ΔC/Δdof	Erest (keV)	EW (eV) or τ	Δχ2/Δdof
UFO $1\&2$	Fe xxvi Lyα em.	7.11	${233}_{-46}^{+46}$	27.19/3	7.13	${97}_{-15}^{+15}$	104.29/3
UFO 1	Fe xxvi Lyα abs.	9.05	${231}_{-59}^{+59}$	16.29/3	9.13	${159}_{-18}^{+18}$	126.36/3
	Fe xxvi Lyβ abs.	10.71	${46}_{-0}^{+0}$	−0.58/1	10.80	${432}_{-11}^{+11}$	1.89/1
	Fe xxvi K edge	12.04	${0.052}_{-0.080}^{+0.042}$	1.36/1	12.15	${0}_{-0}^{+0.012}$	⋯
UFO 2	Fe xxvi Lyα abs.	11.50	${161}_{-76}^{+76}$	9.61/3	11.59	${161}_{-57}^{+57}$	28.68/3
	Fe xxvi Lyβ abs.	13.62	${161}_{-129}^{+129}$	2.98/1	13.72	${161}_{-77}^{+77}$	5.18/1
	Fe xxvi K edge	15.32	${0.021}_{-0.062}^{+0.011}$	0.09/1	15.43	${0}_{-0}^{+0.005}$	⋯
UFO 1&2	${v}_{\mathrm{out}}/c$ em.	$-{0.020}_{-0}^{+0.0229}$			$-{0.023}_{-0.0033}^{+0.0114}$
UFO 1	${v}_{\mathrm{out}}/c$ abs.	$-{0.255}_{-0.019}^{+0.021}$			$-{0.263}_{-0.0050}^{+0.0049}$
UFO 2	${v}_{\mathrm{out}}/c$ abs.	$-{0.463}_{-0.016}^{+0.006}$			$-{0.469}_{-0.0279}^{+0.0321}$
UFO $1\&2$	σ (keV) em.	${0.600}_{-0.071}^{+0}$			${0.449}_{-0.0710}^{+0.0808}$
UFO 1	σ (keV) abs.	${0.472}_{-0.131}^{+0.128}$			${0.282}_{-0.0438}^{+0.0462}$
UFO 2	σ (keV) abs.	${0.100}_{-0}^{+0.060}$			${0.560}_{-0.2465}^{+0.0400}$
C or ${\chi }^{2}$/dof		372.95/353 = 1.057			1137.82/1174 = 0.9692
ΔC or Δχ2/Δdof		6.14/5			12.55/5

Note. Parameters are defined as in Table 3. The last line of the table shows the significance of the fit improvements when the second UFO is considered in addition to the first.

Download table as:  ASCIITypeset image

The combination of the broad emission and the broad and blueshifted absorption Fe xxvi line has a P Cygni-like profile that is associated with the expansion of spherically symmetric stellar winds. Similarly to Nardini et al. (2015), we applied a model from Done et al. (2007), based on the Sobolev approximation with exact integration (SEI; Lamers et al. 1987), to reproduce Fe–K absorption features in this object. The velocity field $w=v/{v}_{\infty }$, i.e., the ratio between the wind velocity v and the terminal velocity ${v}_{\infty }$, is defined in the model by $w={w}_{0}+{(1-{w}_{0})(1-1/x)}^{\gamma }$, with w₀ the velocity at the photosphere and x = r/R₀ the radial distance r normalized to the photospheric radius R₀. The optical depth of the line is described by $\tau {(w)\approx {\tau }_{\mathrm{tot}}{w}^{{\alpha }_{1}}(1-w)}^{{\alpha }_{2}}$, with α₁ and α₂ characterizing the sharpness of the P Cygni profile. We first applied a P Cygni model to replace the Fe xxvi Lyα emission and absorption lines. We fixed the parameters γ and w₀ to 2 and 0.001, respectively (note that these parameters do not have any strong influence on the shape of the profile; by trying different values, we found that they are suitable for the fits on both CN and XN observations). We thus let five parameters free to vary during the fit: the characteristic energy of the profile E₀ (corresponding to the beginning of the absorption feature), the terminal velocity ${v}_{\infty }$, and the parameters τ_tot, α₁, and α₂. We then added another P Cygni model in order to fit the second UFO. We chose the same free parameters, but we tied the energy E₀ of the second P Cygni profile to that of the first P Cygni profile because a single Fe xxvi Lyα emission line was assumed for both UFOs (as also assumed in Section 3.3 and in the analysis from Reeves et al. 2018a, 2018b). This fit with the second UFO was thus made with nine free parameters. We performed the fits on observations CN and XN, and we show the results in Table 5.

Table 5.  Fits Performed above 3 keV with Two P Cygni Profiles

Parameters		Obs. CN	Obs. XN
P Cygni 1&2	E0 (keV)	${6.32}_{-0.05}^{+0.05}$	${6.35}_{-0.04}^{0.04}$
P Cygni 1	${v}_{\infty }/c$	−0.32${}_{-0.027}^{+0.012}$	−0.32${}_{-0.012}^{+0.049}$
	${\tau }_{\mathrm{tot}}$	${0.17}_{-0.03}^{+0.02}$	${0.13}_{-0.016}^{+0.001}$
	${\alpha }_{1}$	${0.93}_{-0.84}^{+1.33}$	${3.60}_{-0.37}^{+1.38}$
	${\alpha }_{2}$	${0.17}_{-0.42}^{+0.77}$	${1.58}_{-0.91}^{+0.91}$
	C or ${\chi }^{2}$/dof	389.6/361 = 1.08	1187.5/1182 = 1.005
P Cygni 2	${v}_{\infty }/c$	−0.52 ${}_{-0.015}^{+0.024}$	−0.53 ${}_{-0.022}^{+0.034}$
	${\tau }_{\mathrm{tot}}$	${0.04}_{-0.004}^{+0.004}$	${0.02}_{-0.005}^{+0.003}$
	${\alpha }_{1}$	${6.17}_{-7.18}^{+3.28}$	${4.13}_{-3.22}^{+1.87}$
	${\alpha }_{2}$	$-{0.42}_{-0.63}^{+0.68}$	$-{0.68}_{-0.17}^{+0.22}$
	C or ${\chi }^{2}$/dof	375.4/357 = 1.05	1170.3/1178 = 0.99
ΔC or Δχ2/Δdof		14.2/4	17.2/4

Download table as:  ASCIITypeset image

The statistics of the fits (last line of Table 5) as well as the residuals (plotted in Figure 5) show that the use of a second P Cygni profile significantly improves the fit in both observations CN (at about 99% confidence level) and XN (at more than 99% confidence level). This supports the hypothesis that a second UFO with a higher velocity than the first one is present in both data sets.

Zoom In Zoom Out Reset image size

After our line identification at high energy (above 3 keV), we added the soft parts of the data (down to 0.4 keV) to our analysis and studied the broad energy band spectra. We initially performed a blind line search and identified individual lines, as described in the Appendix. Together with Sections 3.3 and 3.4, which focused on individual high-energy features, this blind line search aimed at detecting individual lines as signatures of the outflows, particularly at lower energy. Our analysis strategy is to study both individual lines (model-independent analysis) and global models (model-dependent analysis) in order to check the consistency of our results using a dual approach. For this purpose, we implemented a self-consistent photoionization modeling to reproduce the emission and absorption features seen in our broadband spectra of the three data sets CN, XN (from 0.4 to 30 keV), and C (from 0.4 to 8 keV). To do so, we used the XSTAR photoionization code (version 2.39, Kallman et al. 1996, 2004, 2009; Kallman & Bautista 2001). However, instead of using XSTAR tabulated grids for our fits (as has been done in previous works, e.g., Nardini et al. 2015; Matzeu et al. 2017a; Reeves et al. 2018a, 2018b), we employed the analytic XSTAR model warmabs to reproduce the absorption lines and photemis to generate the emission lines. Using these analytic functions requires a larger amount of computation time in comparison to the calculation of XSTAR tabulated grids. However, it allows us to explore the entire range of values for the column density N_H, and the ionization parameter ξ = L_ion/(nr²) (with L_ion the ionizing luminosity, n the hydrogen density, and r the distance to the ionizing source), instead of discrete values that depend on the refinement of the grid. Furthermore, multiple parameters can freely vary during the fitting procedure, including the turbulent velocity.

Similarly to previous works (Matzeu et al. 2016b; Reeves et al. 2018b), we adopted two SEDs that are identical in the UV band. This assumption of an unchanged UV spectrum is supported by recent results from Reeves et al. (2018a). For observation XN, we used the SED reported in Nardini et al. (2015), derived from XMM-Newton EPIC-pn and OM , and from NuSTAR. For observations CN and C, we used an SED that slightly differs from the former only for the X-ray slope (Γ = 2.3 instead of 2.4, as found when fitting the continuum for CN data).

For our photoionization modeling, we added different wind emission and absorption components to the continuum of the three data sets, following the sequence described below and the order shown in the Table 6, which summarizes the results. We first replaced the partial covering wind zpcfabs with a warmabs model convolved with a partcov model (to account for the partial covering absorption required in all observations). To find the best fit, this partial covering absorber was variable in ionization, column density, covering factor, and velocity. Note that the outflowing velocities of the winds modeled with warmabs are derived from the redshift parameters, which were variable in the fitting process. We then included a photoemission component in this baseline model, having variable ionization, normalization, and velocity parameters, because emission lines such as the slightly blueshifted Fe xxvi Lyα have been detected, as described in the previous sections. We used free Doppler-shift parameters for both these models, allowing the emission and absorption features to be blueshifted independently. As seen in the statistics of the best-fit wind emission (ΔC and Δχ² in Table 6), the wind emission is significantly required by XN data and is marginally needed by CN and C data. A low-ionization warm absorber at the systematic Doppler shift is also required by XN data alone, as shown by the significant improvement of the fit when this new component is included in the previous model (see Δχ² value for the warm absorption in Table 6). We then applied a totally covering absorption to the model described above (including the new continuum, the wind emission, and the warm absorber for the XN observation alone), with a variable column density, linking the ionization parameter to that of the photoemission, and linking the blueshift to the partial covering absorption, so we assumed that the partially covering absorber and the fast wind are outflowing at the same velocity, but are independent in ionization parameters and column densities. Adding this first UFO strongly improves the fit for all observations, with ΔC = 21.4 for CN, Δχ² = 406.0 for XN, and ΔC = 19.4 for C, for one degree of freedom, demonstrating that the slowest UFO is significantly required in all our data sets at a confidence level higher than 99.9%. We then included a second total covering absorption in the model with the first UFO, with a varying column density and an independent blueshift (i.e., an independent velocity), and an ionization fixed to the value of the first UFO and the emission component. We see that the inclusion of this higher-velocity UFO in the model improves the fit for observations CN and XN (see Table 6), with ΔC = 8.0 for CN and Δχ² = 49.4 for XN, for two degrees of freedom, showing that the fastest UFO is significantly required at confidence levels of >95% and >99.9%, respectively.

Table 6.  Photoionization Models Applied to Our Data above 0.4 keV, and Corresponding Fit Parameters for Observations CN, XN, and C

Observation	Model
CN	cst∗tbabs∗warmabs[ufo1]∗warmabs[ufo2]∗(photemis+partcov∗warmabs[pc]∗powerlaw+zgauss)
XN	cst∗tbabs∗warmabs[wa]∗warmabs[ufo1]∗warmabs[ufo2]∗(photemis+partcov∗warmabs[pc]∗powerlaw+zgauss)
C	cst∗tbabs∗warmabs[ufo1]∗(photemis+partcov∗warmabs[pc]∗powerlaw+zgauss)

	Parameters	"CN"	"XN"	"C"
Cross-calibration	C	${0.90}_{-0}^{+0.014}$	${0.98}_{-0.007}^{+0.006}$	⋯
Galactic absorption	NHgal (1022 cm−2)	0.24 (f)	(f)	(f)
Power law	Γ	${2.21}_{-0.008}^{+0.012}$	${2.37}_{-0.003}^{+0.003}$	${2.5}_{-0.004}^{+0.001}$
	${F}_{0.4-30\mathrm{keV}}$ (10−12 erg s−1 cm−2)	${11.30}_{-0.10}^{+0.10}$	${15.33}_{-0.018}^{+0.018}$	${17.10}_{-0.17}^{+0.17}$
Soft excess	E (keV)	${0.40}_{-0.02}^{+0.07}$	${0.55}_{-0.007}^{+0.008}$	${0.087}_{-0.010}^{+0.041}$
	σ (keV)	${0.05}_{-0.014}^{+0.008}$	${0.24}_{-0.003}^{+0.003}$	${0.057}_{-0.0002}^{+0.0007}$
	${F}_{0.4-30\mathrm{keV}}$ (10−12 erg s−1 cm−2)	${62.9}_{-15.8}^{+13.8}$	${2.27}_{-0.011}^{+0.011}$	${9.14}_{-3.60}^{+2.98}$
Partial covering	log(${N}_{{\rm{H}}}/{10}^{22}$ cm−2)	${0.52}_{-0.13}^{+0.12}$	${1.22}_{-0.04}^{+0.01}$	${0.83}_{-0.08}^{+0.13}$
	log(ξ/erg cm s−1)	${3.20}_{-0.58}^{+0.18}$	${3.04}_{-0.21}^{+0.02}$	${2.89}_{-0.68}^{+0.08}$
	cf	${0.33}_{-0.054}^{+0.046}$	${0.40}_{-0.003}^{+0.001}$	${0.77}_{-0.04}^{+0.03}$
	vturb (km s−1)	20000 (f)	(f)	(f)
	vout (c)	$-{0.289}_{-0.030}^{+0.025}$	$-{0.268}_{-0.004}^{+0.007}$	$-{0.236}_{-0.064}^{+0.025}$
Wind emission	log(ξ/erg cm s−1)	${6.03}_{-0.28}^{+0.49}$	${6.97}_{-0.51}^{+0.20}$	${5.19}_{-0.09}^{+0.36}$
	norm ($\times {10}^{-2}$)	${8.5}_{-4.2}^{+5.2}$	${5.8}_{-2.2}^{+0.1}$	${5.3}_{-0.028}^{+0.092}$
	vturb (km s−1)	20000 (f)	(f)	(f)
	vout (c)	$-{0.056}_{-0.042}^{+0.006}$	$-{0.089}_{-0.001}^{+0.020}$	$-{0.026}_{-0.021}^{+0.002}$
	ΔC orΔχ2/Δdof	6.1/3	68.1/3	4.8/3
Warm absorption	log(${N}_{{\rm{H}}}/{10}^{22}$ cm−2)	⋯	$-{0.99}_{-0.002}^{+0.009}$	⋯
	log(ξ/erg cm s−1)	⋯	${0.73}_{-0.21}^{+0.05}$	⋯
	vturb (km s−1)	⋯	100 (f)	⋯
	vout (c)	⋯	0 (f)	⋯
	Δχ2 /Δdof	⋯	74.5/2	⋯
Wind absorption 1	log(${N}_{{\rm{H}}}/{10}^{22}$ cm−2)	${1.32}_{-0.22}^{+0.39}$	${1.90}_{-0.46}^{+0.05}$	${0.99}_{-0.24}^{+0.12}$
	log(ξ/erg cm s−1)	${6.03}_{-0.28}^{+0.49}$ (t)	${6.97}_{-0.51}^{+0.20}$(t)	${5.19}_{-0.09}^{+0.36}$ (t)
	vturb (km s−1)	20000 (f)	(f)	(f)
	vout (c)	$-{0.289}_{-0.030}^{+0.025}$ (t)	$-{0.268}_{-0.004}^{+0.007}$ (t)	$-{0.236}_{-0.064}^{+0.025}$ (t)
	C or ${\chi }^{2}$/dof	1660.2/1603 = 1.04	1860.7/1696 = 1.10	1384.2/1325 = 1.04
	ΔC or Δχ2/Δdof	21.4/1	406.0/1	19.4/1
Wind absorption 2	log(${N}_{H}/{10}^{22}$ cm−2)	${1.13}_{-1.31}^{+0.21}$	${1.71}_{-0.38}^{+0.09}$	⋯
	log(ξ/erg cm s−1)	${6.03}_{-0.28}^{+0.49}$ (t)	${6.97}_{-0.51}^{+0.20}$(t)	⋯
	vturb (km s−1)	20000 (f)	(f)	⋯
	vout (c)	$-{0.478}_{-0.094}^{+0.031}$	$-{0.483}_{-0.033}^{+0.001}$	⋯
	C or ${\chi }^{2}$/dof	1652.2/1601 = 1.03	1811.3/1694 = 1.07	⋯
	ΔC or Δχ2/Δdof	8.0/2	49.4/2	⋯

Note. (f) refers to parameters that have been frozen during the fits. (t) refers to parameters that have been tied to other parameters during the fits (such as the ionization parameters, and the velocities of the partially covering absorber and the slowest UFO). Each ΔC or Δχ² value refers to the improvement of the fit when the new component is added to the model including the components listed above in the table.

Download table as:  ASCIITypeset image

For these fits, we allowed the continuum parameters (power law, soft excess, and cross-calibration factor) to vary because the application of the warmabs models slightly changed the overall shape of the broadband spectra. Initially, we varied the turbulent velocity of each XSTAR component. We found that the turbulent velocity reached 20,000 km s⁻¹ for all components in all observations, so we fixed this parameter to this value. This high turbulent velocity is consistent with the large widths of the lines described in the previous sections, and with the values used in previous works with XSTAR grids. All abundances were fixed to solar values.

The fits were performed in ISIS using Cash statistics and the powell method for observations CN and C, and χ² statistics and the mpfit method for observation XN. Using the XSTARDB⁴ library of ISIS S-lang scripts, we list the strongest features predicted by our photoionization models with the fit parameters given in Table 6 for the three observations and for all components. For the UFOs, in addition to the expected Fe xxv Heα, Heβ, Fe xxvi Lyα, and Lyβ absorption lines, the XSTAR model predicts Si xiv and S xvi Lyα absorption lines at lower energy, consistent with the lines identified in the blind line search section (see appendix). Many emission features are also predicted by the photemis model, from Fe xxv Heα and Fe xxvi Lyα, but also from other ions (C vi, O viii, Ne x, Mg xii, Fe xxiv, Si xiv, and S xvi). The lower-ionization partially covering absorber is expected to show absorption features from O vii, O viii, Fe xvii, Fe xviii, Ne ix, Ne x, Mg xi, and Si xiii ions. Some of these lines have also been identified in the blind line search section. Non-blueshifted signatures from O vii and Ne vi ions are predicted for the warm absorber found in the XN observation. We also used a Markov chain Monte Carlo (MCMC) approach to derive the probability distribution of each parameter via Bayesian data analysis. For this purpose, we used the isis_emcee_hammer function⁵ developed in ISIS S-lang by MN. Best-fit values and error bars are reported in Table 6. Figure 6 shows the contour plots characterizing the absorbers in the three observations.

Zoom In Zoom Out Reset image size

We initially modeled the absorption features at high energy (above 3 keV) with blueshifted and highly ionized iron lines to characterize the UFOs in PDS 456. We found that the Fe xxvi Lyα emission and absorption lines, the Fe xxvi Lyβ absorption line, and the Fe K edge corresponding to the fit with only one UFO are significantly or marginally detected in both observations CN and XN. Velocities in emission and in absorption derived from this fit are consistent with results from Nardini et al. (2015; see Table 3). Applying a second set of lines to account for the second possible UFO marginally improved the fit, but not all the lines were detected significantly (see Table 4). This can be explained by the complexity of the model fit to data with an insufficient resolution. Indeed, signatures of the second UFO are located in the energy band that are predominantly covered by NuSTAR. However, the second UFO is required, more convincingly in observation XN than in observation CN, according to the residuals in Figure 4.

We replaced the sets of Gaussian lines with P Cygni profiles using the model from Done et al. (2007), as proposed by Hagino et al. (2015) and Nardini et al. (2015), describing a spherically symmetric outflowing wind. Applying two P Cygni models instead of a single model significantly improved the fits for both observation CN and XN (see the statistics in Table 5 and the residuals in Figure 5). Table 5 shows that the velocities implied by the P Cygni model are slightly higher (by about 15%–20%) than the velocity measured from individual Gaussian lines. The terminal velocity of the wind, ${v}_{\infty }$, is a characteristic parameter of the P Cygni model for a quasi-spherical fully covering outflow. This value is the actual speed of the outflowing gas. However, ${v}_{\infty }$ strongly depends on the P Cygni model used to fit our data, so we instead considered the line-of-sight bulk velocity, which is found to have consistent values when fitting with Gaussian lines and photoionization models, for our analysis and for energetic estimations. We found an optical depth for the first P Cygni profile that is consistent with results from Nardini et al. (2015), the depth of the second P Cygni profile being lower. Parameters α₁ and α₂ are different from values found by Nardini et al. (2015), possibly because they were free to vary independently in our analysis. Their values are also very different between the observations, as well as between the two P Cygni profiles.

For both observations CN and XN, the outflow velocities of both UFOs measured by one method or the other one are consistent between the two different epochs and are also consistent with previous results (e.g., Nardini et al. 2015; Matzeu et al. 2017a; Reeves et al. 2018b).

We applied self-consistent XSTAR photoionization models to our three different epochs of observations (using the entire detector band, i.e., 0.4–30 keV when possible) in order to characterize the physical properties of the absorbers reproducing different absorption lines detected in our blind line search (see Appendix) and in high-energy spectra (see Sections 3.3 and 3.4). The statistics of the fits (ΔC and Δχ² rows in Table 6) demonstrate that including the two UFO components in our modeling significantly improves the fits (at >99.9% confidence level for the slowest UFO, and at > 95 to >99.9% confidence level for the fastest UFO). The high values of column densities N_H and ionization parameters ξ obtained from our modeling are typical of UFOs (e.g., Tombesi et al. 2013) and consistent with previous studies of PDS 456.

As described in Section 3.5, we adopted a fixed value for the turbulent velocity, which helped to constrain the other parameters better because the model is complex with many free parameters. The turbulent velocity of 20,000 km s⁻¹ is consistent with the large widths of individual lines (e.g., see Table 7), and the values reported in previous studies (e.g., Nardini et al. 2015).

The ionization parameter of both UFOs was tied to that of the emission associated with the disk wind to facilitate the fits, as implemented in Nardini et al. (2015) and Reeves et al. (2018b), so the ionization parameters of both UFO are the same, but different from the ionization parameter of the partial covering absorber (even if the partially covering absorber has the same velocity as the slowest UFO). The high-ionization parameters derived from our modeling are consistent with the detection of highly ionized Fe xxvi lines and similar to previous results (e.g., Nardini et al. 2015; Reeves et al. 2018a, 2018b). It is also possible that two different sets of emission and absorption from different stratification layers of a single UFO are associated with the different ionizations. However, the available data would not be able to distinguish the different kinematic components, the overall model being in this case even more complex.

The velocity of the slowest UFO was tied to the velocity of the partially covering absorber. This scenario was adopted to obtain a better constraint. However, we saw in the appendix that when we identified lines that were detected by blind search, the transitions had a slightly different blueshift for lines from the partially covering absorber than from the slowest UFO, in particular in observation CN. It may imply that the actual velocities of both components are slightly different. The high UFO velocities derived from photoionization modeling (about v_out = −0.24 to −0.29 c for the slowest UFO, and v_out = −0.48 c for the fastest UFO) are consistent with previous studies, as well as with results from modeling of high-energy features (see Sections 3.3 and 3.4) and results from the blind line search (see the Appendix). An additional method for the characterization of the kinematic components of the UFOs is also described and applied to PDS 456 data in the Appendix, and the results are also consistent with values found through photoionization modeling.

The MCMC approach gave us the probability distribution of each parameter of the photoionization model, giving contours shown in Figure 6. These contour plots show the behavior of the ionization parameter as a function of the column density, as well as the different velocities in emission and in absorption, for all absorbers detected in the three observations, and the variation of the covering factor as a function of the column density for the partially covering absorber in the three observations. We find that the velocities are consistent between observations CN and XN, with a better constraint for observation XN than for observation CN because of the better signal-to-noise ratio. Despite the poor constraint on the velocity of the UFO detected in observation C (corresponding to the slowest one), we can see hints of a consistency with the other observations, with the best-fit value slightly higher than v_out = −0.2 c and contours coherent with an outflowing velocity of v_out = −0.2 to −0.3 c. Contours of the UFO in observation C are obviously improved when NuSTAR data are combined with Chandra/HETGS data, as in observation CN. However, despite the lack of information at high energy in observation C, the photoionization modeling allowed us to constrain the parameters of the winds at a lower significance level. The contours of all parameters, and in particular, the parameters of the fastest UFO, are more precise for observation XN than for observation CN (regarding velocities, ionization parameters, column densities, and covering factors), as expected because the signal-to-noise ratio in the NuSTAR energy band that we achieved by combining all observations XN1 to XN5 is higher than that from observation CN. The ionization parameters and column densities of the UFOs are relatively stable between the different epochs. The ionization of the partially covering absorber is lower than those of the UFOs, as expected from the detected absorption lines in the soft band, and this parameter is stable between observations, with a modest value (log(ξ) ∼ 3) that is consistent with previous studies (Nardini et al. 2015; Reeves et al. 2018a). The variability of the covering factor (c_f ∼ 0.3–0.8) and of the column density of this absorber (N_H = (3.3–16.6) × 10²² cm⁻²), as seen in the right panel of Figure 6, could be responsible for the continuum shape variability observed in all observations of PDS 456, similar to NGC 5548 (Kaastra et al. 2014; Nardini et al. 2015; Reeves et al. 2018a).

An alternative scenario proposed to explain blueshifted absorption features is that they could be the signatures of reflection on an optically thick plasma that may cover the accretion disk whose inner parts may rotate at extremely high velocities, producing a relativistic blueshifted Fe K-shell feature (see the example of PG 1211+143, Gallo & Fabian 2013, but also the counterargument from Lobban et al. 2016). Nardini et al. (2015) found that such a reflection-dominated scenario underpredicts the strength of the absorption feature by about 9 keV. Behar et al. (2010) found a persistent but small contribution of a reflection component in most observations before the XMM-Newton observation in 2007. A Suzaku observation from 2007 shows a marginally significant hard X-ray excess that could be modeled either by a strong reflection component or by a Compton-thick partially covering absorber (Reeves et al. 2009). Reprocessed and scattered X-ray emission off the surface of an accretion disk wind could explain the low flux and hard X-ray spectrum of the Suzaku observation from 2011, but this scenario is only applicable for this particular state (Reeves et al. 2014). There is no sign of dominating reflection in the spectra of observations CN and XN. Because adding a reflection component to the continuum did not significantly improve the fit (see Section 3.2), we tried to consider neutral reflection in our photoionization modeling (using pexmon) because the ionized absorbers provide a better constraint on the complicated spectral shape around 9 keV than the continuum model. However, adding this pexmon reflection component to the photoionization model with the two UFOs did not improve the fits (ΔC = 1.66 for observation CN and Δχ² = −12.3 for observation XN for three parameters of interest). Furthermore, the parameters of this reflection component, i.e., the reflection factor, the abundance, and the inclination, were poorly constrained during the fits for both observations, even through the MCMC approach.

To investigate the effects of the continuum on the ionization balance and thermal stability of the photoionized gas, we produced thermal stability curves, plotting the temperature of the plasma log(T) as a function of log(ξ/T) (Krolik et al. 1981; Reynolds & Fabian 1995; Krolik & Kriss 2001; Chakravorty et al. 2009; Lee et al. 2013). These curves represent the thermal equilibrium of the gas. On one side of the curve, cooling dominates, while on the other side, at high ionization and low temperature, heating dominates. Generally, a positive gradient of the curve indicates thermal stability for the gas (indeed, a small increase in temperature will increase the cooling, while a perturbation decreasing the temperature will increase the heating). Conversely, a negative gradient is present in regions of instability. We used the two SEDs described in Section 3.5 to produce the thermal stability curves shown in Figure 7. We recall here that these two SEDs are identical in the UV band and slightly differ in X-rays in order to match our three data sets. Note that the stability curve depends on the input assumptions provided to the photoionization code. Its shape is influenced by the ionizing SED (e.g., Lee et al. 2013; Mehdipour et al. 2015) and by the density and chemical composition of the absorber (e.g., Chakravorty et al. 2009). The unusual decrease in temperature at high ionization can be explained by the fact that the SED we used for the calculation of the XSTAR models has a strong soft component (as shown in Figure 5 from Matzeu et al. 2016b), and as the Compton temperature (T_IC = $\langle E\rangle$/4k) depends on the mean energy photon, the value of the Compton temperature is low. When the gas is fully ionized, it reaches this Compton temperature (Kallman & Bautista 2001). However, when the gas is partially ionized, it can reach hotter temperatures, explaining the peak shown in Figure 7. Indeed, in addition to the heating of the gas from the energetic electrons produced by photoionization, these electrons can heat the gas to a higher temperature through secondary collisional ionization of neutral atoms.

Zoom In Zoom Out Reset image size

We overplotted points representing the absorbers, whose ionization parameters have been characterized by the photoionization models in the three observations, on top of their respective stability curves. We can see that the partial covering absorber is in a stable state in all observations (pink points), in the portion of the curve with a positive slope. This is also the case for the non-blueshifted warm absorber (light blue rectangle) detected in observation XN. However, because the UFO components of observations XN, CN, and C are on a negative gradient branch (red points), we might expect that this gas is thermally unstable. However, the cooling time for the UFO gas is much shorter than the outflow time (R/v_out), so the gas should be stable in order for the UFO to be persistent and observable. In fact, the thermal stability of the gas can be determine by the slope of the branches, as described above, only in the case of S shapes. However, for our case, the thermal stability curve, with its unusual negative slope at high ionization, does not follow this rule. Indeed, as shown in Figure 7, at high ionization, cooling dominates above the curve and heating dominates below the curve. On this high-ionization branch, a small perturbation that increases the temperature will therefore bring the gas to a region where cooling dominates, and a small decrease in temperature will increase the heating. The UFOs located on this branch at high ionization are thus stable after all.

Figure 8 shows the distribution of the temperature as a function of the ionization parameter, as well as the distribution of ion fractions of the H-like and He-like ions of the elements identified in the three data sets, as described in previous sections (Fe, Si, S, Ne, and O). This distribution is dependent on the SED used by the photoionization code XSTAR, as we can see in the difference between solid and dashed lines (see the legend of the figure for more details). We overplot the values of the ionization parameters obtained by photoionization modeling, for both partial covering absorbers (in pink) and UFOs (in red), in the three observations (CN: dotted zone, XN: checkered zone, and C: hatched zone). We can see that the ions of Ne and O detected in the three observations can coexist in a single ionization zone at the same velocity, originating from the partial covering absorbers (pink zone). Iron ions are also detected in the three data sets, and can coexist with S and Si ions in a single ionization zone corresponding to the UFOs (red zone). This coexistence is more significant in observation C than in observations XN and CN. S and Si ions could either be produced by the UFOs or by the partially covering absorbers. This uncertain origin is induced by the fact that to simplify the model applied to our data, we linked the velocities of the partially covering absorber and of the slowest UFO together, as well as the ionization parameters between both UFOs. Higher-quality data are required to allow the modeling with untied parameters and hence the precise determination of the origin of the detected S and Si ions. In the present study, despite their dependence on the assumptions made for the fitting, the distributions of ion fractions of different elements show that their coexistence is possible, giving indications of their origins, and thus support the results from the photoionization modeling, as well as the identification of the lines detected by the blind line search.

Zoom In Zoom Out Reset image size

We found that two UFOs are significantly detected in our observations CN and XN, the detection being even more convincing in observation XN. The velocity of the slowest UFO obtained by photoionization modeling on observation XN (${v}_{\mathrm{out}1}=-{0.268}_{-0.004}^{+0.007}c$) is consistent with the average value measured on the five individual XN1 to XN5 observations by Nardini et al. (2015) (${v}_{\mathrm{out}1,\mathrm{average}}=-{0.25}_{-0.01}^{+0.01}c$).

The second UFO was not detected in the individual XN observations, as claimed by Reeves et al. (2018a, 2018b). Indeed, after finding a second, faster absorber in recent XMM-Newton and NuSTAR data from 2017, the authors reanalyzed the XN1 to XN5 observations, but did not find any signature of a second UFO above 10 keV. For our work, we chose to combine all the XN observations together, as explained in Section 3.1, as we wished to achieve the optimum signal-to-noise ratio using all available data. We verified the detection of the fastest UFO in the individual XN observations by applying the fitting procedure described in Section 3.3 on these five data sets, using two sets of Fe lines to look for signatures of both UFOs. Regarding the improvement of the fit statistics when considering two UFOs instead of one, and according to the significant detection of the Fe xxvi Lyα, Lyβ, and K-edge features from the fastest UFO, we found that the second UFO is required by the data XN5 (Δχ²/Δdof = 19.03/5, >99.9% confidence level), XN2, (Δχ²/Δdof = 15.82/5, about 99.5% confidence level) and XN4 (Δχ²/Δdof = 10.16/5, >90% confidence level). However, we did not detect the second UFO in observations XN3 and XN1. Looking at their recent observations from 2017, Reeves et al. (2018b) detected the second UFO with a significance >99.9% (the addition of the Lyα line from the fastest UFO improving the fit by Δχ² = 39.3 for two free parameters). The detection of the lines from the second UFO is less significant in XN5, XN2, and XN4 data than in the more recent data from 2017; this is probably why Reeves et al. (2018a, 2018b) considered that this fastest component was not visible in the individual XN observations. According to the count rate of each observation shown in Table 1 and to the lower flux of the recent data from 2017, there seems to be a hint of an anticorrelation between the statistical confidence level for the detection of the fastest UFO and the flux intensity of PDS 456. This trend is consistent with the proposition from Reeves et al. (2018a), who suggested that PDS 456 has to be in a low state to allow the detection of the second UFO (see further discussion in Section 4.5).

In our observations CN and XN, the column density of the fastest UFO is lower than the density of the slowest UFO (see Table 6), consistent with results from Reeves et al. (2018b).

Figure 9 shows the outflow velocities (top panel) and the ionization parameters (bottom panel) of the UFOs detected in our analysis and in the XMM-Newton/NuSTAR data from 2017 (Reeves et al. 2018b) as a function of the intrinsic flux of PDS 456 between 2 and 10 keV. We can see that the velocities are similar between the different observations. In the study of 12 previous X-ray observations of PDS 456, Matzeu et al. (2017a) depicted a strong correlation between the outflow velocity and the X-ray luminosity, which supports the hypothesis of a radiatively driven wind in PDS 456, boosted by line driving (Hagino et al. 2016). The recent detection of the C iv BAL at the velocity of 0.3 c in HST UV observations may also boost the opacities in UFOs for radiative driving (Hamann et al. 2018). This hypothesis of radiative driving is also strongly supported by the fact that the SMBH of PDS 456 is accreting at a regime near the Eddington limit (King & Pounds 2003; Gofford et al. 2014; Matzeu et al. 2017a). We checked the correlation between velocity and luminosity of the data shown in Figure 9. We only found a possible small correlation between intrinsic flux and velocity of the fastest UFO (Pearson correlation coefficient r = 0.995, p-value = 0.06), and no correlation between the flux and the velocity of the slowest UFO (r = 0.11, p-value = 0.89). The outflowing velocities seem rather to have stable values. We only found a non-significant correlation for the expected positive relation between the ionizing flux and the ionization parameter (r = 0.40, p-value = 0.60), probably because of the less precise constraint on the slowest UFO with Chandra data from observation C compared to the parametric constraints derived from additional high-energy data available in observations CN and XN.

Zoom In Zoom Out Reset image size

The mass outflow rate is

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{out}}\sim {\rm{\Omega }}{N}_{{\rm{H}}}{m}_{p}{v}_{\mathrm{out}}{R}_{\mathrm{in}},\end{eqnarray} \tag{ 3 }$$

where Ω is the solid angle, N_H is the column density, m_p is the proton mass, v_out is the outflowing velocity, and R_in is the starting point of the wind (see details, e.g., in Nardini et al. 2015). We used an approximated value of Ω = 2π for the solid angle, as justified by Nardini et al. (2015; see their discussion in the supplementary material for further details, as well as the argumentation from Reeves et al. 2018a). The column density is directly derived from the photoionization modeling, as well as the outflow velocity (see Table 6). In the case of a radiatively accelerated wind (e.g., Matzeu et al. 2017a), R_in can be approximated as

$$\begin{eqnarray}&&{R}_{\mathrm{in}}\sim 2\left(\alpha \displaystyle \frac{L}{{L}_{\mathrm{Edd}}}-1\right){\left(\displaystyle \frac{{v}_{\infty }}{c}\right)}^{-2},\end{eqnarray} \tag{ 4 }$$

where ${v}_{\infty }$ is the wind terminal velocity and α is a force multiplier factor (Reeves et al. 2018b). PDS 456 is accreting at about the Eddington limit, so we could approximate L/L_Edd = 1. Using a factor of α = 2 (as done by Reeves et al. 2018b) and the outflow velocities derived from photoionization modeling for each of our three observations, we find that the launching radius is ∼5 × 10¹⁵ cm (or 30R_g for M_BH = 10⁹ M_⊙) for the slowest UFO and ∼1 × 10¹⁵ cm (or 9R_g) for the fastest UFO. These values are close to the escape radii ${R}_{\mathrm{escape}}=\tfrac{2{GM}}{{v}_{\mathrm{out}}^{2}}$, i.e., the minimum radii from which winds of a given outflow velocity can be launched. We are probably observing a stratified UFO in PDS 456, with several components from multiple stratification layers having different velocities and being launched from the accretion flow close to the SMBH (e.g., Tombesi et al. 2013; Reeves et al. 2018b). Reeves et al. (2018b) proposed that PDS 456 has to be in a low state for the source to not be extremely luminous (and the iron fully ionized), to allow the detection of the second UFO at such short distances. This trend seems to be verified considering the high confidence level of the detection of the second UFO in individual XN observations as a function of the flux, as explained in Section 4.4. An alternative explanation could be that a partially covering dense gas (as the one already observed in previous studies) is shielding the innermost wind, preventing it from becoming too highly ionized (Matzeu et al. 2016a; Reeves et al. 2018b).

The maximum radial distance of the absorbers can be estimated by considering that ΔR/R < 1, i.e., their thickness cannot exceed their distance from the ionizing source (e.g., Reeves et al. 2003, 2018b). With the definition of the ionization parameter given in Section 3.5, R_max < L_ion/N_Hξ. This gives maximum radii of 7 × 10¹⁷ cm (or 4700R_g) for observation CN, 2 × 10¹⁶ cm (or 130R_g) for observation XN, and 6 × 10¹⁸ cm (or 40,000R_g) for observation C. Thus, the outflows in PDS 456 may extend to the broad line region.

For the slowest UFO, we found a mass outflow rate of 1.2 M_⊙ yr⁻¹ = $0.05{\dot{M}}_{\mathrm{Edd}}$ for observation CN, 4.4 M_⊙ yr⁻¹ = $0.2{\dot{M}}_{\mathrm{Edd}}$ for observation XN, and 0.6 M_⊙ yr⁻¹ = $0.03{\dot{M}}_{\mathrm{Edd}}$ for observation C. For the fastest UFO, we found $0.4\,{M}_{\odot }\,{{yr}}^{-1}\,=0.02{\dot{M}}_{\mathrm{Edd}}$ for observation CN, and 1.6 M_⊙ yr⁻¹ = $0.07{\dot{M}}_{\mathrm{Edd}}$ for observation XN. We found a different value for the mass outflow rate of the slowest UFO for observation XN compared to previous studies (with, e.g., an estimated kinetic power and mass outflow rate of ∼15% and ∼50%, respectively, of Eddington values for Nardini et al. 2015, and 5% and 40% for Gofford et al. 2014) because we obtained slightly different parameters resulting from the photoionization fit and because we used a lower value for the inner radius, calculated in the case of radiatively driven winds, while others used a timing approach to estimate it at a few hundred gravitational radii (Reeves et al. 2009; Nardini et al. 2015; Matzeu et al. 2016a, 2016b). We see that the mass outflow rate of the fastest UFO is lower than the rate from the slowest UFO. We calculated the kinetic power of the winds (${P}_{{kin}}=0.5\dot{M}{v}_{\mathrm{out}}^{2}$) and found kinetic powers of $0.02{L}_{\mathrm{Edd}}$ for observation CN, $0.07\mbox{--}0.08{L}_{\mathrm{Edd}}$ for observation XN, and $0.008{L}_{\mathrm{Edd}}$ for observation C (again slightly lower than in previous works for the same reason as above). Both winds are found to have a similar kinetic energy, which contradicts the prediction of a higher power for the faster UFO, as estimated by Reeves et al. (2018b) using the approximation that ${P}_{\mathrm{kin}}\propto {v}_{\mathrm{out}}^{3}$. The high velocities and high column density characterizing the UFOs detected in our three observations result in a large amount of kinetic power of about 0.8%–8% of the bolometric luminosity, which is sufficient to induce significant AGN feedback according to models of black hole and host galaxy coevolution (King & Pounds 2003; Di Matteo et al. 2005; Hopkins & Elvis 2010).