The response of relativistic outflowing gas to the inner accretion disk of a black hole

TL;DR

This paper reports the detection of a relativistic ultra-fast outflow in an active galactic nucleus, showing rapid variability and a strong link between the outflow and inner accretion disk emission.

Contribution

It presents the first observation of a highly variable, relativistic outflow responding to inner disk emission, revealing a direct connection between accretion processes and outflows in AGNs.

Findings

Detected ultra-fast outflow at 71,000 km/s in AGN IRAS 13224-3809.

Observed wind variability on timescales of hours, much faster than typical quasar winds.

Found the outflow's absorption features anti-correlated with inner disk emission.

Abstract

Active galactic nucleus (AGN) feedback is the process by which supermassive black holes in the centres of galaxies may moderate the growth of their hosts. Gas outflows from supermassive black holes release huge quantities of energy into the interstellar medium, clearing the surrounding gas. The most extreme of these, the ultra-fast outflows (UFOs), are the subset of X-ray detected outflows with velocities higher than 10,000 km/s, believed to originate in relativistic disc winds, a few hundred gravitational radii from the black hole. The absorption features produced by these outflows are variable, but no clear link has been found between the behaviour of the X-ray continuum and the energy or equivalent width of the outflow features due to the long time-scales of quasar variability. Here, we present the detection of multiple absorption lines from an extreme ultra-fast gas flow in the X-ray spectrum of the active galactic nucleus IRAS 13224-3809, at 0.236+/-0.006 times the speed of light (71,000 km/s), where the absorption is strongly anti-correlated with the emission from the inner regions of the accretion disk. If the gas flow is identified as a genuine outflow then it is in the fastest 5 per cent of such winds, and its variability is hundreds of times faster than in other variable winds, allowing us to observe in hours what would take months in a quasar. We find signatures of the wind simultaneously in both low and high energy detectors, which are consistent with a single ionized outflow, linking the two phenomena. The detection of the wind responding to the emission from the inner disk demonstrates a connection between accretion processes occurring on very different scales, with the X-rays from within a few gravitational radii of the black hole ionizing the relativistically outflowing gas as the flux rises.

Tables

Table 1. Table 1: Best fit parameters excluding the absorption feature (model 1) and including it (model 2). Errors are 1 ​ σ 1 𝜎 1\sigma .

Parameter	Model 1	Model 2	Description/unit
$E_{\mathrm{Gauss}}$	-	$8.14\pm 0.03$	Line energy (keV)
${\sigma }_{\mathrm{Gauss}}$	-	${0.1}^{\ast }$	Line width (keV)
$N_{\mathrm{Gauss}}$	-	$0.15\pm 0.03$	Line strength
$q_{\mathrm{in}}$	$>8.5$	${8.3}_{-1.2}^{+0.2}$	Inner emissivity index
$q_{\mathrm{out}}$	$2.4\pm 0.1$	$2.46\pm 0.04$	Outer emissivity index
$r_{\mathrm{break}}$	${2.52}_{-0.07}^{+0.1}$	$2.59\pm 0.04$	Emissivity break radius $\left(r_G\right)$
$a$	$>0.994$	$0.989\pm 0.001$	Black hole spin
$i$	$60\pm 0.1$	$59\pm 1$	Inclination (degrees)
$\Gamma$	${1.99}_{-0.06}^{+0.12}$	${2.06}_{-0.09}^{+0.06}$	Power-law index
$E_{\mathrm{cut}}$	${50}_{-15}^{+17}$	$50\pm 6$	Power-law cutoff energy (keV)
$\log \xi$	${3.25}_{-0.09}^{+0.06}$	$3.34\pm 0.01$	Ionization (erg cm s-1)
$A_{\mathrm{Fe}}$	${4.4}_{-0.5}^{+0.4}$	$3.5\pm 0.2$	Iron abundance (solar)
$R$	$>9.3$	$>9.8$	Reflection fraction
${\chi }^2/$dof	314/285	288/283	Fit statistic

Table 2. Table 2: Best fit physical model parameters for the stacked EPIC-pn spectrum.

	Fe xxv	Fe xxvi
$v$ ($c$)	$0.244\pm 0.09$	$0.210\pm 0.009$
$\log \xi$ (erg cm s-1)	$4.0\pm 0.1$	$>4.5$
$N_H$ ${10}^{22}$cm-2	$8.2\pm 2.5$	$180\pm 50$
${\sigma }_{v1D}$ (km s${}_{}{}^{-}1$)	$1500\pm 1000$	$1500\pm 1000$

Equations

$$\dot{M}=\Omega N_{\mathrm{H}}{m}_{\mathrm{p}}v{R}_{\mathrm{wind}}$$

$$\begin{array}{rl}P& =\frac{1}{2}\dot{M}{v}^2\\ P& =0.006\Omega L_{\mathrm{Edd}}\end{array}$$

Full text

