Radiatively Driven Clumpy X-Ray Absorbers in the NLS1 Galaxy IRAS 13224-3809

We used the 1.5 Ms data obtained by XMM-Newton EPIC-pn (Jansen et al. 2001; Strüder et al. 2001) and OM (Mason et al. 2001) in the summer of 2016, including the unscheduled observations (see Table 1 in Paper I). As our study requires a consecutive data set taken within a short period of time, we did not incorporate any archival data before 2016. In addition to the pn, of which its data processing was described in Paper I, the OM data taken with the UVW1 filter is reduced with the standard SAS routine omichain, and the spectrum is extracted using om2pha.

Since the UFO properties of IRAS 13224-3809 are known to be luminosity dependent (e.g., Parker et al. 2017; Pinto et al. 2018, and Chartas & Canas 2018), we created the intensity-sliced spectra. First, all the 0.3–10.0 keV events were binned into 1 ks intervals, as in previous studies (Parker et al. 2017; Pinto et al. 2018). We grouped the data into 10 intensity levels based on the flux in the 0.3–10 keV range and generated 10 intensity-sliced spectra (see Table 2 in Paper I), which are labeled from A (the brightest) to J (the dimmest). The sliced spectral bins are grouped with 8 and 16 bins in the 0.3–2.0 keV and 2.0–10.0 keV ranges, respectively. These intensity-sliced spectra are shown in the upper panel of Figure 1. Next, the observed spectral ratios are created by taking the ratio of each intensity-sliced spectrum (B–J) to the brightest spectrum A (lower panel of Figure 1).

Zoom In Zoom Out Reset image size

We briefly describe the spectral-ratio model fitting method, the details of which were explained in Paper I. The spectral continuum of IRAS 13224-3809 consists of the power law (PL) and the soft-excess component, and the continuum is absorbed by several intervening gases such as WA, UFO, and clumpy absorbers. WAs do not change their geometry or ionization structure significantly in a timescale of a few weeks, while the partial covering fraction of the clumpy absorbers, whose variation influences the spectral shape below ∼5 keV, significantly changes in this timescale (e.g., Di Gesu et al. 2015; Midooka et al. 2022a). Therefore, taking the spectral ratio below 5 keV enables us to focus on the variability of the clumpy absorbers, canceling out the less time-variable spectral components. It should be noted that a similar spectral-ratio technique has been successfully adopted to study the variation of the high energy cutoff in the AGN spectra (Zhang et al. 2018).

Model spectral ratios in the 0.3–5.0 keV range are calculated using XSPEC (version 12.12.0). We consider a simple model spectrum in which the continuum is covered by an ionized partial absorber (zxipcf; Reeves et al. 2008). This absorber model has three parameters; ionization parameter ξ, hydrogen column density N_H, and partial covering fraction (CF). A table model of the spectral ratio is created by taking the ratio of the absorbed model spectra to the non-absorbed one (CF = 0) with the same continuum parameters (see Figure 2 in Paper I).

The three parameters of the clumpy absorber determine the shape of the spectral ratio, including the characteristic dip structures in the 0.8–1.0 keV range (Figure 2 in Paper I). Fitting the observed and model spectral ratios makes it possible to constrain the partial absorber parameters without being affected by complicated WA spectral features.

One of the purposes of the present study is to constrain the clump outflow velocities. In this section, we show that the flux dependency of the clump velocity can be constrained model independently from the spectral ratios.

We produced nine spectral ratios from the 10 observed intensity-sliced spectra. The upper panel in Figure 2 shows the representative spectral ratios of B/A, F/A, and J/A, each of which has a dip- or cliff-like structure at around 1 keV. We define the dip characteristic energy where the spectral ratio changes from convex upward to convex downward. To do so, first, the observed spectral ratios are fitted with eighth-order polynomial equations (solid lines in the upper panel in Figure 2). Then, the first derivatives of the eighth-order equation models are calculated (lower panel in Figure 2). The energies that give the minimums of the first derivative correspond to the dip characteristic energies defined above (dotted vertical lines in the lower panel of Figure 2). We find that the dip characteristic energy is blueshifted when X-rays become brighter. Note that this result is obtained without assuming any spectral models.

