© The Authors 2023

1. Introduction

Mass outflows are a fundamental physical ingredient of active galactic nuclei (AGN). While mass accretion onto supermassive black holes (SMBHs; with typical black hole masses of M_BH ∼ 10^6 − 10 M_⊙) has been long identified as the main physical mechanism powering AGN (e.g., Rees 1984), only recently have AGN been recognized as able to routinely launch powerful mass outflows, which in turn may profoundly affect the galactic environment (see, e.g., Laha et al. 2021, for a recent review). In particular, accretion disk winds – massive outflows launched on subparsec scales – may account for typical observational signatures of luminous AGN, such as the broad emission and absorption lines in their optical and UV spectra (Murray et al. 1995; Proga et al. 2000; Proga & Kallman 2004). Different from radio jets, which are highly relativistic, highly collimated, and present in a fraction (about 15–20%, Kellermann et al. 1989) of AGN, accretion disk winds are near-relativistic, wide-angle flows of matter whose presence has been inferred to be common in luminous AGN through UV and X-ray spectroscopic studies (Weymann et al. 1991; Gibson et al. 2009; Tombesi et al. 2010; Gofford et al. 2013; Chartas et al. 2021; Matzeu et al. 2023).

The innermost regions around the central SMBHs are a rich gaseous environment (see, e.g., Ramos Almeida & Ricci 2017, for a review), and X-ray and UV spectroscopic observations provide information about the physical conditions of the matter reprocessing the X-ray and UV radiation. In the case of mass outflows, the column density, ionization state, and velocity shift of the matter absorbing photons from the line of sight can be measured, and therefore the properties of accretion disk winds can be inferred. In the UV spectra, the most spectacular evidence of winds are the broad absorption lines (BALs) observed to be blueshifted by ∼0.01 − 0.3c in the C IV ion in about 10–15% of optically selected AGN (e.g., Weymann et al. 1981, 1991; Trump et al. 2006; Gibson et al. 2009). Such winds might be present in most AGN, if the covering fraction C_f of the outflowing matter is less than 1, that is, if the wind does not cover all the solid angle as seen by the continuum source, and also if there is an evolution of the wind properties across cosmic time, depending on the AGN’s physical properties (e.g., Ganguly & Brotherton 2008; Giustini & Proga 2019). BALs are not usually observed in local Seyfert galaxies, but in luminous, high-redshift AGN, which are thus called BAL QSOs (where QSOs stands for quasi-stellar objects; e.g., Vietri et al. 2022). In fact, recent spectroscopic observations of luminous AGN at redshift z ∼ 6 revealed the presence of BALs in half of the sample (Bischetti et al. 2022). The intrinsic fraction of BALs in optically bright AGN is consistent with ∼20% at z ∼ 2 − 4 but increases to almost 50% at z ∼ 6 (Bischetti et al. 2023). The mass outflow rates inferred for high-z BAL QSOs is ∼30 − 400 M_⊙ yr⁻¹ (Fiore et al. 2017). Properties of the wind, such as the mass outflow rate and the geometrical covering fraction – and therefore the energy deposited by the wind in the environment – are expected to depend on the AGN’s spectral energy distribution. The AGN spectral energy distribution depends on fundamental physical properties such as M_BH and the Eddington ratio ṁ = Ṁ_BH/Ṁ_Edd (e.g., Ho 1999; Vasudevan & Fabian 2009; Jin et al. 2012), with Ṁ_Edd being the mass accretion rate corresponding to the Eddington luminosity L_Edd = 4πGm_pM_BHc/σ_T, where m_p is the proton mass and σ_T is the Thomson cross section. Understanding the physics of AGN is therefore likely intimately connected to understanding the physics of their accretion disk winds.

X-ray observations of BAL QSOs are challenging. This is first because of the large cosmological redshift of most of the known sources¹, and then because of the observed X-ray weakness of AGN showing BAL features (e.g., Brandt et al. 2000; Laor & Brandt 2002). This X-ray weakness could be intrinsic (Luo et al. 2014) or caused by the heavy X-ray reprocessing close to the central engine (e.g., Gallagher et al. 2001, 2002; Wang et al. 2022). The X-ray properties of BAL QSOs have therefore usually been inferred by means of statistical studies of samples of sources with a very low number of X-ray photons detected (e.g., Green et al. 2001; Gallagher et al. 2006; Giustini et al. 2008; Fan et al. 2009; Gibson et al. 2009; Saez et al. 2012; Sameer et al. 2019). A few individual sources have been studied in more detail thanks to either gravitational lensing (e.g., Chartas et al. 2002, 2003, 2009) or a low cosmological redshift, as in the case of some AGN with BAL-like features identified via space-based UV observations (e.g., Gallagher et al. 2002; Grupe et al. 2003; Schartel et al. 2005, 2010; Hamann et al. 2018). In these cases, moderate-quality spectroscopy can be performed in order to infer at least the general properties of the absorbing material, such as its column density, covering fraction, ionization state, and, importantly, their variations with time (e.g., Gallagher et al. 2004; Ballo et al. 2008; Giustini et al. 2011; Saez et al. 2021).

After the launch of large effective area X-ray telescopes such as XMM-Newton, ultrafast outflows (UFOs) – X-ray-absorbing winds with very large column densities, ionization states, and velocity shifts – were observed in a growing number of AGN (e.g., Chartas et al. 2002, 2003, 2021; Reeves et al. 2003; Pounds et al. 2003; Cappi 2006; Tombesi et al. 2010; Gofford et al. 2013; Matzeu et al. 2023), and the naturally arising question pertains to the relationship between these powerful X-ray absorbing winds and the winds absorbing the UV photons. The answer might have important implications for the theoretical models for accretion disk winds in AGN.

The connection between the X-ray- and UV-absorbing gas has implications for the theoretical models of accretion disk winds in AGN, as X-ray photons generally inhibit the formation of UV-absorbing winds. A large column density of matter able to absorb X-ray photons is needed in radiation-driven disk wind scenarios, in order to prevent the UV-absorbing wind from becoming overionized (and thus failing). This X-ray absorbing gas has been called ‘shielding gas’ or ‘hitchhiking gas’ (Murray et al. 1995); hydrodynamical simulations have shown that such large columns form naturally in the inner regions of the accretion disk atmosphere, where the gas struggles to escape and forms an inner failed wind (Proga et al. 2000; Proga & Kallman 2004). This inner failed wind in fact acts as a filter of X-ray photons for the outermost UV-absorbing wind, which can then be accelerated farther out (see Giustini & Proga 2021 for the meaning of successful or failed wind). Recently, attempts have been made to link the UFOs to the UV BALs, with promising results (e.g., Mizumoto et al. 2021). However, we caution that these results are based on a specific prescription for the gas opacity, which is likely very simplified (e.g., see the discussion in Sect. 4.2 of Nomura et al. 2020). Future hydrodynamical simulations of AGN winds should take into account a realistic treatment of the opacity in the flow, a task that is at the moment beyond the computational possibilities available.

In a first approximation, the observed properties of accretion disk winds are inferred using 1D photoionization codes and simple spherically symmetric geometries, assuming a constant velocity of expansion (e.g., Tombesi et al. 2012; Gofford et al. 2013). Albeit necessary as a first step of interpretation of the observational results, these 1D spherically symmetric scenarios have implications that are possibly not appropriate for treating realistic accretion disk winds around SMBHs. For example, when multiple absorption troughs of the same ionic species are observed at different velocities, these are usually explained by multiple radial zones of the wind. This might not necessarily be the case if the geometry of the wind is not radial and complex dynamical effects on the wind are taken into account (e.g., Giustini & Proga 2012). The same caveats apply to radial distances estimated using the photoionization approximation and the definition of ionization parameter ξ = L_ion/nR² (Tarter et al. 1969), where L_ion is the ionizing luminosity, n is the gas density, and R is the radial distance of the absorbing gas parcel from the source of ionizing photons: multiple ionization states of the gas might actually co-exist inside a wind at the same radial distances, as demonstrated by hydrodynamical simulations (Waters et al. 2021).

The physical properties of accretion disk winds around SMBHs were first studied with hydrodynamical simulations by Proga et al. (2000) and Proga & Kallman (2004), and these simulations showed similarities between the theoretical absorption line profiles and the observed UV properties of BAL QSOs. Synthetic X-ray spectra based on the output of the hydrodynamical simulations of Proga & Kallman (2004) were first presented by Schurch et al. (2009) using the 1D photoionization code XSTAR (Kallman & Bautista 2001), and later by Sim et al. (2010b) using Monte Carlo radiative transfer methods. Despite being based on the output of hydrodynamical simulations that gave predictions for the UV BALs, these studies also showed remarkable agreement between some of the spectral features predicted and the X-ray UFO features observed in luminous AGN. However, these spectral simulations refer to one single point in the parameter space [M_BH,ṁ] (i.e., the M_BH = 10⁸ M_⊙,ṁ = 0.5 input parameters used in the hydro-simulations of Proga & Kallman 2004), while observations of AGN show that a wide range of parameters are at play.

Radiative transfer calculations based on parameterized (non-hydrodynamical) disk wind models are much faster to compute and have been presented, for example, in Kallman & Bautista (2001), Sim et al. (2008, hereafter S08), Sim et al. (2010a, hereafter S10a), and Luminari et al. (2018). These simulations allow us to probe the parameter space in a physically simplified but computationally much faster way and are needed in order to start constraining the main properties of AGN disk winds, leaving behind the use of simple 1D spherical approximations. Future dedicated hydrodynamical simulations of accretion disk winds will focus on interesting points in the parameter space found through parameterized wind model studies. These are thus fundamental steps toward gaining a more realistic physical picture of the inner accretion and ejection flows around SMBHs.

An analysis of the X-ray spectra of the X-ray luminous AGN PDS 456, I Zw 1, and MCG-03-58-007, modeled within the disk wind scenario of S08 and S10a, demonstrated better agreement between the data and the model compared to the use of 1D photoionization codes (Reeves et al. 2014; Reeves & Braito 2019; Braito et al. 2022). This was also the case for PG 1448+273, where the blueshifted Fe K shell trough was successfully modeled by Laurenti et al. (2021) with the disk wind model WINE of Luminari et al. (2018). In this work, we present the results of a similar experiment on the X-ray-weak AGN PG 1126-041, which is known to host BAL-like features in the UV band and both a partially covering absorber and a UFO in the X-ray band.

The galaxy PG 1126-041 is a low-redshift (z = 0.062, Jones et al. 2009) active galaxy with an optical magnitude of M_B = −22.8 (Schmidt & Green 1983), very close to the threshold (M_B = −23) historically used to divide Seyfert galaxies from the more luminous quasars, or QSOs. PG 1126-041 is “almost” a narrow-line Seyfert 1 galaxy (NLS1) with very strong Fe II and very weak [O III] emission and a full width at half maximum (FWHM) of the Hβ emission line of 2150 km s⁻¹: only slightly larger than the classical FWHM < 2000 km s⁻¹ used to define NLS1s (Osterbrock & Pogge 1985). The black hole mass is estimated to be M_BH = 1.2 × 10⁸ M_⊙ by Dasyra et al. (2007) using velocity dispersion measurements in the CO stellar absorption lines.

