© The Authors 2026

1. Introduction

Active galactic nuclei (AGNs) are primarily powered by accretion onto supermassive black holes (SMBHs) and host some of the most energetic phenomena in the universe. Among the key observational signatures of AGN activity are UFOs, detected as blue-shifted Fe XXV/XXVI absorption lines above 7 keV in X-ray spectra, indicating velocities that can reach a substantial fraction of the speed of light (Pounds et al. 2003; Tombesi et al. 2011). These outflows are thought to originate from the inner regions of the accretion disk, within a few hundred gravitational radii of the black hole (Kraemer et al. 2018). They have been deemed significant for both our understanding of accretion physics and their potential role in AGN feedback (King & Pounds 2015; Gaspari et al. 2020).

Systematic surveys of nearby Seyfert galaxies have revealed that approximately 40% of these AGNs harbor highly ionized UFOs, characterized by velocities exceeding 10 000 km s⁻¹ and averaging between 0.1c and 0.3c (Tombesi et al. 2010; Gofford et al. 2013; Matzeu et al. 2023). Various theoretical models have been proposed to explain the origin and structure of UFOs. These include radiation-driven winds from the accretion disk, magnetically driven outflows, and multiphase cycles such as chaotic cold accretion (CCA), where cold clouds condense from the hot halo and interact with the disk (Gaspari et al. 2013, 2017). Such mechanisms could naturally produce clumpy and transient outflows, rather than smooth and continuous flows (Gaspari & Sądowski 2017; Wittor & Gaspari 2020).

Recent groundbreaking results from the XRISM telescope have revolutionized our understanding of UFOs. High-resolution spectroscopy has revealed that what appeared as single broad absorption features in CCD spectra are, in fact, complex superpositions of multiple narrow absorption lines (Xrism Collaboration 2025; Noda et al. 2025; Mehdipour et al. 2025; Gu et al. 2025; Mizumoto et al. 2025; Reeves et al. 2026). These “forests” of velocity-resolved components indicate that UFOs are highly structured and clumpy, rather than homogeneous in terms of the shape of their flows. However, despite this progress, several key questions remain regarding the launching mechanism, acceleration, and composition of winds, which remain poorly constrained. Detailed time-resolved studies of individual sources are therefore crucial to understand the dynamics and evolution of these powerful phenomena.

2. IRAS 13224-3809

A major advancement in UFO studies was made by Parker et al. (2017) with the identification of a variable, mildly relativistic wind in the bright narrow-line Seyfert 1 (NLS1) galaxy IRAS 13224-3809 (z = 0.0658) based on extensive XMM-Newton observations in 2016 (1.5 Ms). Multiple absorption features were detected, with a significant line at E = 8.6 keV indicating outflow velocities of υ = (0.236 ± 0.006) c. Notably, Fe K absorption features weaken as X-ray luminosity increases, which is interpreted as an overionization of the gas. A direct correlation was observed between UFO velocity and source luminosity, along with flux-dependent trends in equivalent width, ionization state, and column density (Parker et al. 2017; Pinto et al. 2018).

However, Chartas & Canas (2018) raised some important caveats. After reanalyzing with optimized extraction regions to reduce background contamination above 7 keV, they found the flux-dependent trends were artifacts and only the velocity-luminosity correlation survived. This might suggest radiative acceleration (Sądowski & Gaspari 2017), although contrasting results reported in 1H 0707-495 (velocity-luminosity anticorrelation) have complicated the overall interpretation of these results (Xu et al. 2021). The CCD-based conclusions remain highly degenerate and the source’s extreme variability (i.e., a factor of 10 in F_{0.3 − 10 keV} within < 10 ks) demands careful spectral selection.

Despite extensive studies published by Parker et al. (2017), Pinto et al. (2018), Chartas & Canas (2018), Jiang et al. (2018), Jiang et al. (2022), and Midooka et al. (2023), among others, there is still no unified picture of IRAS 13224-3809’s complex X-ray spectrum. Different models have been proposed: the hard spectrum (2–10 keV) has been modeled as either a direct power law or a relativistically reflected continuum; the soft excess (0.3–2 keV) has been attributed to multiple reflection components or to variable partial covering by clumpy absorbers. Meanwhile, the persistent model degeneracy underscores the need for a comprehensive, systematic approach.

This work presents a systematic reanalysis of the full 2016 XMM-Newton (1.5 Ms) and NuSTAR (500 ks) campaign, with careful background treatment and equal-count spectral selections in both flux and time domains. We uniformly applied three alternative spectral models across all intervals, allowing us to track the evolution of UFO parameters robustly. A key novelty of our approach is the ability to measure the wind acceleration in response to X-ray flares, a direct probe of the driving mechanism that has not been systematically explored in previous studies. Here, we focus on the hard X-ray band (E > 3 keV) where UFO features are most prominent and previous claims are based (Parker et al. 2017; Chartas & Canas 2018).

In Section 3, we report on the data reduction. Sections 4 and 5 are dedicated to the spectral modeling, while Sections 7 and 8 present a discussion of the results and our conclusions.

Throughout this paper, unless otherwise specified, errors are at the 68% (∼1σ for Gaussian distributions) confidence level, while the upper and lower limits are at the 90% confidence level. The reported errors on the best-fit parameters were calculated using the likelihood profiling method through the error command in XSPEC (Fig. 1).

Fig. 1. XMM-Newton EPIC-pn light curves of IRAS 13224-3809 in the 3–10 keV band. The plot shows the division of the data into five time intervals (T1–T5, labeled at the top) and five flux intervals (F1–F5, labeled on the left), each selected to contain equal total counts for optimal S/N in subsequent spectral analysis. Vertical dashed lines mark the boundaries of the time intervals, while horizontal dashed lines indicate the flux thresholds. The flux thresholds used are F = 0.037, 0.056, 0.078, 0.114 cts s−1. The bin size of the light curves is 1000 ks.

3. Data reduction and product extractions

We used the full 1.5 Ms (∼17 days) observing campaign with the X-ray Multi-Mirror Mission (XMM-Newton; Strüder et al. (2001)) and 500 ks (almost 6 days) with the Nuclear Spectroscopic Telescope Array (NuSTAR) high-energy X-ray mission (Harrison et al. 2013) taken between July and August 2016. Observation dates, exposure times, and background-subtracted source count rates are listed in Table A.1. All energies are quoted in the source rest frame at z = 0.0658 (Jones et al. 2009). The spectral fitting was performed with the X-ray spectral fitting package XSPEC v12.15.1 (Arnaud 1996). All models included Galactic absorption via TBabs (N_H = 4.76 × 10²⁰ cm⁻²; HI4PI Collaboration 2016). To account for cross-calibration between the XMM-Newton EPIC-pn and the NuSTAR focal plane modules (FPMs), we included a multiplicative constant free to vary between detectors; its value remains within 25% for all fits including NuSTAR data, consistent with calibration uncertainties given the non-simultaneity of the observations and the source variability.

3.1. XMM-Newton data

We limited the analysis to the European Photon Imaging Camera pn (EPIC-pn) data to maximize the S/N in the hard X-ray band, as this instrument has a greater effective area than the two combined MOS detectors between 0.3 and 10 keV. For the reduction, we filtered the EPIC-pn data by selecting events corresponding to instrument PATTERNS in the 0-4 range (single and double pixel events). Several moderate-amplitude background flares were present during several of the XMM-Newton observations and were filtered out by excluding time intervals with rate higher than 0.4 counts/s in the energy range 10–12 keV. To test sensitivity to background non-uniformity, we varied the position of the circular background extraction region (fixed radius of 50 arcseconds) across distances of 100–300 arcseconds from the source on the same CCD chip, finding no differences in spectral shapes or features.

Before proceeding to the flux- and time-resolved extraction of EPIC-pn spectra, we systematically tested different source extraction regions to maximize the S/N in the detector-frame 7–10 keV band, where UFO features are expected. Chartas & Canas (2018) demonstrated that background contamination above ∼7 keV is significant for the larger extraction regions adopted in Parker et al. (2017) and that selecting smaller apertures (250 physical units, or 12.5 arcsec) improves the S/N by a factor of ∼8 by reducing background contribution, while retaining ∼60% of the XMM-Newton point spread function. To reproduce these findings, we tested circular source extraction regions with radii of 30, 25, 20, 15, 12.5, and 10 arcseconds on two representative observations (Obs. ID 0780561601 and 0792180201), comparing the 7–10 keV S/N for each aperture. For the background, we maintained the fixed circular extraction region radius of 50 arcseconds on the same chip as the source. Consistent with Chartas & Canas (2018), we confirmed that a 12.5-arcsecond radius optimizes the hard X-ray S/N and adopted this extraction region for all subsequent analysis. The extraction of the products was performed using version 21.0 of the XMM-Newton Science Analysis System (SAS) and following the XMM-Newton SAS threads. The spectra and light curves were extracted with the evselect task. For each extraction, we generated the ancillary files and response matrices using, respectively, the rmfgen and arfgen tasks.