Zoom In Zoom Out Reset image size

We note that the relativistic reflection model (e.g., Parker et al. 2017; Chartas & Canas 2018; Pinto et al. 2018, and Jiang et al. 2022) can also reproduce these spectral ratios because the individual spectra are fitted. However, changes in the continuum and/or disk reflection normalizations are insufficient to account for the observed blueshift of the dip structure, and the disk geometry and other parameters need to change concordantly to realize such spectral changes. We do not discuss this possibility any further in this paper.

Next, we are going to fit the observed spectral ratios with physical models to constrain absorbers' parameters. Having noticed the blueshifted absorption features around 1 keV, we introduce the clumpy absorber velocity v_out in addition to the three clump parameters (CF, N_H, and ξ). Also, the continuum normalization ratio is a free parameter.

The spectral-ratio table model was created with these four free parameters (CF, ξ, N_H, and v_out; see Table 3 in Paper I), where models are linearly interpolated between the adjacent grid points using the XSPEC built-in function. First, the least-squares fitting was performed for the nine observed spectral ratios allowing all five parameters free (the four parameters above and the normalization). We found that the N_H and ξ were not significantly variable among the nine ratios ($2.7\lt \mathrm{log}\,\xi \lt 3.0$, 23.7 < log N_H <23.9). Thus, we tied these two parameters for all nine spectral ratios and performed simultaneous fitting.

Table 1 shows all the fitting results, and Figure 3 shows examples of the fitting of B/A, D/A, F/A, and J/A. The fitting was successful in all the spectral ratios below ∼2 keV, where the absorbers' parameters are constrained. Note that the model does not fit well above ∼3 keV, especially in dimmer spectral ratios, suggesting a change in the photon index, that is confirmed later in the spectral fitting (Section 3.3). The best-fit parameters are illustrated in Figure 4. We found that the partial covering fraction is smaller and the norm is higher as the X-ray flux becomes higher. The outflow velocity of the partial absorber was 0.2–0.3c, unexpectedly fast, even comparable to that of the UFOs reported in previous studies (Parker et al. 2017; Chartas & Canas 2018; Pinto et al. 2018). Furthermore, the velocity was found to be higher when the X-ray flux is higher.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We aim to confirm that the velocity we have obtained is determined independently, not degenerate with the ionization parameter ξ, since both parameters similarly affect absorption line/edge energies. We perform a Markov Chain Monte Carlo (MCMC) calculation using the emcee package (Foreman-Mackey et al. 2013) to determine the posterior distribution of the best-fit spectral parameters. The initial parameters in the MCMC chains are set to be close to the best-fit parameters obtained by the chi-square fitting. We set a uniform prior distribution over the parameter ranges (Table 3 in Paper I). After the initial 5000 steps are discarded to exclude the burn-in phase, further 10,000 steps are explored by 1000 separate chains (walkers).

We succeeded in the MCMC parameter estimation for all the spectral ratios. The best-fit model ratios are shown in the Appendix (Figure A1). The blue dots in Figure 4 show the mean values obtained by the MCMC calculation. We confirmed that the best-fit parameters and their flux dependency are similar to those by the chi-square fitting (orange dots in Figure 4). Note that ξ is not tied among the different ratios unlike the least chi-square fitting, because the primary aim of this analysis is to ensure whether there is no correlation between v_out and ξ.

Figure 5 shows a corner plot of the two-dimensional correlations between the posterior distributions of the estimated parameter pairs in the spectral ratio F/A. The other corner plots are shown in the Appendix, Figure A2. We found that the velocity is not correlated with other parameters and is determined independently. In contrast, CF and ξ are correlated in most spectral ratios. Still, whereas the CF shows a clear anticorrelation with the X-ray flux, the ξ variation is tiny and will not affect the intensity variations (Figure 4).

Zoom In Zoom Out Reset image size