The UV spectrum of PG 1126-041 shows blueshifted absorption lines in several transitions. An outflow velocity of ∼ − 5000 km s⁻¹ in the C IV and N V species was reported by Wang et al. (1999) based on International Ultraviolet Explorer (IUE) data, while an outflow velocity of ∼ − 2000 km s⁻¹ in the N V, O VI, and P V species was reported by Veilleux et al. (2022) using more recent observations with the Cosmic Origins Spectrograph on board the Hubble Space Telescope (HST-COS). These absorption lines have ionization and velocity similar to BALs but a smaller width (< 2000 km s⁻¹) and are thus called mini-BALs (e.g., Sect. 7.4 of Veilleux et al. 2022 and references therein). In the X-ray band, PG 1126-041 has a relatively large flux and is the second most X-ray-bright PG QSO with known BAL signatures, with a 0.2 − 2 keV observed flux of ∼10⁻¹² erg cm⁻² s⁻¹. There are also clear spectral signatures of strong reprocessing of the nuclear X-ray emission by multiple highly variable ionized gas phases, including a mildly ionized, partially covering absorber and a high-velocity (∼ − 16 500 km s⁻¹), highly ionized component detected with XMM-Newton (Giustini et al. 2011, hereafter G11). We also detected a wind in the Hα and [O III] optical emission lines, observed with the Very Large Telescope (VLT) Multi Unit Spectroscopic Explorer (MUSE) adaptive optics to have a blueshift of a few tens to a few hundred km s⁻¹ on kiloparsec scales (Marasco et al. 2020). The line of sight toward the active nucleus of PG 1126-041 is thus privileged: it offers a view of the nuclear wind on multiple radial scales, in particular of the accretion disk wind originating on UV-emitting and UV-absorbing scales, and also of the disk wind originating on X-ray emitting and absorbing scales. It might therefore hold important clues as to whether such winds are disconnected or to what extent they are connected.

In this work we present an extension of the work of G11, in which we doubled the number of X-ray observations of PG 1126-041, extending the timeline of the study to 11 years (Sect. 2); tested the latest spectral models of accretion disk winds around SMBHs (Sect. 3); and analyzed high-resolution spectroscopic observations of the C IV line profile partially overlapping with the X-ray observations (Sect. 4). We discuss our results in Sect. 5 and present the conclusions in Sect. 6. A cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, q₀ = 0, and Ω_Λ = 0.73 is adopted throughout the paper (Planck Collaboration XIII 2016). The corresponding luminosity distance to PG 1126-041 is 278.4 Mpc. Errors and error bars are at the 1σ level unless otherwise stated.

2. Observations and data reduction

We analyzed eight XMM-Newton pointed observations of PG 1126-041 performed between 2004 and 2015, with exposure times ranging between 10 and 133 ks (principal investigator (PI): M. Giustini for all the observations except for the first one, which has N. Schartel as PI). The first half of these observations was already published by G11. The second half of the data set was taken quasi-simultaneously with observations performed with the HST-COS. Details about the COS observations of PG 1126-041 will be reported in a companion article (Rodríguez Hidalgo et al., in prep.).

All the XMM-Newton observations of PG 1126-041 were taken in full frame mode with the medium optical filter except for the 2004 observation, which was performed with the thin optical filter. The observation data files were processed with the Science Analysis System (SAS) v.18.0.0 using calibration files generated in January 2021. Given the low X-ray count-rate of PG 1126-041, pile-up effects are negligible. The exposure time is also too short for a detection with the high-resolution reflection grating spectrometers (resolving power R ∼ 100 − 500); therefore, we restricted the X-ray analysis to the European Photon Imaging Camera (EPIC) data (R ∼ 10 − 50, see the XMM-Newton Users Handbook for details²). These were reprocessed using standard SAS analysis threads, using the tasks epproc and emproc to concatenate the raw EPIC-pn (pn hereafter) and EPIC-MOS (MOS hereafter) events. A range of thresholds in a count rate of 0.5 − 1.0 (0.35 − 0.5) ct s⁻¹ for high-energy single events (10 keV < E < 12 keV, PATTERN = 0) applied to the light curve of the whole field of view was used to filter out background flares from the pn (MOS) data. Single and double pattern events (PATTERN ≤ 4) or up to quadruple pattern events (PATTERN ≤ 12) with quality flags #XMMEA_EP and #XMMEA_EM, for the pn and the MOS data, respectively, were retained to extract the source and background spectra. Source+background spectra were extracted from circular regions centered on the source coordinates, with optimal radii determined with the SAS task eregionanalyse to be between 27″ and 45″, depending on the signal-to-noise ratio (S/N). Background spectra were extracted from source-free regions of the detector, with circular(annular) shapes in the case of pn(MOS) data. The background extraction regions were always larger than the source ones; the areas were normalized using the backscale SAS task. Ancillary response files and response matrices were generated with the arfgen and rmfgen SAS tasks at the source positions. The MOS spectra and response files were combined with the epicspeccombine SAS task. The XMM-Newton observation identifiers (OBSID), date of observation, exposure time before and after the flaring background filtering, and pn net count rate are reported in Table 1.

3. X-ray spectral analysis

The software Xspec v.12.12.1 (Arnaud et al. 1996) with the python interface pyXspec (Gordon & Arnaud 2021) was used for the X-ray spectral analysis, and the χ² statistic was employed to assess the goodness of fit and estimate measurement errors. The pn and MOS spectra were grouped to a minimum number of 30 counts per energy bin using the task ftgrouppha³ to guarantee a minimum of 20 counts per bin after background subtraction. The energy range considered for the spectral analysis was 0.3 − 10 keV, where the instruments are well-calibrated and their effective area is non-negligible. Bad channels were ignored.

The high-S/N EPIC spectra of the 2009 observation of PG 1126-041 are shown in Fig. 1 unfolded against a power-law model with a photon index of Γ = 2. The black crosses refer to the pn data, while the red open circles refer to the merged data of the two MOS cameras. There is good agreement between the data recorded by the two cameras, and in most of the following plots we will only show the pn data. For each epoch of observation, the pn and MOS spectra were fit to the same model, and uncertainties in cross-calibration between the two instruments were taken into account by adding a multiplicative constant component C_MOS ∈ [0.8, 1.2] between the two spectra of each epoch of observation. A Galactic column density of $N_{\mathrm{H}}^{\mathrm{Gal}}=4.35\times {10}^{20}$ cm⁻² along the line of sight (Kalberla et al. 2005) was applied to all the spectral models using the tbabs model (Wilms et al. 2000).

Fig. 1. EPIC spectra of PG 1126-041 observed in June 2009. Black circles are pn data, while red open circles are MOS data. The spectra are plotted unfolded against a power-law model with Γ = 2.

The 0.3 − 10 keV pn spectra of PG 1126-041 observed by XMM-Newton between 2004 and 2015 are shown in Fig. 2, unfolded against a power-law model (pow in Xspec) with Γ = 2. The X-ray spectra of PG 1126-041 are complex, significantly deviating from the simple power-law model at all the energies probed by the EPIC cameras; this is the signature of strong reprocessing of the intrinsic continuum emission of PG 1126-041 by material along the line of sight.

Fig. 2. XMM-Newton EPIC-pn spectra of PG 1126-041 extracted in the 0.3 − 10 keV band in eight different epochs of observation. Left panel: spectra unfolded against a power law with Γ = 2. Right panel: spectra of each epoch unfolded against a Γ = 2 power-law model are plotted individually over the 2009 spectrum as a reference.

A fit to a phenomenological power law gives a very flat photon index, ⟨Γ⟩∼0.7, compared to both the expected theoretical value (1.5 < Γ < 2.5, Haardt et al. 1994) and the typical value observed in AGN (⟨Γ⟩∼1.8 − 2; see, e.g., Piconcelli et al. 2005). The spectral residuals of the eight epochs of observation for the power-law model are shown in the left column of Fig. 3.

Fig. 3. EPIC-pn spectral residuals of PG 1126-041 to the power-law model [pow] (left column) and to the baseline model [(partcov*xstar500)*(pow)] (right column). The top panel reports the residuals of all eight epochs together, while the smaller individual panels correspond to different epochs of observation, as listed in Table 1.

The X-ray spectral complexity of PG 1126-041 can be reproduced to first order by the addition, along the line of sight, of a layer of ionized gas that is only partially covering the source of X-ray continuum emission (G11). We modeled this gas with the code XSTAR, which computes the physical conditions of a geometrically and optically thin shell of gas illuminated by a point source continuum in the 0.1 − 20 keV energy range, assuming photoionization equilibrium (Kallman & Bautista 2001). We used XSTAR v2.54a to generate a grid of spectra assuming a power-law ionizing continuum with Γ = 2 and a luminosity 10⁴⁴ erg s⁻¹, ionizing a gas shell with a density of n_e = 10¹² cm⁻³ and an intrinsic turbulent velocity of⁴υ_turb = 500 km s⁻¹. The gas column density and ionization parameter were logarithmically sampled at 15 points between N_H = 3 × 10²² cm⁻² and N_H = 4 × 10²³ cm⁻², and at five points between log ξ = 1.5 and log ξ = 3, respectively. From the resulting grid, we generated a multiplicative table to be read into Xspec with the task xstar2xspec, which we named xstar500. In order to account for the absorber only partially covering the source, the xstar500 table was convolved with the partcov model, parameterized by the covering fraction C_f. This is the fraction of the X-ray emission source covered by the absorber, leaving the remaining (1 − C_f) X-ray flux unabsorbed.

The model (partcov*xstar500)*pow, hereafter the baseline model, gives a fit statistic of χ²/ν = 2779/1923, where ν is the number of degrees of freedom, for a joint fit with all the parameters free to vary among the eight epochs of observation. The average photon index is ⟨Γ⟩ = 2.04. The column density along the line of sight ranges from a minimum of $N_{\mathrm{H}}=6.6_{-0.1}^{+0.3}\times {10}^{22}$ cm⁻² during 2008B to a maximum of $N_{\mathrm{H}}=21.0_{-0.6}^{+0.4}\times {10}^{22}$ cm⁻² during 2014A, with an average of ⟨N_H⟩ = 1.5 × 10²³ cm⁻². The average ionization parameter is ⟨log ξ⟩ = 1.94, and the covering fraction ⟨C_f⟩ = 96.3. Fit results are reported in Table 2, and spectral residuals are shown in the right column of Fig. 3. Negative residuals at E > 6 keV are visible in all the spectra except for 2014A, while positive residuals between 4 and 6 keV are visible in half of the spectra.