The response of relativistic outflowing gas to the inner accretion disk of a black hole

Michael L. Parker

Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK

Ciro Pinto

Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK

Andrew C. Fabian

Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK

Anne Lohfink

Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK

Douglas J. K. Buisson

Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK

William Alston

Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK

Erin Kara

Department of Astronomy, University of Maryland, College Park, MD 20742-2421, USA

Edward M. Cackett

Department of Physics and Astronomy, Wayne State University, Detroit, MI 48201, USA

Chia-Ying Chiang

Department of Physics and Astronomy, Wayne State University, Detroit, MI 48201, USA

Thomas Dauser

Remeis Observatory and ECAP, Universitat Erlangen-Nurnberg, Sternwartstr. 7, D-96049, Bamberg, Germany

Barbara De Marco

Max-Planck-Institut fur extraterrestrische Physik, Giessenbachstrasse, D-85748 Garching, Germany

Luigi C. Gallo

Department of Astronomy and Physics, Saint Mary’s University, 923 Robie Street, Halifax, Nova Scotia B3H 3C3, Canada

Javier Garcia

Space Radiation Laboratory, California Institute of Technology, 1200 E California Blvd, MC 249-17, Pasadena, CA 91125, USA

Fiona A. Harrison

Space Radiation Laboratory, California Institute of Technology, 1200 E California Blvd, MC 249-17, Pasadena, CA 91125, USA

Ashley L. King

Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA 94305, USA

Matthew J. Middleton

School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK

Jon M. Miller

Department of Astronomy, University of Michigan, 1085 South University Avenue, Ann Arbor, MI 48109, USA

Giovanni Miniutti

Centro de Astrobiologia (CSIC-INTA), Dep. de Astrofísica, ESAC, PO Box 78, 28691 Villanueva de la Cañada, Madrid, Spain

Christopher S. Reynolds

Department of Astronomy, University of Maryland, College Park, MD 20742-2421, USA

Phil Uttley

Astronmical Institute ‘Anton Pannekoek’, University of Amsterdam, Postbus 94249, NL-1090 GE Amsterdam, the Netherlands

Ranjan Vasudevan

Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK

Dominic J. Walton

Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK

Daniel R. Wilkins

Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA 94305, USA

Abderahmen Zoghbi

Department of Astronomy, University of Michigan, 1085 South University Avenue, Ann Arbor, MI 48109, USA

Active galactic nucleus (AGN) feedback is the process by which supermassive black holes in the centres of galaxies may moderate the growth of their hosts[1]. Gas outflows from supermassive black holes release huge quantities of energy into the interstellar medium[2], clearing the surrounding gas. The most extreme of these, the ultra-fast outflows (UFOs), are the subset of X-ray detected outflows with velocities higher than 10,000 km s-1, believed to originate in relativistic disc winds, a few hundred gravitational radii from the black hole[3]. The absorption features produced by these outflows are variable[6], but no clear link has been found between the behaviour of the X-ray continuum and the energy or equivalent width of the outflow features due to the long time-scales of quasar variability. Here, we present the detection of multiple absorption lines from an extreme ultra-fast gas flow in the X-ray spectrum of the active galactic nucleus IRAS 13224–3809, at $0.236\pm 0.006$ times the speed of light (71,000 km s*-1*), where the absorption is strongly anti-correlated with the emission from the inner regions of the accretion disk. If the gas flow is identified as a genuine outflow then it is in the fastest 5% of such winds, and its variability is hundreds of times faster than in other variable winds, allowing us to observe in hours what would take months in a quasar. We find signatures of the wind simultaneously in both low and high energy detectors, which are consistent with a single ionized outflow, linking the two phenomena. The detection of the wind responding to the emission from the inner disk demonstrates a connection between accretion processes occurring on very different scales, with the X-rays from within a few gravitational radii of the black hole ionizing the relativistically outflowing gas as the flux rises. **

IRAS 13224–3809 is a low-redshift ($z=0.0658$) NLS1. These galaxies are characterised by having low mass, high accretion rate black holes[7], some of which have evidence for relativistic outflows [8]. IRAS 13224–3809 is the most variable AGN in the X-ray band[9], changing by up to two orders of magnitude in flux in only a few hours. This source was the target of a recent 1.5 Ms continuous observing campaign with the X-ray Multi-Mirror Mission (XMM-Newton[10]) and 500 ks with NuSTAR[11] intended to study the rapid variability, from minutes to weeks. IRAS 13224–3809 is well known as a source with strong relativistic reflection from the inner disk [12], where fluorescent emission, produced by X-ray illumination of the inner accretion disk, is blurred and skewed by the strong relativistic effects close to the event horizon. This produces a characteristic broad iron line at 6–7 keV and a Compton scattered ‘hump’ at high energies. In this case, reflection spectroscopy revealed near maximal spin, and a spectrum dominated by reflection, implying a high degree of light-bending close to the event horizon. The iron K and L lines can also be modelled as a relavitistic absorption edge [13, 14], however the measurement of X-ray reverberation, which measures the delay between the continuum and reflected emission[15], demonstrated that the line is produced by reflection and that we have an unobstructed view of the inner accretion disk.