To obtain an average spectrum with a high S/N, we stacked all EPIC-pn exposures using the epicspeccombine task. In this case, S/N was not a limiting factor, and we binned the spectrum to not oversample the instrumental energy resolution by more than a factor of three. With this criterion alone, we found that the average spectrum already had more than 25 counts per bin in the 3–10 keV band.

To complement the X-ray dataset, we extracted UV data from the Optical Monitor (OM) onboard XMM-Newton, which observed IRAS 13224-3809 in imaging mode using the UVW1 filter. We processed the OM data with the omichain SAS task and filtered the source list to retain only calibrated photometry of the target. The average UVW1 magnitude for the whole 2016 campaign is m_UVW1 = 15.36 with a standard deviation of 0.03 mag. The peak-to-peak variability amplitude is ∼0.17 mag; hence, the source has a stable UV emission varying by only < 2% during the whole XMM-Newton campaign.

3.2. X-ray spectral analysis: time- and flux-resolved approaches

To investigate spectral variability with a high statistical significance, we adopted a strategy based on equal-count selections in both time and flux domains. For the time-resolved analysis, we split the background-subtracted 7–10 keV light curve into five intervals, each containing an equal number of counts. This approach ensured that each time-resolved spectrum had comparable S/N in the iron Kα region, optimizing our ability to study the Fe XXV/XXVI features at high significance. For the flux-resolved analysis, we divided the 3–10 keV light curve into five bands, again using thresholds that yield five spectra with equal total counts. This choice allowed us to explore flux-dependent spectral changes across the broader continuum, minimizing potential biases from absorption features that may affect the higher energy band. We also tested finer selections, including a ten-flux division following Pinto et al. (2018) and a 20-interval time-resolved selection. While these approaches are consistent with our main results, the lower S/N in each spectrum limited the robustness of the conclusions. Therefore, we focused on the five time and five flux selections, which provided a balance between the temporal and flux resolution and the statistical quality, enabling stronger constraints to be placed on spectral variability. To mitigate the impact of statistical noise while preserving sensitivity to narrow spectral features, the spectra were rebinned using the specgroup task within the SAS software package. The binning was performed such that no spectral channel oversamples the instrumental energy resolution by more than a factor of three. In our baseline analysis, we adopt a minimum of five counts per bin. This choice offers an improved sampling of the Fe Kα region and the associated outflow signatures, which would otherwise be excessively smoothed by more aggressive binning. The spectral fitting was carried out using the W statistic in XSPEC, which is appropriate for Poisson-distributed data and remains valid across the full range of count regimes explored here. We verified robustness by testing unbinned spectra and groupings of 5 and 25 counts per bin (maintaining bins ≥1/3 resolution element) and comparing W and χ² statistics. The results remain consistent within uncertainties. We adopted a five-count binning with resolution-based grouping and W statistic, optimally balancing statistical robustness and spectral resolution.

3.3. NuSTAR data

NuSTAR observed IRAS 13224-3809 for a total of ∼600 ks (July-August 2016). The source is very soft and NuSTAR has a smaller effective area than EPIC-pn, so the S/N above ∼10 keV is intrinsically low. To optimize the S/N in the 7–15 keV band we tested source and background apertures (r_a = 40, 35, 30, 25, 20 arcsec) and found that a 25″ source radius (50″ background on the same chip) gave the best results; this radius encloses only an encircled energy fraction of ∼40% (An et al. 2014). The ancillary response file produced by nuproducts includes an encircled energy fraction correction, so the aperture-dependent effective area loss is accounted for in our modeling. Only by stacking all NuSTAR exposures could we reliably recover the source up to ∼20 keV; for the average spectrum we applied coarse binning (8 channels per bin) to ensure ≳25 counts per bin. Given the source variability (factor of ∼20), the time-averaged NuSTAR spectrum is not suitable for flux-resolved fits, so NuSTAR was excluded from interval-resolved analyses and retained only for the stacked average spectrum (with FPM A and B fitted separately and plotted as their average; Figure 3).

4. X-ray average spectrum

We first analyzed the average spectrum from stacked XMM-Newton EPIC-pn and NuSTAR FPM observations (0.3–20 keV), focusing on E > 3 keV where the prominent soft excess is less dominant, enabling clearer assessment of continuum and absorption features. A simple absorbed power-law model (Γ = 2.56 ± 0.02) provides a poor fit to the data, with reduced χ_ν² ∼ 4.4, indicating significant spectral complexity. To systematically identify additional features, we performed a blind Gaussian line scan (Appendix B), which revealed both broad emission and absorption structures across the spectrum (Fig. 2). We tested three distinct spectral models, building in complexity to explore the origin of these features. First, we introduced photoionized absorbers using XSTARv2.58 table models Kallman (1999); see Appendix C for further details). The photoionized gas in the XSTAR simulations is assumed to be optically thin and to have a single velocity component. One absorber layer improved the fit (χ_ν² = 1.53); residuals in the Fe K region motivated a second component, yielding χ_ν² = 1.2. This absorption-only model (model A) is

$$\begin{array}{c}\hfill \mathtt{TBabs}\times \mathtt{XSTAR\_1}\times \mathtt{XSTAR\_}2\times \mathtt{zpowerlaw}.\end{array}$$

Fig. 2. Result of a blind Gaussian line search performed on the average XMM EPIC-pn and NuSTAR FPM spectra, plotted in the 3–12 keV range in the IRAS 132224-3809 rest frame. The line width is fixed to σline = 100 eV. The ordinate axis shows the line normalization. Negative values indicate absorption lines. The purple and orange solid contours represent the 3σ and 5σ confidence levels corresponding to Δχ2 = 11.8 and Δχ2 = 28.7 for the addition of two degrees of freedom, respectively.

Second, we added a broad Gaussian emission line to account for the 6–7 keV excess, yielding χ_ν² = 0.99 (Δχ² = 46.5, 3 d.o.f.) with an equivalent width (EW) of $1.{42}_{-0.38}^{+0.40}$ keV (model B) via

$$\begin{array}{c}\hfill \mathtt{TBabs}\times \mathtt{XSTAR\_1}\times \mathtt{XSTAR\_}2\times (\mathrm{zpowerlaw}+\mathrm{zgauss}).\end{array}$$

Finally, the best overall description was achieved using relxill v2.6 (Dauser et al. 2020), which self-consistently computes intrinsic power-law continuum and relativistic reflection from an ionized Kerr disk. The relxill fit adopts a radial emissivity profile with index α = 3, as expected for a coplanar point source in flat spacetime, and an accretion disk extending from the innermost stable circular orbit to R_out = 400 R_g. We fixed E_cut = 250 keV (typical and unconstrained by NuSTAR). A single absorber layer sufficed (Δχ² = 56.5, 4 d.o.f.; χ_ν² = 0.95), yielding (model C) via

$$\begin{array}{c}\hfill \mathtt{TBabs}\times \mathtt{XSTAR}\times \mathtt{relxill}.\end{array}$$

The model indicates reflection dominance, high iron abundance (A_Fe ∼ 3), and a highly ionized inner disk (log ξ > 3.5 erg cm s⁻¹). We also tested whether a low cutoff (< 100–200 keV) could mimic the spectral downturn, finding no improvement. Best-fit parameters are listed in Table 1. The fitted spectra are given in Fig. 3.

Fig. 3. Average XMM-Newton EPIC-pn (black) and NuSTAR FPM (gray) spectra of IRAS 13224-3809 fitted with model A (absorption only), model B (absorption + Gaussian emission), and model C (relativistic reflection) in the 3–20 keV range. The three panels show the data-to-model ratio for each model, where we removed the absorption components to better visualize their effect on the continuum. The vertical dashed lines indicate the Fe K absorption centroids of each layer. model A prefers lower-ionization absorbers to account for the curvature observed around 7 keV. See Table 1 for details.

5. Flux- and time-resolved spectral analysis

To investigate spectral variability, we applied the same three models tested on the average spectrum to all flux- and time-resolved intervals. As a first step, we performed a blind Gaussian line scan across the spectra to identify significant emission and absorption features, guiding our choice of model components (Figures B.1 and B.2). We also report results from a purely phenomenological model comprising an absorbed power-law continuum plus Gaussian absorption lines with free line widths (Tables D.1 and D.5). Unlike Parker et al. (2017) and Chartas & Canas (2018), we left the line width free to vary to monitor possible changes in UFO broadening. We systematically find broadening values larger than the σ = 0.1 keV assumed in previous works, which inflates the equivalent width (EW) measurements.