A blind line search for a narrow (σ = 10 eV) Gaussian line with free normalization and centroid energy applied to the baseline model was performed in the Fe K band between 5 − 11 keV (rest-frame) for each epoch of observation, using a uniform step in energy ΔE = 25 eV (e.g., Miniutti & Fabian 2006). Results are shown in Fig. 4, where the Δχ² contours correspond to 68%, 90%, 99%, and 99.9% confidence level (from the outermost to the innermost, corresponding to 1, 1.6, 2.6, and 3.3 standard deviations σ) in the centroid energy-normalization parameter space in the rest frame of PG 1126-041. The three dashed vertical lines mark the rest-frame energy of Fe I, Fe XXV, and Fe XXVI Kα transitions.

Fig. 4. Results of the scan of χ2 statistical spaces between 5 and 11 keV with a Gaussian line. Top panels: confidence contours at (from the outermost to the innermost) 68%, 90%, 99%, and 99.9% significance levels for the centroid energy and normalization of a Gaussian emission or absorption line applied to the baseline model [(partcov*xstar500)*pow]. The contours are filled with a color intensity proportional to the Δχ2 represented in the color map to the right of each panel. The three dashed vertical lines mark the rest-frame energy of Fe I, Fe XXV, and Fe XXVI Kα transitions. We note the different y-axes in the 2008A and 2008B panels. Bottom panels: observed spectra (filled circles) and corresponding background (shaded area).

Absorption features with a significance level of > 99% are visible in most of the spectra, with energies either between E = 7 − 7.5 keV or at E > 9 keV. The small panels below each contour plot show the observed pn spectra as filled circles and the background as filled areas. During the 2004, 2009, and 2015 observations, there was a strong background emission line at 9 keV (likely due to the Ni, Cu, and Zn in the detector), which is causing a spurious absorption feature once the background is subtracted from the source spectrum. Another emission line at ∼10 keV is present in the 2009 background spectrum. We conclude that the highest-energy absorption features shown in the Δχ² contour plots during the three epochs of observation 2004, 2009, and 2015 are not intrinsic to PG 1126-041 but are an effect of background subtraction. On the contrary, the absorption features between 7 − 7.5 keV are confirmed to be intrinsic to the source at a > 99% confidence level in five out of eight observations. Emission features at a > 99% confidence level are also present in half of the observations.

While the emission features’ centroid energy is compatible with rest-frame neutral iron emission or even redshifted emission, most of the absorption features are blueshifted, indicating that outflowing matter along the line of sight is present in most of the observations of PG 1126-041. The absorption features are too blueshifted to be associated with lowly ionized Fe, and the most conservative identification in terms of derived outflowing velocity is with the highly ionized Fe XXV or Fe XXVI Kα or Kβ transitions (see, e.g., the discussion in Sect. 4.1 of Tombesi et al. 2010). We attempted to model these residuals with two models: the phenomenological model in Sect. 3.1, and the physical accretion disk wind model in Sect. 3.2.

3.1. Phenomenological model

Residuals larger than 3σ above and below the baseline model are present in the spectra of PG 1126-041 at E > 4 keV. The negative residuals indicate the presence of absorbing gas of an even higher ionization state than the log ξ ∼ 2 partially covering absorber; a result already found by G11. We generated a second absorption table with XSTAR, expanding the parameter space toward larger column densities and ionization parameters: the first was sampled at five points between N_H = 5 × 10²² and N_H = 10²⁴ cm⁻², while the latter was sampled at five points between log ξ = 2.5 and log ξ = 5. Based on a fit of the high S/N 2009 data set with warmabs, the turbulent velocity for this grid was set to 5000 km s⁻¹ and the name of the grid to xstar5000.

The baseline model with the addition of the highly ionized absorber [(partcov*xstar500)*(xstar5000*pow)] was fit with the xstar5000 ionization parameter tied between epochs, the column density and the velocity shift free to vary, and all the baseline components free to vary between epochs. The fit statistic is χ²/ν = 2370/1906. The ionization parameter is log ξ ∼ 3.5, the column density is in the range between 3 − 5 × 10²³ cm⁻² and the velocity shift is −(0.045 − 0.11)c. The velocity measured with xstar5000 is a lower limit on the actual velocity of the wind because we could only measure the velocity component projected along our line of sight.

Despite the improvement of the fit statistic by Δχ²/Δν = 409/17 with respect to the baseline model, the positive residuals larger than 3σ at E > 4 keV are not reproduced by the model. We added to the model a Gaussian emission line at the redshift of the source (zgauss in Xspec), with the centroid energy fixed to 6.4 keV, corresponding to the Fe I Kα transition. The width of the emission line was kept constant between epochs, while the normalization of the line was left free to vary. We tested the two scenarios where the emission line is affected (or not) by the partially covering absorber and found better fit statistics in the first case (Δχ²/Δν = 78/9 compared to Δχ²/Δν = 56/9). The line is narrow ($\sigma ={150}_{-50}^{+70}$ eV) and has an average equivalent width (EW) of ⟨EW⟩ ∼ 200 eV, with a minimum during 2008A, 2008B, and 2014B (when only upper limits could be placed) and a maximum during 2014A (EW ∼ 500 eV). The fit statistic for the model [(partcov*xstar500)*(xstar5000*pow + zgauss1)] is χ²/ν = 2292/1898.

In order to account for the residuals observed at energies lower than 6.4 keV, we added a second Gaussian emission line with centroid energy and the width left free to vary, but we kept it linked between epochs. The fit statistic improves by Δχ²/Δν = 116/10 (F-test probability > 99.999%) for a broad (σ = 1100 ± 100 eV) emission line centered at E = 5.3 ± 0.1 keV. The line has an average ⟨EW⟩ ∼ 580 eV. We considered this component not necessarily a physical emission line but a way to model a continuum more complex than a power law. The fit statistic for this model [(partcov*xstar500)*(xstar5000*pow + zgauss1 + zgauss2)] is χ²/ν = 2176/1888. The intrinsic power-law emission has an average ⟨Γ⟩ = 1.98 and carries an average 2 − 10 keV luminosity (corrected for absorption) ⟨L_2 − 10⟩∼1.9 × 10⁴³ erg s⁻¹, with a minimum value of ∼7 × 10⁴² erg s⁻¹ during 2014A and a maximum value of ∼3 × 10⁴³ erg s⁻¹ during 2008B. Spectral parameters are reported in Table 3, while spectral residuals for this phenomenological model are shown in the left column of Fig. 5.

Fig. 5. EPIC-pn spectral residuals of PG 1126-041 to the phenomenological model [(partcov*xstar500)*(xstar5000*pow + zgauss1 + zgauss2)] (left column) and to the disk wind model [(partcov*xstar500)*(fast32*pow)] (right column). The top panel reports the residuals of all eight epochs together, while the smaller individual panels correspond to different epochs of observation, as listed in Table 1.

3.2. Disk wind model

The positive and negative residuals with respect to the baseline model visible in the spectra of PG 1126-041 (see Fig. 4) might be explained by the scattering and absorption of photons in the accretion disk wind, a scenario introduced by S08 and S10a. We tested this scenario using the extended grid fast32, which is a large collection of spectral simulations of the S10a disk wind model extended by Matzeu et al. (2022). The wind is assumed to be smooth and stationary, with a biconical axisymmetric geometry and a wind opening angle of 45° with respect to the polar axis (Fig. 1 of Matzeu et al. 2022). The X-ray source of continuum emission is located at the origin of the coordinate system and has a size of 6r_g (r_g ≡ GM_BH/c² is the gravitational radius, where G is the gravitational constant and c is the speed of light), and the wind inner launching radius is R_min = 32r_g. While in the fast32 grid the slope of the ionizing continuum is a variable parameter, in the extended fast32 grid used in this work it has been fixed, to save computational time, to a power law with Γ = 2. The velocity structure of the wind is assumed to follow a simple β−law, υ(R)∝υ_∞(1 − R_min/R)^β with β = 1, and special relativity effects are taken into account. The ionization state of the wind is computed self-consistently, where K-, L-, and M-shell transitions of the most abundant cosmic ions are taken into account, and absorption, scattering, and reflection of photons into the wind are computed through Monte Carlo radiative transfer methods (see Matzeu et al. 2022 for details). The extended fast32 disk wind model free parameters are the mass outflow rate normalized to the Eddington value ${\dot{\mathcal{M}}}_{\text{w}}={\dot{M}}_{\mathrm{w}}/{\dot{M}}_{\mathrm{Edd}}$; the ratio of 2 − 10 keV luminosity over the Eddington luminosity L_X/L_Edd; the cosine μ of the inclination angle θ between the line of sight and the polar axis; and the ratio between the terminal velocity on the wind streamline and the escape velocity, f_υ = υ_∞/υ_esc, thus from equating the observed velocity with the escape velocity at a radius of 32r_g, for f_υ = 1 and then υ_∞ = 0.25c.

We added the grid fast32 to the baseline model (fit statistics χ²/ν = 2779/1923), resulting in the model [(partcov*xstar500)*(fast32*pow)]. The fast32 parameters ${\dot{\mathcal{M}}}_{\text{w}}$ and μ were kept constant during the different epochs, while L_X/L_Edd was allowed to scale proportionally to the 2 − 10 keV flux (corrected for absorption) in each epoch. We tested three different thicknesses of the wind, with an outer launching radius of R_max/R_min = 1.5, 3, 5, and found the best representation of the data for the case R_max/R_min = 3. In general, a larger wind thickness implies a smaller density for a given column density, and this is compensated by a larger ${\dot{\mathcal{M}}}_{\text{w}}$ (Matzeu et al. 2022). We found a fit statistic of χ²/ν = 2327/1912 for a large inclination angle of our line of sight with respect to the biconical wind polar axis, θ ∼ 82°, meaning that we are looking through the base of the wind. The wind terminal velocity takes into account the geometry of the system and is found to be ⟨υ_∞⟩= − 0.22c, a factor almost 4× larger than the velocity along the line of sight measured with the 1D XSTAR model in the previous section. The mass outflow rate normalized to Eddington is ${\dot{\mathcal{M}}}_{\text{w}}=0.24$, while the 2 − 10 keV ionizing luminosity ranges from 0.1% to 0.47% of Eddington.

The partially covering absorber parameters are consistent with those found with the phenomenological model in the previous Sect. 3.1. In fact, modeling the spectra with the fast32 wind component takes away the need for both the xstar5000 component and the broad emission line at E ∼ 5.3 keV and does not strongly change the covering fraction, column density, or ionization state of the partially covering absorber. The same conclusion holds when the fast32 model is replaced altogether by the relativistically blurred reflection model relxill (García et al. 2014), confirming the ability of the baseline model to reproduce the broadband X-ray spectral shape of PG 1126-041. In this case, the excess of photons in the 4 − 6 keV energy range is well modeled, but the negative residuals in the Fe K band are not taken into account, giving much worse fit statistics than the fast32 model (χ²/ν = 3120/1910) overall.