The source is spectrally very soft, hence the very high energy band ($>8$ keV) has not been well studied previously, due to a shortage of photons in any given observation. By combining all the data from the 2016 observing campaign into a single stacked spectrum and constraining the continuum flux using NuSTAR, we are able to reveal the presence of a structured absorption feature at $\sim 8.6$ keV in the rest frame of the source, significant at $>4\sigma$ (Fig. 1, top), along with a possible secondary line at $\sim 9.2$ keV. The 8.6 keV line is also present, at $>3\sigma$, in the stacked spectrum of the archival data (Fig. 3), meaning that it must have been present for at least 5 years. Including the absorption line does not affect the results of reflection fitting, the parameters returned are entirely consistent with those found by previous authors[12, 16] (near maximal spin, inclination of $\sim 60$ degrees, and a steep emissivity profile).

Additional evidence for the presence of a UFO is given by the simultaneous detection of three more blueshifted absorption lines in the high spectral resolution Reflection Grating Spectrometer (RGS[23]) data, which is capable of resolving narrow features at low energies. These features have a combined significance of $>5\sigma$, and are found at the same blueshift as the 8.6 keV line. Three broad features are detected at $\sim$9.5, 13.0, and 16 Å in the low-flux spectrum. They have less significance or are undetected in the higher flux spectra with possible evidence of slight velocity change. Their wavelengths match with the strongest lines predicted by the photoionized UFO model in the RGS energy band: 10.0 $\pm$ 0.5 Å (Ne X + Fe XVIII-XXII blend), 13.2 $\pm$ 0.5 Å (O VIII K$\beta$ + Fe XVIII) and 15.8 $\pm$ 0.5 (O VIII K$\alpha$).

Depending on the ionization of the source, the 8.6 keV feature could be due to either Fe XXV or Fe XXVI, with corresponding outflow velocities of 0.24$c$ or $0.21c$. However, a joint fit to the EPIC-pn and RGS data with a single absorber gives a best fit velocity of $0.236\pm 0.006c$, and robustly prefers the faster Fe xxv solution. The typical outflow velocity of radio quiet AGN is $\sim 0.1c$ [4], this line identification puts this outflow in the fastest $\sim$5% of sources. The power output is also substantial. Assuming a solid angle $\Omega=2\pi$, the kinetic power of the wind is $\sim 4$% of the Eddington luminosity (see methods), comparable to the 15% found in the quasar PDS 456[3], which is 2–3 orders of magnitude higher in mass ($\sim 10^{9}M_{\odot}$, compared to $\sim 6\times 10^{6}M_{\odot}$ for IRAS 13224-3809). This indicates that the wind in IRAS 13224-3809 may be capable of driving feedback. Similarly, the mass outflow rate is 0.4 $\dot{M}_{\mathrm{Edd}}$, implying that a large fraction of the mass is lost in the wind and that the accretion rate outside the wind launching radius may be super-Eddington. This is very similar to the wind found in the ultraluminous X-ray source NGC 1313 X-1[17, 18], which also responds to changes in the continuum[19], implying that super-Eddington relativistic winds may be a common phenomenon in high accretion rate black holes[20, 21, 22].

IRAS 13224-3809 varies dramatically in flux on short time scales, giving us a unique opportunity to probe the fast variability of the wind. We split the large dataset up into three different flux levels such that the total number of counts in each of the three spectra is equal (Fig. 1). In the lowest flux state, the line is extremely significant ($7\sigma$), with an equivalent width of 0.24 keV. However, when the flux is higher the line disappears completely, and is not required in the medium and high flux spectra. This variability of the 8.6 keV line strength must be occurring on relatively short timescales ($<10$ks), as the source variability is so extreme that any delay in the absorption responding to the continuum flux would remove the observed correlation between the two. The average time between the end of a low flux interval and the start of a high flux interval is 4–5 ks, so the line must be responding within this time. Spectrally, the changes in the equivalent width can be explained by changes in either the column density or the ionization of the absorber. However, a causal relationship between the flux of the source and the ionization is expected, and both the flux and the model ionization increase by a factor of $\sim 4$ between the low and high flux spectra, consistent with the ionization interpretation. Because of the strong relativistic reflection observed in this source, we can be confident that we are seeing the inner regions of the AGN, and the gas is responding directly to that emission.

The strength of the RGS absorption lines also anti-correlates with the flux, and therefore strongly argues in favour of the EPIC and RGS absorbers being part of the same extreme wind. In Fig. 2 we show the significance of the spectral features obtained adopting the (RGS fitted) spline continuum. We computed the confidence level of the three main absorption lines in the low-flux spectrum, where they are significantly detected, as with the EPIC-pn. We obtain 2.1$\sigma$, 2.9$\sigma$, and 3.4$\sigma$ for the 9.5 Å, 13 Å, and 16 Å absorption lines, respectively, which, since they have the same velocity shift and broadening, gives a cumulative 5.1$\sigma$ detection.

