Relativistic components of the Ultra-Fast Outflow in the Quasar PDS 456 from Chandra/HETGS, NuSTAR and XMM-Newton observations

TL;DR

This study analyzes multi-epoch X-ray observations of quasar PDS 456, revealing a persistent ultra-fast outflow with relativistic velocities and a new faster component, providing insights into AGN feedback mechanisms.

Contribution

First detection of a relativistic outflow component at -0.48c in PDS 456 across multiple epochs, with detailed photoionization modeling of the outflow properties.

Findings

Confirmed persistent ultra-fast outflow at -0.24 to -0.29c

Detected an additional faster outflow component at -0.48c

Outflows are powerful enough to influence host galaxy evolution

Abstract

We present the spectral analysis of Chandra/HETGS and NuSTAR observations of the quasar PDS 456 from 2015, and XMM-Newton and NuSTAR archival data from 2013-2014, together with Chandra/HETGS data from 2003. We analyzed these three different epochs in a consistent way, looking for absorption features corresponding to highly ionized blueshifted absorption lines from H-like and He-like ions of iron (and nickel), as well as of other elements (O, Ne, Si, and S) in the soft band. We confirm the presence of a persistent ultra-fast outflow (UFO) with a velocity of v_out=-0.24 - -0.29c, previously detected. We also report the detection of an additional faster component of the UFO with a relativistic velocity of v_out=-0.48c. We implemented photoionization modeling, using XSTAR analytic model warmabs, to characterize the physical properties of the different kinematic components of the ultra-fast outflow and of the partial covering absorber detected in PDS 456. These two relativistic components of the ultra-fast outflow observed in the three epochs analyzed in this paper are powerful enough to impact the host galaxy of PDS 456 through AGN feedback.

Tables

Table 1. Table 1 Observation log

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

Table 2. Table 2 Continuum parameters