The narrow emission line at 6.4 keV is instead not reproduced by fast32, and it can either be modeled with a phenomenological Gaussian emission line or with a self-consistent reflection model. In the former case, the improvement in the fit statistic is Δχ²/Δν = 33/9 (F-test probability > 99.7%), with the line detected during the 2009, 2014A, and 2014C epochs. In the latter case, we used the xillver reflection model (García & Kallman 2010; García et al. 2013) fixing the photon index Γ = 2, the density of the reflecting material n = 10¹⁵ cm⁻³, and the metal abundance to the solar value. The reflection component needs to be absorbed by the partially covering gas, and we obtain a Δχ²/Δν = 22/9 (F-test probability > 94%) for an inclination angle < 45° and a low ionization parameter of log ξ < 0.15. In both the Gaussian and the reflection scenarios, the partially covering absorber and the disk wind model parameters are not affected by the inclusion of the Fe K emission components.

The intrinsic power-law continuum has a photon index of ⟨Γ⟩ = 1.85. The power-law normalization best-fit value, which depends on the number of photons lost from our line of sight due to absorption or scattering, is much larger when fast32 is used compared to the phenomenological modeling. In fact, fast32 includes the effect of electron scattering, and for a given optical depth τ = σN_H takes into account the decrease in amplitude of the observed emission, which is ∝e^−τ. This effect is not included in the simple xstar5000 model, and, as a consequence, the measured intrinsic luminosity of the source is underestimated. The average 2 − 10 keV luminosity corrected for absorption computed with the disk wind model is ⟨L_2 − 10⟩ = 5.7 × 10⁴³ erg s⁻¹, with a minimum of ∼2 × 10⁴³ erg s⁻¹ during 2014A and a maximum of ∼9 × 10⁴³ erg s⁻¹ during 2008B, a factor almost 3× higher compared to the phenomenological modeling.

We untied the value of the mass outflow rate between epochs, but we did not obtain a significant improvement of the fit statistic. Spectral parameters of the [(partcov*xstar500)*(fast32*pow)] model, along with their errors and the associated fluxes and luminosities, are reported in Table 4. Spectral residuals are shown in the right column of Fig. 5.

3.3. Partially covering X-ray absorber variability: A Bayesian approach

Independent of the modeling adopted to reproduce the broadband 0.3 − 10 keV spectral shape of PG 1126-041, variability of the partially covering absorber was observed between all epochs of observation. In order to constrain these variations taking into account possible degeneracies between spectral parameters, and thus test the robustness of our spectral fit results obtained with the χ² statistic, we explored the parameter space of both the phenomenological model and the disk wind model using a Bayesian approach. We used the Bayesian X-ray Analysis (BXA) v.4.0.6 (Buchner et al. 2014), which is an interface between Xspec and the Bayesian inference package Ultranest, that uses the most advanced nested sampling algorithm in terms of robustness, correctness, and speed (Buchner 2021). The Bayesian approach allows us to scan the parameter space without losing information due to, for example, spectral binning, without having to make assumptions on the connections between different model parameters, and without having to worry ending up in local statistical minima. The whole parameter space is scanned at once with nested sampling methods, with the likelihood function being marginalized with probability weights given by the prior probability density. We used a Poisson log-likelihood function and analyzed both the grouped and ungrouped spectra, finding consistent results.

We assumed no prior knowledge and fit first the highest S/N, 2009 epoch, and uninformed priors on all the parameters. The power-law photon index was allowed to range between Γ = [1.5 − 2.3] and its normalization between N_1 keV = [3 × 10⁻⁴ − 3 × 10⁻²]; and the partial covering column density, ionization state, and covering fraction between N_H = [3 − 30 × 10²²] cm⁻², log ξ = [1.5 − 2.5], and C_f = [0.8 − 1.0]. In the phenomenological model, the xstar5000 column density, ionization state, and outflow velocity were allowed to range between N_H = [2 − 10 × 10²³] cm⁻², log ξ = [2.5 − 4], and |υ_out|=[0.001 − 0.2]c. To speed up the computational time, when analyzing the data with the phenomenological model we did not include the emission lines (modeled with Gaussians in Sect. 3.1). The presence or absence of these high-energy (E > 5 keV) components has a negligible effect on the estimate of the properties of the partially covering X-ray absorber, which affects lower energy photons. In the disk wind model, the fast32 mass outflow rate, cosine of the inclination angle, terminal velocity, and X-ray ionizing luminosity were allowed to range from ${\dot{\mathcal{M}}}_{\text{w}}=[0.02-1.25]$, μ = [0.05 − 0.9] (i.e., θ between 25° and 87°) to υ_∞/c = [0 − 0.25]c and L_X/L_Edd = [0.026 − 2.5] (the full range of parameters included in the grid).

Figure 2 shows that there are no dramatic variations either in flux or in the spectral shape of PG 1126-041 between the different epochs of XMM-Newton observations. The median values of the posterior probability distributions of the parameters found for the 2009 epoch were therefore used as informed Gaussian priors when fitting the other epochs. All the epochs were fit independently of each other.

Results are reported in the appendix for the 2009 data, where we plot the 1D and 2D histograms of the marginal posterior probability distribution for each spectral parameter (corner plot) using the phenomenological model (Fig. A.1) and the disk wind model (Fig. A.2). In the top right corner of these figures, we plot the posterior probability distributions of the theoretical model (top) and of the model convolved with the instrumental response, with the EPIC-pn data overplotted (bottom). Every solid line is a representation of the model that gives a solution drawn from the posterior probability distribution; therefore, the darker and thicker the line, the more probable the solution.

Most of the statistical solutions for the model parameters found with the χ² minimization are close to the median of the posterior probability distribution found with the Bayesian analysis. The Bayesian parameter estimation scans the whole parameter space at once and is able to clearly show any interdependency between parameters as well as the presence of multiple minima in the χ² space. The parameter interdependencies are visible in the corner plots as 2D histograms that strongly deviate from being symmetric, such as the one for the power-law normalization N_1 keV versus the column density N_H of the highly ionized absorber in the phenomenological model, or the one for the wind inclination angle and terminal velocity in the disk wind modeling. In this latter case, a double minimum in the statistical space is revealed: the solution found with the χ² minimization is not unique but is accompanied by another with lower inclination (θ ∼ 60°) and smaller terminal velocity (υ_∞ ∼ 0.13c).

The partially covering X-ray-absorber parameters, on the other hand, are well-constrained in all the epochs of observation. Figure 6 shows the median values of the posterior probability distribution of the column density, the covering fraction, and the ionization parameter of the X-ray partially covering the absorber measured at different epochs. Here, we plot the values referring to the phenomenological model with open squares and those referring to the disk wind model with filled circles; each error bar represents the 2σ equivalent probability, as derived from the 2D histograms of the posterior probability distribution. The Bayesian analysis confirmed that there are variations in column density of > 10% between every consecutive observation, independent of the model adopted. The average value ⟨N_H⟩ = 13.5 × 10²² cm⁻² is much larger than the values measured for the X-ray warm absorbers and at the lower range of the columns measured in high velocity X-ray UFOs (e.g., Laha et al. 2021). The moderate spectral resolution of the EPIC cameras does not allow us to measure the velocity of the partially covering absorber. A minimum column density is measured in the 2008B observation, N_H = 6.1 ± 0.5 × 10²² cm⁻², and a maximum is measured in the 2014A observation, with N_H = 1.9 ± 0.1 × 10²² cm⁻². The total variation in column density is by a factor larger than 3×.

Fig. 6. Column density of partially covering X-ray absorber versus its covering fraction (top panel) and ionization parameter (bottom panel), derived from the Bayesian analysis of the EPIC-pn spectra in different epochs assuming the phenomenological model (empty squares) or the disk wind model (filled circles). Error bars represent 2σ deviations from the median posterior probability distribution of the parameters.

Remarkably, there are variations of N_H on timescales as short as the separation between the 2014A, 2014B, and 2014C observations (panels d, e, and f in Fig. 2 and green, cyan, and blue points in Fig. 6). The column density decreases by > 20% (ΔN_H ∼ 4 × 10²² cm⁻²) during the 11 days elapsed between 2014A and 2014B and then increases again by ∼20% during the 16 days elapsed between 2014B and 2014C (ΔN_H ∼ 3 × 10²² cm⁻²). Finally, there is a further drop in N_H by ∼30% between 2014C and 2015 (ΔN_H ∼ 6 × 10²² cm⁻²).

4. UV C IV absorption variability

High-resolution (R ∼ 20 000, dispersion 12.23 mÅ/pixel) UV observations were taken very close in time to the XMM-Newton observations using the HST-COS during the 2014-2015 period: HST programs 13429 and 13836 in cycles 21 and 22, respectively (PI: M. Giustini). The average S/N per resolution element of the COS data in the C IV region is ∼7.3. Table 5 shows the gratings and central wavelengths used, the exposure times, and the time difference between the starting time of the HST and XMM-Newton observations. All the HST observations in 2015 were taken during the 18 ks XMM-Newton observation. The 2014 HST observations are very close in time to the corresponding XMM-Newton observations, with the minimum time separation being 0.05 days and the maximum being 0.30 days. Here, we focus on the analysis of C IV absorption, while the full details and analysis of the COS spectrum will be reported in a companion article (Rodríguez Hidalgo et al., in prep.).

Figure 7 (top panel) shows the region around the C IV emission line in all epochs. The spectra in this figure have been smoothed with a boxcar filter 31 pixels wide. In this top figure, the four spectra were only matched in the 1595 − 1605 Å wavelength region, so any observed variability is due to both changes in C IV emission and C IV absorption and not to changes in the continuum emission level. C IV absorption is present from ∼ 1605 − 1612 Å, ∼ 1616 − 1622 Å, and ∼ 1625 − 1637 Å. In order to quantify the strength of these absorption features, we normalized the four spectra taking into account the emission and continuum around the absorption features. Similarly to the work described in Rodríguez Hidalgo et al. (2013), we used second-order polynomial functions to mimic the slope of the blue side of the C IV emission line. We fit the functions to four regions where absorption is not present in either spectra; the same regions were used for all epochs: 1602 − 1605 Å, 1612.5 − 1615.5 Å, 1623 − 1624 Å, and 1637.8 − 1638.1 Å, but different polynomial functions were used for each epoch (see Appendix B). We then measured continuous absorption troughs from 1605 − 1637 Å, defining absorption features as troughs present below 0.9× of the continuum level for > 200 km s⁻¹. The strong absorption line at 1608 Å is likely due to intervening Galactic Fe II absorption and is not variable; thus, it was removed from the spectra by fitting it with a single Gaussian. The bottom three panels of Fig. 7 show a zoomed-in view of the region where outflowing C IV was detected, comparing different epochs in pairs.

