Mapping the X-ray corona evolution of IRAS 13224-3809 with the power spectral density

We develop the power spectral density (PSD) model to explain the nature of the X-ray variability in IRAS 13224-3809, including the full effects of the X-ray reverberation due to the lamp-post source. We utilize 16 XMM-Newton observations individually as well as group them into three different luminosity bins: low, medium and high. The soft (0.3-1 keV) and hard (1.2-5 keV) PSD spectra are extracted and simultaneously fitted with the model. We find that the corona height changes from h ~ 3 $\ r_{\rm g}$ during the lowest luminosity state to ~ $25 \ r_{\rm g}$ during the highest luminosity state. This provides further evidence that the source height from the reverberation data is significantly larger than what constrained by the spectral analysis. Furthermore, as the corona height increases, the energy spectrum tends to be softer while the observed fractional excess variance, $F_{\rm var}$, reduces. We find that the PSD normalization is strongly correlated with $F_{\rm var}$, and moderately correlated with the PSD bending index. Therefore, the normalization is dependent on accretion rate that controls the intrinsic shape of the PSD. While the intrinsic variability of the disk is manifested by the reverberation signals, the disk and corona may evolve independently. Our results suggest that, during the source height increases, the disk itself generates less overall variability power but more high-frequency variability resulting in the PSD spectrum that flattens out (i.e. the inner disk becomes more active). Using the luminosity-bin data, the hint of Lorentzian component is seen, with the peak appearing at lower frequencies with increasing luminosity.