For the model selection, we used the Akaike information criterion (AIC; Burnham & Anderson 2002; Liddle 2007). To demonstrate the underlying physical consistency across models, we report the best-fit parameters and statistical improvement of the primary absorber across all intervals. Additional UFO components (e.g., the second absorber in model A) were retained only when supported by strong statistical evidence, requiring ΔAIC > 10 (evidence ratio P₀ = exp(−ΔAIC/2) < 0.006).

For model A (absorption only), we began with a simple absorbed power-law continuum and iteratively added photoionized absorber components using XSTAR table models. Up to two absorbers were required in lower flux intervals, while a single absorber sufficed at higher fluxes (Tables D.2 and D.6).

For models B and C, we tested the addition of either a broad Gaussian emission line or a relativistic reflection component to the absorbed power-law continuum. For model B, we allowed the Gaussian centroid, width, and normalization to vary but imposed a limit of σ ≤ 1.5 keV to maintain fit stability, as broader widths led to unphysical solutions and strong degeneracies with the continuum. This broad Gaussian is intended as a phenomenological proxy for wind emission: Parker et al. (2022) have shown that the 3–10 keV emission feature in IRAS 13224-3809 can be reproduced not only by relativistic reflection but also by P-Cygni profiles from a multicomponent wind. We verified that varying the upper bound over the range σ = 1.0–1.5 keV does not significantly alter the inferred UFO parameters (velocity, ionization, column density). For the relativistic reflection model, we fixed the spin parameter to a = 0.3, the disk inclination to i = 56 degrees, and the iron abundance to A_Fe = 3 as found in the average spectrum fit (1), to reduce parameter degeneracy and improve fit stability since we do not expect these parameters to vary on the timescales of two months. However, the spin parameter is unconstrained in the average spectrum fit, with the error range spanning from 0 to 1; a = 0.3 is simply the face value from the fit, and any choice of spin within this interval does not strongly affect the results. Within 1σ, these parameters are consistent with those reported by Parker et al. (2017) (a = 0.989 ± 0.001, i = 58 ± 1°, and A_Fe = 3.5 ± 0.2). In higher flux intervals for models B and C, spectra were well described by continuum and emission alone. In lower flux intervals, one additional UFO component is still required (see Tables D.3, D.7, D.4, D.8). The absence of further absorbers is reasonable, as the spectral drop above 8 keV is already accounted for by the emission line blue wing or the reflection component.

No single model was statistically preferred across all intervals, with fit statistics remaining comparable (Figure 4). The best-fit models for flux- and time-resolved spectra are shown in Figures 5 and 6, demonstrating that all three approaches successfully describe the spectral variability across different intervals.

Fig. 4. Histograms showing the distribution of reduced W statistic (Wν) values for models A, B, and C across all flux- and time-resolved spectra. Each model demonstrates comparable fit quality, with no single model consistently outperforming the others across the different intervals. This indicates that all three modeling approaches are viable for describing the spectral variability of IRAS 13224-3809 within the statistical uncertainties of the data.

Fig. 5. Best-fit models for the flux-resolved XMM-Newton EPIC-pn spectra of IRAS 13224-3809 in the 3–11 keV range. Data points are rebinned at 3σ for easier visualization. Each panel corresponds to a different flux interval, labeled F1 through F5, with F1 being the lowest flux state and F5 the highest. The vertical dashed lines highlight the evolution of the centroid energy of the primary absorption feature associated with the UFO as phenomenologically modeled with a Gaussian.

Fig. 6. Same as Fig. 5 but for the time-resolved spectra labeled T1 through T5, with T1 being the earliest time interval and T5 the latest.

All the column densities reported in this work for the UFO components have been corrected for special relativistic effects, following the prescription in Luminari et al. (2024) (Eq. 14):

$$\begin{array}{c}\hfill N_{\mathrm{H}}=\frac{1+v/c}{1-v/c}{N}_{\mathrm{H}}^{\mathrm{obs}}.\end{array}$$(1)

This relation allows us to derive the intrinsic N_H from the observed value under the assumption of a radial, uniform outflow along the line of sight, with a single velocity per component, optically thin absorption, and no transverse velocity components or relativistic ray-tracing effects.

6. Model comparison

The results from the different modeling approaches (models A, B, and C) provide complementary insights, although they also highlight some discrepancies and challenges in interpretation. Figure D.3 summarizes the evolution of the best-fit parameters derived from the spectral analysis of all flux- and time-resolved spectra. We adopted the three different models and a purely phenomenological modeling (power law plus Gaussian absorption lines only). In this latter model, we derive outflow velocities from the line centroid energy, assuming it associates with Fe XXVI Kα, as most of our best-fit photoionization modeling indicates (see Appendix C). Model A (absorption only) consistently requires an absorber with high significance (≳3σ) across all time- and flux-resolved intervals, yielding clear trends in the UFO parameters. While the statistical significance of the UFO in models B and C can decrease at higher fluxes where the emission components account for more of the spectral curvature, the inferred physical parameters and specifically the outflow velocities, remain remarkably consistent with those of model A (see Fig. D.3 in the appendix). This indicates that the detection of the wind is not an artifact of the continuum choice, but a robust feature whose significance is simply modulated by the underlying modeling of the Fe K emission. In flux-resolved spectra, although models B and C do not strictly require an absorber in the highest flux bins (ΔAIC < 6), the inclusion of a UFO component in those fits yields velocity estimates that are fully compatible with the clear acceleration trend seen in model A. Importantly, we constructed all time- and flux-resolved intervals to have comparable S/N. This ensured that variations in absorber detection reflected genuine spectral variability rather than statistical fluctuations. Finally, all models recovered the softer-when-brighter behavior of the continuum. In model A, this trend is unambiguous (Γ increasing from ∼1.5 to ∼2.5 with flux; Fig. 7, panel A), while models B and C show a flatter trend, as the emission component favors steeper power-law slopes (Fig. D.3, panel a, b). Due to the small sample size, we assess significance using an exact permutation test, which enumerates all 5! = 120 rank orderings and requires no distributional assumptions. Model A yields p = 0.008 (the minimum attainable value for n = 5), while models B and C give p = 0.06 and p = 0.025, respectively. To test whether the relation would remain non-flat when parameter uncertainties are included, we performed a bootstrap linear regression. Following Laurenti et al. (2026), we considered uncertainties in both the x and y directions. Recomputing the slope over 10 000 iterations, with each point drawn from a uniform distribution within its 1σ bounds, gives a positive slope in more than 99.9% of realizations, strongly disfavoring a flat trend. This behavior likely reflects intrinsic changes in the corona or accretion flow, with higher flux states corresponding to a softer, more efficient Comptonization spectrum (e.g., Sobolewska & Papadakis 2009; Serafinelli et al. 2017).

Fig. 7. Correlations among key spectral and wind properties of IRAS 13224-3809. Left: Photon index Γ versus 3–10 keV luminosity for the flux-resolved spectra (F1–F5), illustrating the softer-when-brighter trend. Center: Outflow velocity of the UFO as a function of 3–10 keV luminosity for the time-resolved spectra (T1–T5), highlighting their coupled evolution. Right: Relation between UFO outflow velocity and Fe Kα line width, tentatively indicating that faster winds are associated with greater velocity dispersion. Shaded regions represent 1σ uncertainties on the best-fit correlations, derived from 10 000 bootstrap regressions sampling the error boxes. Each panel includes the p-value from an exact permutation test assessing null hypotheses of no correlation and the best-fit relation.

Regarding the specific parameters of the emission components, in model C, the fits find an iron abundance A_Fe ∼ 3 times solar and a high reflection fraction, often exceeding 30. The spin parameter remains unconstrained (0 < a < 1), so we fixed a = 0.3 from the average fit to reduce degeneracies. No significant trend with flux or time is evident for the disk ionization, which remains high (log ξ > 3.5 erg cm s⁻¹) across all intervals. For model B, the broad Gaussian emission line requires full widths at half maximum (FWHM) exceeding 2 keV, with the equivalent width remaining high across flux states and the centroid energy showing a modest luminosity dependence.

7. Discussion

7.1. Physical plausibility of the models

The physical plausibility of the emission models remains questionable. Both the phenomenological broad Gaussian and the relativistic reflection model require extreme parameters, such as very high equivalent widths for the emission lines and high reflection fractions. Our modeling confirms previous claims in the literature that relativistic reflection can reproduce the observed spectra. However, we find this interpretation incompatible with the clear UFO parameter trends reported elsewhere (e.g., Parker et al. 2017; Pinto et al. 2018; Jiang et al. 2018). Indeed, both flux- and time-resolved analyses show that when we include relativistic reflection (model C), we detect the UFO only in the lowest flux or earliest time intervals. At higher fluxes or later times, we can adequately model the spectra without an absorber. Therefore, we cannot simultaneously claim a reflection-dominated spectrum and maintain the clear UFO parameter trends.