Fig. 7. HST-COS spectra o C IV region of PG 1126-041 taken between 2014 and 2015. Top panel: four spectra are plotted overlapped after being normalized and matched at 1595 − 1605 Å to show both changes in emission and absorption. Bottom three panels: zoomed-in view of the 1600 − 1640 Å region, comparing different couples of epochs in velocity scale, after normalization using a second-order polynomial fit. Three velocity systems are marked, as well as three likely C IV doublets in system II. The strong absorption line at 1608 Å is due to Galactic Fe II and was removed prior to the measurements, so it does not appear in the bottom three panels. All spectra in this figure have been smoothed with a boxcar filter with a width of 31 pixels (corresponding to less than 0.4 Å or 75 km s−1) for visualization purposes.

For each absorption feature, we measured the maximum and minimum velocity (υ_max and υ_min, respectively, where zero velocity lies at the quasar redshift z = 0.062 and velocity limits are defined when the flux goes back up to > 0.9× the normalized flux), the EW, and the maximum depth of the C IV absorption troughs. The results are reported in Table 6. Errors in the measurements derive mostly from the systematic uncertainty in the placement of the continuum fit. We followed a similar procedure to Rodríguez Hidalgo et al. (2011).

The C IV absorption in PG 1126-041 displays three distinct absorption systems. System I is the system with the shortest wavelength (1605 Å ≲ λ ≲ 1612 Å), highest velocity absorption. It shows the shallowest depth and the largest variability between all epochs of the three absorption systems. It was not detected in 2015 and was only marginally detected during 2014B. Its maximum velocity is υ_max ∼ −6700 km s⁻¹ in epochs 2014A and 2014C and is slightly lower during 2014B (υ_max ∼ −6530 km s⁻¹). Its depth is maximum during 2014A and significantly decreases during the 2014B epoch, then increases again during 2014C; the EW variations follow a similar trend. System II shows absorption from 1617 − 1623 Å with a maximum velocity of ∼ − 4700 km s⁻¹. It shows minimal variability in velocity but shows significant changes in EW and depth, it being stronger during 2014A and 2014C and weaker during 2014B and 2015. System III is the system with the longest wavelength and lowest velocity absorption (υ_max ∼ −3400) km s⁻¹. It shows the largest EW and depth and almost no variability between any two epochs.

The most extreme variability in the C IV absorption features occurs between observations 2014A and 2014B (second panel of Fig. 7), which are just separated by 11 days (10 days in the rest frame of PG 1126-041). System III overlaps most between observations. However, systems I and II become shallower between these two observations; indeed, system I is present in 2014A at 1606 − 1607 Å but almost disappears in 2014B. Fifteen days later, during the 2014C observation (dark blue in the third panel of Fig. 7), system II returns to a strength and depth similar to the absorption in the 2014A spectrum, but system I remains as weak as in 2014B. In 2015 (pink in the fourth panel of Fig. 7), the absorption profiles are overall the weakest and resemble the ones observed in 2014B.

System II offers a remarkable pattern where the absorption returns to very similar depths at different observations: after weakening in 2014B, the absorption profile in the 2014C spectrum is very similar to the one in 2014A (shown in the third panel of Fig. 7). Similarly, its absorption profile observed in 2015 is very similar to the one in 2014B (see the fourth panel of Fig. 7). System I and the highest velocity part of the absorption complex of system II (υ_max ∼ −4700 km s⁻¹) show the largest variability in EW and depth. System I is also the weakest of the three absorption complexes. This resembles what we observe in high-redshift BAL QSOs, where the strongest variability of the UV absorption troughs is observed in the weakest BALs and in those outflowing at the highest velocities (Capellupo et al. 2011; Aromal et al. 2023).

5. Summary and discussion

5.1. The X-ray spectral properties of PG 1126-041

The nucleus of PG 1126-041 displays significant X-ray spectral variability, remarkably also between the three observations of 2014, which are only separated by about 10–15 days. The observed X-ray flux does not vary dramatically between the eight different epochs of observation except for 2008B and 2014A, which stand out with the largest and lowest observed X-ray flux. The average observed 0.3 − 10 keV flux is 1.5 ± 0.2 × 10⁻¹² erg cm⁻² s⁻¹; it is a factor 2× higher during 2008B and a factor 3× lower during 2014A.

The spectral features that characterize the 0.3 − 10 keV spectra of PG 1126-041 are as follows. First, there is a broadband spectral curvature best reproduced by ionized absorption partially covering the X-ray continuum emission source, which is modeled with a power law with a photon index of Γ ∼ 1.9 (the “baseline model”). Second, there are complexities at E ∼ 4 − 10 keV, with neither X-ray emission nor absorption features taken into account by the baseline model.

The ionized, partially covering absorber gas is detected in all the epochs of observation, while emission and absorption features are detected in about a half of them. The ionization parameter of the partially covering absorber gives the maximum opacity to the continuum photons at E ≲ 3 keV, where a substantial spectral curvature is predicted along with a deep and broad absorption trough from 0.6 − 1 keV. The spectral curvature is a deviation of the observed X-ray photon flux from the power-law continuum emission model. It is due to the photoelectric cutoff and resonant absorption lines and moves to higher energies for larger column densities (e.g., Kallman & Bautista 2001). The broad absorption trough from 0.6 − 1 keV is due to a large number of absorption lines and edges, and its depth is diluted by the presence of the non-negligible fraction of the power-law emission that escapes unaffected from the partially covering absorber. In particular, the decrease in covering fraction increases the soft X-ray flux at E < 2 keV, and this flux dilutes the broad absorption trough.

The different ways in which variations of column density or covering fraction affect the X-ray spectral shape can be seen in Fig. 8. Here, a power-law emission with Γ = 2 is partially covered by a layer of gas with log ξ = 2 and different column densities and covering fractions. The top panel shows the effect of variations of the covering fraction for a fixed column density N_H = 10²³ cm⁻²: at the highest C_f value, the broad absorption trough from ∼0.6 − 1 keV is very deep, while it becomes shallower as C_f decreases. The bottom panel shows the effects of variations in column density for a fixed C_f = 95%: the X-ray flux is absorbed at higher energies for larger column densities compared to lower column densities. These theoretical expectations can be observed in the average flux state spectra (2009), compared to the lowest (2014A, panel d in Fig. 2) and to the maximum flux state spectra (2008B, panel c in Fig. 2). During 2008B, the covering fraction is maximum and the observed absorption trough from ∼0.6 − 1 keV is the deepest; during 2014B, the covering fraction is minimum and the absorption trough is the shallowest. The spectral deviation from a simple Γ = 2 power law happens at E ∼ 3 keV during 2014A, when the column density is maximum, and at E ∼ 1.5 keV during 2008B, when the column density is minimum; the 2009 case is between the two. In summary, in the X-ray spectrum of PG 1126-041, the variations of the partially covering absorber covering fraction can be most appreciated in the soft X-ray band at E ≲ 1 keV, while variations in its column density leave their signature at harder energies, E ≳ 1 keV.

Fig. 8. Theoretical models showing effects of varying covering fraction Cf and column density NH of a log ξ = 2 absorber partially covering a power-law emission with Γ = 2. Top panel: column density is fixed to NH = 1023 cm−2, and the covering fraction increases from Cf = 90% (green), to 95% (red), to 99% (black); the solid lines represent the total models, while the dashed lines represent the (1 − Cf) fractions of unabsorbed power-law emission. Bottom panel: covering fraction is fixed to Cf = 95% (the 5% fraction of the unabsorbed power law is plotted with a dashed line), and the column density increases from NH = 5 × 1022 cm−2 (green), to NH = 1023 cm−2 (red), to NH = 2 × 1023 cm−2 (black). The normalization of the y-axis is arbitrary.

The spectral complexities with respect to the baseline model are visible in the right column of Fig. 3. The results of a blind search for a Gaussian line either in emission or in absorption in the Fe K band of the eight epochs are shown in Fig. 4. Residuals in absorption significant at > 99.9% confidence are present during 2004, 2009, and 2015. During 2009 and 2014C, there was a prominent and broad emission feature extending redward of 6 keV (rest-frame). The negative residuals can be reproduced by a highly ionized outflowing absorber, the positive residuals by a phenomenological Gaussian emission line (Sect. 3.1), or all together by a disk wind model that takes into account both the X-ray photons absorbed along observer’s line of sight and those scattered back into it (Sect. 3.2).

Modeling the highly ionized absorber with the 1D photoionization code XSTAR gives a velocity blueshift ⟨υ_out⟩ of about −0.06c. This is larger than the 10 000 km s⁻¹ threshold used to define ultra-fast outflows (UFOs, Tombesi et al. 2010); therefore, we refer to this component as a UFO in the following. The UFO velocity observed in PG 1126-041 is at the lower end of the UFO velocity distribution observed in local Seyferts, while the column density ⟨N_H⟩∼5 × 10²³ cm⁻² and the ionization parameter log ξ ∼ 3.5 are near the average (Tombesi et al. 2011; Gofford et al. 2013).

The Gaussian emission line used to model the residuals in excess of the baseline model is very broad (σ ∼ 1 keV) and centered at E ∼ 5.3 keV in the source rest frame; therefore, it is unlikely that it corresponds to a physical individual emission feature. An excess of counts at E ∼ 4 − 6 keV is often observed in the X-ray spectra of local AGN, and it is often interpreted within a relativistic reflection scenario or a complex partial covering absorption scenario (e.g., Fabian et al. 2002; Mizumoto et al. 2014). The two scenarios often give statistically equivalent fits to the data; for example, in Mrk 335 (e.g., Gallo et al. 2013, 2015; Grupe et al. 2008), in 1H 0419-577 (e.g., Fabian et al. 2005; Turner et al. 2009; Di Gesu et al. 2014), in 1H 0707-495 (e.g., Gallo et al. 2004; Tanaka et al. 2004; Dauser et al. 2012), and in PG 1535+547 (Ballo et al. 2008). In some cases, both relativistic reflection and complex absorption may contribute to shaping the X-ray spectra of AGN (e.g., Risaliti et al. 2009a; Patrick et al. 2012; Parker et al. 2021).

Both the photons in emission at E ∼ 4 − 6 keV and the photons absorbed at E > 7 keV might also be produced by scattering and absorption in an accretion disk wind (Sim et al. 2008, 2010a,b). It was recently demonstrated by Parker et al. (2022) that the emission expected from the accretion disk wind models is in fact spectrally degenerate, with the relativistic reflection emission expected close to the SMBH, but the latter does not include the effects of absorption along the line of sight. The accretion disk wind scenario has been successfully applied to reproduce AGN spectra at E > 2 keV (Tatum et al. 2012; Hagino et al. 2015, 2016) as well as the broadband X-ray spectra of individual AGN, namely PDS 456 (Reeves et al. 2014), I Zw 1 (Reeves & Braito 2019), and MCG-03-58-007 (Braito et al. 2022).