INTRODUCTION
The active galactic nucleus (AGN) of the narrow-line Seyfert 1 galaxy IRAS 13224-3809 contains a maximally spinning supermassive black hole Chiang et al. 2015;Jiang et al. 2018) with the mass of ∼ 2×10 6 M (Alston et al. 2020). It is one of the AGNs with a complex energy spectrum and complex X-ray variability over a broad range of timescales (e.g. Alston et al. 2019). The energy-integrated spectrum during the long 2011 XMM-Newton observation could be explained by a patchy disk model that produced two reflection components from separate ionized elements, with an overabundance of iron . The disk black-body emission and the narrow emission line at ∼ 6.4 keV were required. Furthermore, the time delays between X-ray variability in the reflection and continuum dominated energy bands (referred to as reverberation lags) were measured by Kara et al. (2013). They found that during low-flux stages, changes in the observed reverberation lags supported a compact coronal geometry. Chiang et al. (2015) discovered that the reflection component of IRAS 13224-3809 is substantially less variable than the power-law emission, supporting the light-bending framework (e.g. Miniutti & Fabian 2004). The spectral fitting revealed a constant emitting region that produced a soft thermal emission following the expected L ∝ T 4 blackbody relation (Chiang et al. 2015;Jiang et al. 2018), which could have an accretion disk origin (Ponti et al. 2010). Chainakun et al. (2016) performed simultaneous modelling of the mean and lag-energy spectra taken into account the full dilution and ionization effects, and found that the soft excess was dominated by the narrow components which could be produced by the distant reflection from the cold torus, rather than dominated by the broad features from the inner disk reflection.
Furthermore, Parker et al. (2017) identified a series of variable peaks in the long-term X-ray variability spectra that could be interpreted as the strong absorption lines from an ultra-fast outflow (UFO). Alston et al. (2019) found that the power spectrum density (PSD) of IRAS 13224-3809 observed in 2002 (∼ 64 ks), 2011 (∼ 500 ks), and 2016 (∼ 1.5 Ms) contained non-stationary multiple-peaked components whose normalization increased while the low-frequency peak moved to higher frequencies when the source flux decreased. The intrinsic PSD shape is mostly determined by the mass accretion rate fluctuations that are propagated along the accretion disk (Lyubarskii 1997;Churazov et al. 2001;Arévalo & Uttley 2006;Ingram & Done 2011). The accretion rate of IRAS 13224-3809 implied from the PSD spectra was comparable to that of black hole X-ray binaries in very-high/intermediate states (Alston et al. 2019).
Recently, Jiang et al. (2022) utilized the high-density disk model to fit the time-averaged spectra of IRAS 13224-3809. They considered a broken power-law emissivity and a free reflection fraction parameter. Based on ionization and density parameters assuming the geometry is a lamppost, the source heights were calculated to be ∼ 3-6 gravitational radii (r g = GM/c 2 ; M is the central balck hole mass, G is the gravitational constant and c is the speed of light). By considering the lag-frequency spectra in multiple observations during 2002-2016, the source height was found to increase with increasing luminosity, from ∼ 6 r g to 20 r g (Alston et al. 2020). Furthermore, Caballero-García et al. (2020) analyzed the combined spectral-timing data simultaneously in different flux periods. For the maximally spinning case, they also found a tendency of rising source height, from ∼ 3 r g to 10 +10 −1 r g , with luminosity. Nevertheless, the PSD profiles can be imprinted with X-ray reverberation patterns (Emmanoulopoulos et al. 2016;Papadakis et al. 2016;Chainakun 2019), providing an independent way to probe the coronal geometry. Due to the gravitational light-bending effects, the reflection component is less variable compared to the power-law continuum. The reverberation signals then reduce the fractional excess variance (F var ) so producing the dip in the PSD profiles that is more prominent in the more reflection-dominated band. The models for the X-ray reverberation signatures for the lag-spectra of AGN have been investigated extensively (e.g. Wilkins & Fabian 2013;Cackett et al. 2014;Emmanoulopoulos et al. 2014;Chainakun et al. 2016;Epitropakis et al. 2016), especially in IRAS 13224-3809 (Alston et al. 2020;Caballero-García et al. 2020). We then choose to investigate the X-ray variability power using the PSD model that includes the full effects of the X-ray reverberation. We focus on the energy-dependence of the PSD shapes for 16 observations of IRAS 13224-3809 under the lamp-post assumption. We explore two energy bands: 0.3-1 keV and 1.2-5 keV as the representative of the reflection-dominated and continuum-dominated bands, respectively. This allows us to check for consistency the implied framework and lamp-post geometry using different timing profiles.
In Section 2, we present the data reduction and how the PSD spectra are produced. The theoretical PSD models including X-ray reverberation effects are described in Section 3. Section 4 explains the model grid created in this study as well as the fitting procedure. The best-fit results are presented in Section 5. Discussion and conclusion is given in Section 6 and 7, respectively.