We investigate this in more detail by further subdividing the data into 10 flux-resolved spectra, again with the same number of counts in each. Fitting these data with the same model and allowing the normalization of the Gaussian line to vary between them with the velocity fixed gives a very strong correlation between the equivalent width of the line and the 3–10 keV flux (Fig. 3). The significance of this correlation is $>99.99$%. This rapid variability clearly demonstrates, for the first time, that the outflow is responding to the source flux. Variability of outflows has been reported before[6], but has never shown a clear correlation with the continuum behaviour. The correlation we find can potentially be explained as the absorption lines disappearing due to the increased ionization when the flux is high. Alternatively, this could be a geometric effect, where the line of sight from the X-ray corona intersects the wind at low source heights and misses it when the source is high or vertically extended, or the line could be produced in an ionized layer on the surface of the inner disk, absorbing the reflected emission only[24]. Since the reflection fraction scales inversely with flux [16], this could produce the observed correlation.

The discovery of a rapidly variable outflow, analogous to a quasar wind, in a low-mass, low-redshift AGN gives us a new window onto AGN feedback. The simultaneous detection of a relativistic outflow in the RGS and EPIC-pn spectra means that we can accurately characterise the physical properties of the wind, demonstrating that both sets of features can be produced by the same physical wind. The rapid variability, in a source where we can see clearly into the inner accretion disk, allows us to measure the wind responding to the X-rays originating a few gravitational radii away from the black hole, linking the physical processes occurring on these different scales. This discovery may allow us to uncover how relativistic winds are launched and accelerated.

Acknowledgements

M.L.P., C.P., A.C.F., and A.L. acknowledge support from the European Research Council through Advanced Grant on Feedback 340492. W.N.A. and G.M. acknowledge support from the European Union Seventh Framework Programme (FP7/2013-2017) under grant agreement n.312789, StrongGravity. D.J.K.B. acknowledges support from the Science and Technology Facilities Council (STFC). This work is based on observations with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. D.R.W. is supported by NASA through Einstein Postdoctoral Fellowship grant number PF6-170160, awarded by the Chandra X-ray Center, operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. This work made use of data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by NASA. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center and the California Institute of Technology.

Author contributions

M.L.P. wrote the manuscript with comments from all authors and performed the flux-resolved EPIC-pn analysis and line detections. C.P. analysed the RGS data and did the physical modelling. A.C.F. lead the XMM-Newton proposal. All authors were involved with the proposal at various stages.

Author Information Statement

Reprints and permissions information is available at www.nature.com/reprints. Correspondence and requests for materials should be addressed to mlparker@ast.cam.ac.uk. The authors declare no competing financial interests.

Data availability statement

All data used in this work is publicly available. The XMM-Newton observations can be accessed from the XMM-Newton science archive (nxsa.esac.esa.int/nxsa-web/) and the NuSTAR data from the HEASARC archive (heasarc.gsfc.nasa.gov/docs/archive.html). Figure data is available from the authors.

Methods

Data reduction

We use all the available XMM-Newton data, both from our recent observing campaign (PI Fabian) and from the archive. The EPIC-pn data are reduced using XMM-Newton Science Analysis System (SAS) version 15.0.0 epproc tool. The EPIC observations were made in large window mode. We extract source counts from a 30*′′* diameter circular region centred on the source coordinates, and background counts from a $\sim 60^{\prime\prime}$ circular region nearby on the same chip, avoiding contaminating sources, chip edges, and the region where the internal background due to copper is high, and filter the data for background flares. We create separate stacked spectra of the archival and new data using the addspec ftool. We extract full band (0.3–10 keV) lightcurves for each spectrum, shown in Extended Data Figure 1, and divide the lightcurve into low, medium, and high flux intervals such that each flux band contains the same total number of counts (thus the low flux intervals are much longer than the high flux intervals). We then extract spectra corresponding to each flux level from each observation, and combine them using addspec. We bin all the EPIC-pn spectra to a signal-to-noise ratio of 6, after background subtraction, and to oversample the spectral resolution by a factor of 3.

The RGS camera consists of two similar detectors, which have high effective area and high spectral resolution between 7 and 38 Å. The second order spectra cover the 7–18 Å wavelength range and provide double the spectral resolution. We correct for contamination from soft-proton flares following the XMM-SAS standard procedures. For each exposure, we extracted the first- and second-order RGS spectra in a cross-dispersion region of 1′ width, centred on IRAS sky coordinates. We have extracted background spectra by selecting photons beyond the 98% of source point-spread-function. The background spectra were consistent with those from blank field observations. Using the SAS task rgscombine, we stack all RGS 1 and 2 spectra obtaining two high-quality spectra for both the first and the second order with a total, clean, exposure of 1.529 Ms each. We group the RGS spectra in channels equal to 1/3 of the PSF, and use C-statistics, because it provides the optimal spectral binning and avoids over-sampling. RGS spectral fitting is performed using the SPEX package, with contributions from XSPEC, in particular for reflection models. Flux resolved spectra are extracted using the same good time interval (GTI) files as used for the EPIC-pn analysis.