The accretion disk wind model applied in this work is an extension of the spectral grids presented in Matzeu et al. (2022). The data require a very large inclination angle of the line of sight with respect to the biconical wind polar axis, θ ≡ arccosμ ∼ 80°. In this case, our line of sight goes through the wind base, and the observed spectrum is dominated by emission reprocessed by the wind. Since the wind is assumed to be biconical and the inclination angle high, the terminal velocity of the wind is much larger than the observed projected velocity, υ_out ≪ υ_∞ ∼ −0.2c. While in principle a wind launched at larger radii (and therefore with a lower terminal velocity) could lead to the same observed projected velocity for a lower inclination angle, the main constraint for such a large inclination angle comes from the very large depth of the Fe K absorption trough. This is shown in Fig. 9, with a simplified sketch of two geometries for the biconical accretion disk wind in PG 1126-041: the one corresponding to our best-fit scenario on the left in blue, and a wind launched at larger radii on the right, which has a terminal velocity close to the projected velocity observed in PG 1126-041, in red. The middle panel shows the theoretical prediction for the Fe K band in the two scenarios: the wind launched farther out allows for a lower inclination angle (about 65°), almost parallel to the wind streamline, giving the same blueshift as the wind launched closer in and observed at a larger inclination angle (about 80°), but with much narrower absorption line profiles. The bottom panel shows the high S/N 2009 EPIC-pn data compared to the two models; the slower wind predicts a shallower absorption trough in the Fe K band compared to the faster wind, which is therefore preferred by the data. Future X-ray microcalorimeter observations of the Fe K band of PG 1126-041, with, for example, Resolve on board XRISM (XRISM Science Team 2020) or X-IFU aboard Athena (Barret et al. 2018, 2023), should allow a definitive distinction between the two scenarios.

Fig. 9. Two different geometries for the interpretation of the blueshifted absorption troughs in the Fe K band of PG 1126-041. A wind with a smaller Rin (corresponding to υ∞ = −0.25c) and observed at a larger inclination angle (θ = 80°) is shown in blue, while a wind with a larger Rin (corresponding to υ∞ = −0.0625c) observed at a smaller inclination angle (θ = 65°) is shown in red. Top panel: sketch of the geometry of the disk wind in the two cases (only Rin has been sketched for simplicity), with streamlines represented by arrows with different thicknesses, proportional to the terminal velocity. Central panel: two theoretical models. Bottom panel: 2009 pn data (in gray) overplotted on the two models folded with the instrumental response.

Such a large inclination angle is so close to the base of the wind (or the atmosphere of the accretion disk) that very large column densities -highly Compton-thick- are expected by more realistic numerical simulations (Sim et al. 2010b). For such equatorial lines of sight, the dusty torus on parsec scales might also be intercepted, if present (e.g., Ramos Almeida & Ricci 2017); however, the geometry of the torus is expected to vary with the evolution of an AGN (e.g., Hopkins et al. 2012). Furthermore, the biconical geometry of the wind considered here is simple, and the physics ignores the effects of, for example, gas pressure and magnetic fields. Therefore, the results that we present should be taken as indicative and not as absolute measurements. It is possible that future hydrodynamical simulations could predict such large velocities and absorption depths for winds launched at larger radii than those considered here, or observed at lower inclination angles. In any case, the line of sight toward the nucleus of PG 1126-041 is likely intercepting the wind at every epoch of observation.

The intrinsic power-law photon index is found to be Γ ∼ 1.8, slightly flatter than the Γ = 2 assumed in both the XSTAR and fast32 model calculations. The main effect of increasing the input Γ in the calculation of the fast32 model is to produce deeper absorption troughs in the Fe K band (Matzeu et al. 2022), thanks to a smaller number of hard X-ray photons able to over-ionize the iron atoms of the wind. Therefore, by assuming a slightly higher Γ in our calculations, we might have underestimated the amount of matter necessary to produce the deep absorption trough, which in the model is parameterized by ${\dot{\mathcal{M}}}_{\text{w}}$. The current uncertainties on the value of Γ should be significantly reduced in the near future, thanks to approved broadband X-ray observations with XMM-Newton + NuSTAR (PI: J. N. Reeves) that should reveal the intrinsic continuum slope of PG 1126-041. Overall, the disk wind model reproduces the X-ray spectra of PG 1126-041 well, if highly variable; massive clumps, which are represented here by the partially covering X-ray absorber, are included in the modeling.

5.2. The X-ray/UV connection

In the X-ray band, we detected two absorbers: one highly ionized ultra-fast outflow absorbing mainly in the Fe K band (E ∼ 7 − 10 keV), and one partially covering the X-ray continuum emission source and at a lower ionization state, with the largest opacity between E ∼ 0.5 − 2 keV. While both absorbers are found to be variable with time, the spectral variability observed in PG 1126-041 is dominated by variations in the column density of the partially covering absorber (Sect. 3.3). In the UV band, we focused on the absorption in the C IV region, which shows three systems of blueshifted absorption lines with maximum velocities between −3300 and −6700 km s⁻¹. These systems show variability in a coordinated way between different epochs, which is strongest during 2014A and 2014C and weakest during 2014B and 2015.

The variability observed in the C IV-absorbing wind appears to be coordinated with the column density variations of the partially covering X-ray absorber. The 2014A and 2014C epochs are the ones with the largest X-ray-absorbing column density of the partially covering gas (N_H ∼ 1.8 − 1.9 × 10²³ cm⁻²). The partially covering X-ray absorber column density is a factor of 20% and 40% smaller in 2014B and 2015, respectively, than in 2014A and 2014C. Concurrently, the C IV absorption line profiles are much more similar and weaker in 2014B and 2015 compared to 2014A and 2014C. In Fig. 10, we plot the EW of each C IV absorption system against the median of the posterior probability of the partially covering X-ray absorber column density for each epoch of observation, computed with the [(partcov*xstar500)*(fast32*pow)] model. This coordinated variability strongly suggests that the partially covering X-ray absorber and the highly blueshifted C IV absorbers are connected. The UV absorption troughs observed in BAL QSOs are often inferred to be shaped by partially covering absorption of the UV-continuum source, and the covering fractions are comparable to those obtained for the partially covering X-ray absorber in PG 1126-041 (C_f ∼ 80 − 90%; see, e.g., Arav et al. 1999; Rodríguez Hidalgo et al. 2011). This might suggest similar geometries for the UV-absorbing gas and the X-ray-absorbing gas relative to their sources of continuum emission.

Fig. 10. EW of C IV absorption systems against the column density of the partially covering X-ray absorber in each epoch of coordinated XMM-Newton/HST observations of PG 1126-041. From left to right: system I, system II, and system III (system I was not detected in 2015).

It would be interesting to compare the velocities of the outflowing absorbers. Given the limited spectral resolution of the EPIC detectors, the multitude of closely spaced spectral transitions produced in the log ξ ∼ 2 partially covering absorber are not individually resolved in our X-ray observations, and a velocity shift cannot be measured. On the contrary, the UFO with log ξ ∼ 3.5 only imprints a few strong transitions in the Fe K band, which can then be used as a velocity shift marker even at moderate spectral resolution. The velocity projected along the line of sight of the UFO has a blueshift of ∼0.06c (Table 3), a factor of three larger than the maximum outflow velocity projected along the line of sight for the UV absorber. When the UFO is modeled within the accretion disk wind scenario, the biconical geometry assumed for the wind gives an even larger deprojected velocity of ∼ − 0.2c. In order to compare this value with the UV maximum velocity, a geometry for the UV-absorbing wind should also be assumed.

One possibility is that, from farther out to closer toward the central SMBH of PG 1126-041, our line of sight goes through the UV-emitting and -absorbing region, then encounters the partially covering X-ray absorber, then the X-ray UFO, and finally the intrinsic continuum source. We might therefore be observing an innermost wind, where the highly ionized (log ξ ∼ 3.5) X-ray absorption is produced and the matter is accelerated (up to υ_∞ ∼ −0.2c), together with the portion of the outflow situated farther out from the acceleration zone, which is thermally unstable and therefore clumpy and produces the partially covering X-ray absorption (Waters et al. 2022). The partially covering X-ray absorber then acts as a “filter” (a variable patchy screen) of the incoming ionizing soft X-ray photons for the UV-absorbing gas (Misawa et al. 2007). The UV absorber is exposed to a number of X-ray photons proportional to (1 − C_f)+C_f e^−τ, where the first term is the unobscured photon flux and the second term is the photon flux making it through any shielding gas. Generally speaking, the closer to the SMBH the wind launching point is, the faster its terminal velocity must be. Assuming that the observed velocities are proportional to the wind terminal velocity, the following scenario might explain the observations: when more massive clumps along the line of sight are covering the X-ray-continuum source (e.g., during 2014A and 2014C), a lower number of X-ray photons are reaching the UV-absorbing wind. When the X-ray absorption along the line of sight is reduced (e.g., during 2014B and 2015), then the X-ray flux reaching the C IV-producing region is larger, thus over-ionizing the UV-absorbing wind.

Hints as of the nature of the partially covering X-ray absorber come from the existence of an anti-correlation between its column density and covering fraction in different epochs, as shown in the top panel of Fig. 6. One possibility is that the 1 − 10% of X-ray light that escapes unaffected by the partially covering absorber is in fact scattered light. As the column density increases, so does the electron scattering optical depth, which allows more photons from the background (continuum) source to scatter off the electrons of the absorber without changing their spectrum. Scattering might take place within the partially covering absorber itself and in the accretion disk wind. One way to differentiate between these two scenarios would be to use X-ray polarimetric observations. In fact, the flux from ∼0.6 − 1 keV is expected to be polarized if it is produced by electron scattering, while this is not the case for directly transmitted flux.

Similarly to PG 1126-041, the blueshifted absorption line system in the UV spectrum of the AGN HS 1603+3820 was shown to display coordinated variability on timescales of years (Misawa et al. 2007), which could be explained either by electron scattering or by the effects of a variable clumpy “screen” between the mini-BAL gas and the continuum source. The latter hypothesis was supported by UV spectropolarimetric observations (Misawa et al. 2010). The physical characteristics of the X-ray absorbers in sources with intrinsic blueshifted UV absorption are not clear. For example, Hamann et al. (2013) conducted a study of eight AGN with high-velocity (outflow velocities ∼0.1 − 0.2c) UV-absorbing winds, finding no significant X-ray absorption. They suggest that the UV absorber ionization is low, not due to shielding gas but due to high densities likely resulting from magnetic confinement. However, only cold (neutral) X-ray absorption was considered, and the limits in neutral column density found (N_H ∼ 0.3 − 5 × 10²² cm⁻²) might be consistent with larger column densities of gas along the line of sight for more complicated scenarios such as ionized or partially covering X-ray absorption. On the basis of the observed spectral complexities, Giustini (2016) showed that the X-ray spectra of AGN with BAL and mini-BAL features are consistent with large column densities N_H = 10^{22 − 24} cm⁻² of gas along the line of sight, if the ionization state is significantly different from neutral, or if the gas is only partially covering the X-ray continuum source.