OBSERVATIONS AND DATA REDUCTION
The X-ray data of IRAS 13224-3809 used in this work were previously observed by XMM-Newton observatory (Jansen et al. 2001) and were obtained from XMM-Newton Science Archive. 1 Since one purpose of this study is to analyze each observational data individually, to get high signal to noise data, we selected only the observations which have the total exposure time of 100 ks. The XMM-Newton observational data used here are tabulated in Table 1. To avoid the combination of data from pn and MOS detectors which were observed with different time resolutions, we chose to consider only the pn data which have higher time resolution and effective area. 2 We performed the data reduction using Science Analysis Software (SAS) version 19.1.0 with the latest version of the calibration files (CCF). 3 The pn observation data files were reprocessed using the SAS task epproc with the default parameter values. Then, the observational periods which were affected by high background flaring activity were also removed by the SAS task 1 http://nxsa.esac.esa.int 2 https://xmm-tools.cosmos.esa.int/external/xmm user support/documentation/uhb/epic.html 3 https://www.cosmos.esa.int/web/xmm-newton/download-and-install-sas Note. a Good exposure time after data cleaning. b The background-subtracted source count rate in the 1-4 keV energy band and c the luminosity bin in which the data belong to, regarding to its count rate (see text).
espfilt using the method histogram with the parameter allowsigma of 2.5 (default value). The remaining exposure time after removing background flaring for each observation is shown in the fourth column of Table 1. For all observations, we consider the lightcurves of IRAS 13224-3809 extracted in two energy bands which are 0.3-1 keV (reflection dominated) and 1.2-5 keV (continuum dominated) referred to as the soft and hard energy bands, respectively. The background-subtracted light curve in each energy band was extracted from the events flagged with PATTERN 4 and #XMMEA EP using the SAS task evselect and epiclccorr; the source extraction region was defined as a circular area centred at the source position with the radius of 20 arcsec, while the background region was defined as a 60 arcsec radius circle, located on a source-free area that is still on the same CCD chip with that of the source region. Then, the produced light curves were converted into the power spectral density (PSD) utilising the ftools task powspec; 4 concisely, each light curve was divided into a number of segments with the length and the time-bin resolution of 20 ks and 179 s, respectively. Then all segments were converted into the PSDs, in which their output frequency (f ) is corresponding to 0.05 -2.8 mHz, and averaged over; the Poisson noise was also subtracted during this step. Finally, the obtained PSD was rebinned logarithmically, in which the width of the next, higher frequency bin is larger by a factor of 1.06 (f → 1.06f ), to get the PSD for analysing further.
Moreover, we investigated the variability of PSD as a function of the luminosity when increasing their signal to noise. To do this, we categorised the observational data into three groups -low, medium and high luminositiesbased on the source's luminosity; here we used the background-subtracted, instrumental count rate, i.e. that of the pn detector, of the source in the 1-4 keV energy band as a proxy of the luminosity since this band should be dominated by the primary X-ray emission (Caballero-García et al. 2020). The low, medium and high luminosity observations were defined as the observations that have the count rate < 0.19 counts s −1 , 0.19 counts s −1 count rate 0.28 counts s −1 , and count rate > 0.28 counts s −1 , respectively. The count rate and the group that each observation belongs to are shown in column 5 and 6 of Table 1. The light curves in each luminosity bin were grouped together and then converted in the single PSD of each luminosity bin following the method explained above. The PSDs obtained from individual observations and grouped observations were then used in our analysis.
Note that the common 1-4 keV energy band was used only for grouping the observational data into 3 different luminosity states: low, medium, and high luminosity. When we calculated the PSD in the soft and hard bands, we selected to follow Alston et al. (2020) where the 0.3-1 keV and 1.2-5 keV bands were used to represent the soft (reflection dominated) and the hard (intrinsic emission dominated) bands, respectively.

VARIABILITY POWER MODEL
The PSD of AGN can be generally explained by a broadband power-law with one or two bend frequencies where the profile changes its slope (Papadakis 2004;González-Martín & Vaughan 2012). We use the PSD model in the form of a bending power-law (e.g. Emmanoulopoulos et al. 2016): where the PSD has a low-frequency slope −1 that is bent gradually to the high-frequency slope s above the bending frequency f b . A is the normalization factor of the variability power. The reverberation signals imprint the oscillatory structures on the intrinsic PSD profile Chainakun 2019). We use the ray-tracing simulations to generate the response functions of the disk reflection under the lamp-post scenario, by tracing photon paths along the Kerr geodesics from the source to the disk and to the observer's sky (e.g. Karas et al. 1992;Fanton et al. 1997;Wilkins & Fabian 2013;Cackett et al. 2014;Chainakun et al. 2016;Emmanoulopoulos et al. 2014;Epitropakis et al. 2016;Caballero-García et al. 2018. The observer is stationary at 1000 r g from the black hole. The X-ray reprocessing by the disk is calculated via the reflionx model (George & Fabian 1991;Ross et al. 1999;Ross & Fabian 2005), where the photon index and the iron abundance are fixed at Γ = 2.4 and A Fe = 15, respectively (Chainakun et al. 2016). The disk extends from the inner-most stable circular orbit (ISCO) to 400 r g . The black hole spin is fixed at a = 0.998 Chiang et al. 2015;Jiang et al. 2018), so the ISCO is at ∼ 1.23 r g . The disk response function is produced by collecting the flux of the reflected photons as a function of time at the observer's sky.
Let us assume that Ψ(f, E j ) is the response function for the X-ray reverberation scheme, the observed PSD including reverberation effects can be computed via Uttley et al. (2014); Papadakis et al. (2016); Chainakun (2019); Chainakun et al. (2021b) as where P 0 (f ) is the intrinsic PSD profile expressed in eq. 1. Therefore, P 0 (f ) acts as a driving signal that is filtered by the response function, causing the observable dip and oscillatory structure in the P rev (f ). These echo features are more prominent in a more reflection-dominated band. We normalize the area under the response function to 1 and employ the reflection fraction R = (reflection flux)/(continuum + reflection) to reduce the effect of the energy band on the PSD profiles (e.g. Chainakun 2019), which is otherwise responsible for an energy-dependent amount of dilution applied to the reverberation calculations. This also accounts for the uncertainty in the variations of, e.g., continuum photon index and iron abundance among different individual observations which are not constrained here. Additional narrow features in the broadband PSD profiles were previously suggested in IRAS 13224-3809 (Alston et al. 2019), to resemble the observed PSD shape of very-high/intermediate state black hole X-ray binaries (Remillard & McClintock 2006). We model these narrow features using the Lorentzian function given in the form of: where N is a normalization factor, σ lor is the FWHM of the Lorentzian line and f lor is the centroid frequency of the line. The PSD data are fitted using the P rev model, or the (P rev + P lor ) model if the additional Lorentzian component is required by the data.