The NuSTAR data are reduced using the NuSTAR data analysis software (NuSTARDAS) version 1.6.0 and CALDB version 20160731. We extract source counts from a 30*′′* diameter circular region, centred on the source, and background counts from a large circular region on the same chip. We combine all the NuSTAR data into a single spectrum, as the count rate is very low due to the extremely soft spectrum, and bin to a signal-to-noise ratio of 6 and oversampling of 3.

Pileup

In the high flux intervals, the source flux is above the nominal EPIC-pn large window mode pile-up limit of 3 counts s*-1*[2], hitting $\sim 9$ counts s*-1* at times. This risks distorting the spectrum and potentially affecting the detection of the UFO. However, the count rate of IRAS 13224-3809 is dominated by photons from below 1 keV (the count rate from 0.5–1 keV is an order of magnitude higher than the 2–3 keV count rate), because it is an extremely soft source. This means that the effects of pile-up are strongest below 2–3 keV. We test this by extracting the same high flux spectrum using an annular region, instead of a circle, with an excised core of $7^{\prime\prime}$ which encircles the central 4 piled up pixels. Above 2 keV, we find no difference in spectral shape between the two spectra, so we conclude that our analysis (restricted to $E>3$ keV) is robust to this effect. The absorption feature is still present in both the mean spectrum and the low flux spectrum when an annular extraction region is used. We also repeat this test using only single events, and again find no difference. The low flux spectrum and the RGS spectra are not affected by pile-up.

Copper contamination

One potential cause of a false detection of a UFO around 8 keV is the complex of emission lines, dominated by Cu K$\alpha$, in the instrumental background[3]. Over-subtracting these features would result in an artificial absorption feature at the corresponding energy, which would depend on the source to background flux ratio, giving an anticorrelation between the equivalent width of the line and the source flux. The copper background is only high in the outer regions of the detector, outside the central  $300^{\prime\prime}$, leaving a central ‘hole’ where contamination is minimal. We are careful to avoid the region where the copper background is high when selecting background regions, which should prevent contamination (see Extended Data Figure 2). The easiest way to show that the UFO line is not an artefact of background over-subtraction is simply to not subtract the background and check the line remains. While this is not always optimal (it may remove genuine but weak lines, or introduce new features), strong absorption features should remain in the spectrum. In Extended Data Figure 2, we show the low flux spectrum with no background subtraction, fit with a power-law. The iron line is weaker, due to the additional high energy contribution from the background, but the UFO line is clearly still present. If the observed line were produced by over-subtraction of the background, the (negative) flux of the line should be constant, the equivalent of an additional constant (positive) line in the background. This is trivial to test, by measuring the strength of an additive line with flux, rather than the multiplicative line we use elsewhere. We find clear variability, and an anti-correlation between the line flux and source flux (Fig. 3, inset), which is impossible if the line is a background feature. Finally, we note that the lines seen in the RGS spectrum are independent of this effect. We conclude that the line is genuine, and produced by absorption in the AGN spectrum.

The potential secondary feature at $\sim 8.7$ keV (observer’s frame, 9.2 keV source frame) is coincidental with the Zn K$\alpha$ line, and appears as an emission feature when the background is not subtracted. We cannot therefore robustly determine whether it is a genuine spectral feature, a statistical fluctuation, or due to the background.

Spectral modelling

EPIC-pn/NuSTAR Stacked Spectra

We fit the stacked 2016 EPIC-pn spectrum from 3–10 keV (outside the band where pileup effects are present, and where the spectrum is relatively simple and unambiguous), and the stacked NuSTAR spectrum from 3–40 keV. We model the spectrum with the relxill relativistic reflection model [4]. The relativisitc blurring parameters are consistent with those found by previous authors[12] (see Table 1), but a strong absorption feature remains at $\sim 8.6$ keV. When we include an additional Gaussian absorption line (modelled with gabs, with $\sigma_{\mathrm{Gauss}}$ fixed at 0.1 keV), the fit improves by $\Delta\chi^{2}=26$, for 2 additional free parameters. Parameters for both these models are given in Extended Data Table 1. We also test allowing $\sigma_{\mathrm{Gauss}}$ to vary, but we find no significant difference in the fit statistic and no impact on the other fit parameters.