We performed spectral fitting in the 0.3–10 keV band to constrain the UFO velocity with the partial absorption model based on the ratio model fitting. The spectral data are binned to have at least 20 counts per new energy bin. The XSPEC spectral model adopted is represented as

$$\begin{eqnarray}&&\begin{array}{l}{\mathtt{phabs}}\ast {\mathtt{zxipcf}}\ast ({\mathtt{powerlaw}}+{\mathtt{diskbb}})\\ \ast \,{\mathtt{const}}\ast {\mathtt{kabs}},\end{array}\end{eqnarray} \tag{ 1 }$$

where each spectral component is explained below.

The column density of phabs for the interstellar absorption was fixed at the Galactic value N_H = 5 × 10²⁰ cm⁻² (Bekhti et al. 2016), where we set the cosmic abundances to the wilms values (Wilms et al. 2000). In addition to the model for the ratio fitting, the kabs model (Ueda et al. 2004, updated by Tomaru et al. 2020) is incorporated to explain the Fe xxv and/or Fe xxvi UFO absorption features. The kabs model calculates the Voigt profile absorption lines (Kα and Kβ) for Fe xxv and Fe xxvi for the given column density of the ion N_atom and the turbulent velocity v_turb. The redshift z due to the line-of-sight velocity is also a free parameter. Although Pinto et al. (2018) reported the presence of an O viii UFO absorption line in the RGS spectra, we do not include the absorption line in the spectral model for simplicity. Also, disk reflection components are not required in our spectral model.

The initial parameters of the spectral fitting are set to the values constrained by the ratio model fitting. Both the continuum parameters and the partial absorption parameters were free to vary for the fitting. It is difficult to place limits on the ionization parameter and the turbulent velocity v_turb simultaneously from the absorption line profiles above 8 keV. Therefore, v_turb was fixed at 10,000 km s⁻¹, and the fitting was performed assuming that the UFO absorption lines are produced either by He-like Fe or H-like Fe.

In the spectral fitting, we include Kα and Kβ lines as a pair. Three Kα lines (accompanied by three Kβ lines) at maximum are necessary to fit the UFO feature depending on the spectrum. This indicates that more than one velocity component of the UFO was observed. The fitting result of a typical intensity-sliced spectrum F is shown in Figure 6, while the others are shown in Figure A1. Best-fit parameters are given in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Best-fit Parameters Determined by the Spectral Fitting in the 0.3–10.0 keV Band

Components	Parameters	A	B	C	D	E
zxipcf	NH (1022cm−2)	⋯	${36}_{-7}^{+13}$	${46}_{-6}^{+9}$	${51}_{-5}^{+4}$	${56}_{-5}^{+5}$
	log ξ	⋯	${2.72}_{-0.04}^{+0.04}$	${2.73}_{-0.05}^{+0.03}$	${2.73}_{-0.11}^{+0.05}$	${2.72}_{-0.03}^{+0.03}$
	z	⋯	$-{0.33}_{-0.02}^{+0.01}$	$-{0.33}_{-0.01}^{+0.01}$	$-{0.32}_{-0.01}^{+0.01}$	$-{0.32}_{-0.01}^{+0.01}$
powerlaw	Γ	${2.67}_{-0.03}^{+0.03}$	${2.71}_{-0.03}^{+0.03}$	${2.60}_{-0.03}^{+0.03}$	${2.60}_{-0.04}^{+0.03}$	${2.54}_{-0.03}^{+0.03}$
	norm (10−3)	${2.2}_{-0.1}^{+0.1}$	${2.3}_{-0.1}^{+0.1}$	${2.2}_{-0.1}^{+0.1}$	${2.3}_{-0.2}^{+0.1}$	${2.4}_{-0.1}^{+0.1}$
diskbb	kTin (eV)	${182}_{-1}^{+1}$	${178}_{-2}^{+2}$	${177}_{-2}^{+2}$	${179}_{-3}^{+2}$	${177}_{-2}^{+2}$
	norm	${731}_{-33}^{+34}$	${757}_{-62}^{+45}$	${845}_{-55}^{+46}$	${781}_{-57}^{+89}$	${831}_{-59}^{+72}$