Previous studies investigating wind evolution have typically not applied a fully systematic relativistic reflection model. Reflection was usually tested only on selected spectra, while the broader wind analysis relied on phenomenological or photoionized absorption models. For example, Pinto et al. (2018) modeled the continuum as a power law plus a Gaussian convolved with a relativistic profile and a photoionized absorber. Chartas & Canas (2018) used a simple absorbed power law with Gaussian lines for the systematic UFO analysis. Similarly, Midooka et al. (2023) used a power law with up to three Voigt absorption lines. We further note that our spectral fitting covers only E > 3 keV. Relativistic reflection imprints diagnostic features across the full X-ray band. Restricting the bandpass reduces sensitivity to disk spin, inner emissivity, and ionization profile, making these parameters harder to constrain and the reflection scenario harder to test conclusively. However, extending the bandpass would not straightforwardly improve constraints on the reflection geometry. Previous broadband analyses of IRAS 13224-3809 (e.g., Pinto et al. 2018; Jiang et al. 2018, 2022) were unable to model the full spectrum with a single reflection component: they required either an additional phenomenological Gaussian or a second independent reflection component to account for the strong soft excess below 2 keV. This additional complexity introduces further degeneracies that limit the constraining power even with broader energy coverage. Model C faces several physical challenges:

• The high iron abundance (A_Fe ∼ 3) is difficult to justify for a local Seyfert galaxy like IRAS 13224-3809, particularly for iron alone rather than for all elements.

• The reflection fraction is often extremely high (greater than 30), which is hard to explain within standard reflection scenarios. To achieve such high reflection fractions, one would need to invoke extreme light-bending effects or highly anisotropic emission geometries, which are not typically expected, given also that polarimetric studies favor an extended corona rather than a very compact one (e.g., Gianolli et al. 2023). Notably, such high reflection fractions are physically inconsistent with the baseline emissivity index α = 3, as they would require much more compact illumination. We tested steeper emissivity profiles (up to α = 9), but found that this does not solve the parameter degeneracies, with the black hole spin remaining unconstrained and the reflection fraction still extremely high.

While previous studies reported tighter spin constraints (e.g., Parker et al. 2017), our systematic analysis across all flux states reveals that CCD data do not robustly determine this parameter (0 < a < 1). The shallow, intrinsically broad line profiles, limited spectral resolution, and the restriction of the modeling to energies E > 3 keV all contribute to the lack of meaningful spin constraints.

Jiang et al. (2022) modeled IRAS 13224-3809 using reflection models with enhanced disk densities, but the same physical challenges persist. Absorption requires an iron abundance of ∼3 solar for iron alone, a second independent reflection component with a different ionization state is needed for the soft excess, and even with enhanced disk densities the spin is constrained to near-maximal values (a > 0.98).

Our application of the relxill model thus shows that adopting a relativistic reflection scenario requires extreme parameters and precludes a coherent characterization of the UFO evolution, consistent with previous reports (Parker et al. 2017; Jiang et al. 2018, 2022).

Model B (including a broad Gaussian emission line) inherits some of the problems of model C concerning the absorption features. In the flux-resolved spectra, the presence of strong emission erodes the need for absorbers, making it difficult to track their behavior. However, in contrast to model C, time-resolved analysis always requires an outflowing component with physical parameters consistent with those derived from model A. As noted in Section 5, the FWHM values for the emission component exceed 2 keV. However, to maintain fit stability, we had to limit the line width to σ ≤ 1.5 keV, since broader widths resulted in strong degeneracies with the continuum, where the power-law normalization becomes negligible. This limit is still large enough to encompass a broad, structured emission component as expected from wind emission. Parker et al. (2022) already demonstrated on IRAS 13224-3809 that the 3–10 keV emission is consistent with P-Cygni profiles from a multicomponent wind, where reproducing the feature using dedicated radiative-transfer simulations required even two wind components. We should therefore read the broad Gaussian in model B as a phenomenological stand-in for this wind emission. The inferred UFO parameters are robust against this choice, as verified over the range σ = 1.0–1.5 keV. A self-consistent treatment using the WINE wind model (Luminari et al. 2018, 2024) will be presented in a forthcoming paper. Despite their nature, models B and C confirm that the data are compatible with an emission structure that accompanies the absorption, as previously shown by Parker et al. (2022).

7.2. Trends for UFO parameters: flux and time dependence

A key result of our time-resolved analysis is the direct measurement of UFO acceleration in response to X-ray flares. The wind rapidly accelerates during luminosity increases, reaching velocities up to ∼0.38c (see top panel of Fig. 9), providing the first direct observational constraint on the wind acceleration timescale and its coupling to the central engine.

Regarding the ionized outflowing absorber, model A consistently shows that the outflow velocity increases with higher fluxes, indicating a direct and dynamic response of the UFO to variations in the central engine. This velocity–luminosity correlation is robustly confirmed in both flux- and time-resolved analyses, demonstrating it is intrinsic to the source and not an artifact of the selection procedure (see top panel of Fig. 9, center panel of Fig. 7). Model B independently confirms the velocity–luminosity trend and the wind acceleration. It recovers a significant absorber (ΔAIC > 6) in all five time-resolved intervals and in the three lowest flux intervals (F1–F3), with outflow velocities comparable with the ones inferred from model A (Tables D.7 and D.3). Where model C detects the absorber (T1–T3 and F1–F3), its inferred velocities agree with models A and B within 3σ, consistently tracing the same velocity increase with luminosity. In the highest flux states (T4, T5, F4, F5), the reflection component in model C introduces curvature in the iron K-band that absorbs the absorption residuals, which is a limitation at high flux rather than a contradiction of the velocity–luminosity trend seen robustly in the other two models. Although the statistical significance of the UFO in models B and C drops below the ΔAIC = 6 ≈ 2σ threshold in some of these high flux intervals, the physical outflow velocities inferred from these models remain consistent with model A within their statistical uncertainties. This mutual agreement across all three models (Fig. D.3, panels c, d) directly mitigates concerns of model degeneracy, demonstrating that our main physical conclusions about the observation of wind acceleration and the velocity-luminosity relation, are robust regardless of the baseline continuum model. Previous studies (Pinto et al. 2018; Chartas & Canas 2018) reported the velocity-luminosity correlation, and our analysis strengthens this evidence by confirming it simultaneously in both flux and time domains, across independent spectral prescriptions, and by directly measuring the wind acceleration. Applying the same exact permutation test described above, we obtain p = 0.025 for model A and p = 0.008 for model B in the time-resolved intervals, supporting the correlation across independent spectral prescriptions. Furthermore, bootstrap linear regression that includes the measured parameter uncertainties gives a positive slope in more than 99.9% of realizations, yielding for Model A the relation log(v_out/c) = (0.9 ± 0.2)log(L_3 − 10)− (39 ± 9).

Regarding the anticorrelation between EW and 3–10 keV luminosity reported by Parker et al. (2017) and Pinto et al. (2018) but not confirmed by Chartas & Canas (2018), we find no significant trend; the EW is consistent with being constant across all flux and time intervals. This result is consistent with the line scans shown in Fig. B.1 and Fig. B.2, where the line normalization is plotted as a function of line energy. In lower flux states, the absorption lines appear more tightly localized, whereas in higher flux states the line detection region becomes more extended. As a consequence, the EW remains approximately constant across the intervals, although the detection of the line is more robust in lower flux states. Allowing the line width, σ, to vary confirms that the UFO features are systematically broader than 0.1 keV, in agreement with the turbulent velocities inferred from the XSTAR modeling, which are often in excess of 10 000 km s⁻¹. In the appendix (Fig. D.1), we directly demonstrate the absence of a correlation between EW and luminosity in the time-resolved spectra, and compare our results with the best-fit correlation reported by Parker et al. (2017). The limited spectral resolution of the EPIC-pn data and the large observed equivalent widths limit our ability to fully disentangle the UFO physical properties, particularly between the column density and ionization. However, this degeneracy is not absolute: the ionization parameter is reasonably constrained across intervals (log ξ ≈ 3.3–4.4), and for log ξ ≳ 3.6 erg cm/s the 3–10 keV opacity is dominated by H-like iron (Fe XXVI), which reduces ambiguities in line identification (see Appendix C).