5.3. PG 1126-041 in the grand scheme of AGN

The AGN PG 1126-041 is low-redshift with powerful nuclear winds and is thus similar to BAL QSOs, which usually have a much lower X-ray flux and therefore fewer detected photons at the detectors, which prevents complex spectral models being applied. The general observed spectral properties of PG 1126-041 and similar sources might therefore be used to infer the general properties of large samples of sources with a much lower number of X-ray photons available, thus paving the way for future studies with larger X-ray telescopes.

The partial obscuration of the central X-ray continuum emission source by large column densities of ionized gas is common in AGN, and, in fact, partially covering absorbing gas might be an AGN ingredient on all scales. Partially covering X-ray absorption by cold, Compton-thick clouds transiting on torus or broad line region scales (∼10³ − 10⁶r_g) has been detected in several AGN, notably NGC 1365 (Risaliti et al. 2005, 2009b; Maiolino et al. 2010; Rivers et al. 2015), Mrk 766 (Risaliti et al. 2011), H0557-385 (Longinotti et al. 2009). Thanks to long temporal baseline observations of a large sample of AGN, the occurrence of such occultation events is relatively common (Markowitz et al. 2014). Following deep monitoring studies, several historically unabsorbed AGN have recently been discovered that are affected by transient obscuration events, which depress the observed X-ray flux and where blueshifted UV absorption lines are observed to emerge (e.g., Kaastra et al. 2014; Ebrero et al. 2016; Mehdipour et al. 2017, 2021; Kriss et al. 2019). Partially covering X-ray-absorbing gas has also been detected very close to the central SMBH, as inferred by the rapid variability, in the luminous AGN PDS 456 (at tens to hundreds of r_g; e.g., Nardini et al. 2015; Matzeu et al. 2016; Reeves et al. 2018a).

In the case of PG 1126-041, we can use the time variability constraints between the closely spaced 2014 observations to obtain an estimate of the distance of the partially covering X-ray absorber from the continuum source. The velocity of the partially covering absorber cannot be constrained with our data, but we can consider the range in velocities given by the most blueshifted UV absorber (−0.02c) and the deprojected X-ray UFO (−0.2c). If one assumes that the transverse velocity is comparable to the radial velocity, so that ΔR is the distance the absorber moves during time Δt, then the variations observed on timescales of about ten days would correspond to a thickness ΔR ∼ 5 × 10^{14 − 15} cm. Given n ∼ N_H/ΔR, one would obtain n ∼ 2.6 × 10^7 − 8 cm⁻³ for ⟨N_H⟩∼1.3 × 10²³ cm^−2. Inserting these density estimates into the ionization parameter definition ξ = L/nR² and using L = 10⁴⁴ erg s⁻¹ and log ξ = 2, one would obtain a distance for the partially covering absorber R_PC ∼ (6.2 − 20)×10¹⁶ cm. This distance corresponds to about 3500 − 11 000r_g for a black hole mass estimate of M_BH = 1.2 × 10⁸ M_⊙ (Dasyra et al. 2007).

Assuming that the variations in column density are due to a single cloud moving across the line of sight, another estimate of the location of the partially covering X-ray absorber can be placed by calculating the absorber transverse velocity necessary to move across the X-ray emitting region during the shortest separation between consecutive observations and comparing it to the orbital (Keplerian) distance. The size of the X-ray-emitting region in PG 1126-041 was constrained to be < 13r_g by Giustini et al. (2011), using the duration of a continuum flare of 8 ks. An absorber with a transverse velocity of υ_K ∼ 0.008c would cross such a region during the 11 days elapsed between 2014A and 2014B, and this velocity would correspond to orbits at ∼15 000r_g from the central SMBH. An absorber with a higher velocity would of course cross the same region in a shorter time, for example in about 4.4 days at the velocity of the UV absorber and in about 10.5 h at the velocity of the X-ray UFO. The partially covering absorber would therefore be on the broad line region scales, consistently with the observed variations in UV-absorbing gas coordinated with variations in X-ray absorbing column density.

For comparison, the location of the UFO can be estimated in an analogous way using the variations observed on timescales of about 8 ks during 2009 (Giustini et al. 2011), an average column density 5 × 10²³ cm⁻² and column density variations ΔN_H ∼ 2.5 × 10²³ cm⁻², obtaining a distance estimate of R_UFO ∼ 30 − 40r_g. The distance of the UFO to the X-ray continuum source might therefore be comparable to the size of the latter; this would have implications for the detailed theoretical modeling of the X-ray UFO. As in this case, the important assumption of a point-source continuum as seen by the absorber may no longer be valid.

In any case, our distance estimates should be considered simple approximations, because in our calculations we neglected the presence of the wind itself, with all the physical consequences this might have. These are, notably, a different (than spherical) geometry and a highly dynamical environment, where the effects of pressure and temperature gradients across the flow are expected to generate much more complex observational scenarios (e.g., Giustini & Proga 2012; Ganguly et al. 2021; Waters et al. 2022). What can be said with our distance estimates is that the partially covering X-ray absorber either lies between the UFO and the UV-absorbing gas, or it is co-spatial (within 1.5 light hours, the minimum time separation between XMM-Newton and HST-COS observations of PG 1126-041) with the latter. Dynamical thermal instability has been recently demonstrated to be a mechanism that can make AGN winds clumpy beyond the acceleration zone (Waters et al. 2022).

The X-ray obscurer detected in NGC 5548 since 2012 has a comparable covering fraction, column density, and ionization state to the partially covering X-ray absorber in PG 1126-041 (Kaastra et al. 2014; Mehdipour et al. 2016). Multi-year monitoring of NGC 5548 in the X-ray and UV band between 2013 and 2022 showed coordinated variability of the C IV EW and the amount of X-ray obscuration, as in PG 1126-041 (Mehdipour et al. 2022). The column density and the ionization state of the partially covering absorber in PG 1126-041 are also similar to those of the obscurer detected in NGC 3783 (Mehdipour et al. 2017; De Marco et al. 2020; Costanzo et al. 2022), although the covering fraction in the latter AGN is much lower (by about a half) than the covering fraction in PG 1126-041. A good match in covering fraction and column density is found between the partially covering X-ray absorber in PG 1126-041 and the X-ray obscurer in Mrk 817 (Kara et al. 2021), which shows variability of column density on timescales of days and weeks and coordinated intensity of the phosphorus ion P V absorption troughs similar to the one observed in the C IV absorption troughs in PG 1126-041. Therefore, the partially covering X-ray absorber in PG 1126-041 has the characteristics of a massive, highly ionized X-ray obscurer.

The general X-ray spectral properties of PG 1126-041 are also remarkably similar to those of PDS 456, the most well-studied luminous QSO hosting an X-ray-absorbing accretion disk wind (Reeves et al. 2003). PDS 456 is absorbed by variable ionized gas and has one or more UFOs with observed velocities υ_out = −(0.25 − 0.3)c (Reeves et al. 2009; Hamann et al. 2018). The observed velocity shift in PG 1126-041 is much smaller, ⟨υ_out⟩∼ − 0.06c; however, υ_out is the projected velocity along the line of sight, and therefore a lower limit on the actual wind velocity. When taking into account the geometry of the accretion disk wind, a terminal velocity ⟨υ_∞⟩∼ − 0.2c was recovered for the accretion disk wind of PG 1126-041. Therefore, the differences in observed velocity shift in PDS 456 and PG 1126-041 might be mainly due to a different inclination angle of our line of sight with respect to the wind: in PDS 456 the inclination angle is smaller, of the order of 50° (Matzeu et al. 2022), compared to the 80° derived for PG 1126-041. Recently, an even larger component with υ_out = −0.45c was unveiled in the X-ray spectrum of PDS 456 using a broadband (joint XMM-Newton + NuSTAR) observation (Reeves et al. 2018b). Future hard X-ray observations of PG 1126-041 should help clarify whether there is a higher velocity UFO phase in this AGN as well. In any case, it is likely that the accretion disk wind in PG 1126-041 is a scaled-down version of the extremely powerful wind of PDS 456 in terms of energetics. The black hole mass estimate for PDS 456 is about one order of magnitude larger than for PG 1126-041; therefore, all the timescales of variability in PG 1126-041 should be more rapid by a factor of ∼10. Given the observed variability of the partially covering X-ray absorber in PDS 456 on timescales of 100 ks (Matzeu et al. 2016), it could be worth investigating whether such a component varies on timescales of ∼10 ks in PG 1126-041.

6. Conclusions

The results of the analysis of eight XMM-Newton observations of PG 1126-041 between 2004 and 2015, the last four of which were taken quasi-simultaneously with HST-COS exposures, are as follows.

• The spectral shape of PG 1126-041 from 0.3 − 10 keV is complex, showing strong reprocessing of the intrinsic X-ray emission. A fit to a phenomenological power law gives a very flat photon index Γ ∼ 0.6; however, once accounting for complex absorption along the line of sight, a more typical Γ ∼ 1.9 is measured.

• A massive partially covering X-ray absorber is detected in all epochs, independently of the underlying broadband continuum modeling. The ionization parameter is log ξ ∼ 2, the covering fraction C_f ∼ 95%, and the column density is in the range of N_H ∼ (5 − 20)×10²² cm⁻².

• X-ray spectral variability of PG 1126-041 is observed among every epoch probed, with time separations as short as 11 days; the spectral variability is driven by the variable column density of the partially covering absorber.

• The column density of the partially covering X-ray absorber shows coordinated variability with the highest-velocity components of the C IV absorber observed in the UV with the COS during the last four epochs of observation. In particular, during 2014A and 2014C the X-ray absorbing column density was at its maximum and the C IV absorber had the maximum EW; during 2014B and 2015, the column density was smaller by 20% and 40%, respectively, and the C IV absorber had the minimum EW.

• In addition to the partially covering absorber, a highly ionized outflowing absorber (UFO, detected with a statistical confidence > 99.9% in 3/8 epochs) and a broad emission feature modeled with a phenomenological Gaussian emission line centered at E = 5.45 keV and with a σ ∼ 1 keV (present in half of the epochs) are needed to reproduce the spectral shape of PG 1126-041 in the Fe K band. Both components can be accounted for by a biconical accretion disk wind model observed along an equatorial line of sight.

• When the UFO is reproduced by 1D photoionization models, a column density of ∼6 × 10²³ cm⁻², an ionization parameter of log ξ ∼ 3.5, and a velocity projected along the line of sight υ_out ∼ −0.06c are measured. When the UFO is modeled within the accretion disk wind scenario, the biconical geometry assumed for the wind gives a deprojected terminal velocity of υ_∞ ∼ −0.2c.

• While the accretion disk wind model used ignores the effects of gas pressure in computing the physical structure of the outflow and assumes a very simple geometry of a thin cone, this is a necessary step toward a more realistic treatment of AGN accretion disk winds. For example, it takes into account both scattering and absorption of X-ray photons, and the measured intrinsic X-ray luminosity is therefore larger (by a factor of ∼4×) compared to the one measured using 1D photoionization codes, which only consider the effects of absorption of photons along the line of sight.