There are some differences between this result and those found by previous authors, which likely stem from the different energy range used. In particular, the photon index, high energy cut-off, and iron abundance are different. The continuum parameters are not of great importance to this work, so long as the continuum is adequately described. The steeper $\Gamma$ in archival results ($\sim 2.7$ [16]) is likely due to the inclusion of the soft excess, which past authors have fit with a two-component reflection model. This requires a steep power-law to produce enough soft photons to fit the soft excess. This model is not unique, as the soft excess generally has limited spectral features due to the lower resolution of the EPIC-pn at these energies, and other factors, such as density of the disk, may alter the parameters from such a fit[5]. A visual comparison of the archival data and the new data (Extended Data Figure 3 a) does not show any major changes in the structure of the iron line or UFO absorption.

Similarly, the iron abundance is largely determined by the relative strengths of the iron line and soft excess or Compton hump. Given the steep power-law in the dual-reflection model, a high iron abundance is required to produce enough flux in the iron line. This is not required here, as we do not fit the soft excess and the Compton hump is only weakly constrained. This is important, as the iron abundance is potentially degenerate with the strength of the 8.6 keV absorption feature: an increased iron abundance produces a larger iron absorption edge in the reflection spectrum. We can be confident that this is not having a significant effect on our results, because the iron abundance is free to vary in all our fits, including the fits without the absorption modelled, and the feature still remains. We have explicitly searched for degeneracies using MCMC, and find no degeneracy between the strength of the line and the iron abundance.

Following on from this, we performed a blind line scan over the 6–10 keV band, stepping an unresolved Gaussian line ($\sigma_{\mathrm{Gauss}}=0.01$, allowed to be positive or negative) across the energy band, varying the normalization, and recording the $\Delta\chi^{2}$ at each point on this grid (Fig. 1, top). We use the same underlying relxill model, allowing the same parameters to vary. We calculate the significance of this by taking the probability of the maximum $\Delta\chi^{2}$ for two additional free parameters, and correcting by the number of trials (i.e. the number of resolution elements from 6.7–10 keV). This gives a final chance probability of $1.5\times 10^{-5}$, which corresponds to a $4.3\sigma$ detection. No other features are significant above $\sim 1\sigma$.

We also fit the absorption with a series of physical models - warmabs in Xspec (shown in Extended Data Figure 3 b), which uses grids of xstar photoionization models, and xabs and pion in spex. The three models give consistent results, with a degeneracy between two possible solutions: $v=0.210\pm 0.009$ and $v=0.244\pm 0.09$, corresponding to Fe xxv and Fe xxvi. These solutions have different column densities and ionizations, which are summarized in Extended Data Table 2. The velocity broadening is not strongly constrained, but does not appear to impact any of the other wind parameters. We test this by fixing the broadening to lower and higher values, and find no change in the column density, velocity, or ionization of the fit.

RGS Stacked Spectrum

The RGS spectrum is complex, showing several broad emission-like features at 15 Å and 18 Å. This spectrometer is the most sensitive to narrow ($\lesssim$1000 Å) features, but higher effective area and broader energy range EPIC detectors are more efficient to determine the spectral continuum. We therefore perform an independent analysis of the RGS spectra using either a phenomenological spline continuum model fitted to the RGS spectrum or the physical reflection model provided by the best-fit reflection (relxill) model of the EPIC-pn stacked spectrum. When fitting the RGS spectrum, the spline is corrected for redshift and Galactic interstellar (ISM) absorption. We search for features in the RGS spectrum following an advanced procedure[17]. We include a Gaussian spanning the 7–38 Å wavelength range in increments of 0.05 Å, and assume linewidth of 1000 km s*-1* (comparable to the RGS resolution). This broadening will also tend to strengthen the detection of any warm-absorber and UFO lines with respect to interstellar absorption lines, since the latter are typically narrower[6] ($\leq$200 km s*-1*). We take into account the absorption edges of neutral neon (14.3 Å), iron (17.5 Å), and oxygen (23.0 Å), but we exclude the corresponding 1s–2p absorption lines in order to detect and compare any spectral feature intrinsic to IRAS 13224-3809 or to the ISM. The strongest non-Galactic absorption feature detected is a broad depression around 16 Å, which is also clear in the RGS stacked spectrum (see Extended Data Figure 4). The other two putative, weaker, absorption-like features appear at 10 Å and 13 Å. Interestingly, the photoionization model of the EPIC spectrum predicts three broad ($\sim$1500 km s*-1*) UFO absorption lines that match the three RGS absorption features. We have tested different line widths (from 100 to 5000 km s*-1*) without finding a major effect on their detection. The significance of the rest-frame absorption lines of Galactic O VII and O VIII instead increase up to 5$\sigma$ for narrower widths ($\leq$200 km s*-1*) confirming the results obtained with the grating spectra of the brightest X-ray binaries[6].