In light of recent XRISM micro-calorimeter results (e.g., Xrism Collaboration 2025), broad absorption troughs seen with CCDs can be resolved into multiple narrow, kinematically distinct components. Therefore, our single-zone fits should be read as effective, ensemble descriptions: the measured v_out acts as a centroid-weighted velocity tracing the global wind kinematics, while the fitted log ξ and N_H represent averaged phases of a structured, multiphase absorber rather than a unique homogeneous medium. In contrast, N_H measurements carry large uncertainties. If we fix one parameter to its average value (e.g., log ξ = 3.4 erg cm/s or log N_H = 23.7 cm⁻²), the other shows a clearer trend with flux (Figure 8); when both are left free, the large statistical errors make the data consistent with a weak trend or with no trend. Several N_H and v_turb estimates reach the boundaries of our model grids, as indicated by limits in the relevant tables, reducing the interpretability of trends for these variables, although the outflow velocity remains robustly constrained. Importantly, v_out is robustly measured since it primarily sets the line centroid, and does not display the same degeneracy. If N_H and log ξ were fully degenerate, we would expect a correspondingly larger uncertainty in v_out, which we do not observe. Thus, while our CCD-based errors are large, our results are broadly consistent with previously reported trends (e.g., Pinto et al. 2018) but do not permit equally tight independent constraints on both N_H and log ξ.

Fig. 8. Ionization parameter (log ξ) of the UFO and column density (NH) corrected for special relativistic effects for the flux-resolved spectra (F1-F5) of IRAS 13224-3809. Here we are directly visualizing the trends when fixing one of the two parameters to its average value across the intervals.

This complex behavior finds a qualitative physical interpretation within the framework of chaotic cold accretion (CCA; e.g., Gaspari et al. 2013, 2017). In this scenario, the multiphase rain and cloud structure is inherently stochastic and variable, which naturally explains the (i) large scatter and non-monotonicity observed in ξ and N_H and (ii) the episodic appearance of the absorption features, as clouds condense and dissolve or pass transiently through the line of sight. Multiphase wind structure is also predicted by complementary physical mechanisms. Radiation-hydrodynamic simulations show that super-Eddington flows fragment into clumpy structures via Rayleigh-Taylor instabilities on hour-to-day timescales (Takeuchi et al. 2013), consistent with our observed variability. Similarly, thermal instability rapidly condenses initially smooth outflows into multiphase clumps when plasma cools below the equilibrium curve, with growth times comparable to our observed acceleration timescale (Dannen et al. 2020; Waters et al. 2021). These mechanisms naturally produce the clumpy substrate upon which magnetic acceleration operates, potentially explaining both the structured absorption features and the rapid velocity response to X-ray flares we observe.

7.3. Degeneracies and modeling challenges

The application of model A sometimes results in the detection of a second UFO component (v_out ∼ 0.4c) in low-flux intervals (either flux- or time-selected). While multiple components are plausible (e.g., MCG-03-58-007 Braito et al. 2021 or PDS 456; Xrism Collaboration 2025), this is model-dependent: in models B and C, steeper continua allow Fe Kβ or emission wings to explain these features without additional absorbers. Our data cannot robustly constrain multiple components. When considering the average spectrum, model A gives a poorer fit (χ_ν² ∼ 1.2) than models B and C (χ_ν² ∼ 0.95–0.99), suggesting that the stacked spectrum contains genuine broad emission features near 6–7 keV that cannot be reproduced by absorption models alone. In model A, the fit compensates by lowering the ionization parameter (log ξ ∼ 3.1 erg cm/s in the average, versus ∼3.6 erg cm/s in the resolved intervals), enhancing lower-ionization transitions that broaden the continuum curvature. A similar approach was adopted by Midooka et al. (2023), though it does not adequately reproduce the emission-like feature near 6–7 keV.

To check whether stacking spectra with different absorption properties could artificially produce emission features, we simulated and stacked synthetic spectra using the observed trends in continuum slope and UFO velocity. The resulting stacked spectrum showed only smeared absorption with no significant emission, confirming that the emission seen in the average spectrum is intrinsic (see Figure D.2). We cannot rule out a contribution from variability below 4 keV, where the soft spectrum is not well understood.

Model differences appear above 10 keV: models A and B predict a flatter spectrum; model C shows a Compton hump at 20–30 keV. However, NuSTAR lacks sufficient S/N above 10 keV, and the Swift Burst Alert Telescope is insensitive to the flux predicted by model A (F_{14 − 195 keV} ∼ 6 × 10⁻¹³ erg s⁻¹ cm⁻², well below its sensitivity limit). Current hard X-ray data cannot discriminate.

In summary: model A (absorption only) consistently reveals clear trends in the UFO parameters across all intervals, but provides a poorer fit to the average spectrum and cannot reproduce the observed broad emission features. Model B (absorption plus broad Gaussian emission) also tracks the UFO parameters well, but requires an extremely broad emission line (FWHM ≳ 2 keV). Model C (relativistic reflection) fits the average spectrum best, but requires extreme parameters such as very high iron abundance, reflection fractions, and broad emission lines, and in high-flux states, often does not require significant absorption, weakening constraints on UFO variability. We caution that these conclusions rest on fitting the E > 3 keV band alone; broadband coverage would provide stronger leverage on the reflection geometry and on the degeneracy between reflection and absorption. This highlights the need for more advanced models (e.g., those incorporating clumpy, multiphase outflows or wind emission features) and higher-sensitivity broad-band data from future missions to robustly disentangle absorption and emission components in AGN spectra.

7.4. Width of UFO absorption features and implications for UFO structure

A persistent feature in our analysis, regardless of the specific model employed, is the significant width of the absorption features. When left free to vary, the turbulent velocity (v_turb) consistently converges to values exceeding 10 000 km s⁻¹, and often hits the upper limit of our model grid (30 000 km s⁻¹). This is consistent with previous studies of this source, which required similarly large broadening to fit the data. In the context of recent high-resolution X-ray spectroscopy, such extreme broadening requires careful interpretation. Groundbreaking results from the XRISM satellite, starting with PDS 456 (Xrism Collaboration 2025; Noda et al. 2025; Mizumoto et al. 2025), have revealed that broad absorption troughs observed in CCD spectra may not be due to single, intrinsically broad lines. Instead, they are resolved into complex “forests” of narrower lines (σ ∼ 1000 km s⁻¹) originating from multiple clumps with distinct line-of-sight velocities. It is highly probable that the “broad” UFO we observe in IRAS 13224-3809 follows this pattern. The large effective width likely represents the velocity dispersion of a clumpy, multiphase outflow rather than turbulence within a homogeneous gas. If we consider the velocity range implied by our best-fit widths, the absorption trough spans from approximately 0.2c to 0.4c. This suggests we are observing a superposition of kinematic components ranging from slower, perhaps newly launched gas, to faster, fully accelerated clumps. This interpretation aligns with the rapid variability we observe: as the central source flares, we may be seeing the sudden acceleration of specific clumps or the ionization of new velocity components within this broad kinematic structure. High-resolution micro-calorimeter spectroscopy with XRISM and NewAthena will be crucial for resolving the wind substructure.

We note a tentative positive correlation between the outflow velocity and the line width (see panel C in Fig. 7), where higher velocities correspond to broader features. Its exact permutation significance is low (p = 0.1), and therefore this trend should be treated with caution. A bootstrap linear regression that includes the measured uncertainties yields a positive slope in most realizations, with log(v_out/c) = (0.5 ± 0.2)log σ_Fe + (0.0 ± 0.1). This is consistent with the idea that faster winds are more structured, with a larger velocity dispersion among the absorbing material. A similar but even weaker trend links velocity dispersion and X-ray continuum luminosity. A positive correlation between outflow velocity and line width, as hinted in our time-resolved spectra, is found in larger AGN samples. Laurenti et al. (2026) recently report that UFOs with higher velocities tend to show broader absorption features, based on a systematic analysis of over 100 AGN outflows. This trend may indicate that faster winds are more structured or clumpy, with a larger velocity dispersion among the absorbing material. Our results for IRAS 13224-3809 are consistent with this broader picture, but given the limited number of data points and the low spectral resolution of EPIC-pn, this trend should be considered as tentative.

7.5. Bolometric luminosity

To estimate the total radiation output, we complemented the X-ray luminosity with a UV measurement that is robust against variability. Using a typical X-ray luminosity (log L_X ∼ 42.8 erg s⁻¹) and the population-averaged bolometric correction (K_X ∼ 15; Duras et al. 2020), we obtained log(L_bol,X/erg s⁻¹) ≈ 44.1. We then used the XMM-Newton OM UVW1 filter (λ_eff ≈ 2910 Å) photometry (Vagnetti et al. 2010; Laurenti et al. 2024). We corrected the flux density for Galactic extinction (Fitzpatrick 1999) and then extrapolated the monochromatic luminosity to rest-frame 3000 Å assuming a canonical UV spectral slope of α_ν = −0.5 (f_ν ∝ ν^α) (Richards et al. 2006). The resulting isotropic luminosity (L_iso) was calculated using the log-linear bolometric correction provided by Runnoe et al. (2012) and expressed as