Parameters	Obs. CN	Obs. XN	Obs. C
Powerlaw:
Photon Index	${2.30}_{-0.012}^{+0.0048}$	${2.33}_{-0.008}^{+0.009}$	${2.47}_{-0.01}^{+0.01}$
$F_{7-30keV}$	${1.28}_{-0.011}^{+0.011}$	${1.24}_{-0.001}^{+0.001}$	${1.27}_{-0.013}^{+0.013}$
$F_{0.4-30keV}$	${12.27}_{-0.10}^{+0.10}$	${13.31}_{-0.0015}^{+0.0015}$	${16.49}_{-0.16}^{+0.16}$
Partial covering:
$N_H$	${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}$
$\sigma$	${0.10}_{-0.0017}^{+0.0017}$	${0.19}_{-0.004}^{+0.004}$	${0.06}_{-0.006}^{+0.003}$
$F_{0.4-30keV}$	${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

Table 3. Table 3 Fits performed above 3 keV with one set of Fe lines

Parameters		Obs. CN			Obs. XN
Line		$E_{rest}$ (keV)	EW (eV) or $\tau$	$\Delta$C / $\Delta$dof	$E_{rest}$ (keV)	EW (eV) or $\tau$	$\Delta {\chi }^2$ / $\Delta$dof
Fe xxvi Ly$\alpha$	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$\beta$ 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_{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}$
$\sigma$ (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

Table 4. Table 4 Fits performed above 3 keV with two sets of Fe lines

Parameters		Obs. CN			Obs. XN
Line		$E_{rest}$ (keV)	EW (eV) or $\tau$	$\Delta$C / $\Delta$dof	$E_{rest}$ (keV)	EW (eV) or $\tau$	$\Delta {\chi }^2$ / $\Delta$dof
UFO 12	Fe xxvi Ly$\alpha$ em.	$7.11$	${233}_{-46}^{+46}$	27.19 / 3	$7.13$	${97}_{-15}^{+15}$	104.29 / 3
UFO 1	Fe xxvi Ly$\alpha$ abs.	$9.05$	${231}_{-59}^{+59}$	16.29 / 3	$9.13$	${159}_{-18}^{+18}$	126.36 / 3
	Fe xxvi Ly$\beta$ 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$\alpha$ abs.	$11.50$	${161}_{-76}^{+76}$	9.61 / 3	$11.59$	${161}_{-57}^{+57}$	28.68 / 3
	Fe xxvi Ly$\beta$ 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 12	$v_{out}/c$ em.	$-{0.020}_{-0}^{+0.0229}$			$-{0.023}_{-0.0033}^{+0.0114}$
UFO 1	$v_{out}/c$ abs.	$-{0.255}_{-0.019}^{+0.021}$			$-{0.263}_{-0.0050}^{+0.0049}$
UFO 2	$v_{out}/c$ abs.	$-{0.463}_{-0.016}^{+0.006}$			$-{0.469}_{-0.0279}^{+0.0321}$
UFO 12	$\sigma$ (keV) em.	${0.600}_{-0.071}^{+0}$			${0.449}_{-0.0710}^{+0.0808}$
UFO 1	$\sigma$ (keV) abs.	${0.472}_{-0.131}^{+0.128}$			${0.282}_{-0.0438}^{+0.0462}$
UFO 2	$\sigma$ (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
$\Delta$C or $\Delta {\chi }^2$ / $\Delta$dof		6.14 / 5			12.55 / 5

Table 5. Table 5 Fits performed above 3 keV with two P Cygni profiles

Parameters		Obs. CN	Obs. XN
P Cygni 12	$E_0$ (keV)	${6.32}_{-0.05}^{+0.05}$	${6.35}_{-0.04}^{0.04}$
P Cygni 1	$v_{\infty }/c$	-0.32${}_{}{}^{+0.012}{}_{-0.027}{}^{}$	-0.32${}_{}{}^{+0.049}{}_{-0.012}{}^{}$
	${\tau }_{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.024}{}_{-0.015}{}^{}$	-0.53 ${}_{}{}^{+0.034}{}_{-0.022}{}^{}$
	${\tau }_{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
$\Delta$C or $\Delta {\chi }^2$ / $\Delta$dof		14.2 / 4	17.2 / 4

Table 6. 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)

Table 7. Table 7 Blind line search, performed above 0.4 keV, on observation “CN”.

${\text{E}}_{\text{rest}}$	$\sigma$	EW (eV)	$\Delta$C	ID	${\text{E}}_{\text{lab}}$	${\text{z}}_{\text{out}}$
${0.69}_{-0.002}^{+0.001}$	${0.005}_{-0.0019}^{+0.0059}$	$220.00\pm 164.00$	30.39	O viii Ly$\alpha$ ($v_{em}$)	0.65	$-{0.058}_{-0.005}^{+0.009}$
${0.94}_{-0.004}^{+0.004}$	${0.007}_{-0.0032}^{+0.0043}$	$15.90\pm 6.99$	9.43	O viii edge/rrc ($v_{em}$)	0.87	$-{0.074}_{-0.007}^{+0.008}$
${1.04}_{-0.089}^{+0}$	${0.001}_{-0}^{+0.001}$	$6.52\pm 2.80$	8.10
${1.06}_{-0}^{+0.001}$	${0.001}_{-0.0004}^{+0.0017}$	$5.57\pm 2.78$	6.71
${1.10}_{-0.001}^{+0.004}$	${0.001}_{-0}^{+0.0006}$	$-5.95\pm 1.55$	9.01	O viii Ly$\beta$ ($v_{abs1}$, pc)	0.77	$-{0.300}_{-0.001}^{+0.003}$
${1.19}_{-0.001}^{+0.001}$	${0.001}_{-0.0001}^{+0.0013}$	$-2.95\pm 0.98$	7.08	Fe xvii 2p-3d ($v_{abs1}$, pc)	0.83	$-{0.303}_{-0.001}^{+0.001}$
${1.31}_{-0}^{+0.003}$	${0.005}_{-0.0021}^{+0.0021}$	$-5.34\pm 1.09$	7.19	Ne ix He$\alpha$ ($v_{abs1}$, pc)	0.92	$-{0.298}_{-0.001}^{+0.003}$
${1.37}_{-0.001}^{+0.001}$	${0.002}_{-0.001}^{+0.0013}$	$-2.78\pm 0.90$	6.48	Ne x Ly$\alpha$ ($v_{abs1}$, pc)	1.02	$-{0.255}_{-0.001}^{+0.001}$
${1.76}_{-0.001}^{+0.001}$	${0.001}_{-0}^{+0.0008}$	$1.93\pm 0.77$	6.93
${1.84}_{-0.001}^{+0.001}$	${0.001}_{-0.0004}^{+0.0009}$	$3.08\pm 1.23$	12.18
${2.93}_{-0.002}^{+0.003}$	${0.001}_{-0}^{+0.0031}$	$7.24\pm 3.15$	6.70	Si xiv edge/rrc ($v_{em}$)	2.67	$-{0.089}_{-0.001}^{+0.002}$
${3.59}_{-0.005}^{+0.003}$	${0.002}_{-0.0006}^{+0.0014}$	$-9.77\pm 3.15$	7.52	S xvi Ly$\alpha$ ($v_{abs1}$)	2.62	$-{0.270}_{-0.001}^{+0.001}$
${5.37}_{-0.002}^{+0.008}$	${0.001}_{-0.0002}^{+0.0535}$	$17.50\pm 5.75$	8.26
${5.80}_{-0.056}^{+0.032}$	${0.071}_{-0.0426}^{+0.0297}$	$49.60\pm 20.90$	10.71
${5.96}_{-0.014}^{+0.009}$	${0.006}_{-0.0025}^{+0.0087}$	$-42.10\pm 0.11$	11.52
${6.91}_{-0.337}^{+0.003}$	${0.600}_{-0.0528}^{+0}$	$238.00\pm 22.40$	92.71	Fe xxv He$\alpha$ ($v_{em}$)	6.70	$-{0.030}_{-0.058}^{+0.000}$
				+Fe xxvi Ly$\alpha$ ($v_{em}$)	6.97	${0.009}_{-0.060}^{+0.000}$
${9.11}_{-0.425}^{+0.238}$	${0.434}_{-0.426}^{+0.1607}$	$-169.00\pm 0.00$	7.74	Fe xxv He$\alpha$ ($v_{abs1}$)	6.70	$-{0.265}_{-0.076}^{+0.031}$
				+Fe xxvi Ly$\alpha$ ($v_{abs1}$)	6.97	$-{0.235}_{-0.079}^{+0.032}$
${11.34}_{-0.082}^{+0.075}$	${0.016}_{-0.0065}^{+0.146}$	$-112.00\pm 0.00$	8.94	Fe xxv He$\alpha$ ($v_{abs2}$)	6.70	$-{0.409}_{-0.005}^{+0.011}$
				+Fe xxvi Ly$\alpha$ ($v_{abs2}$)	6.97	$-{0.385}_{-0.005}^{+0.012}$
				+Fe xxvi Ly$\beta$ ($v_{abs1}$)	7.88	$-{0.305}_{-0.005}^{+0.013}$
				+Ni xxviii Ly$\alpha$ ($v_{abs1}$)	8.11	$-{0.285}_{-0.006}^{+0.014}$
${13.76}_{-0.1}^{+0.061}$	${0.024}_{-0.0158}^{+0.0783}$	$-161.00\pm 0.62$	6.88	Fe xxvi Ly$\beta$ ($v_{abs2}$)	7.88	$-{0.427}_{-0.005}^{+0.006}$
				+Ni xxviii Ly$\alpha$ ($v_{abs2}$)	8.11	$-{0.411}_{-0.005}^{+0.006}$
${15.18}_{-0.633}^{+0.127}$	${0.349}_{-0.0781}^{+0.2506}$	$-28.30\pm 16.10$	7.18

Table 8. Table 8 Blind line search, performed above 5 keV, on observation “XN”.

${\text{E}}_{\text{rest}}$	$\sigma$	EW (eV)	$\Delta {\chi }^2$	ID	${\text{E}}_{\text{lab}}$	${\text{z}}_{\text{out}}$
${7.19}_{-0.068}^{+0.069}$	${0.549}_{-0.055}^{+0.050}$	124.64 $\pm$10.84	124.90	Fe xxv He$\alpha$ ($v_{em}$)	6.70	$-{0.068}_{-0.016}^{+0.015}$
				+ Fe xxvi Ly$\alpha$ ($v_{em}$)	6.97	$-{0.031}_{-0.017}^{+0.016}$
${9.22}_{-0.082}^{+0.066}$	${0.317}_{-0.023}^{+0.031}$	-174.48 $\pm$7.62	208.38	Fe xxv He$\alpha$ ($v_{abs1}$)	6.70	$-{0.273}_{-0.008}^{+0.008}$
				+Fe xxvi Ly$\alpha$ ($v_{abs1}$)	6.97	$-{0.244}_{-0.009}^{+0.008}$
${11.24}_{-0.131}^{+0.069}$	${0.593}_{-0.091}^{+0.007}$	-173.50 $\pm$25.12	58.42	Fe xxvi Ly$\beta$ ($v_{abs1}$)	7.88	$-{0.299}_{-0.014}^{+0.005}$
				+ Ni xxviii  Ly$\alpha$ ($v_{abs1}$)	8.11	$-{0.278}_{-0.015}^{+0.005}$
				+Fe xxv He$\alpha$ ($v_{abs2}$)	6.70	$-{0.404}_{-0.012}^{+0.004}$
				+Fe xxvi Ly$\alpha$ ($v_{abs2}$)	6.97	$-{0.449}_{-0.012}^{+0.004}$
${12.64}_{-0.154}^{+0.062}$	${0.005}_{-0.003}^{+0.003}$	-45.17 $\pm$23.88	6.79
${13.62}_{-0.065}^{+0.063}$	${0.009}_{-0.002}^{+0.002}$	-106.03 $\pm$25.86	18.15	Fe xxvi Ly$\beta$ ($v_{abs2}$)	7.88	$-{0.422}_{-0.003}^{+0.003}$
				+ Ni xxviii Ly$\alpha$ ($v_{abs2}$)	8.11	$-{0.405}_{-0.003}^{+0.003}$

Table 9. Table 9 Blind line search, performed above 0.4 keV, on observation “C”.

${\text{E}}_{\text{rest}}$	$\sigma$	EW (eV)	$\Delta$C	ID	${\text{E}}_{\text{lab}}$	${\text{z}}_{\text{out}}$
${0.76}_{-0.005}^{+0.005}$	${0.008}_{-0.0069}^{+0.0074}$	$-16.8\pm 37.7$	8.88	O vii He$\alpha$ ($v_{abs1}$, pc)	0.57	$-{0.250}_{-0.012}^{+0.012}$
${0.87}_{-0.003}^{+0.001}$	${0.001}_{-0}^{+0.0005}$	$-5.34\pm 0.96$	6.98	O viii Ly$\alpha$ ($v_{abs1}$, pc)	0.65	$-{0.253}_{-0.003}^{+0.001}$
${0.88}_{-0.002}^{+0.002}$	${0.004}_{-0.0012}^{+0.0021}$	$9.97\pm 4.09$	11.01	O viii edge/rrc ($v_{em}$)	0.87	$-{0.011}_{-0.004}^{+0.005}$
${1.25}_{-0.022}^{+0.042}$	${0.102}_{-0.0258}^{+0.0302}$	$-23.90\pm 0.00$	32.53	Ne ix He$\alpha$ ($v_{abs1}$, pc)	0.92	$-{0.264}_{-0.029}^{+0.040}$
				+Ne x Ly$\alpha$ ($v_{abs1}$, pc)	1.02	$-{0.184}_{-0.032}^{+0.045}$
${1.31}_{-0.001}^{+0}$	${0.001}_{-0}^{+0}$	$3.80\pm 0.02$	8.54
${1.39}_{-0.001}^{+0.001}$	${0.002}_{-0.0006}^{+0.0009}$	$3.64\pm 1.64$	7.68	Ne x edge/rrc ($v_{em}$)	1.36	$-{0.022}_{-0.001}^{+0.001}$
${1.47}_{-0.001}^{+0.001}$	${0.002}_{-0.001}^{+0.0013}$	$-3.66\pm 1.28$	8.50
${1.93}_{-0.001}^{+0.001}$	${0.001}_{-0}^{+0.0011}$	$-2.67\pm 0.82$	7.18
${2.22}_{-0.064}^{+0.215}$	${0.204}_{-0.1517}^{+0.086}$	$14.30\pm 8.07$	9.55
${2.45}_{-0.003}^{+0.003}$	${0.004}_{-0.0029}^{+0.0036}$	$-5.98\pm 6.81$	6.43
${2.55}_{-0.001}^{+0.001}$	${0.001}_{-0}^{+0}$	$7.03\pm 2.76$	8.35
${2.60}_{-0.001}^{+0.002}$	${0.001}_{-0}^{+0.0004}$	$-6.46\pm 1.91$	8.97	Si xiv Ly$\alpha$ ($v_{abs1}$)	2.01	$-{0.227}_{-0.0003}^{+0.0007}$
${2.65}_{-0.001}^{+0.001}$	${0.001}_{-0}^{+0.0023}$	$10.30\pm 4.00$	9.51
${3.08}_{-0.002}^{+0.002}$	${0.001}_{-0}^{+0}$	$-8.23\pm 0.17$	6.64
${3.21}_{-0.004}^{+0.004}$	${0.001}_{-0}^{+0.0112}$	$11.30\pm 4.63$	8.51
${3.47}_{-0.007}^{+0.192}$	${0.033}_{-0}^{+0.1709}$	$-22.60\pm 32.80$	9.84	S xvi Ly$\alpha$ ($v_{abs1}$)	2.62	$-{0.245}_{-0.002}^{+0.071}$
${3.55}_{-0.328}^{+0.003}$	${0.001}_{-0.0003}^{+0.0004}$	$-9.47\pm 1.06$	6.26
${3.92}_{-0.004}^{+0.004}$	${0.005}_{-0.0039}^{+0.0063}$	$-11.70\pm 25.40$	10.76
${5.69}_{-0.005}^{+0.005}$	${0.003}_{-0.001}^{+0.0035}$	$-23.40\pm 0.00$	7.91
${7.03}_{-0}^{+0}$	${0.002}_{-0.0007}^{+0.0225}$	$54.90\pm 39.00$	7.14	Fe xxv He$\alpha$ ($v_{em}$)	6.70	$-{0.047}_{-0.0001}^{+0.003}$
${7.35}_{-0.01}^{+0.017}$	${0.013}_{-0.0125}^{+0.027}$	$90.30\pm 80.00$	11.03	Fe xxvi Ly$\alpha$ ($v_{em}$)	6.97	$-{0.052}_{-0.0029}^{+0.006}$

Equations

$${\text{HR}}_1=\frac{C{R}_{\text{ 5-10 keV}}^{\text{ XMM}}-C{R}_{\text{ 0.5-4 keV}}^{\text{ XMM}}}{C{R}_{\text{ 5-10 keV}}^{\text{ XMM}}+C{R}_{\text{ 0.5-4 keV}}^{\text{ XMM}}}$$

$${\text{HR}}_2=\frac{C{R}_{\text{ 15-30 keV}}^{\text{ NuSTAR}}-C{R}_{\text{ 3-8 keV}}^{\text{ XMM}}}{C{R}_{\text{ 15-30 keV}}^{\text{ NuSTAR}}+C{R}_{\text{ 3-8 keV}}^{\text{ XMM}}}$$

$${\dot{M}}_{out}\sim \Omega N_H{m}_p{v}_{out}{R}_{in}$$

$$R_{in}\sim 2(\alpha \frac{L}{L_{Edd}}-1)(\frac{v_{\infty }}{c}{)}^{-2}$$

Full text

1. Introduction

Outflows are commonly detected in Active Galactic Nuclei (AGN), through absorption lines visible in X-rays and in UV (Crenshaw et al., 2003), with moderate velocity of hundreds to several thousands km/s. The so-called warm absorbers are found in more than 50% of AGN (Blustin et al., 2005; Piconcelli et al., 2005; McKernan et al., 2007). They are thought to result from efficient accretion onto supermassive black holes, if the radiative energy exceeds the binding energy of the gas (King & Pounds, 2003; King, 2010). Higher velocity ($v_{out}=-0.1$-$-0.4c$) and dense outflows ($N_{H}\sim 10^{23}\text{ cm}^{-2}$) have been also observed through blueshifted iron absorption lines above 7 keV (Chartas et al., 2002; Reeves et al., 2009; Tombesi et al., 2010) in about 40% of AGN (Tombesi et al., 2010; Gofford et al., 2013), and in particular in bright and distant quasars (Chartas et al., 2003; Reeves et al., 2003; Pounds et al., 2003; Lanzuisi et al., 2012; Chartas et al., 2014; Vignali et al., 2015; Dadina et al., 2018). These are often called Ultra-Fast Outflows (UFOs). Their frequent detection suggests that they are characterized by a wide angle (e.g. Nardini et al. 2015; Reeves et al. 2018a). Their outflowing rates can reach several solar masses per year (up to $\sim 10^{3}M_{\odot}\text{yr}^{-1}$), for a kinetic power of $10^{45}-10^{46}$ erg s*-1* (Reeves et al., 2009; Tombesi et al., 2011, 2013). The link between galaxy parameters and the growth of the central supermassive black hole (SMBH), such as the black-hole mass – velocity dispersion $M-\sigma$ relation (Ferrarese & Merritt, 2000; Gebhardt et al., 2000), is thought to be regulated by AGN feedback processes, via these powerful winds. The kinetic power of UFOs, launched from the accretion disk at a few gravitational radii ($R_{g}$) from the SMBH, can be 0.5-5% of the bolometric luminosity, and affect the AGN host galaxy, by sweeping away the galaxy’s reservoir of gas and hence quench the star formation (Silk & Rees, 1998; King & Pounds, 2003; King, 2010; Di Matteo et al., 2005; Hopkins & Elvis, 2010; Alexander & Hickox, 2012; Fabian, 2012; Kormendy & Ho, 2013; Tombesi et al., 2015; Fiore et al., 2017). Several mechanisms have been suggested to play a role in the acceleration of these fast winds: UFOs can be radiatively driven winds (e.g. Proga & Kallman 2004; Sim et al. 2010; King 2010; King & Pounds 2015; Reeves et al. 2014; Hagino et al. 2015; Nomura & Ohsuga 2017), thermal winds (e.g. Begelman et al. 1983), or magneto-hydrodynamic (MHD) flows (e.g. Blandford & Payne 1982; Fukumura et al. 2010; Chakravorty et al. 2016; Fukumura et al. 2015, 2018a; Kraemer et al. 2018; Fukumura et al. 2018b).

PDS 456, which was first identified by Torres et al. (1997) in the Pico dos Dias survey (PDS), is the most luminous radio-quiet quasar in the local universe (z=0.184), with a bolometric luminosity of $L_{Bol}=10^{47}$ erg s*-1* (Reeves et al., 2000) and a black-hole mass of $M_{BH}=10^{9}M_{\odot}$ (Reeves et al., 2009). Its initial observations taken with RXTE, ASCA, BeppoSAX and then XMM-Newton (Reeves et al., 2000; Vignali et al., 2000; Reeves et al., 2002) have shown a steep absorption feature around 9 keV, as well as a rapid variability, suggesting the presence of an absorbing highly ionized outflow in the line-of-sight. The deep absorption trough detected above 7 keV has been associated with the blueshifted K-shell transition of highly ionized iron (Fe xxvi ). Broad absorption features in XMM-Newton/RGS data near 1 keV are additional signatures of a dense and highly ionized UFO with an extreme velocity of about $v_{out}=-0.1$-$-0.2c$ (Reeves et al., 2003; Behar et al., 2010).

This persistent UFO was also identified in Suzaku observations, via an absorption feature near 9 keV corresponding to a Compton-thick and clumpy wind of velocity $v_{out}=-0.25$-$-0.30c$ (Reeves et al., 2009). Rapid variability of the high velocity iron K-shell absorption lines may result from wind clumpiness (Gofford et al., 2014; Reeves et al., 2016). The spectral variability of PDS 456 can be explained by the variations in the partial covering absorber observed in most of the data (e.g. Reeves et al. 2014; Nardini et al. 2015; Matzeu et al. 2016a; Reeves et al. 2016; see also the Principal Component Analysis, PCA, performed by Parker et al. 2018 on XMM-Newton and Suzaku data). In this model, the clouds responsible for the partial covering have a size of 20 $R_{g}$ and would be the denser clumps of the inhomogeneous accretion disk wind (Matzeu et al., 2016b; Reeves et al., 2018a). Luminari et al. (2018) recently determined that the $0.23c$-velocity wind may be a wide-angle outflow, according to the covering factor of 0.7 and the opening angle of 71 deg resulting from the application of an AGN wind emission model (WINE) to XMM-Newton and NuSTAR data.

Nardini et al. (2015) analyzed five simultaneous XMM-Newton and NuSTAR observations from 2013-2014 and found in all these observations a broadened Fe K emission line in addition to an absorption trough that they first fitted with individual lines and then with a P Cygni profile characteristic of a spherically symmetric expanding gas. Their analysis revealed for the first time several blueshifted absorption lines (Fe xxvi Ly$\alpha$, Ly$\beta$ and K-edge) with the same outflow velocity $v_{out}=-0.25c$. Nardini et al. (2015) performed photoionization modeling of PDS 456, using the photoionization program XSTAR (Kallman & Bautista, 2001), constrained with its Spectral Energy Distribution (SED) based on data from UV to hard X-rays (from 2 eV to 30 keV) collected by NuSTAR and XMM-Newton EPIC-pn and OM observations from 2013-2014, approximated by a three-segment broken powerlaw (Nardini et al., 2015). The hard X-ray powerlaw has a photon index of $\Gamma=2.4$ according to XMM-Newton/EPIC and NuSTAR data. The OM photometric data are well described by a slope of $\Gamma=-0.7$, the connecting slope between 10 and 500 eV being $\Gamma=3.3$. A more complex description of this SED with a multi-temperature Comptonized accretion disk model was applied by Matzeu et al. (2016b) to Suzaku data. Analysis of the RGS data was performed by Reeves et al. (2016), who found Broad Absorption Line (BAL) profiles around 1 keV, identified as He- and H-like Neon and L-shell Iron lines blueshifted by a velocity of $v_{out}=-0.1$-$-0.2c$, which could be the signature of a lower ionization and clumpy phase of the accretion disk wind responsible for the absorption trough around 9 keV. The emission lines detected by Nardini et al. (2015), Reeves et al. (2016) and Matzeu et al. (2017a) in Suzaku data are likely associated with the re-emission from the outflow in PDS 456. Reflection on ionized material has also been considered in several studies, but the scenario of an ultra-fast outflow absorbing the hot corona emission was often preferred (Reeves et al., 2009; Behar et al., 2010; Reeves et al., 2014; Nardini et al., 2015; Chiang et al., 2017).

BAL profiles are expected in the UV spectrum of PDS 456 as signatures of the fast outflowing gas detected in X-rays. O’Brien et al. (2005) indeed detected a Ly$\alpha$ BAL in the UV spectrum of the Hubble Space Telescope (HST), blueshifted with a velocity of $v_{out}=-0.05$-$-0.08c$, that could be the signature of a decelerating cooling outflow. However, a recent analysis of the same data by Hamann et al. (2018) identified this UV BAL as a C iv line blueshifted with a velocity of $v_{out}=-0.3c$, similar to the one measured in X-ray data. Such BAL could come from dense and low-ionization clumps embedded in the X-ray UFO.

Considering the fact that the supermassive black hole in PDS 456 is accreting at about the Eddington limit, its detected UFO should be radiatively driven (Matzeu et al., 2017b; Nardini et al., 2015), supported by a correlation between the velocity of the outflow and the ionizing luminosity, with UV line driving strongly contributing to the wind acceleration (Hagino et al., 2015; Hamann et al., 2018; Reeves et al., 2018a). However, a MHD-wind model has also been successfully applied recently to XMM-Newton and NuSTAR data of PDS 456 by Fukumura et al. (2018b).

Despite the fact that PDS 456 is a radio-quiet quasar, it has been observed recently by the European VLBI Network (EVN) at 5 GHz (Yang et al., 2019). This observation revealed a radio structure made of two components, faint, diffuse and separated by about 20 pc, that could either be the signature of a recent jet, or the radio-emission of an outflow launched in the vicinity of the central supermassive black-hole. In the latest hypothesis, the radio-emission could originate from shocks produced by the interaction of the known powerful mildly-relativistic X-ray outflow and the surrounding material.

Recent observations with XMM-Newton and NuSTAR in 2017 have shown the presence of an additional faster UFO with a velocity of $v_{out}=-0.46c$, identified from a deep absorption line around 11 keV, that may be part of the multiple velocity components of the ultra-fast outflow around the accretion disk, launched from a radius of about 10$R_{g}$, and may be visible only during the lowest states (Reeves et al., 2018a, b). In addition to the two fast components of the UFO, a less ionized soft absorber with a variable covering factor was detected, probably originating from denser clumps further out along the stratified outflow (Reeves et al., 2018a). Simultaneous UV data from HST did not show significant absorption signature as the X-ray absorbers might be too ionized to be visible in this energy band.

In the present paper, we aim to analyze multiple epochs of PDS 456, namely Chandra/HETGS data (from 2003 and 2015), NuSTAR data (from 2015), and XMM-Newton and NuSTAR data from 2013-2014, to look for signatures of the persistent UFO previously detected and of the faster UFO recently identified, and constrain their physical properties. We describe the data, their selection and their reduction in section 2. We then present our dual-approach spectral analysis. We first explain the combination and binning of the data in section 3.1 and the continuum fitting in section 3.2. We next present a model-independent analysis, involving the modeling of the absorption features at high energy with gaussian lines (section 3.3) and with P Cygni profiles (section 3.4). We finally present in section 3.5 the model-dependent approach that consists in photoionization modeling, before discussing the results in section 4 and conclude in section 5.

2. Data

2.1. Data selection

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 datasets, as presented on Table 1. 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 with NuSTAR data, therefore XMM-Newton and NuSTAR data from 2013-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 performed consistently along the different datasets allowed us to constrain the winds 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 re-analyze RGS data, already studied in details by Reeves et al. (2016), as we focus on high velocity and high ionization outflows, whose signatures are expected to be outside 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.

2.2. Data reduction

PDS 456 was observed with the High Energy Transmission Grating Spectrometer (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 May 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 a full-width at half-maximum (FWHM) resolution of 0.023 Å, and the high energy gratings (HEGs), having a 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=\pm 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 (RMF) and effective area files (ARF) were generated with mkgrmf and mkgarf.

PDS 456 was simultaneously observed by NuSTAR during the observation “CN” of Chandra, from 2015 July 21 to July 23 (74 ks exposure) and from July 24 to July 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 stray light, 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 (CMEs) 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.

3. Spectral analysis

We performed our dual-approach spectral analysis of the data listed in Table 1 using the Interactive Spectral Interpretation System (ISIS; version 1.6.2-40, Houck 2002). We analyzed the simultaneous data from Chandra/HETGS and NuSTAR from 2015 (observation “CN”), the simultaneous time-averaged XMM-Newton and NuSTAR observation (combined data “XN1” to “XN5”) and the archival Chandra/HETGS observation from 2003 (observation “C”) in a similar way, adapting the procedure to the available datasets, as described in the following. After describing the combination and binning of our data in section 3.1, and the continuum modeling in section 3.2, we performed a model-independent analysis of the selected observations, looking at individual features at high energy (sections 3.3 and 3.4) and in the broadband spectra (blind line search described in appendix). We then performed a model-dependent analysis using photoionization modeling (section 3.5). This dual-method involving both model-dependent and model-independent analyses is an accurate way to check that the results found via different procedures are consistent and thus reliable. In this paper, all errors are quoted at 1$\sigma$ confidence level. All figures are presented in rest frame energy. Figure 1 shows the spectra of the three observations (coarsely rebinned for presentation purpose), and illustrates the changes in fluxes and shapes between the different epochs, as already noticed in previous studies (e.g. Matzeu et al. 2017b).

3.1. Combination and binning

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_datasets function in ISIS. For the observation “CN” (see Table 1), we put the HEG data on the same grid as MEG and combined all HETGS spectra (HEG-1, HEG+1, MEG-1, MEG+1), and we also combined NuSTAR FPMA and FPMB data of both observations 90101008002 and 90101008004. We used a combination of binning by signal-to-noise and channel binning, in order to get an adequate signal-to-noise even below 1 keV (as discussed by Danehkar et al. 2018). We have thus binned Chandra and NuSTAR data with a minimum of 1 channel per bin (binned to Half Width Half Maximum of MEG resolution) and a minimum signal-to-noise ratio of 2. We took into account the background. The same binning of 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. For the estimation of the error bars and the evaluation of the significance of our results, we used a $\Delta$C value 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 $\Delta\chi^{2}$ for the case of Gaussian statistics (see for example 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 dataset, we combined the pn spectra of the 5 observations “XN1” to “XN5”, as well as the FPMA and FPMB spectra. We refer to this combined dataset as “XN”. Despite a significant variability of the shape of the spectra 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 Fig. 2), using count rates $CR$ in different energy bands:

[TABLE]

and

[TABLE]

in order to study the possible impact of spectral changes among the observations. By their definitions, $\text{HR}_{1}$ gives information about spectral variability regarding the whole EPIC/pn energy band, while $\text{HR}_{2}$ 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, 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 including and excluding this second observation, we didn’t notice any significant change in the spectral shape. In the bottom panel of Figure 2, we see that the hardness ratio representative of the variability around the absorption features is consistent among the five “XN” observations. The hardness ratios presented in Fig. 2 show that all “XN” observations are consistent spectrally and can be combined for our analysis. The stability of $\text{HR}_{2}$ 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, and to get an average shape of the P Cygni-like profile detected in individual observations, to be compared with the “CN” observation. We binned NuSTAR data to a minimum of 2 channels per bin and a minimum signal-to-noise of 2, and XMM-Newton data to a minimum of 1 channel per bin and a minimum signal-to-noise of 5 (during the analysis described in the following sections, we checked that re-binning the data to a minimum of 5 channels per bin does not change the results on the absorption lines). The binning of these data is sufficient to get enough counts per channel to allow us to use $\chi^{2}$ statistics for the analysis.

3.2. Continuum modeling

To perform a precise analysis of the combined datasets 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 fitted the hard energy band above 3 keV with a simple powerlaw, taking into account Galactic absorption (with $N_{H}^{\text{gal}}$ fixed to $2.4\times 10^{21}\text{cm}^{-2}$ 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 to fit simultaneously data from different satellites, however the photon index of the powerlaw was the same. The fit resulted in a C statistic of 472.8 for 364 dof. In order to account for the spectral curvature as found in previous studies between 2 and 5 keV (e.g. Behar et al. 2010; Reeves et al. 2009; Turner & Miller 2009; Reeves et al. 2014, 2018a), we added a partial covering absorption to the powerlaw, using the zpcfabs model, which improved the fit significantly, with a $\Delta$C of 15.4 (for two parameters of interest) and a ftest F-value of 4.05 and 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 powerlaw with a pexmon model. However, the improvement on the fits is less significant, with a $\Delta$C of 11.1 (for three free parameters), a F-value of 2.97 but a p-value of 0.03. Furthermore, the spectra of our datasets do not show an excess around 30 keV or a fluorescence Fe K$\alpha$ line at 6.4 keV. The spectral complexity of our datasets above 3 keV is well illustrated by the data-to-model ratios of observations “CN” and “XN” on Figure 3.

When adding the soft part of the spectra and extrapolating this continuum model fitted in hard X-rays, 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 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) (as we were re-analyzing the same data), and because using a more usual black-body component instead did not improve the fit. The resulting continuum model is:

constanttbabszpcfabs*(powerlaw+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.

3.3. Fe K emission and absorption features

After establishing the continuum model, we looked for absorption features above 7 keV. We fitted the spectra above 3 keV, adding to the base continuum model Gaussian and edge models, 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$\alpha$ emission and absorption lines (at $E_{lab}=$6.97 keV), the Fe xxvi Ly$\beta$ 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 was also free to vary, independently from the ones of the absorption features. This way we were able to measure the outflow velocities of the UFO in emission and in absorption. We did this fit for both “CN” and “XN” observations, because the NuSTAR data helped to define the continuum so well. The significance of each line is given by $\Delta$C for observation “CN” and $\Delta\chi^{2}$ for observation “XN”, adding first to the continuum the emission line (with three parameters of interest), then the Ly$\alpha$ absorption line (three free parameters), then the Ly$\beta$ absorption line (linked in shift and width, so only one free parameter), and finally the edge (linked in shift, leaving one parameter of interest). Parameters of this fit are shown in Table 3. Our results show that the Fe xxvi Ly$\alpha$ emission and absorption line are significantly detected, in addition to the Fe xxvi Ly$\beta$ absorption line and the Fe xxvi K-edge (at lower significance), with a slight emission blueshift of $-0.01$-$-0.02c$ and a large absorption blueshift of $-0.27c$. Note that we assessed the statistical significance of the lines using $\Delta\chi^{2}$ and $\Delta$C, however this method only gives an approximated significance according to Protassov et al. (2002).

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. Parameters of this fit are shown in Table 4. The results show that the Fe xxvi Ly$\alpha$ emission and absorption lines for both UFOs ($v_{out1}$=-0.26c and $v_{out2}$=-0.47c) are detected significantly, the Fe xxvi Ly$\beta$ absorption line from the second UFO being detected at a lower confidence level. However, the Fe xxvi Ly$\beta$ 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 the observation “CN”, since it can only be established from NuSTAR data that have a smaller spectral resolution compared to 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$\sigma$ confidence level).

3.4. P Cygni-like profiles

The combination of the broad emission and the broad and blueshifted absorption Fe xxvi line has a P Cygni-like profile 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_{0}$ the velocity at the photosphere and $x=r/R_{0}$ the radial distance $r$ normalized to the photospheric radius $R_{0}$. The optical depth of the line is described by $\tau(w)\approx\tau_{tot}w^{\alpha_{1}}(1-w)^{\alpha_{2}}$, with $\alpha_{1}$ and $\alpha_{2}$ characterizing the sharpness of the P Cygni profile. We first applied a P Cygni model to replace the Fe xxvi Ly$\alpha$ emission and absorption lines. We fixed the parameters $\gamma$ and $w_{0}$ 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 those 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_{0}$ (corresponding to the beginning of the absorption feature), the terminal velocity $v_{\infty}$, and the parameters $\tau_{tot}$, $\alpha_{1}$ and $\alpha_{2}$. We then added another P Cygni model in order to fit the second UFO. We chose the same free parameters, however we tied the energy $E_{0}$ of the second P Cygni profile to that of the first P Cygni profile, because of the assumption of a single Fe xxvi Ly$\alpha$ emission line for both UFOs (as also supposed in section 3.3 and in the analysis from Reeves et al. (2018a, b). This fit with the second UFO was thus done with nine free parameters. We performed the fits on observations “CN” and “XN”, showing the results in Table 5.

The statistics of the fits (last line of Table 5) as well as the residuals (plotted on 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), supporting the hypothesis of the presence of a second UFO having a higher velocity than the first one in both datasets.

3.5. Photoionization modeling

After our line identification at high energy (above 3 keV), we then added to our analysis the soft parts of the data (down to 0.4 keV) 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 previous sections 3.3 and 3.4 focusing on high energy individual 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 datasets “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; Kallman & Bautista 2001; Kallman et al. 2004, 2009). However, instead of using XSTAR tabulated grids for our fits (as it has been done in previous works, e.g. Nardini et al. 2015; Matzeu et al. 2017a; Reeves et al. 2018a, b), we employed the analytic XSTAR models warmabs for reproducing the absorption lines, and photemis for generating the emission lines. Using such analytic functions requires a larger amount of computation time in comparison to the calculation of XSTAR tabulated grids. However, it allows to explore the entire range of values for the column density $N_{H}$, and the ionization parameter $\xi=L_{ion}/(nr^{2})$ (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 a SED that slightly differs from the former one only for the X-ray slope ($\Gamma=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 datasets, following the sequence described below and the order shown in the Table 6 that 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, since emission lines such as slightly blueshifted Fe xxvi Ly$\alpha$ 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 ($\Delta$C and $\Delta\chi^{2}$ in Table 6), the wind emission is significantly required by “XN” data, and marginally needed by “CN” and “C” data. A low-ionization warm absorber at the systematic Doppler shift is also required by “XN” data only, as shown by the significant improvement of the fit when including this new component in the previous model (see $\Delta\chi^{2}$ 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 only), 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 partial 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 $\Delta$C=21.4 for “CN”, $\Delta\chi^{2}$=406.0 for “XN”, and $\Delta$C=19.4 for “C”, for one degree of freedom, demonstrating that the slowest UFO is significantly required in all our datasets at a confidence level larger 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 (so 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 larger velocity UFO in the model improves the fit for the observations “CN” and “XN” (see Table 6), with $\Delta$C=8.0 for “CN” and $\Delta\chi^{2}$=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.

For these fits, we allowed the continuum parameters (powerlaw, soft-excess and cross-calibration factor) to vary, since the application of 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 20000 km/s for all components in all observations, so we fixed this parameter to this value. Such a large 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 have been fixed to solar values.

The fits have been performed in ISIS using Cash statistics and the powell method for observations “CN” and “C”, and $\chi^{2}$ statistics and the mpfit method for observation “XN”. Using the XSTARDB111http://space.mit.edu/cxc/analysis/xstardb/index.html library of ISIS S-lang scripts, we listed the strongest features predicted by our photoionization models with the fitted parameters given in Table 6, for the three observations and for all components. For the UFOs, in addition to the expected Fe xxv He$\alpha$, He$\beta$, Fe xxvi Ly$\alpha$ and Ly$\beta$ absorption lines, the XSTAR model predicts Si xiv and S xvi Ly$\alpha$ 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$\alpha$ and Fe xxvi Ly$\alpha$, but also from other ions (C vi , O viii , Ne x , Mg xii , Fe xxiv , Si xiv and S xvi ). The lower-ionization partial 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 get the probability distribution of each parameter via Bayesian data analysis. For this purpose, we used the isis_emcee_hammer function222https://www.sternwarte.uni-erlangen.de/wiki/index.php/Emcee 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.

4. Discussion

4.1. Modeling the highly-ionized features

We initially modeled the absorption features at high energy (above 3 keV) with blueshifted and highly-ionized iron lines to characterize the ultra fast outflows in PDS 456. We found that the Fe xxvi Ly$\alpha$ emission and absorption lines, the Fe xxvi Ly$\beta$ 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 fitted to data whose resolution is not sufficient. Indeed, signatures of the second UFO are located in the energy band covered predominantly by NuSTAR. However, the second UFO is required, more convincingly in the observation “XN” than in the 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 one significantly improved the fits for both observation “CN” and “XN” (see the statistics in Table 5 and the residuals on Figure 5). Table 5 shows that the velocities implied by the P Cygni model are slightly higher (of about 15-20$\%$) than the one 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, that 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 one of the second P Cygni profile being smaller. Parameters $\alpha_{1}$ and $\alpha_{2}$ 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 consistent also with previous results (e.g. Nardini et al. 2015; Matzeu et al. 2017a; Reeves et al. 2018b).

4.2. Photoionization results

We applied self-consistent XSTAR photoionization models to our three different epochs of observations (using the entire detectors 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 ($\Delta$C and $\Delta\chi^{2}$ 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->99.9% confidence level for the fastest UFO). The high values of column densities $N_{H}$ and ionization parameters $\xi$ obtained from our modeling are typical of ultra-fast outflows (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, that helped to constrain the other parameters better, as the model is complex with a lot of free parameters. The turbulent velocity of 20000 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 has been 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 one of the partial covering absorber (even if the partial 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, b). 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 partial covering absorber. This scenario was adopted to obtain a better constraint. However, we saw in appendix that when we identified lines detected by blind search, the transitions had a slightly different blueshift for lines from the partial covering absorber than from the slowest UFO, in particular in the observation “CN”. It may imply that the actual velocities of both components are slightly different. The large UFO velocities derived from photoionization modeling (about $v_{out}=-0.24$ - $-0.29c$ for the slowest UFO, and $v_{out}=-0.48c$ 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 Appendix). An additional method for the characterization of the kinematic components of the UFOs is also described and applied to PDS 456 data in Appendix, and the results are also consistent with values found through photoionization modeling.

MCMC approach gave us the probability distribution of each parameter of the photoionization model, giving contours shown on 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 partial covering absorber along the three observations. We find that the velocities are consistent between the observations “CN” and “XN”, with a better constraint for the observation “XN” compared to the observation “CN” because of the better signal-to-noise. Despite the poor constraint on the velocity of the UFO detected in the 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.2c$ and contours coherent with an outflowing velocity of $v_{out}=-0.2$ - $-0.3c$. Contours of the UFO in the observation “C” are obviously improved when NuSTAR data are combined with Chandra/HETGS data as in the observation “CN”. However, despite the lack of information at high energy in the observation “C”, the photoionization modeling allowed us to constrain the parameters of the winds at a lower significance level. The contours of the all parameters, and in particular parameters of the fastest UFO, are more precise for the observation “XN” compared to the observation “CN” (regarding velocities, ionization parameters, column densities and covering factors), as expected because the signal-to-noise in the NuSTAR energy band brought by combining all the observations “XN1” to “XN5” is higher than that from the observation “CN”. The ionization parameters and column densities of the UFOs are relatively stable between the different epochs. The ionization of the partial 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(\xi)\sim 3$) consistent with previous studies (Nardini et al., 2015; Reeves et al., 2018a). Variability of the covering factor ($c_{f}\sim 0.3-0.8$) and of the column density of this absorber ($N_{H}=3.3-16.6\times 10^{22}\text{ cm}^{-2}$), 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 example of PG 1211+143, Gallo & Fabian 2013, but counter-argument from Lobban et al. 2016). Nardini et al. (2015) found that such a reflection-dominated scenario under-predicts the strength of the absorption feature around 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 partial 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”. As 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 compared to the continuum model. However, adding this pexmon reflection component to the photoionization model with the two UFOs did not improve the fits ($\Delta$C=1.66 for observation “CN”, and $\Delta\chi^{2}$=-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.

4.3. Thermal stability

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(\xi/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 of 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 remind here that these two SEDs are identical in the UV band, and slightly differ in X-rays in order to match our three datasets. 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 of 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}=<E>/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.

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), on a 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 the observations “XN”, “CN” and “C” are on a negative gradient branch (red points), we might expect that such gas should be 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 on Figure 7, at high ionization the cooling dominates above the curve and the heating dominates below the curve. So on this high ionization branch, a small perturbation increasing the temperature will bring the gas to a region where the cooling dominates, and a small decrease of 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 datasets, 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 the difference between solid and dashed lines (see legend of the figure for more details). We overplotted 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, “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 datasets, and can coexist with S and Si ions in a single ionization zone corresponding to the UFOs (red zone). This co-existence 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 partial covering absorbers. This uncertain origin is induced by the fact that, to simplify the model applied to our data, we linked together the velocities of the partial covering absorber and of the slowest UFO, 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 co-existence 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.

4.4. Compare results from different epochs

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_{out1}=-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_{out1,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, b). Indeed, after finding a second, faster absorber in recent XMM-Newton and NuSTAR data from 2017, the authors re-analyzed 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 wanted to get optimum signal-to-noise using all available data. We performed a simple check of the detection of the fastest UFO in the individual “XN” observations, by applying the fitting procedure described in Section 3.3 on these five datasets, using two sets of Fe lines to look for signatures of both UFOs. Regarding the improvement of the fits statistics when considering two UFOs instead of one, and according to the significant detection of the Fe xxvi Ly$\alpha$, Ly$\beta$ and K-edge features from the fastest UFO, we found that the second UFO is required by the data “XN5” ($\Delta\chi^{2}/\Delta\text{dof}$=19.03/5, >99.9% confidence level), “XN2” ($\Delta\chi^{2}/\Delta\text{dof}$=15.82/5, about 99.5% confidence level) and “XN4” ($\Delta\chi^{2}/\Delta\text{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$\alpha$ from the fastest UFO improving the fit by $\Delta\chi^{2}$=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, b) 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 anti-correlation between the statistical confidence level for the detection of the fastest UFO and the flux intensity of PDS 456. Such a 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 smaller than the one of the slowest UFO (see Table 6), consistent with results from Reeves et al. (2018b).

Figure 9 shows the outflow velocities (top figure) and the ionization parameters (bottom figure) 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 twelve 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.3c 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 supermassive black hole 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 the observation “C” compared to the parametric constraints derived from additional high energy data available in the observations “CN” and “XN”.

4.5. Winds mass outflow rates and energetics

The mass outflow rate is:

[TABLE]

where $\Omega$ is the solid angle, $N_{H}$ is the column density, $m_{p}$ is the mass of the proton, $v_{out}$ is the outflowing velocity, and $R_{in}$ is the starting point of the wind (see details for example in Nardini et al. 2015). We used an approximated value of $\Omega=2\pi$ for the solid angle, as justified by Nardini et al. (2015) (see their discussion in 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:

[TABLE]

where $v_{\infty}$ is the wind terminal velocity and $\alpha$ 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 $\alpha=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 $\sim 5\times 10^{15}$ cm (or 30$R_{g}$ for $M_{BH}=10^{9}M_{\odot}$) for the slowest UFO, and $\sim 1\times 10^{15}$ cm (or 9$R_{g}$) for the fastest UFO. These values are close to the escape radii $R_{escape}=\frac{2GM}{v_{out}^{2}}$, i.e. the minimum radii from which winds of a given outflow velocity can be launched. We are probably observing a stratified ultra-fast outflow in PDS 456, with several components from multiple stratification layers having different velocities and being launched from the accretion flow close to the supermassive black hole (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 small 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 getting too highly ionized (Matzeu et al., 2016a; Reeves et al., 2018b).

The maximum radial distance of the absorbers can be estimated by considering that $\Delta R/R<1$, i.e. their thickness cannot exceed their distance from the ionizing source (e.g. Reeves et al. 2003, 2018b). Using the definition of the ionization parameter given in section 3.5, $R_{max}<L_{ion}/N_{H}\xi$. This gives maximum radii of $7\times 10^{17}$ cm (or 4700$R_{g}$) for observation “CN”, $2\times 10^{16}$ cm (or 130$R_{g}$) for observation “XN”, and $6\times 10^{18}$ cm (or 40000$R_{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.2M_{\odot}/yr=0.05\dot{M}_{Edd}$ for observation “CN”, $4.4M_{\odot}/yr=0.2\dot{M}_{Edd}$ for observation “XN”, and $0.6M_{\odot}/yr=0.03\dot{M}_{Edd}$ for observation “C”. For the fastest UFO, we found $0.4M_{\odot}/yr=0.02\dot{M}_{Edd}$ for observation “CN”, and $1.6M_{\odot}/yr=0.07\dot{M}_{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 for example estimated kinetic power and mass outflow rate of $\sim$ 15% and $\sim$ 50% of Eddington values respectively for Nardini et al. 2015, 5% and 40% for Gofford et al. 2014) because we obtained slightly different parameters resulting from the photoionization fit and because we used a smaller 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 hundreds of gravitational radii (Reeves et al., 2009; Nardini et al., 2015; Matzeu et al., 2016a, b). We see that the mass outflow rate of the fastest UFO is smaller than the one from the slowest UFO. We calculated the kinetic power of the winds ($P_{kin}=0.5\dot{M}v_{out}^{2}$) and found kinetic powers of $0.02L_{Edd}$ for observation “CN”, $0.07-0.08L_{Edd}$ for observation “XN” and $0.008L_{Edd}$ for observation “C” (again a bit smaller than in previous works, for the same reason as above). Both winds are found to have a similar kinetic energy that contradicts the prediction of a larger power for the faster UFO, estimated by Reeves et al. (2018b) using the approximation that $P_{kin}\propto v_{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, that is sufficient to induce significant AGN feedback according to models of black hole and host galaxy co-evolution (King & Pounds, 2003; Di Matteo et al., 2005; Hopkins & Elvis, 2010).

5. Conclusion

We presented the analysis of simultaneous Chandra/HETGS and NuSTAR observations of PDS 456 from 2015 (“CN”), simultaneous XMM-Newton and NuSTAR data from 2013-2014 (“XN”), and Chandra/HETGS data from 2003 (“C”). We performed a dual-approach study of these selected observations of the quasar, analyzing data from the three different epochs in a consistent way, using both model-independent and model-dependent techniques.

We confirmed the presence of the persistent ultra-fast outflow at velocity of $v_{out}=-0.24$-$-0.29c$, that was observed in previous studies (e.g. Reeves et al. 2009; Nardini et al. 2015). We also detected a faster UFO ($v_{out}=-0.48c$) in the “CN” and “XN” observations, that was reported previously only in very recent observations (Reeves et al., 2018a, b). In the model-independent approach, we observed their signatures via deep absorption troughs at about 9 and 11 keV, corresponding to blueshifted, highly ionized iron K-shell transitions, that form P Cygni profiles when considering the associated blueshifted emission. We also identified other lines in the HETGS spectra, blueshifted at the same extreme velocities at lower energies, e.g. Si xiv Ly$\alpha$, S xvi Ly$\alpha$, O viii Ly$\alpha$, Ne x Ly$\alpha$, and possible contribution from nickel at high energy.

In the model-dependent approach, we performed photoionization modeling to characterize both UFOs ($log(\xi)\sim 6-7\text{ erg cm s}^{-1}$, $N_{H}\sim 1-8\times 10^{23}\text{ cm}^{-2}$) as well as the partial covering absorber ($log(\xi)\sim 3\text{ erg cm s}^{-1}$, $N_{H}\sim 3\times 10^{22}-2\times 10^{23}\text{ cm}^{-2}$, $c_{f}\sim 0.3-0.8$). We found that all the winds detected in the three datasets are thermally stable and can coexist.

The outflow of PDS 456 is probably composed of several components from multiple layers having different velocities and ionizations, launched from the accretion flow close to the supermassive black hole, and certainly radiatively driven. Both relativistic components of the outflow are powerful enough to play a role in the evolution of the host galaxy, with mass outflow rates of 2-20% and kinetic powers of 0.8-8% of the Eddington values.

We performed an analysis using different methods that led to consistent results. However, we made some assumptions in order to be able to constrain the winds characteristics, because of the complexity of the models applied to our data. Further simultaneous and high signal data are required in order to test our assumptions. Future high resolution instruments such as ARCUS and Athena will be useful to determine more precisely the structure of the high velocity winds in PDS 456.

We thank very much the anonymous referee for useful comments and corrections which helped to improve this paper. We gratefully acknowledge Claude Canizares for the Chandra/HETGS GTO time to observe PDS 456. We sincerely thank Tim Kallman for our private communication about stability curves for photoionized gas. RBM acknowledge Emanuele Nardini and Chris Done for providing the P Cygni profile model used in this paper. We thank Fiona Harrison, PI of NuSTAR, for DDT time. Support for this work was provided in part by the National Aeronautics and Space Administration (NASA) through the Smithsonian Astrophysical Observatory (SAO) contract SV3-73016 to MIT for support of the Chandra X-Ray Center (CXC), which is operated by SAO for and on behalf of NASA under contract NAS8-03060.

Blind line detection

We performed a blind line search in the broad energy band for the observations “CN” and “C”, and the hard energy band for the observation “XN”. To perform this blind line search, we added 50 gaussian lines on top of our continuum (described in section 3.2) that we allowed to be either in emission or in absorption. We used gabs models in order to take into account eventual saturation of the lines. We first fixed the continuum, added one gaussian model (with free energy, width and strength, i.e. three parameters of interest), fitted the data, varied the continuum parameters and fitted again the data, and then fixed everything before running again this sequence for 50 iterations. We kept only statistically significant lines for which we got a $\Delta$C or $\Delta\chi^{2}$ larger than 6.25, corresponding to 90% confidence level for three degrees of freedom. Note that these estimations of the significance of the lines might be slightly overestimated according to Protassov et al. (2002). The results of this blind line search and line identification are presented in Table 7 and Figure 11 for observation “CN”, Table 8 and Figure 12 for observation “XN”, and Table 9 and Figure 13 for observation “C”. The tables give the rest energy of each line, together with its width, equivalent width and significance. We tentatively identified some of the detected emission and absorption lines (reported in the tables), considering the transitions of the strongest lines resulting from the photoionization modeling described in section 3.5. We also considered weaker contributions from higher Z elements, such as the H-like transition of Nickel. Even if Ni is 20 times less abundant than Fe (Grevesse & Sauval, 1998), and thus negligible in XSTAR photoionization modeling, such a transition could possibly exist in gas with such high column densities and high ionizations as found in Table 6.

For the observation “CN”, some lines detected blindly have been tentatively identified as blueshifted lines with three different ranges of values (see Figure 11), as shown by the $z_{out}$ column in Table 7. Taking into account errors on the energy of the line and its width, and taking also into account that the actual shape of the line may be different from a Gaussian curve (like for radiative recombination continuum), we found a velocity in emission $v_{em}$ that is consistent with the value and error bars found when fitting the data with XSTAR models, i.e. $v_{em}=-0.098$ - $-0.050c$. This is also the case for the velocity in absorption of the second wind $v_{abs2}=-0.512$ - $-0.436c$, with absorption lines from Fe xxv He$\alpha$, Fe xxvi Ly$\alpha$, Fe xxvi Ly$\beta$ and Ni xxviii Ly$\alpha$, and the lower significance Si xiv Ly$\alpha$ line. The velocity of the first wind determined using photoionization models, $v_{abs1}=-0.319$ - $-0.264c$, is consistent with the identification of absorption lines Fe xxv He$\alpha$, Fe xxvi Ly$\alpha$, S xvi Ly$\alpha$, Ne x Ly$\alpha$ and Ni xxviii Ly$\alpha$. However, the lines tentatively identified as Fe xxvi Ly$\beta$, Fe xvii 2p-3d, O viii Ly$\alpha$ (detected only at lower significance), O viii Ly$\beta$ and Ne ix He$\alpha$ show a slightly higher velocity of $-0.355$ - $-0.332c$. This is consistent with the fact that the Fe xxvi Ly$\beta$ absorption line of the first wind was not detected significantly when fitting the data with sets of Fe K lines (see section 3.3). However, this Fe xxvi Ly$\beta$ transition may contribute to the broad absorption feature around 11 keV. The O viii and Ne ix He$\alpha$ absorption lines modeled by our broadband fitting with warmabs come mostly from the partial covering absorber. In our model, we chose to tie the velocity of this partial covering component to the one of the slowest UFO, consistently with Nardini et al. (2015), to improve its constraint when fitting the data. Furthermore, we also linked the turbulent velocities for the same reason. The photoionization modeling thus tended to make these O viii and Ne ix lines broader than when detected during the blind line search, it can thus explain the slight difference of blueshift for these lines, and it suggests that the partial covering component may have a slightly different velocity than the slowest UFO.

For the observation “XN”, some lines have been tentatively identified in the spectrum above 5 keV, consistent with the identification for the observation “CN” (see Table 8 and Figure 12). Due to the large widths of the detected lines, several transitions can be attributed to the same line, as they contribute to the absorption feature. Considering the uncertainties, the blueshift values resulting from this identification are close to the values resulting from the photoionization modeling (see Table 6). The slight difference could be due to the fact that the blind line search has been done only above 5keV, while the photoionization model takes the broadband spectrum into account, determination of the velocities are thus more precise in the later case.

For the observation “C”, only some emission lines and absorption lines from one UFO have been identified. The blueshifted Fe xxv He$\alpha$ and Fe xxvi Ly$\alpha$ absorption features detected in the other observations were not detected but some other lines (Si xiv , S xvi ) coming from the slowest UFO (and from the partial covering wind that is set to have the same velocity in our photoionization model) have been identified, as shown in Table 9 and Figure 13. Some emission lines have also been tentatively identified. The velocities derived from these identifications are globally consistent with the velocity values resulting from photoionization modeling (see Table 6).

Additional kinematic diagnostic

An additional diagnostic for identifying unknown kinematics of the absorbers, as proposed by Danehkar et al. (2018) in PG 1211+143, is to take the 1st-order Chandra/HETGS counts, and use the lines tentatively identified in Tables 7 and 9 as a set of reference wavelengths, to transform the grid repeatedly to velocities, resulting in Figure 10 for the observations “CN” and “C”. We can see that velocity values resulting from photoionization modeling (orange zone: emission, red zone: absorption from the slowest UFO, green zone: absorption from the fastest UFO for the observation “CN” only) are consistent (considering uncertainties) with the peaks of emission and absorption at certain velocities (black vertical lines) for both observations “CN” and “C”. This analysis provides an additional confirmation of the kinematic components of the UFO detected in PDS 456.