A full description of the RGS spectral modelling and the corresponding flux-resolved high-resolution X-ray spectroscopy will be discussed in a forthcoming paper. Here we provide the main result obtained with the overall spectrum and some first interpretation of the wind variability. We have modelled the RGS stacked spectrum with both a spline and a reflection continuum in order to constrain the characteristics of the UFO. The interstellar medium is modelled following the detailed multi-phase gas model constrained with the low-mass X-ray binaries[6]. We model the UFO absorption features in the RGS spectrum with an outflowing gas in photoionization equilibrium (xabs model in SPEX 3.02). The best fit of the RGS stacked spectrum provides the column density $N_{\mathrm{H}}=9.5\pm 0.5\times 10^{22}$ cm*-2* (90% error) the ionization parameter $\log{\xi}=4.0\pm 0.1$ erg cm s*-1* and the line width $\sigma_{v}=2000\pm 1000$ km s*-1*. The RGS velocity shift $v=-0.231\pm 0.007c$ fits between the EPIC Fe XXV ($-0.244\pm 0.009c$) and the Fe XXVI ($-0.210\pm 0.009c$) solutions, and does not fully constrain which solution is most likely, but slightly prefers the Fe XXV ($-0.244c$) solution, which is consistent within the 90% confidence level.

Joint fit

We investigate a combined fit to both EPIC-pn and RGS spectra, fitting with the same absorption model but different continuum models for each spectrum (the physical reflection model for the pn, and a spline for the RGS). We also include a photoionized emission component, modelled with photemis. The soft and hard absorption features are consistent with being from the same absorber (freeing the parameters between the two results in an improvement to the fit of only $\Delta\chi^{2}=3$, for 4 additional free parameters). The joint fit clearly prefers the Fe xxv solution, with final best fit parameters of $v=0.236\pm 0.006c$, $\sigma_{v}=4000\pm 1000$, $\log(\xi)=4.14\pm 0.13$ and $N_{\mathrm{H}}=2.2^{+0.8}_{-1.6}\times 10^{23}$ cm*-2*. The increased broadening with respect to the individual spectrum fits may be due to a small offset between the EPIC-pn and RGS spectra, which could be caused by gain shift in the EPIC-pn. However, it is consistent at the 90% level with that found from the RGS alone.

The inclusion of the emission component improves the fit significantly ($\Delta\chi^{2}=21$, for 2 additional free parameters), accounting for the residuals at $\sim 8.3$ keV and other possible features. The velocity of this component is $0.213\pm 0.015c$, and the luminosity $1.1\pm 0.5\times 10^{41}$ erg s*-1*. If this component is genuine, it is made up of scattered emission from the wind, and can in principle be used to determine the wind geometry. However, it is likely that much of the P-Cygni profile, including any redshifted emission, is obscured by the relativistic iron line, which is very strong in this source. One possible approach to take here would be searching for the emission component of the P-Cygni profile in the lag spectra, as the reverberation timescale should be much longer than for the relativistic reflection component, due to the greater distance from the source.

EPIC-pn Flux-resolved Spectra

We also fit the three flux-resolved spectra, tying parameters we expect to be constant (such as $a$ and $i$) between the different spectra. The model parameters are consistent with those given in Table 1, with the reflection fraction inversely proportional to flux. We perform the same line scan over these spectra simultaneously, stepping the line across in energy then recording the $\Delta\chi^{2}$ for each spectrum individually. The line is only significantly detected in the low flux spectrum, with a maximum $\Delta\chi^{2}$ of 59.7, for 2 additional free parameters. This gives a corrected probability of $1.96\times 10^{-12}$, and a significance of $7.0\sigma$.

We also check the robustness of the low-flux line detection using a Monte-Carlo test. We draw parameters from an MCMC chain, used to evaluate the errors and degeneracies in the best fit parameters, and use them to simulate 10,000 fake spectra. We then fit these spectra using the same procedure. None of the simulated spectra have higher significance features, in either emission or absorption, setting a lower limit of $P>99.99\%$ on the significance. Given the expected fraction of $4.6\times 10^{-12}$, it is unfeasible to test sufficient spectra to establish the true significance with this method.

We perform a similar analysis with 10 flux-resolved spectra, again with the same number of counts in each. We fit the spectra simultaneously, using relxill and a Gaussian absorption line, allowing the reflection fraction, power-law index, and normalizations of the reflection and Gaussian components to vary between each spectrum. This gives a reasonable fit ($\chi^{2}/$d.o.f. $=966/857=1.13$). We then record the equivalent width and flux of the absorption line in each spectrum, and the 3–10 keV flux. These are plotted against each other in Fig. 3, showing a strong correlation. We use Bayesian regression to perform a linear fit which incorporates the upper limits, and draw samples from the posterior distribution to calculate the uncertainty. We calculate the probability of a stronger correlation being found from a constant absorption feature by simulating 10,000 sets of points with the same errors, assuming that the line strength is constant in each case, and performing the same analysis. In no case do we find a stronger correlation.

RGS Flux-resolved Spectra