Assuming Fe xxv
kabs	NFexxv (1018cm−2)	⋯	${2}_{-1}^{+2}$	${4}_{-1}^{+2}$	${4}_{-1}^{+2}$	${5}_{-2}^{+2}$
	z	⋯	$-{0.26}_{-0.10}^{+0.04}$	$-{0.29}_{-0.02}^{+0.03}$	$-{0.24}_{-0.02}^{+0.01}$	$-{0.25}_{-0.02}^{+0.02}$
kabs	NFexxv (1018cm−2)	⋯	⋯	⋯	${8}_{-3}^{+6}$	${6}_{-2}^{+5}$
	z	⋯	⋯	⋯	$-{0.33}_{-0.01}^{+0.01}$	$-{0.32}_{-0.02}^{+0.02}$
kabs	NFexxv (1018cm−2)	⋯	⋯	⋯	⋯	⋯
	z	⋯	⋯	⋯	⋯	⋯
Assuming Fe xxvi
kabs	NFexxvi (1018cm−2)	⋯	${4}_{-3}^{+4}$	${7}_{-3}^{+5}$	${8}_{-3}^{+4}$	${10}_{-4}^{+5}$
	z	⋯	$-{0.23}_{-0.11}^{+0.04}$	$-{0.26}_{-0.02}^{+0.03}$	$-{0.21}_{-0.02}^{+0.01}$	$-{0.22}_{-0.02}^{+0.02}$
kabs	NFexxvi (1018cm−2)	⋯	⋯	⋯	${17}_{-7}^{+12}$	${12}_{-5}^{+9}$
	z	⋯	⋯	⋯	$-{0.30}_{-0.01}^{+0.01}$	$-{0.29}_{-0.01}^{+0.02}$
kabs	NFexxvi (1018cm−2)	⋯	⋯	⋯	⋯	⋯
	z	⋯	⋯	⋯	⋯	⋯
χ2/dof		760.5/677	731.3/628	715.4/639	616.7/608	747.9/619
Components	Parameters	F	G	H	I	J
zxipcf	NH (1022cm−2)	${52}_{-5}^{+6}$	${56}_{-5}^{+5}$	${54}_{-4}^{+4}$	${57}_{-4}^{+4}$	${60}_{-3}^{+3}$
	log ξ	${2.76}_{-0.16}^{+0.04}$	${2.78}_{-0.03}^{+0.03}$	${2.81}_{-0.03}^{+0.03}$	${2.85}_{-0.02}^{+0.02}$	${2.91}_{-0.02}^{+0.02}$
	z	$-{0.31}_{-0.01}^{+0.01}$	$-{0.31}_{-0.01}^{+0.01}$	$-{0.30}_{-0.01}^{+0.01}$	$-{0.29}_{-0.01}^{+0.01}$	$-{0.28}_{-0.01}^{+0.00}$
powerlaw	Γ	${2.47}_{-0.04}^{+0.04}$	${2.40}_{-0.04}^{+0.04}$	${2.34}_{-0.04}^{+0.04}$	${2.24}_{-0.05}^{+0.05}$	${1.95}_{-0.05}^{+0.05}$
	norm (10−3)	${2.2}_{-0.2}^{+0.3}$	${2.2}_{-0.2}^{+0.2}$	${2.0}_{-0.2}^{+0.2}$	${1.8}_{-0.2}^{+0.2}$	${1.2}_{-0.1}^{+0.1}$
diskbb	kTin (eV)	${173}_{-2}^{+2}$	${170}_{-2}^{+2}$	${168}_{-1}^{+2}$	${165}_{-1}^{+1}$	${165}_{-1}^{+1}$
	norm	${986}_{-99}^{+85}$	${1085}_{-89}^{+85}$	${1238}_{-83}^{+83}$	${1406}_{-78}^{+77}$	${1589}_{-62}^{+63}$