MODEL GRID AND FITTING PROCEDURE
The response functions corresponding to the source height h between 2-40 r g are simulated with equal spacing of the model grid of 2 r g . We vary the bending power-law index s in between 1-4 with the grid spacing of 0. For each observation, the source height h is tied between the soft (0.3-1 keV) and hard (1.2-5 keV) band data sets, while s, f b , R, and A are allowed to be free. Our parameters then consist of the source height, h, which is tied between Table 2. Best-fit parameters from the Prev model. The errors show 90% confidence intervals (∆χ 2 = 2.71). The uncertainty with '−' is quoted when the upper or lower limit cannot be estimated due to the finite model-grid extension. two energy bands, the PSD bending index (s s and s h ), bending frequency (f b,s and f b,h ), reflection fraction (R s and R h ), and the PSD normalization (A s and A h ) in the soft and hard bands for all the parameters. Different combinations of these parameter values represent different grid cells of the model. We perform the simultaneous fits of the soft and hard PSD data in isis (Houck & Denicola 2000) by stepping through the model grid cell. The χ 2 statistics are calculated using the subplex optimization method. The best-fit parameters are those that match to the grid cell with the lowest χ 2 value. Finer local grids are created if necessary, and the fitting is repeated to acquire the new best-fit parameter values.

RESULTS
Firstly, we used the reverberation PSD model (P rev ) to fit the individual PSD data of IRAS 13224-3809 without including the Lorentzian component. The reverberation effect imprinted in the PSD profiles can plausibly explain the oscillatory features seen in the data. Examples of the fitting results for certain individual observations are presented in Fig. 1. Some observations show clues of narrow features, but due to the poor quality of the data (low signal-to-noise), these narrow features cannot be robustly identified. In fact, even without the additional Lorentzian component, the model can provide good fits for the majority of the data, with adequate fits in some observations. We then select to exclude the Lorentzian component when we analyze the individual observations, and include it later when the data are combined using the luminosity bin.
The best-fit parameters for 16 observations of IRAS 13224-3809 are presented in Table 2. Note that, for each observation, the PSD profiles extracted in the soft and hard bands are simultaneously fitted, with the source height (h) being tied together. The errors correspond to 90% confidence intervals (∆χ 2 = 2.71) around the best-fitting values and, where necessary, are estimated by linear interpolation between the consecutive model-grid spacing for each parameter. The obtained reflection fraction in the soft band is larger than those in the hard band, which is consistent with the X-ray reflection framework, in which the reflection flux is more dominant in the soft band. We find the average values of f b,s ∼ 6.9 × 10 4 Hz, f b,h ∼ 9.2 × 10 4 Hz, s s ∼ 1.91 and s h ∼ 1.63. The source height found in these 16 observations varies between h ∼ 3-25 r g .
In Fig. 2, the source height is plotted against the count rate, as well as the fractional excess variance observed in the soft band (F var,s ). To elaborate more on the parameter correlations, we employ the photon index (Γ) from the average spectral model fits of Alston et al. (2020). The Γ values are represented by the size and color of the data points. The  data correlations and visualizations are analyzed using the Orange data mining platform in Python (Demšar et al. 2013). We find that the Pearson and Spearman correlation coefficients between h and count rate are r p = 0.23 and r s = 0.55, with the p-values of 0.39 and 0.03, respectively. This suggests that monotonic relationships between h and count rate are not statistically significant (p > 0.01). This is probably due to the small sample and large scatter. On the other hand, the variables h and F var,s are moderately anti-correlated in a non-linear way, with r s = −0.61 (p = 0.01). Fig. 3 represents the distribution of the parameters h, Γ, F var,s and F var,h where the observational data are divided into three different luminosity groups (i.e., low, medium, and high luminosity), as specified in Table 1. A layered kernel density estimate (KDE) is used to create the sample density distribution for each group. The observations displaying larger L seem to relate with larger source h and Γ, and smaller F var . Obviously, F var,s is correlated with F var,h . Furthermore, we perform both univariate and multivariate nonparametric two-sample test with bootstrap probabilities using Cramer-Test implemented in CRAN package (Feigelson & Babu 2012) in order to verify if the visual differences of the plots between 3 luminosity groups are significant. We produce 1,000 bootstrap-replicates with the normal Monte-Carlo-bootstrap method and set the confidence level of test to be 95%. Our hypothesis is that the samples in one luminosity group can be distributed as those in another group, which is accepted or rejected based on the estimated p-value. We find that the sample distributions in the F var,s -Γ space are clearly and significantly separated. The patterns of associations between h and F var in both energy bands also show a hint of separation of sources in different luminosity bins, but the differences are not significant. Fig. 4 shows the correlations between the constrained A and F var for both energy bands. We observe a high correlation between A and F var , with r p = 0.79 and 0.89 for the soft and hard bands, respectively, and r s = 0.90 for both energy bands. Furthermore, Fig. 5 shows the pairwise relationships between the PSD parameters of the soft band with the KDE density distribution, where the sources are divided into two groups: low and high A s . The multivariate two-sample test with bootstrap probabilities, as mentioned earlier, is also performed to identify whether the visual differences are significant. We can see clearly that the group of low A s most likely corresponds to a source with smaller s s and F var,s (blue profiles along the diagonal plots significantly shift to the left). However, the association of A s with Figure 3. Pairwise relationships between the variable h, Γ, Fvar,s and F var,h , with the data being divided into 3 groups: low, medium and high luminosity as noted in Table 1. The distributions are derived using 16 observations of IRAS 13224-3809, with a layered kernel density estimate (KDE). The diagonal plots represent a univariate marginal distribution of the data in each column. The labels AS, PS, and NS mean that the visual differences of all 3 grouped samples are significant, only pair of them is significant, and none of them is significant, respectively. See text for more details. f b,s cannot be clearly resolved (i.e. the base line of blue and orange profiles along the diagonal plots almost covers the same parameter space).
Furthermore, we illustrate how the overall parameters of IRAS 13224-3809 change as a heat map in Fig. 6, where the obtained values for each parameter are transformed to a scale from 0 to 1. This clearly reveals that the observation with higher count rate seems to have higher h and Γ, but smaller s, f b , A, and F var . This implies an increase of source height occurring with a decrease of the observed F var , with the flatter PSD slope. Scatter plots representing the overall correlations of the model parameters with respect to the count rate are also shown in Fig. 7. A hint of moderate-to-strong monotonic relationships between the count rate and h, Γ, F var,s , and A s (| r s |> 0.5) can be observed.
Due to the small sample, we also examine the distribution of correlation coefficients with paired-bootstrap resampling in order to construct the 95% confidence interval of the paired-sample statistics. This is to ensure that the correlation coefficients that show p < 0.01 in Fig. 7 are certain. The random sample is taken with replacement from our 16 samples to form a paired-bootstrap sample similar size to the original sample. This process is repeated 1,000 times to produce the paired-bootstrap distribution. For the count rate and Γ, we find that 0.34 ≤ r s ≤ 0.93, with 95% confidence. We also find that −0.93 ≤ r s ≤ −0.38 for the count rate and F var,s , and −0.85 ≤ r s ≤ −0.28 for the count rate and A s . The correlation coefficients then seem to be varied, but the trend in which each parameter is correlated or anti-correlated with count rate is quite certain despite of our limited sample. Now, we fit the PSD of IRAS 13224-3809 when the observational data are combined into three groups: low, medium and high luminosity. In this case the signal-to-noise ratio is relatively high compared to when we fit each individual  PSD data. Therefore, the additional Lorentzian component is included when necessary to explain the narrow features appearing in the profiles. The line's centroid frequency (f lor ) and FWHM (σ lor ) are allowed to be free. The best-fit results are presented in Fig. 8 and the parameters are shown in Table 3. The model can provide good fits for the low and medium luminosity spectra, with adequate fits for the high luminosity spectrum. Both energy bands show a hint of the narrow peak moving towards lower frequencies as the source luminosity increases, especially when considering the medium and high luminosity bins. The FWHM of the narrow line is σ lor ∼ 10 −4 Hz for both energy bands and all luminosity bin spectra, with the exception of the hard band of low luminosity data where the narrow feature is not significantly required.
The results from the luminosity bin spectra also show that the PSD shape flattens as energy increases (s s > s h ), but there is a less clear evolution of the bending PSD index with luminosity. Despite of this, the results still support the trend of increasing source height with increasing luminosity, from ∼ 3 r g to 20 r g , as predicted by the individual fitting. Figure 5. Pairwise relationships between the obtained source height (h), count rate, bending index (ss), bending frequency (f b,s ) and fractional excess variance (Fvar,s). The density distribution is produced via KDE, with two PSD normalization intervals: low and high As of < 20 and ≥ 20, respectively. The labels AS and NS mean that the visual differences are and are not significant, respectively . We investigate the corona evolution of IRAS 13224-3809 using 16 XMM-Newton observations solely through the PSD analysis. We take into account the reverberation effects caused by the lamp-post source, which represent an echo filter to the intrinsic variability signals (e.g. Papadakis et al. 2016;Chainakun 2019). There might be the systematic uncertainties due to the choices of Γ and i. While Alston et al. (2020) derived Γ from the reflection fits using the low-density disk, Jiang et al. (2022) considered a high-density model and found Γ is consistently higher. Note that Γ could change the flux contribution in each energy band, which should be compensated with allowing the reflection fraction to be free. Jiang et al. (2022) also reported i ∼ 60 • -70 • , while in this work we use the intermediate value of i = 45 • as in Caballero-García et al. (2020). These uncertainties should not have significant effects in our key results. Nevertheless, we notice that the reverberation signatures on the PSD slightly change with inclinations .