We perform high-resolution flux-resolved X-ray spectroscopy with the RGS data, consistent with that performed with EPIC: we extract RGS 1 and 2 first- and second-order spectra with the good time intervals defined according to the EPIC flux-prescriptions. We stack the RGS 1 and 2 spectra for each flux-range obtaining three high-quality RGS spectra with comparable statistics. As previously seen for the overall stacked spectrum, there are some non-interstellar absorption-like features (9.5 Å, 13 Å, 16 Å) which show evidence of variability, being both stronger and possibly bluer in the low-flux spectrum. In order to probe the strength of the features in each spectrum, we apply the same technique used for the stacked spectrum by fitting a Gaussian over the 7–38 Å wavelength range in increments of 0.05 Å. In Fig. 2, we show the significance of the spectral features obtained adopting the (RGS fitted) spline continuum. The three broad features are still detected at 9.5, 13.0, and 16 Å in the low-flux spectrum. They have less significance or are undetected in the higher flux spectra with possible evidence of slight velocity change. Their wavelengths match with the strongest lines predicted by the 0.24c UFO model in the RGS energy band: $10.0\pm 0.5$ Å (Ne X + Fe XVIII-XXII blend), $13.2\pm 0.5$ Å (O VIII K$\beta$ + Fe XVIII) and $15.8\pm 0.5$ (O VIII K$\alpha$). The strength of the absorption lines anti-correlate with the flux in agreement with the EPIC result and therefore strongly argues in favour of a connection between the EPIC and RGS absorbers as being part of the same extreme wind. We computed the confidence level of the three main absorption lines in the low-flux spectrum, where they are significantly detected as in EPIC. Accounting for the number of trials due to bins of 0.05 Å and an outflow-velocity range from 0-to-0.3c, we obtain 2.1$\sigma$, 2.9$\sigma$, and 3.4$\sigma$ for the 9.5 Å, 13 Å, and 16 Å absorption lines, respectively, which – since they have the same velocity shift – gives a cumulative 5.1$\sigma$ detection.

Energetics

We can estimate the mass outflow rate by combining the velocity and column densities[3], and the mass of $6\times 10^{6}M_{\odot}$[7] (estimated using the empirical reverberation relation[8]):

[TABLE]

where $\Omega$ is the solid angle of the wind and $R_{\mathrm{wind}}$ the radius.

We cannot be confident of the value of $\Omega$, as the emission from the P-Cygni profile, if there is any, is obscured by the blue horn of the iron line, which is extremely strong in this source. However, given that the absorption line is found in the stacked archival data (most of which is from 2011), this implies that the feature has been present and roughly constant (as a function of flux) for at least 5 years, which would argue for a reasonably large covering fraction, otherwise any clump along the line of sight would likely have moved away. Similarly, we do not know the radius of the wind. However, we know that it must be variable on timescales $\lesssim 5$ ks, which corresponds to 170 $R_{\mathrm{G}}$. Assuming a radius of 100 $R_{\mathrm{G}}$, we find $\dot{M}=2\times 10^{23}\Omega$ g s*-1* ($0.03\Omega\ M_{\odot}$ year*-1*) for Fe xxv, while the Eddington accretion rate for a black hole of this mass is $2.7\times 10^{24}$ g s*-1*, assuming an efficiency of 0.3 for near-maximal spin. In either case, a large fraction of the matter accreted by the disk is lost to the wind, possibly implying super-Eddington accretion at large radii.

We can then calculate the power in the wind:

[TABLE]

For $\Omega=2\pi$, this gives a power of 4% of the Eddington luminosity, implying that a significant fraction of the accretion power must be lost into the wind. For the same assumed covering fraction, the power of PDS 456 is 15% of the Eddington luminosity[3].

Alternative Interpretations

The prevailing interpretation for highly blueshifted absorption features in the X-ray spectra of AGN is that they are due to outflowing gas. However, it is possible that some of these features may instead be due to absorption by a diffuse absorbing surface layer on the approaching side of the accretion disc, which naturally gives relativistic velocities [9]. The absorption line then appears in the reflection component. Aberration means that the blue side is brighter than the red side. For a disc inclination of 60 degrees the absorption layer needs to extend from $\sim 5$ to 10 $R_{\mathrm{G}}$ to give an observed line at 8.2 keV. If the brighter parts of the light curve are associated with the corona rising above 10 $R_{\mathrm{G}}$, then reduced light bending and irradiation of the inner disc weakens both the reflection component and the absorption, consistent with observation.

It is also possible that the variability is produced by a geometric effect. Previous authors have suggested that the relatively constant spectrum of the relativistic reflection component can be produced by changes in the height or extent of the X-ray corona above the disk[10], and the covering fraction of a low scale-height wind could similarly depend on the size or position of the compact X-ray source.

Code availability

All the code used for the data reduction is available from the respective websites. xspec and spex are freely available online. Code used for generating figures, calculating flux-resolved extraction intervals, and calculating line significance, is available upon request to the lead author.