$$\begin{array}{c}\hfill log(L_{\mathrm{iso}}/\text{erg}{\text{s}}^{-1})=1.85+0.97log(\lambda L_{\lambda ,3000}/\text{erg}{\text{s}}^{-1}).\end{array}$$(2)

Finally, we applied the recommended anisotropy correction factor of 0.75 (Runnoe et al. 2012) to account for the geometric bias of the accretion disk viewing angle, yielding our final, sustained bolometric luminosity estimate of log(L_bol/erg s⁻¹) ≈ 44.38, or L_bol ∼ 2.4 × 10⁴⁴ erg s⁻¹. This value is adopted for all subsequent outflow momentum and power calculations. We note that, within the typical bolometric correction uncertainty of ∼0.37 dex (Duras et al. 2020), the UV-based bolometric luminosity is consistent with the X-ray-derived estimate.

7.6. UFO acceleration and energetics

Using the average values from the time-resolved analysis with model A, we estimate the energetics of the UFO in IRAS 13224-3809. For these calculations, we adopt a black hole mass of M_BH = 10⁷ M_⊙ as estimated by Emmanoulopoulos et al. (2014) from X-ray reverberation measurements. Adopting typical parameters from our analysis (v_out ∼ 0.3c, N_H ∼ 10²⁴ cm⁻², log ξ ∼ 3.7) and the relation in Tombesi et al. (2013), we compute the minimum launching radius assuming v_out = v_esc,

$$\begin{array}{c}\hfill r_{\min }=2G{M}_{\mathrm{BH}}/v_{\mathrm{out}}^2\sim 3.3\times {10}^{13}\text{cm}\sim 20r_g.\end{array}$$

Then the mass outflow rate is

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{out}}=4\pi C_f\mu m_p{r}_{\mathrm{wind}}{N}_{\mathrm{H}}{v}_{\mathrm{out}}\sim 0.42C_f{M}_{\odot }/\text{yr}\end{array}$$

where C_f is the covering factor of the wind, μ ∼ 1.4 is the mean atomic mass per proton, and m_p is the proton mass and we used the minimum radius r_min as an estimate for r_wind. Since both Ṁ_out and Ė_kin are linear in C_f and r_wind, adopting r_min and a covering factor of C_f = 0.3, as representative of UFO population studies (Tombesi et al. 2011; Matzeu et al. 2023), yields robust estimates for the energetic quantities. With these choices, we have Ṁ_out ∼ 0.13 M_⊙/yr and a kinetic power ${\dot{E}}_{\mathrm{kin}}=\frac{1}{2}{\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}^2\sim 1.1\times {10}^{44}$ erg/s. We note that in some time-evolving studies where the absorber distance can be measured (e.g., Serafinelli et al. 2025), the wind velocity falls below the local escape velocity, indicating a failed wind destined to fall back. The time evolution of the UFO velocity in IRAS 13224-3809 also provides valuable insight into the acceleration mechanism at play. Using the time-resolved spectra, we estimate the wind acceleration by computing the velocity difference between consecutive time intervals divided by the time separation between their mean observing times. The outflow velocities are taken from model A. Model B recovers the absorber with ΔAIC > 6 in all cases and the inferred outflow velocity is consistent with model A within the uncertainties (Table D.7). The acceleration measurement is therefore not contingent on a flat power law continuum assumption of model A.

The outflow responds promptly to changes in the X-ray continuum, with a velocity increase of Δv ≈ 0.1c from interval T3 to T4 on a timescale of about 6 days (Δt ≈ 590 ks), yielding an acceleration of ∼50 m s⁻² (Fig. 9, middle panel); model B yields a lower estimate of ∼30 m s⁻², still a significant increase. This rapid coupling to the X-ray luminosity suggests the wind is highly sensitive to the physical processes powering the corona.

Fig. 9. Three-panel summary of the time-resolved properties of the UFO in IRAS 13224-3809, fitted with model A (absorption only). Top: 3–10 keV power-law luminosity and UFO outflow velocity as a function of time, showing their coupled evolution. Middle: UFO acceleration, derived from velocity changes between intervals, highlighting rapid increases during bursts. Bottom: UFO kinetic power expressed as a percentage of the adopted bolometric luminosity; the right axis shows the wind momentum rate normalized to the X-ray radiation momentum rate. The left axis traces how energetic the wind is, while the right axis traces how difficult it is for the instantaneous X-ray radiation field alone to account for that wind momentum. In all panels, the background displays the 3–10 keV light curve for context.

To quantify the contribution of radiation pressure to the UFO acceleration, we estimated the acceleration imparted to a column of gas according to Gu et al. (2025):

$$\begin{array}{c}\hfill a_{\mathrm{rad}}=\frac{\int F(E)[1-T(E)]dE}{c\mu m_p{N}_{\mathrm{H}}},\end{array}$$(3)

where F(E) is the incident X-ray flux at the wind location, T(E) is the transmission as a function of energy, c is the speed of light, μ is the mean atomic mass per proton, m_p is the proton mass, and N_H is the column density. We consider the observed 3–10 keV X-ray luminosity and, using the average N_H and ξ values (accounting for their uncertainties), extract the multiplicative opacity spectrum from our XSTAR absorption model. By integrating this spectrum over the 3–10 keV band, we obtain the total band-averaged transmission, yielding T = 0.89. This value is representative of all time intervals given the uncertainties and degeneracy between ξ and N_H (log ξ ≃ 3.7 erg cm/s, log N_H ≃ 24 cm⁻²). The choice of energy range is not restrictive, as the wind imprints are confined to the 3–10 keV band, with negligible impact outside this interval. Moreover, while the X-ray luminosity shows strong variability, the bolometric luminosity, dominated by the optical emission, exhibits no significant changes on these timescales. This indicates that the UFO responds specifically to variations in the inner disk/corona rather than to global changes in the accretion disk luminosity. We place the wind at the minimum launching radius r_min, representing the most favorable configuration for radiation driving. Under these assumptions, we obtain a_rad ∼ 10 m s⁻², which is a factor of ∼3 − 5 below the observed peak acceleration. This estimate should be regarded as an optimistic upper limit, as it assumes high column densities and the smallest possible wind radius. Even in this extreme scenario, radiation pressure alone cannot account for the prompt and strong acceleration of the UFO. Instead, the immediate response and high acceleration favor magnetic processes, such as magnetic reconnection or magnetically driven winds, as the dominant mechanism, as recently proposed for NGC 3783 and NGC 4051 using XRISM Resolve observations (Gu et al. 2025; Reeves et al. 2026). In this scenario, reconnection events in the inner accretion disk and/or corona can rapidly release stored magnetic energy, simultaneously producing X-ray flares and accelerating UFOs. The reconnection timescale is set by the Alfvén speed, which can approach relativistic values in AGN coronae, allowing for nearly instantaneous wind acceleration. This naturally explains the tight correlation between continuum flares and UFO velocity variations observed in IRAS 13224-3809, in analogy with solar coronal mass ejections. Our results therefore favor a magnetic origin for the UFO acceleration.

To assess the energetics and feedback potential of the UFO in IRAS 13224−3809, we computed its momentum rate and kinetic power using our estimates. To minimize the impact of parameter uncertainties, we again adopted representative average values for the ionization parameter and column density, which dominate the systematic uncertainties in the energetics. The UFO momentum rate and kinetic power were computed as ṗ_wind = Ṁ_outv_out and ${\dot{E}}_{\mathrm{kin}}=\frac{1}{2}{\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}^2$, respectively. Adopting a covering factor of C_f = 0.3 and assuming the UFO is launched at the escape radius r_min, the UFO momentum rate exceeds the radiation momentum rate ṗ_rad = L_X/c, with ${\dot{p}}_{\mathrm{wind}}^{\min }\gtrsim 60{\dot{p}}_{\mathrm{rad}}$ (bottom panel of Fig. 9). Under these assumptions, this momentum excess disfavors single-scattering radiation pressure as the dominant driving mechanism and instead points to either highly efficient multiple scattering or an additional driving agent, such as magnetic stresses. The detailed trends in the bottom panel are nevertheless model dependent, since they rely on the conservative choice R = r_min, which should be regarded as providing lower-bound energetics rather than a unique wind geometry. For the momentum ratio we use the 3–10 keV X-ray luminosity to compute ṗ_rad = L_X/c, since the UV-based bolometric luminosity remains constant throughout the campaign and cannot track the observed variability. The UFO responds directly to changes in the X-ray continuum, so this choice allows us to follow the evolution of ṗ_rad,X and compare it meaningfully to the UFO momentum rate. By contrast, the kinetic-power fraction in the bottom-left panel is normalized by the constant UV-based bolometric luminosity adopted in Section 7.5. Using the constant UV-based bolometric luminosity would yield only a single, static value ṗ_wind ≃ 3 ṗ_rad,UV; even then, the UFO momentum rate is significantly larger than the radiation momentum rate. The kinetic power of the outflow reaches Ė_kin ∼ 1.1 × 10⁴⁴ erg s⁻¹, corresponding to ∼45% of the bolometric luminosity, as shown in the bottom panel of Fig. 9. This value comfortably exceeds the canonical ∼0.5–5% L_bol threshold required for effective AGN feedback (King & Pounds 2015), implying that the UFO in IRAS 13224−3809 is capable of exerting substantial influence on its host galaxy (Gaspari et al. 2020).