Assuming Fe xxv
kabs	NFexxv (1018cm−2)	${3}_{-1}^{+1}$	${2}_{-1}^{+1}$	${8}_{-2}^{+3}$	${3}_{-1}^{+1}$	${7}_{-1}^{+2}$
	z	$-{0.18}_{-0.02}^{+0.01}$	$-{0.14}_{-0.02}^{+0.02}$	$-{0.24}_{-0.01}^{+0.01}$	$-{0.20}_{-0.01}^{+0.01}$	$-{0.22}_{-0.01}^{+0.01}$
kabs	NFexxv (1018cm−2)	${5}_{-1}^{+2}$	${10}_{-2}^{+3}$	${6}_{-2}^{+4}$	${7}_{-2}^{+3}$	${7}_{-2}^{+4}$
	z	$-{0.25}_{-0.01}^{+0.01}$	$-{0.23}_{-0.01}^{+0.01}$	$-{0.32}_{-0.02}^{+0.01}$	$-{0.26}_{-0.01}^{+0.01}$	$-{0.31}_{-0.01}^{+0.01}$
kabs	NFexxv (1018cm−2)	${9}_{-3}^{+26}$	${5}_{-2}^{+3}$	⋯	${3}_{-1}^{+2}$	⋯
	z	$-{0.34}_{-0.04}^{+0.01}$	$-{0.31}_{-0.02}^{+0.02}$	⋯	$-{0.33}_{-0.02}^{+0.02}$	⋯
Assuming Fe xxvi
kabs	NFexxvi (1018cm−2)	${6}_{-2}^{+3}$	${5}_{-2}^{+2}$	${16}_{-4}^{+7}$	${7}_{-3}^{+3}$	${16}_{-3}^{+5}$
	z	$-{0.14}_{-0.02}^{+0.02}$	$-{0.10}_{-0.02}^{+0.02}$	$-{0.21}_{-0.01}^{+0.01}$	$-{0.16}_{-0.01}^{+0.01}$	$-{0.19}_{-0.01}^{+0.01}$
kabs	NFexxvi (1018cm−2)	${11}_{-3}^{+5}$	${20}_{-4}^{+7}$	${13}_{-5}^{+9}$	${15}_{-5}^{+7}$	${16}_{-4}^{+8}$
	z	$-{0.22}_{-0.01}^{+0.01}$	$-{0.20}_{-0.01}^{+0.01}$	$-{0.29}_{-0.02}^{+0.01}$	$-{0.22}_{-0.01}^{+0.01}$	$-{0.29}_{-0.01}^{+0.01}$
kabs	NFexxvi (1018cm−2)	${18}_{-6}^{+45}$	${10}_{-4}^{+7}$	⋯	${6}_{-3}^{+4}$	⋯
	z	$-{0.31}_{-0.03}^{+0.01}$	$-{0.28}_{-0.02}^{+0.02}$	⋯	$-{0.30}_{-0.02}^{+0.02}$	⋯
χ2/dof		676.4/625	712.9/624	766.9/629	665.8/612	779.4/642

Notes. The CF of zxipcf and const are fixed to the values determined by the chi-square ratio fitting (upper panel of Table 1). Fitting was made assuming that the UFO absorption lines are produced either by He-like Fe or H-like Fe. The parameters except kabs are determined assuming the He-like Fe. The fitting results are shown in Figure A3 and A4. Hereafter the error ranges correspond to the 90% confidence level.

Download table as:  ASCIITypeset image

Note that our previous study based on a similar partial covering model required a strong absorption edge around 1 keV to explain the residual feature (Yamasaki et al. 2016). Now, the fit is successful without such an additional edge, which suggests that the residual feature was caused by the energy shift of the outflowing clumpy absorbers.

In this section, we model the broadband spectral energy distribution (SED) with the pn and OM data to estimate the bolometric luminosity, affecting the disk winds and ionization degree in the ambient matter. Since we did not find significant flux variations in the OM data compared to the X-ray variation, we use the time-averaged OM spectrum.