DISCUSSION
By fitting the soft and hard PSD simultaneously, we find that the source height h seems to be correlated with luminosity. This is consistent with Alston et al. (2020) and Caballero-García et al. (2020) who utilized different spectral and temporal profiles and reported a tendency of increasing source height with luminosity. Furthermore, Chainakun et al. (2022) investigated the correlations between the reverberating AGN parameters and reported an anticorrelation between the Fe-K lag and F var . Given that the amplitude of the lag (either in Fe-K or Fe-L band) increases with h, the hint of the anticorrelation between h and F var found here is well justified. This demonstrates that a consistent framework can be inferred even when diverse timing data and analyses are used. The source height in this work is varying between ∼ 3-25 r g . The spectral analysis by Jiang et al. (2022) using almost the same set of observations revealed that the source height h varied between 0.43-1.71 f AD /f INF r g , where f AD /f INF is the ratio between the coronal flux that reaches the accretion disk and the flux at infinity, which is geometry dependent. By assuming a lamp-post geometry and using the best-fit ionization and density parameters, they approximated the source location to be of 3-6 r g , which is much smaller than our analysis. This provides further evidence for the inconsistency in the parameters obtained from the timing (reverberation) data and the time-averaged spectral analysis, which was recently reported in several AGN such as NGC 5506 (Zoghbi et al. 2020) andMCG-5-23-16 (Zoghbi et al. 2021).
Meanwhile, the spectral components of IRAS 13224-3809 have been found to be complicated, particularly in the soft excess below 2 keV (e.g. Fabian et al. 2013;Chainakun et al. 2016). The time-averaged spectrum, for example, can be fitted with either the standard disk model with an additional soft-excess blackbody component (Jiang et al. 2018) or the reflection from a high density disk model (Jiang et al. 2022). Different spectral models result in different amount of dilutions applied to the reverberation lags (e.g. increasing the disk density results in an increase of reflection flux especially in the soft band of < 2 keV). Here, we do not determine which model spectrum is preferable over the others; instead, we choose to let the reflection fraction R to be free. There is also the change in flux contribution if we consider the returning radiation (Wilkins et al. 2020), or take into account the flux from the ultra-fast outflows (Parker et al. 2017(Parker et al. , 2021 that may be present. The parameter R then takes into account the impacts of the energy band chosen as well as incorporates the spectral complexity which is not explicitly modelled here. Nevertheless, the best-fit model prefers R s R h , which is consistent with the reflection framework that the reflection flux is mostly contributing in the soft band. The intrinsic shape of the PSD represents the variability power that tightly depends on the physical properties of the accretion flow in the framework of the propagation of mass accretion (e.g. Arévalo & Uttley 2006;Ingram & Done 2011;Mahmoud et al. 2019;Chainakun et al. 2021a;Ashton & Middleton 2022). An increase in the index γ of the disk emissivity, (r) = r −γ , and the disk parameter, α(H/R) 2 , where α is the viscosity parameter and H/R is the disk scale-height ratio, can produce more high-frequency variability power at higher energies. Therefore, varying γ and α(H/R) 2 is analogous to changing the bending power-law baseline model (e.g., bending index s and f b ).
Furthermore, our results show that the PSD shape flattens out as the energy increases (s s > s h ), indicating that harder X-rays vary more at higher frequencies. This is well consistent with what is commonly observed in AGN (Nandra & Papadakis 2001;Vaughan et al. 2003;McHardy et al. 2004McHardy et al. , 2005Papadakis 2004;González-Martín & Vaughan 2012;Ashton & Middleton 2022). We find that the PSD normalization in both energy bands (A s and A h ) are strongly correlated with F var (see Fig. 4). This is expected since the fractional excess variance is the integrated area under the PSD curve. We also find that the normalization is moderately correlated with the bending index s. The PSD shape then flattens not only with an increase in energies, but also with a decrease in A. Since s is also related to the accretion parameters such as α(H/R) 2 , this supports the framework that the PSD amplitude governed by A is dependent on the properties of the accretion rate (Georgakakis et al. 2021).
Piecing these results together (e.g., Fig. 6-7), we can infer that as the corona height increases, the observed luminosity increases, and the disk itself also evolves in such a way that produces more high-frequency variability power (PSD spectrum is more flat) but produces less overall variability power (smaller F var ). It is not straightforward how an increase in corona height, with increasing observed luminosity, induces more variability power at high frequencies. Perhaps, the corona and the disk variability evolve separately. Note that the correlation does not always imply causation. Therefore, it might be more intuitive to think that, when the corona height increases, the accretion disk itself becomes relatively more active at the inner regions so producing more power at high frequencies, rather than what it would be directly caused by the corona evolution. High luminosity 1.2-5 keV Figure 8. Data, model and residuals from fitting the (Prev + P lor ) model (red) to the data (blue) extracted in the 0.3-1 keV band (left panels) and 1.2-5 keV band (right panels) for different luminosity bins: low (top panels), medium (middle panels) and high (lower panels) luminosity. The obtained parameters are listed in Table 3. Emmanoulopoulos et al. (2016) fitted the reverberation echoes in the PSD data of AGN. By extracting IRAS 13224-3809 PSD in the 0.5-1 keV band using the combined XMM-Newton light curve, they found the bending frequency of ∼ 0.7 × 10 −4 Hz and the bending index of ∼ 2.26. Although the soft bands used here are slightly different, our average f b,s is comparable to Emmanoulopoulos et al. (2016), whilst the average s s is smaller. In fact, our f b,s falls within the bending frequencies reported by Alston et al. (2019), who used two bending power-law components to fit the PSD data of IRAS 13224-3809. The results, as well, agree with Alston et al. (2019) that a higher bending frequency is required for the harder energy band.
In addition, the constrained f b in both energy bands has a relatively weak correlation with the corona height. This suggests that the characteristic bend times-scales are less dependent on the coronal geometry. By assuming the IRAS 13224-3809 mass to be of 2 × 10 6 M , the f b obtained here is comparable (i.e., within the same order of magnitude) to the observation-based bending frequency that scales with the black hole mass (Papadakis 2004;González-Martín & Vaughan 2012).
Last but not least, the luminosity bin data show a hint of the Lorentzian peak shifting towards lower frequencies for higher luminosity, consistent with Alston et al. (2019). The FWHM of the line, when present, is ∼ 10 −4 Hz. The Lorentzian features may be present in some individual observations, but they cannot be easily distinguished from the oscillatory structure that can be described by reverberation. Chainakun et al. (2021b) found that the reverberation features in the AGN PSD profiles can potentially be retrieved using the machine learning (ML) technique, allowing the source height to be inferred accurately. Applying ML techniques to track changes in geometry of all current reverberating AGN in a systematic way using the PSD data is a subject of future work (Mankatwit et al., in preparatory).

CONCLUSION
By utilizing the PSD analysis of IRAS 13224-3809, we report the consistent trend of increasing lamp-post source height with increasing luminosity, and compared it with previous literature. We find that the PSD model that includes the effects of X-ray reverberation can reasonably explain the oscillatory structures seen in the PSD profiles. However, while our source height (∼ 3-25 r g ) is comparable to those observed in the time-lag data (Alston et al. 2020), it is substantially larger when compared with the source locations indirectly implied using the energy-integrated spectra (Jiang et al. 2022). This provides further evidence that the lamp-post parameters inferred from reverberation data are not consistent with those inferred from the time-averaged spectral analysis.
The model shows that when the corona height increases, the source luminosity increases as well, while the source spectrum tends to be softer. In addition, the observed fractional excess variance reduces. For the model to explain the data, a drop in F var requires a smaller PSD normalization (A) and a smaller PSD bending index (s). Furthermore, both A and s are certainly linked to the accretion phenomena. This means that as the lamp-post source moves away from the black hole, the accretion disk evolves, probably independently, in a manner that produces less X-ray variability power making PSD spectrum to flatten out towards high frequencies.
Also, the trend of the PSD slopes that flatten out as energy increases is clearly observed. The Lorentzian features in individual observations cannot be clearly distinguished after accounting for reverberation effects. We can only see a hint of shifting the Lorentzian peak to low frequencies for high luminosity using the luminosity-bin data. Interestingly, the FWHM of the line remains almost constant regardless of the luminosity bin. High quality data (i.e. high signal-to-noise ratio) and long observations will be required to place robust constraints on the evolution of the PSD parameters.