8. Conclusions

We present a comprehensive time- and flux-resolved spectral analysis of the narrow-line Seyfert 1 galaxy IRAS 13224–3809, based on the 2016 XMM-Newton (1.5 Ms) and NuSTAR (500 ks) monitoring campaign. We systematically explored three spectral frameworks: photoionized absorption (model A), phenomenological emission (model B), and relativistic reflection (model C). Based on this work, we have drawn the following key conclusions.

1. Model comparison and physical interpretation: All three spectral models provide statistically comparable fits and independently confirm the presence of a highly ionized, variable UFO. However, they differ substantially in physical plausibility. Model A (absorption-dominated) yields the most coherent and physically interpretable trends for the wind properties across both flux- and time-resolved spectra. Model B tracks wind parameters with similar robustness. The inclusion of a broad Gaussian emission component improves the fit to the average spectrum, serving as a phenomenological proxy for UFO emission. However, in flux- and time-resolved spectra, this component requires an extreme width that becomes degenerate with the continuum slope. This degeneracy challenges the physical interpretation of the model and suggests that a single Gaussian does not adequately capture the complex emission structure of the source. In contrast, model C (relativistic reflection) requires extreme and challenging parameters (reflection fractions of R ≳ 30 and iron abundances of A_Fe ≳ 3 A_Fe, ⊙), and the dominance of reflection features suppresses the detectability of the UFO at high fluxes.

2. Velocity-luminosity correlation: We robustly confirm a positive correlation between the UFO velocity (v_out) and the X-ray luminosity (L_X). The outflow responds dynamically to changes in the source brightness, with velocities increasing from ∼0.25c in low-flux states to ∼0.37c during flaring episodes.

3. Line width and equivalent width behavior: No significant anticorrelation is found between the EW of the absorption features and the source flux; the EW remains consistent with being constant across all intervals. The UFO absorption features are intrinsically broad (σ ≳ 300 eV), implying turbulent velocities of v_turb ≳ 10⁴ km s⁻¹. We find tentative evidence for a positive correlation between outflow velocity and line width, suggesting that faster winds are increasingly kinematically complex, possibly reflecting the velocity dispersion within a clumpy, multiphase outflow consistent with CCA predictions.

4. Acceleration mechanism: The UFO reacts rapidly to X-ray flares, reaching accelerations of up to ∼50 m s⁻². Radiation pressure alone cannot account for this behavior, as the inferred force exceeds the expected radiative force by a factor of ∼5. This strongly favors a magnetic driving mechanism, potentially associated with magnetic reconnection events in the inner accretion disk, analogous to solar coronal mass ejections. Importantly, this conclusion is robust against the choice of spectral model used to infer the outflow velocity.

5. Energetics and feedback implications: Adopting a representative covering factor (C_f = 0.3) and the minimum launching radius (r = r_min), the kinetic luminosity reaches Ė_kin ∼ 10⁴⁴ erg s⁻¹, corresponding to nearly half of the bolometric luminosity. The UFO momentum rate significantly exceeds the radiation momentum rate (ṗ_wind ≳ 60 ṗ_rad), demonstrating that the UFO in IRAS 13224–3809 is energetically capable of driving substantial AGN feedback onto its host galaxy.

The combination of near-Eddington kinetic power, large momentum flux, and compact launch radius places this source among the most powerful known AGN winds, comparable to those observed in luminous quasars such as PDS 456 (e.g., Xrism Collaboration 2025). Future works will be aimed at resolving the remaining spectral degeneracies by applying the WINE model, which self-consistently accounts for both emission and absorption from the outflow. This approach promises a more unified picture of the wind structure and dynamics.

Acknowledgments

FT, EP, and ML acknowledge funding from the European Union – Next Generation EU, PRIN/MUR 2022 (2022K9N5B4). MG acknowledges support from the ERC Consolidator Grant BlackHoleWeather (101086804). RM acknowledges financial support from the INAF Scientific Directorate. RS acknowledges funding from the CAS-ANID grant No. CAS220016. LZ acknowledges financial support from the Bando Ricerca Fondamentale INAF 2022 Large Grant “Toward an holistic view of the Titans: multiband observations of z > 6 QSOs powered by greedy supermassive black holes” and from the European Union – Next Generation EU, PRIN/MUR 2022 2022TKPB2P – BIG-z. PC thanks Flor Arevalo Gonzalez and Vicente Madurga Favieres for useful discussions. AI tools (Gemini and Claude) were employed for language refinement and technical code assistance. The integrity of the research, including all analysis and results, was ensured through independent verification by the authors. Figures in this work were produced using the Pubplotlib Python package.

References

Appendix A: Observations and data availability

The data used in this work are publicly available archival data from the XMM-Newton and NuSTAR observatories. The list of observations used in this work is shown in Table A.1. In figure A.1 we show the XMM-Newton Optical Monitor (OM) UVW1 light curve of IRAS 13224-3809 during the 2016 campaign. All XMM-Newton EPIC-pn observations were performed in Large Window Mode.

Fig. A.1. XMM-Newton Optical Monitor (OM) UVW1 light curve of IRAS 13224-3809 during the 2016 campaign. The source shows a stable UV emission with a peak-to-peak variability amplitude of only ∼0.17 mag (< 2%) over the entire observation period. The average magnitude is mUVW1 = 15.36 with a standard deviation of 0.03 mag.

Appendix B: Gaussian line scan and significance estimation

To identify line features above the continuum across a broad energy range, we conducted a blind search using an absorbed power-law continuum model combined with a Gaussian line component, which could be either in emission or absorption. Line energies greater than 3 keV were scanned in the source rest frame, with 10 eV increments. At each step, the normalization of the line was varied, and any resulting statistical improvement in the fit was noted. The line width was fixed to σ_line = 100 eV since we know from previous analyses that we are dealing with broad features and from the fact we are hitting the instrumental resolution limit. To draw the significance contours of the detected features, we used the likelihood ratio test. Alternatively, significance can be estimated using the F-test, as commonly done in previous UFO studies, or information criteria such as AIC; we tested all these approaches and all yielded consistent results. To remain test-independent, we generally report an improvement in the fit.

Fig. B.1. Gaussian line scans performed on the flux-resolved spectra of IRAS 13224–3809 in the 3–10 keV range in the source rest frame. The line width is fixed to σline = 100 eV, and the ordinate axis shows the best-fitting line normalization, with negative values indicating absorption features. The colored contours indicate thresholds in ΔAIC relative to strong and very strong evidence for the inclusion of an additional line, corresponding to ΔAIC = 10 (purple) and ΔAIC = 20 (orange), respectively. These conservative ΔAIC thresholds are adopted to robustly highlight only the most significant structures; the same regions are also recovered when using likelihood-based Δχ2 (ΔW) contours commonly adopted in previous analyses.

Fig. B.2. Same as Fig. B.1 but for the time-resolved spectra.

Appendix C: Photoionization model

We used the XSTARv2.58 photoionization code (Kallman 1999) to generate a grid of photoionization absorption spectra, aiming to model the UFO in IRAS 13224-3809. The code requires an ionizing spectral energy distribution (SED) and parameters describing the physical properties of the absorbing gas. For our SED, we adopted a simple power law with a photon index of Γ = 2 to remain generic, since we could favor any of the models A, B, or C based on statistical grounds. To test the impact of a steeper SED, as expected in high accretion regimes, we also generated a grid with Γ = 2.5. We found no significant differences in the resulting absorption features, except for a systematic shift towards higher ionization parameter values (ξ), as anticipated. This shift arises because a steeper power law requires a higher overall luminosity to maintain the same ionizing flux in the hard X-ray band, thereby increasing the ionization parameter. The ionic abundances of Fe XXIV, Fe XXV, and Fe XXVI, those most relevant to the observed absorption, are similarly affected, with the steeper SED causing only a systematic shift in ξ that in our case is on the order of an increase of log ξ = 1.1 (see Luminari et al. (2024) for a deeper discussion about the effect of the SED on X-ray winds that confirms our conclusion for the high ionization potential features in UFOs). In our model, the column density, N_H, ranges logarithmically from 10²² to 10²⁵ cm⁻², distributed over 10 steps, while log ξ was varied between 3.0 and 6.0 with steps of 0.2 dex, and the turbulent ranged from 1000 to 30000 km/s with 8 steps. The default solar abundances from XSTAR were applied for the elements. Around log ξ ≈ 3.5, which corresponds to the average value measured for the UFO, the iron ions are distributed with approximately 40% as Fe XXV and 40% as Fe XXVI, while the rest are fully ionized. Below this threshold, the opacity is dominated by Fe XXV; above it, it is dominated by Fe XXVI. At log ξ ∼ 4.6, the gas is almost completely ionized.