Since IRAS 13224-3809 is a super-Eddington object and the AGN component is significant, we ignore the host galaxy contamination. In this study, we assume that spectrum A is not partially absorbed (CF = 0), so we model SED of the spectrum A with the agnslim model (Kubota & Done 2019), multiplied by the reddening (redden; Cardelli et al. 1989) and the interstellar absorption (phabs). Here agnslim is a broadband spectral model for super-Eddington AGN based on the slim disk emissivity. The accretion flow is radially stratified; the inner hot Comptonization region (R_in to R_hot) with the electron temperature kT_hot and photon index Γ_hot, the intermediate warm Comptonization region to produce the soft X-ray excess (R_warm to R_hot) with the electron temperature kT_warm and photon index Γ_warm, and the outer standard disk (R_warm to R_out).

For a highly super-Eddington slim disk, the inner radius of the disk (R_in) is determined not only by the black hole spin but also by the gas pressure, so that R_in can be smaller than the innermost stable circular orbit (Watarai et al. 2000; Kubota & Done 2019). We assume the spin parameter a* = 0, and set R_in to the radius calculated by Kubota & Done (2019) (see Equation (1) in their paper). We also fixed the column density of the interstellar absorption N_H to 5 × 10²⁰ cm⁻² and E(B − V) = 1.7 × N_H/10²² cm⁻² (Bohlin et al. 1978).

The averaged OM spectrum and pn spectrum A are successfully fitted with the assumed model. Figure 7 shows the ν f_ν plot of the broadband SED in the 3 eV–10 keV range. The gray lines indicate the three representative intensity-sliced spectra A, E, and J, and the blue line represents the best-fit SED model for spectrum A. Best-fit parameters are given in Table 3.

Zoom In Zoom Out Reset image size

Table 3. Best-fit Parameters Determined by the SED Fitting Obtained by OM and pn

Components	Parameters	Best-fit values
phabs	NH (1020cm−2)	5 (fixed)
redden	E(B − V)	0.086 (fixed)
agnslim	mass (M⊙)	5.4 × 106 (fixed)
	$\dot{M}$ (${\dot{M}}_{\mathrm{Edd}}$)	6.30 ± 0.01
	a*	0 (fixed)
	cos i	0.5 (fixed)
	kThot (keV)	100 (fixed)
	kTwarm (keV)	0.166 ± 0.005
	Γhot	2.67 ± 0.04
	Γwarm	${1.66}_{-0.10}^{+0.09}$
	Rhot (Rg)	8.3 ± 0.2
	Rwarm (Rg)	10.6 ± 0.1

Download table as:  ASCIITypeset image

In Section 3, we constrained the outflow velocity, v_out, of the partial clumpy absorbers using the unique spectral-ratio fitting method. We also performed the conventional spectral fitting to determine the v_out of the UFO. A comparison of the outflow velocities of the UFOs and partial absorbers is shown in Figure 8, where these velocities are plotted against the 3.0–10.0 keV X-ray flux of the intensity-sliced spectra. When a UFO is composed of multiple absorption lines with different velocity components, the average velocity weighted by the column density of each absorption line (N_atom) is plotted.

Zoom In Zoom Out Reset image size

We fitted the velocities of the UFOs and clumpy absorbers as a function of the X-ray flux with linear functions. The orange and blue dotted lines in Figure 8 show the best-fit models with 1σ uncertainties. Although the slopes of the two best-fit functions are slightly different, we see similar increasing trends with the X-ray flux of the clump velocities and the UFO velocities either assuming the He-like or H-like absorption lines. In particular, the orange dots in Figure 8 assuming the He-like Fe absorption show similar velocity values with those of the clumpy absorbers. This finding is consistent with the hot inner and clumpy outer wind scenario by Mizumoto et al. (2019), where the outer clumpy absorbers are the consequence of the inner hot UFO winds.