• Assuming a launching radius for the wind inversely proportional to the observed velocity, a possible scenario to interpret the XMM-Newton/HST-COS observations of PG 1126-041 is the following: a wind with an inner ultrafast component launched at a few tens r_g from the central SMBH (the UFO) fragments beyond its acceleration zone due to thermal instability; the subsequent clumps that form in the flow might be the partially covering X-ray absorber, which might be co-located (or at distances shorter than 10 light days) with the UV-absorbing wind. The partially covering X-ray absorber clumps act as a variable screen between the X-ray photons transmitted through and scattered off the wind and the UV absorber. Therefore, in epochs of low X-ray column density, the highest velocity, innermost UV-absorbing wind would be destroyed by the large flux of X-ray photons able to reach it.

• We will need to account for the presence of massive clumps in the structure of accretion disk wind models at some point in the future.

This work adds PG 1126-041 to the growing number of AGN studied with parameterized accretion disk wind models based on radiative transfer calculations. These studies can allow us to estimate physical parameters of the inner accretion and ejection flow around SMBHs and are a necessary step toward the development of full-hydrodynamical disk wind models that might explain all the observed complexities. Coordinated multiwavelength spectroscopic analyses of AGN are very promising in constraining the dynamical properties of their accretion disk winds. The line of sight toward PG 1126-041 offers a privileged view through a highly dynamical accretion disk wind, and it is worth exploring with future studies in order to shed light on the general connection between the accretion and the ejection flow around SMBHs.

Acknowledgments

We thank the referee for the comprehensive and constructive review that helped improve the presentation of the results. MG is supported by the “Programa de Atracción de Talento” of the Comunidad de Madrid, grant number 2018-T1/TIC-11733. P.R.H. acknowledges support from start-up funding provided by UW Bothell, from start-up funding from the College of Natural Resources and Sciences, and from the Sponsored Programs Foundation at Humboldt State University (now Cal Poly Humboldt) through a Research, Scholarly, & Creative Activity grant. GP acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 865637). JNR acknowledges support from NASA cooperative grant 80NSSC22K0474. VB acknowledges support from NASA cooperative grant 80NSSC22K0220. GM was partly supported by grant n. PID2020-115325GB-C31 funded by MCIN/AEI/10.13039/501100011033. We thank Stuart S. Sim for making his disk wind model available to help us interpret our observations, Jeremy Sanders for writing Veusz, Craig Gordon and Keith Arnaud for developing pyXspec, and Johannes Buchner for developing UltraNest and BXA and for kind and helpful discussions. Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. We acknowledge financial support from NASA/XMM-Newton via the grant NNX10AE11G. This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with programs 13429 and 13836. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.

References

Appendix A: Bayesian analysis results

We report an example of the results of the Bayesian analysis described in Sect. 3.3 referring to the 2009 data. In Fig. A.1, we report the posterior probability distributions for the model [(partcov*xstar500)*(xstar5000*pow)]. In Fig. A.2, we report the posterior probability distributions for the model [(partcov*xstar500)*(fast32*pow)].

Fig. A.1. Bayesian analysis results for the parameters of the model [(partcov*xstar500)*(xstar5000*pow)] given the 2009 data. The corner plot reports the 1D (histograms) and 2D (credible contours) posterior probability distributions for the model parameters. The top right panels report the posterior probability distributions of the theoretical model (top) and of the model convolved with the instrumental response, with the binned EPIC-pn data overplotted (bottom).

Fig. A.2. Bayesian analysis results for the parameters of the model [(partcov*xstar500)*(fast32*pow)] given the 2009 data. The corner plot reports the 1D (histograms) and 2D (credible contours) posterior probability distribution for the model parameters. The top right panels report the posterior probability distributions of the theoretical model (top) and of the model convolved with the instrumental response, with the binned EPIC-pn data overplotted (bottom).

Appendix B: Normalization of the HST-COS spectra

We show the normalization procedure for the four HST-COS epochs of observation of PG 1126-041 described in Sect. 4. The left panel of Fig. B.1 shows the four HST-COS original spectra before the normalization process. All of them are very similar in their continuum. The right panel of Fig. B.1 shows the second-order polynomial fits to the region of interest where the absorption is present. The black horizontal lines show the regions selected to anchor the fit, which were the same for all spectra: 1602 − 1605 Å, 1612.5 − 1615.5 Å, 1623 − 1624 Å, and 1637.8 − 1638.1 Å. All the spectra in both figures have been smoothed with a boxcar filter with a width of 31 pixels (corresponding to less than 0.4 Å or 75 km s⁻¹) for visualization purposes.

Fig. B.1. Normalization procedure for the C IV region of the HST-COS spectra of PG 1126-041 taken between 2014 and 2015. Left panel: Observed spectra before normalization. Right panel: Second-order polynomial fits to the region of interest where the absorption is present. The black horizontal lines show the regions selected to anchor the fit.

All Tables

All Figures

	Fig. 1. EPIC spectra of PG 1126-041 observed in June 2009. Black circles are pn data, while red open circles are MOS data. The spectra are plotted unfolded against a power-law model with Γ = 2.
In the text

	Fig. 2. XMM-Newton EPIC-pn spectra of PG 1126-041 extracted in the 0.3 − 10 keV band in eight different epochs of observation. Left panel: spectra unfolded against a power law with Γ = 2. Right panel: spectra of each epoch unfolded against a Γ = 2 power-law model are plotted individually over the 2009 spectrum as a reference.
In the text

	Fig. 3. EPIC-pn spectral residuals of PG 1126-041 to the power-law model [pow] (left column) and to the baseline model [(partcov*xstar500)*(pow)] (right column). The top panel reports the residuals of all eight epochs together, while the smaller individual panels correspond to different epochs of observation, as listed in Table 1.
In the text

Fig. 4. Results of the scan of χ2 statistical spaces between 5 and 11 keV with a Gaussian line. Top panels: confidence contours at (from the outermost to the innermost) 68%, 90%, 99%, and 99.9% significance levels for the centroid energy and normalization of a Gaussian emission or absorption line applied to the baseline model [(partcov*xstar500)*pow]. The contours are filled with a color intensity proportional to the Δχ2 represented in the color map to the right of each panel. The three dashed vertical lines mark the rest-frame energy of Fe I, Fe XXV, and Fe XXVI Kα transitions. We note the different y-axes in the 2008A and 2008B panels. Bottom panels: observed spectra (filled circles) and corresponding background (shaded area).

In the text

	Fig. 5. EPIC-pn spectral residuals of PG 1126-041 to the phenomenological model [(partcov*xstar500)*(xstar5000*pow + zgauss1 + zgauss2)] (left column) and to the disk wind model [(partcov*xstar500)*(fast32*pow)] (right column). The top panel reports the residuals of all eight epochs together, while the smaller individual panels correspond to different epochs of observation, as listed in Table 1.
In the text

Fig. 6. Column density of partially covering X-ray absorber versus its covering fraction (top panel) and ionization parameter (bottom panel), derived from the Bayesian analysis of the EPIC-pn spectra in different epochs assuming the phenomenological model (empty squares) or the disk wind model (filled circles). Error bars represent 2σ deviations from the median posterior probability distribution of the parameters.

In the text

Fig. 7. HST-COS spectra o C IV region of PG 1126-041 taken between 2014 and 2015. Top panel: four spectra are plotted overlapped after being normalized and matched at 1595 − 1605 Å to show both changes in emission and absorption. Bottom three panels: zoomed-in view of the 1600 − 1640 Å region, comparing different couples of epochs in velocity scale, after normalization using a second-order polynomial fit. Three velocity systems are marked, as well as three likely C IV doublets in system II. The strong absorption line at 1608 Å is due to Galactic Fe II and was removed prior to the measurements, so it does not appear in the bottom three panels. All spectra in this figure have been smoothed with a boxcar filter with a width of 31 pixels (corresponding to less than 0.4 Å or 75 km s−1) for visualization purposes.

In the text

Fig. 8. Theoretical models showing effects of varying covering fraction Cf and column density NH of a log ξ = 2 absorber partially covering a power-law emission with Γ = 2. Top panel: column density is fixed to NH = 1023 cm−2, and the covering fraction increases from Cf = 90% (green), to 95% (red), to 99% (black); the solid lines represent the total models, while the dashed lines represent the (1 − Cf) fractions of unabsorbed power-law emission. Bottom panel: covering fraction is fixed to Cf = 95% (the 5% fraction of the unabsorbed power law is plotted with a dashed line), and the column density increases from NH = 5 × 1022 cm−2 (green), to NH = 1023 cm−2 (red), to NH = 2 × 1023 cm−2 (black). The normalization of the y-axis is arbitrary.

In the text

Fig. 9. Two different geometries for the interpretation of the blueshifted absorption troughs in the Fe K band of PG 1126-041. A wind with a smaller Rin (corresponding to υ∞ = −0.25c) and observed at a larger inclination angle (θ = 80°) is shown in blue, while a wind with a larger Rin (corresponding to υ∞ = −0.0625c) observed at a smaller inclination angle (θ = 65°) is shown in red. Top panel: sketch of the geometry of the disk wind in the two cases (only Rin has been sketched for simplicity), with streamlines represented by arrows with different thicknesses, proportional to the terminal velocity. Central panel: two theoretical models. Bottom panel: 2009 pn data (in gray) overplotted on the two models folded with the instrumental response.

In the text

	Fig. 10. EW of C IV absorption systems against the column density of the partially covering X-ray absorber in each epoch of coordinated XMM-Newton/HST observations of PG 1126-041. From left to right: system I, system II, and system III (system I was not detected in 2015).
In the text

Fig. A.1. Bayesian analysis results for the parameters of the model [(partcov*xstar500)*(xstar5000*pow)] given the 2009 data. The corner plot reports the 1D (histograms) and 2D (credible contours) posterior probability distributions for the model parameters. The top right panels report the posterior probability distributions of the theoretical model (top) and of the model convolved with the instrumental response, with the binned EPIC-pn data overplotted (bottom).

In the text

Fig. A.2. Bayesian analysis results for the parameters of the model [(partcov*xstar500)*(fast32*pow)] given the 2009 data. The corner plot reports the 1D (histograms) and 2D (credible contours) posterior probability distribution for the model parameters. The top right panels report the posterior probability distributions of the theoretical model (top) and of the model convolved with the instrumental response, with the binned EPIC-pn data overplotted (bottom).

In the text

	Fig. B.1. Normalization procedure for the C IV region of the HST-COS spectra of PG 1126-041 taken between 2014 and 2015. Left panel: Observed spectra before normalization. Right panel: Second-order polynomial fits to the region of interest where the absorption is present. The black horizontal lines show the regions selected to anchor the fit.
In the text