All the reported values of column densities in this work are corrected for special relativistic effects using the prescription described in Luminari et al. (2024).

Appendix D: Additional materials

Fig. D.1. EW of the UFO as a function of the 3-10 keV observed flux for the time-resolved spectra (T1-T5) of IRAS 13224-3809. The lack of correlation between the two quantities is directly visualized here, in contrast to the work of Parker et al. (2017), Pinto et al. (2018). The dashed line represents the best-fit anticorrelation reported in Parker et al. (2017) for comparison.

Fig. D.2. Simulated average spectrum obtained by stacking synthetic spectra generated for each flux-resolved interval, incorporating observed trends in continuum slope, outflow velocity and the tentative ones in ionization and column density. The resulting stacked spectrum (black) is compared to absorbed power laws for each simulated flux interval (dashed gray). The absence of significant emission features in the stacked spectrum indicates that the emission observed in the actual average spectrum is unlikely to be an artifact of stacking variable flat spectra.

Fig. D.3. Visualization of the key findings in the flux- and time-resolved analysis. Missing values are reported as × and upper–lower limits as triangles. Panels c) and d) show the outflow velocities per each flux- and time-resolved spectrum, respectively showing clearly the consistency of vout across the models. For better visibility, values from each model are offset by a small amount in the x-axis.

All Tables

All Figures

Fig. 1. XMM-Newton EPIC-pn light curves of IRAS 13224-3809 in the 3–10 keV band. The plot shows the division of the data into five time intervals (T1–T5, labeled at the top) and five flux intervals (F1–F5, labeled on the left), each selected to contain equal total counts for optimal S/N in subsequent spectral analysis. Vertical dashed lines mark the boundaries of the time intervals, while horizontal dashed lines indicate the flux thresholds. The flux thresholds used are F = 0.037, 0.056, 0.078, 0.114 cts s−1. The bin size of the light curves is 1000 ks.

In the text

Fig. 2. Result of a blind Gaussian line search performed on the average XMM EPIC-pn and NuSTAR FPM spectra, plotted in the 3–12 keV range in the IRAS 132224-3809 rest frame. The line width is fixed to σline = 100 eV. The ordinate axis shows the line normalization. Negative values indicate absorption lines. The purple and orange solid contours represent the 3σ and 5σ confidence levels corresponding to Δχ2 = 11.8 and Δχ2 = 28.7 for the addition of two degrees of freedom, respectively.

In the text

Fig. 3. Average XMM-Newton EPIC-pn (black) and NuSTAR FPM (gray) spectra of IRAS 13224-3809 fitted with model A (absorption only), model B (absorption + Gaussian emission), and model C (relativistic reflection) in the 3–20 keV range. The three panels show the data-to-model ratio for each model, where we removed the absorption components to better visualize their effect on the continuum. The vertical dashed lines indicate the Fe K absorption centroids of each layer. model A prefers lower-ionization absorbers to account for the curvature observed around 7 keV. See Table 1 for details.

In the text

Fig. 4. Histograms showing the distribution of reduced W statistic (Wν) values for models A, B, and C across all flux- and time-resolved spectra. Each model demonstrates comparable fit quality, with no single model consistently outperforming the others across the different intervals. This indicates that all three modeling approaches are viable for describing the spectral variability of IRAS 13224-3809 within the statistical uncertainties of the data.

In the text

Fig. 5. Best-fit models for the flux-resolved XMM-Newton EPIC-pn spectra of IRAS 13224-3809 in the 3–11 keV range. Data points are rebinned at 3σ for easier visualization. Each panel corresponds to a different flux interval, labeled F1 through F5, with F1 being the lowest flux state and F5 the highest. The vertical dashed lines highlight the evolution of the centroid energy of the primary absorption feature associated with the UFO as phenomenologically modeled with a Gaussian.

In the text

	Fig. 6. Same as Fig. 5 but for the time-resolved spectra labeled T1 through T5, with T1 being the earliest time interval and T5 the latest.
In the text

Fig. 7. Correlations among key spectral and wind properties of IRAS 13224-3809. Left: Photon index Γ versus 3–10 keV luminosity for the flux-resolved spectra (F1–F5), illustrating the softer-when-brighter trend. Center: Outflow velocity of the UFO as a function of 3–10 keV luminosity for the time-resolved spectra (T1–T5), highlighting their coupled evolution. Right: Relation between UFO outflow velocity and Fe Kα line width, tentatively indicating that faster winds are associated with greater velocity dispersion. Shaded regions represent 1σ uncertainties on the best-fit correlations, derived from 10 000 bootstrap regressions sampling the error boxes. Each panel includes the p-value from an exact permutation test assessing null hypotheses of no correlation and the best-fit relation.

In the text

	Fig. 8. Ionization parameter (log ξ) of the UFO and column density (NH) corrected for special relativistic effects for the flux-resolved spectra (F1-F5) of IRAS 13224-3809. Here we are directly visualizing the trends when fixing one of the two parameters to its average value across the intervals.
In the text

Fig. 9. Three-panel summary of the time-resolved properties of the UFO in IRAS 13224-3809, fitted with model A (absorption only). Top: 3–10 keV power-law luminosity and UFO outflow velocity as a function of time, showing their coupled evolution. Middle: UFO acceleration, derived from velocity changes between intervals, highlighting rapid increases during bursts. Bottom: UFO kinetic power expressed as a percentage of the adopted bolometric luminosity; the right axis shows the wind momentum rate normalized to the X-ray radiation momentum rate. The left axis traces how energetic the wind is, while the right axis traces how difficult it is for the instantaneous X-ray radiation field alone to account for that wind momentum. In all panels, the background displays the 3–10 keV light curve for context.

In the text

	Fig. A.1. XMM-Newton Optical Monitor (OM) UVW1 light curve of IRAS 13224-3809 during the 2016 campaign. The source shows a stable UV emission with a peak-to-peak variability amplitude of only ∼0.17 mag (< 2%) over the entire observation period. The average magnitude is mUVW1 = 15.36 with a standard deviation of 0.03 mag.
In the text

Fig. B.1. Gaussian line scans performed on the flux-resolved spectra of IRAS 13224–3809 in the 3–10 keV range in the source rest frame. The line width is fixed to σline = 100 eV, and the ordinate axis shows the best-fitting line normalization, with negative values indicating absorption features. The colored contours indicate thresholds in ΔAIC relative to strong and very strong evidence for the inclusion of an additional line, corresponding to ΔAIC = 10 (purple) and ΔAIC = 20 (orange), respectively. These conservative ΔAIC thresholds are adopted to robustly highlight only the most significant structures; the same regions are also recovered when using likelihood-based Δχ2 (ΔW) contours commonly adopted in previous analyses.

In the text

	Fig. B.2. Same as Fig. B.1 but for the time-resolved spectra.
In the text

	Fig. D.1. EW of the UFO as a function of the 3-10 keV observed flux for the time-resolved spectra (T1-T5) of IRAS 13224-3809. The lack of correlation between the two quantities is directly visualized here, in contrast to the work of Parker et al. (2017), Pinto et al. (2018). The dashed line represents the best-fit anticorrelation reported in Parker et al. (2017) for comparison.
In the text

Fig. D.2. Simulated average spectrum obtained by stacking synthetic spectra generated for each flux-resolved interval, incorporating observed trends in continuum slope, outflow velocity and the tentative ones in ionization and column density. The resulting stacked spectrum (black) is compared to absorbed power laws for each simulated flux interval (dashed gray). The absence of significant emission features in the stacked spectrum indicates that the emission observed in the actual average spectrum is unlikely to be an artifact of stacking variable flat spectra.

In the text

Fig. D.3. Visualization of the key findings in the flux- and time-resolved analysis. Missing values are reported as × and upper–lower limits as triangles. Panels c) and d) show the outflow velocities per each flux- and time-resolved spectrum, respectively showing clearly the consistency of vout across the models. For better visibility, values from each model are offset by a small amount in the x-axis.

In the text