If the proposed scenario is valid, there should be a correlation between the clumps and UFOs along the line of sight. To investigate this possibility, we plotted the clump covering fraction (see Table 1) and the equivalent width of the UFO absorption lines (Figure 9). The equivalent width was calculated as the sum of the individual equivalent widths of multiple absorption lines. Figure 9 indicates a clear correlation between the two variables, in particular, a surprising linearity in the B–G range. This finding suggests a simple fact that more material in the line of sight will result in more prominent UFOs and partial covering, providing further evidence of the hot inner and clumpy outer wind scenario.

Zoom In Zoom Out Reset image size

These results have also been suggested in PDS 456, which is another super-Eddington source. Matzeu et al. (2016) searched for the outflow velocity of the clumpy absorber to give the best spectral fit, and reported the global minimum solution at a velocity of 0.25c. Matzeu et al. (2017) also found in the same target that the centroid energy of Fe K UFO absorption increases with the ionizing flux. Reeves et al. (2021) similarly explained the spectrum of PDS 456 in terms of the partial absorption; Figure 16 in their paper clearly shows that the UFO absorption lines are more prominent when the soft X-rays are more significantly absorbed by the partial clumps.

As shown in Figure 1, the observed 0.3–10 keV flux variation ranges over more than one order of magnitude. On the other hand, the ratio fitting results show that the intrinsic X-ray flux is varied by only a factor of ∼4, indicating that a considerable amount of the apparent flux variation is caused by changes in the partial covering fraction.

The orange dashed line in Figure 7 shows the deabsorbed UV/X-ray spectrum for spectrum A in which the covering fraction is assumed to be null. From this deabsorbed model, the intrinsic UV flux in the 13.6 eV–0.1 keV range and the intrinsic X-ray flux in the 0.1–10.0 keV range are estimated respectively to be 1.8 × 10⁻¹⁰ and 6.1 × 10⁻¹¹ erg cm⁻² s⁻¹, using the flux command in XSPEC. The proportion of the intrinsic UV flux to X-ray flux is approximately 3: 1. Since the intrinsic X-ray flux varies by a factor of ∼4 while the UV flux is invariable, the total flux variation of UV and X-rays can be estimated to be about a factor of ∼16/13, assuming that the obtained ratio (UV flux: X-ray flux = 3:1) is valid throughout the observation period.

In case the disk wind is driven by continuum radiation, the total flux variation is expected to contribute directly to the increase in the kinetic energy of the disk wind. Since the outflow velocity is proportional to the square root of the wind kinetic energy, the velocity change is estimated to be about a factor of $\sqrt{16/13}\sim$ 1.1.

By extrapolating the best-fit model in Figure 8, the velocities of the UFO and clumpy absorbers in spectrum A (1.4 × 10⁻¹² erg cm⁻² s⁻¹ in the 3–10 keV range) are estimated to be 0.31c and 0.37c. As a result, the variations of the UFO and clump velocities from the dimmest spectrum (J) to the brightest one (A) are about a factor of 1.1 and 1.4, respectively. These values are consistent with the simple estimation above, thus both the clumpy absorbers and the UFOs are considered to be driven by the same UV-dominant continuum radiation.

Ionization parameter (ξ = L/nr²) values of the clumpy absorbers determined from the model fitting are hardly variable regardless of the X-ray flux (Tables 1 and 2). This fact suggests that the gas density of the clump, n, increases or that the location where the clump is formed, r, becomes more distant as the X-ray luminosity increases.

Takeuchi et al. (2014) suggest that the clump is formed by a sort of radiation hydrodynamic instability in addition to the Rayleigh–Taylor instability. The radiation instability should work only at a certain opacity, i.e., a certain ξ. If the instability is activated and the clumps are produced at this specific distance, which increases with the luminosity, the result that ξ being constant is reasonable. This may be a common feature of the super-Eddington objects with continuum-driven outflows. It will be intriguing to investigate the flux dependence of ξ in other super-Eddington sources in the future.

IRAS 13224-3809 has been observed many times and extensively investigated to date. This section will discuss how our model can explain these observed properties of IRAS 13224-3809.

In our model, spectral variability is explained by the variation of the PL continuum component and the partial absorption component. On the other hand, Parker et al. (2017) used the principle component analysis (PCA) to decompose the spectral variability, and argued that the PL emission dominates the spectral variability, with a minor contribution from the reflection and UFO component. The PCA is a powerful data science tool to reduce the dimensionality of the input data, but one must be careful with its applicability in X-ray spectral analysis. Pike et al. (2017) performed a spectral component decomposition of the simulated X-ray spectral variations of NLS1 galaxies using the PCA and the non-negative matrix factorization (NMF). Despite the PCA's good performance in determining the number of independent spectral components, the decomposed components derived by the PCA include the negative values as residuals from the average spectrum, being very different from the original components. The other method, NMF, has the advantage of being composed of non-negative values, making it easier to interpret the spectral components, in particular when the variable components are additive. Still, NMF is not very useful when variable spectral components are multiplicative, just like our model. We expect that forthcoming analyses using new techniques will reveal the main cause of the spectral variation without prior assumptions.

X-ray observations provide not only spectroscopic information but also timing information. Characteristic power spectral density (PSD) and soft/hard lags have been reported in IRAS 13224-3809. The PSD peaks at timescales of about 30 ks (e.g., Alston et al. 2019). The radiation-magnetohydrodynamic simulations suggest the clumpy absorbers exist at a few hundred R_s from the center (Takeuchi et al. 2013). Assuming the clump exists at 500 R_s, the Keplerian rotation period is calculated to be about 200 ks. Since the typical clump size is around 10 R_s (Takeuchi et al. 2013), a variation of 30 ks is reasonable with our model in terms of changes in the partial covering fraction.

The lag-frequency plot can be derived from the light curves of the primary X-ray emission (typically 1.0–4.0 keV) and the soft excess (below ∼1.0 keV). In this plot, a hard lag is observed in the low-frequency band below 10⁻⁴ Hz, while a soft lag is observed in the high-frequency band above 10⁻³ Hz (e.g., De Marco et al. 2013 and Alston et al. 2020). The low-frequency hard lag has a PL-like shape and decreases with frequency. This is considered to be due to the inward propagation of the fluctuations through the accretion flow and corona (e.g., Arévalo & Uttley 2006).

Since the mechanism of the soft X-ray excess has not yet been clearly understood, two major hypotheses exist as to the cause of the high-frequency soft lag. One is reverberation due to reflection by the inner region of the accretion disk around ∼10 R_s (e.g., Fabian et al. 2009 and Zoghbi et al. 2012), where the reprocessed emission is delayed by the light travel time from the source. The other hypothesis says that the lag is a solely mathematical artifact due to phase wrapping (Miller et al. 2010b; Legg et al. 2012; Mizumoto et al. 2019), which is also favored by our model.

When examining the energy dependence of the lag in the high-frequency range where the soft lag is observed, a significant peak structure is often found in the Fe K band. This is thought to be caused by the primary X-rays scattering off the surrounding material, generating the Fe emission lines with a reverberation lag (e.g., Zoghbi et al. 2012 and Kara et al. 2013). If this lag is literally taken as the light propagation time due to the reflection in the inner accretion disk, the size is about 10 R_s. On the other hand, if we take into account the dilution effect due to the weaker iron emission line compared to the continuum, the lag can be attributed to the reflections from the materials located as far away as ∼100 R_s (e.g., Miller et al. 2010a). We speculate that the scattering by the disk wind in our model at ∼100 R_s can explain the observed reverberation lag (Mizumoto et al. 2018, 2019).

X-ray light curves of almost any kind of accreting object, including AGNs, exhibit a linear relationship between the rms amplitude and the flux variations (e.g., Uttley & McHardy 2001 and Alston et al. 2019). The fact that this relationship has been observed for any objects with accretion disks suggests that it is due to intrinsic X-ray variability, not due to apparent absorption effects. We suppose that the intrinsic PL variability in our model is responsible for this relation, while quantitative estimation is a future subject.