Coronal height constraint in IRAS 13224–3809 and 1H 0707–495 by the random forest regressor

X-ray reverberation in active galactic nucleus (AGN) occurs when the primary X-ray continuum from the corona illuminates an accretion disc, producing the reflected X-rays that are observed with a time delay after the direct continuum (see Uttley et al., 2014; Cackett, Bentz, & Kara, 2021, for a review). The reverberation lags due to the delay between the variations of the reflection and the direct continuum components were detected via the Fourier-based technique (e.g. Fabian et al., 2009; De Marco et al., 2013; Kara et al., 2016). The lag amplitude depends on the light travel distance between the corona and the disc, providing valuable information about the geometry of the system. Most previous work investigated X-ray reverberation signatures in the lag-frequency spectra using the standard lamp-post geometry where the corona is an isotropic point-like source located on the rotational axis of the black hole (Cackett et al., 2014; Emmanoulopoulos et al., 2014; Chainakun et al., 2016; Epitropakis et al., 2016; Caballero-García et al., 2018), while some investigated it in an extended corona environment (Wilkins et al., 2016; Chainakun & Young, 2017; Chainakun et al., 2019).

The AGN emission fluctuates rapidly on a wide range of timescales. The amount of the X-ray variability power that the AGN emits can be quantified using the power spectral density (PSD), whose intrinsic shape can also be imprinted with X-ray reverberation effects (Papadakis et al., 2016; Chainakun, 2019). The observed PSD profiles tend to be relatively flat at low Fourier frequencies, with a steep drop-off at higher frequencies (González-Martín & Vaughan, 2012). The intrinsic shape of the PSD depends on the specific characteristics of the accretion disc (Arévalo & Uttley, 2006; Ashton & Middleton, 2022). An attempt at probing the reverberation signature (appeared as an oscillatory structure) in the PSD data of AGN to constrain the coronal height were carried out by Emmanoulopoulos et al. (2016); Chainakun et al. (2022). In this way, the coronal geometry can be constrained, but high accuracy prediction may still require the intrinsic PSD shape to be known in advance.

Furthermore, Chainakun et al. (2021) developed a machine learning (ML) model using dictionary-learning and support vector machine to predict the discrete values of the coronal height. Interestingly, while the reverberation-imprinted and power-law based PSD profiles were used to train the model, it was suggested that the accuracy of predicting the source height was still high even when the model was tested with a bending power-law PSD. Here, we aim to further develop the model so that it can predict a continuous value of the coronal height and can directly be applied to real observational data. We use a random forest regressor (RFR) which is an ensemble ML algorithm that is composed of many individual decision trees being trained on subsets of the data. Each tree makes its individual prediction, and the final prediction is achieved by statistical averaging the predictions of all individual decision trees.

We select to apply our developed RFR model to the AGN IRAS 13224–3809 and 1H 0707–495, both of which are of particular interest due to their unique X-ray spectral properties and extreme variability. IRAS 13224–3809 is a type 1 AGN that is known to have a supermassive black hole at its center, with a mass of $\sim 2\times 10^{6}\ M_{\odot}$ (Alston et al., 2020). It has been thoroughly studied via both spectral and timing analysis (e.g. Fabian et al., 2013; Kara et al., 2013a; Chiang et al., 2015; Parker et al., 2017; Jiang et al., 2018). Along the line of these studies, the changing of the coronal height with luminosity was evidenced (Alston et al., 2020; Caballero-García et al., 2020; Chainakun et al., 2022), where the height increases from $\sim 3$–$5~{}r_{\rm g}$ to $\sim 10$–$25~{}r_{\rm g}$ ($r_{\rm g}=GM/c^{2}$, $G$ is the gravitational constant, $M$ is the black hole mass and $c$ is the speed of light) when the source changes from lowest to highest flux state. The model that simultaneously fitted the energy and lag spectra in various flux states also favoured a maximum spin black hole (Caballero-García et al., 2020). Meanwhile, the 1H 0707–495 corona geometry is very uncertain and under debate. While an analysis using the emissivity profiles suggested an extended corona for 1H 0707–495 (Wilkins & Fabian, 2012; Wilkins et al., 2014), a much more compact corona was obtained using a newly developed spectral model accounting for the spatial extent and the rotation of the corona (Szanecki et al., 2020). Both intrinsic and environmental absorption were also found when analysing the fractional excess variance of 1H 0707–495 (Parker et al., 2021).

This paper is organised as follows. The observations and data reduction are described in Section 2. How the PSD profiles are produced is explained in Section 3. In Section 4 we present the RFR algorithms used in this work and explain how the RFR model is developed. The results and analysis towards the performance of the model is given in Section 5, as well as the results when we apply the model to predict the coronal height in IRAS 13224–3809 and 1H 0707–495. The discussion and conclusion is provided in Section 6.

The X-ray data of the AGN IRAS 13224–3809 and 1H 0707–495 used in this study were previously observed by XMM-Newton observatory (Jansen et al., 2001) and were obtained from the XMM-Newton Science Archive.¹¹1Available online at http://nxsa.esac.esa.int/ All observational data are shown in Table 1. Here, only pn data are used to obtain high quality data. To analyse the data, we created the PSDs from these observations following the method and criteria explained in Chainakun et al. (2022). In brief, the pn data were reprocessed using the Science Analysis Software (SAS) following the procedure suggested in XMM-Newton Science Analysis System’s user guide.²²2https://xmm-tools.cosmos.esa.int/external/xmm_user_support/documentation/sas_usg/USG.pdf The observational periods which were affected by high background flaring activity were also removed during these steps. The remaining useful exposure time of each observation after removing the background flaring events is shown in column 4 of Table 1. Note that, in this step, we only considered the observations with a remaining exposure time of $\ga$30 ks so that the signal to noise of their PSDs is sufficiently high.

We then extracted the background-subtracted light curves of the AGN in the soft (0.3–1 keV) and hard (1.2–5 keV) energy bands – defined as reflection dominated and continuum dominated energy bands, respectively – from the observational data flagged with PATTERN $\leqslant$ 4 and #XMMEA_EP for each observation. The corresponding PSDs were derived from the light curves which were divided into the length of $\sim$20 ks and re-binned to have a time resolution of 179 s. To improve the data signal to noise, we also re-binned the PSDs in which every bin width is equal in logarithmic space of base 1.06. These PSDs were then used as a representative of the real AGN sample in the RFR analysis in Section 5.2.

The variability power of the AGN on different timescales can be represented using the PSD whose intrinsic shape follows a broadband power-law profile usually simplified in the form of

where $f$ is the Fourier frequency and $\gamma$ is the power spectral index.

However, some characteristic features such as the dip and oscillatory structures caused by the X-ray reverberation can also be imprinted on the PSD profiles. We can assume that the total light curve includes both direct continuum and reverberation flux, and the total response function of the system can be written in the form of (e.g. Papadakis et al., 2016)

where $\delta(t)$ is a delta function representing the continuum flux response and $\psi(t)$ is the disc response function representing the amount of the reflection flux detected by a distant observer, as a function of time, after the primary continuum flux. The precise $\psi(t)$ is obtained via the ray-tracing simulations (e.g. Cackett et al., 2014; Chainakun et al., 2016; Wilkins et al., 2016; Epitropakis et al., 2016). Here, we use the kynxilrev model Caballero-García et al. (2018, 2020) to compute the disc response functions under the lamp-post geometry where an accretion disc is optically thick and geometrically thin and is illuminated by a point-like corona located on the rotational axis of the central black hole. The area under the disc response function is normalized to 1, so that the reflection fraction $R_{F}$ determines the ratio of (reflection flux / continuum flux).

Based on the convolution theorem, the PSD including the X-ray reverberation effects can be obtained by Uttley et al. (2014); Papadakis et al. (2016); Chainakun et al. (2022)

where $\Psi(f)$ is the Fourier form of the total response function. The intrinsic $P_{0}(f)$ is then an input signal which is processed by the total response function resulting in a reverberation feature appearing as an oscillatory structure in the $P_{\text{obs}}(f)$ (Papadakis et al., 2016; Chainakun, 2019).

By using the kynxilrev model (Caballero-García et al., 2018, 2020) to produce the disc response functions, the continuum spectrum is represented by a cut-off power-law: $F(E)\propto E^{-\Gamma}E^{-E/E_{c}}$, where $E_{c}=300$  keV. The disc extends from the innermost stable circular orbit (ISCO) to 1000 $r_{\rm g}$. We initially fix the black hole spin to be $a=1$. Other parameters in the kynxilrev model, if not specified, are set to their default values.

Examples of the simulated PSD profiles are presented in Fig. 1. In this illustration, we assume the black hole mass of $10^{6}~{}M_{\odot}$ and the inclination $i=45^{\circ}$. The orange lines in the left panels represent the PSD that already include the reverberation effects due to the lamp-post corona. The main reverberation dip, as expected, shifts towards higher frequencies (shorter timescales) for lower coronal height. The corresponding noisy PSD profiles (blue lines) are simulated using the Timmer & Koenig (1995) method, then are binned into $\sim 30$ frequency bins within the frequency range of $\sim 3\times 10^{-5}$–$3\times 10^{-3}$ Hz (red points in the right panel), similar to what is normally extracted from the XMM-Newton observations of IRAS 13224–3809. After being binned, these noisy PSD data are then used to train the machine, as their shapes and frequency bins are quite close to the observational ones.

The machine learning (ML) algorithm used here is the random forest regressor (RFR). It aims to probe the characteristic features of the reverberation signatures that appear on the PSD profiles in order to predict the disc corona geometry, such as the coronal height and the disc inclination. We follow the standard procedure by splitting the PSD data into 80/20 for training/testing. All ML algorithms and related functions for the RFR model development are adopted from scikit-learn³³3https://scikit-learn.org/ (Pedregosa et al., 2012).

Firstly, we investigate the effectiveness of the RFR model in estimating three physical parameters of AGN including the source height, inclination, and black hole mass. Our results reveal that the accuracy for forecasting the inclination is poor ($R^{2}<0.5$) because either the reverberation feature is indistinct to be extracted or the model requires a larger number of PSDs varying with inclination to be trained. Similarly, the accuracy in determining the black hole mass simultaneously with the coronal height is poor as the characteristics of PSDs varying with the mass can resemble those changing with the source height (Papadakis et al., 2016). For example, the data with smaller mass and larger height can also be reasonably explained by the models with larger mass but smaller height. The degeneracy between the mass and the source height is also found in other timing profiles such as the lag-frequency spectrum (e.g. Cackett et al., 2014; Caballero-García et al., 2018, 2020), and still cannot be easily broken by the model.

In the case of predicting the source height alone, the obtained accuracy is excellent ($R^{2}\gtrsim 0.9$), particularly when the reflection fraction $R_{F}\geqslant 1$. This is because the reverberation feature in the PSD, such as the dip, is more pronounced. As a result, we focus on the PSD data extracted in the energy band of $0.3$–1 keV where the reflection flux is usually dominated and choose to predict the source height by fixing the black hole mass and the inclination angle.

Chainakun et al. (2021) proved that the ML technique can be used to extract the reverberation signatures on the PSD profiles to predict the source height. However, for the ML approach, there is for a classification problem where the source heights are categorized into several distinct classes (i.e. discrete values), each of which represents a range of source heights. Here, we use the RFR algorithms so that the source height can be predicted as continuous values. The simulated PSD profiles being used to train, cross-validate and test the machine are also specifically produced using a similar frequency range and a similar number of frequency bins to each observational data. Higher Poisson noise can lead to higher standard deviation of the measured lag, without affecting the lag value (Kara et al., 2013b). Rather than quantifying the Poisson noise level in each observation and involving the errors from the PSD measurement, we select to directly remove the Poisson-noise dominated bins instead. This is acceptable as long as the training and test data are processed in the same way. The prediction accuracy of $R^{2}>0.9$ can still be achieved if the number of Poisson-noise dominated bins being removed is less than 10, which is likely the case for our observational data.

Once the developed RFR model is applied to observational data, we find the IRAS 13224–3809 corona varies between $\sim 5-18~{}r_{\rm g}$. This is comparable with the constrained coronal height of $\sim 6-20~{}r_{\rm g}$ using the lag-frequency spectra (Alston et al., 2020). The combined spectral-timing analysis by Caballero-García et al. (2020) also reveals the coronal height varying between $\sim 3-10^{+10}_{-1}~{}r_{\rm g}$, assuming a maximally spinning black hole. Recently, tracing the corona evolution from the reverberation signatures that appear in the PSD of IRAS 13224–3809, but with conventional fitting technique, suggests the varying corona located at the height of $\sim 3-25~{}r_{\rm g}$ (Chainakun et al., 2022). We also find a moderate monotonic correlation showing the trend of increasing source height with the count rate, $r_{S}\sim 0.55$ ($p\sim 0.03$), similar to Alston et al. (2020); Caballero-García et al. (2020); Chainakun et al. (2022). This is because the photons emitted from the corona located further away from the black hole can escape the strong gravitational field more easily than those emitted from the corona closer to the black hole, leading to higher source luminosity. The corona evolution in IRAS 13224–3809 then is quite consistent among these studies, even though different methods are used.

Here, we fix the IRAS 13224–3809 mass to be $2\times 10^{6}~{}M_{\odot}$ (Alston et al., 2020) and employ the intermediate inclination of $45^{\circ}$ (Caballero-García et al., 2020). The result suggests that high accuracy can still be obtained as long as the fixed mass is $\lesssim 10$ per cent different from the true value. The spectral fits by Jiang et al. (2022) using a high-density disc model suggested the inclination of $60^{\circ}$–$70^{\circ}$. The choice of inclination, however, should not significantly affect the accuracy in predicting the coronal height since the reverberation features in timing profiles are less dependent with the inclination (e.g. Cackett et al., 2014; Papadakis et al., 2016). Recently, Hancock, Young, & Chainakun (2023) fitted an extended corona model, simplified using dual lamp-post sources, to the lag-frequency spectra of IRAS 13224–3809. By adopting the inclination $i=64^{\circ}$ as reported by Fabian et al. (2013) and placing the first, lower source at $2~{}r_{\rm g}$, they found that the second, upper source could vary between $\sim 3-20~{}r_{\rm g}$. Although different coronal geometry and inclination are assumed, the disc-corona distance constrained here is comparable with the upper source height suggested in Hancock, Young, & Chainakun (2023). Developing the RFR model to probe the spatial extent of the corona is possible, but it is a subject for future research.

Recently, Chainakun et al. (2023) employed the Granger-causality test to probe the intrinsic reverberation lags in IRAS 13224–3809 to be $\sim 200$–500 s. They directly converted the intrinsic lags to the true light-travel distance of $\sim 20$–50 $r_{\rm g}$, and by assuming a face-on disc the coronal height could be roughly approximated to be $\sim 10$–25 $r_{\rm g}$. Furthermore, they reported the evidence of the corona variability within some individual observations (e.g. the corona can change from $\sim 10$–25 $r_{\rm g}$ to $\sim 55$ $r_{\rm g}$, or vice versa, from the beginning to the end of the observation). The method outlined in Chainakun et al. (2023) is a statistical test and the large coronal height (e.g. $>25$ $r_{\rm g}$) is seen only in some partial segments of the light curves. The different lags from different segments can lead to an insignificance of the lags at a particular value when analysing the full light curve, or the multiple lags can be averaged out. However, Chainakun et al. (2023) still found, for the majority of the IRAS 13224–3809 observations, that the corona situated at $\sim 10$–25 $r_{\rm g}$, so with the measurement error of $\sim 5$ $r_{\rm g}$ suggested by their method, the coronal height constrained there is mostly comparable to what obtained in this work.

For 1H 0707–495, the RFR model places the corona at $h\sim 5-18~{}r_{\rm g}$ as well. The hint of increasing source height with count rate is also observed. Moreover, Chainakun et al. (2021) suggested that even if the ML is trained using the power-law PSDs but the observed data have a bending power-law shape, the coronal height can still be forecasted with the accuracy of $\gtrsim 0.75$. Then, the obtained accuracy reported here can be lower if the intrinsic shape of the PSD, which is not known prior, has a bending power-law form. Perhaps, the best way might be to train the machine with the PSD with various intrinsic shapes as much as possible, but this will progressively increase computational time consuming and is beyond the scope of this work. Also, when the range and number of the photon indices are changed, the regression results show slight modifications (see Test A and Test B in Fig. 8). This effect, however, should not change the overall trend of the correlation if the range of the photon indices are large enough.

Nevertheless, the nature of 1H 0707–495 is quite uncertain. Some suggested that it should have a low-spin black hole (Done & Jin, 2016), while the mass can possibly be $\sim 5-10\times 10^{6}~{}M_{\odot}$ (e.g. Zoghbi, Uttley, & Fabian, 2011; Pan et al., 2016), which is larger than the mass assumed in, e.g., this work and Szanecki et al. (2020) by a factor of $\sim 2$–5. Szanecki et al. (2020) analysed the relativistic reflection spectra from the disc around a maximally spinning black hole and found the corona is very compact, with an extended size no larger than $\sim 1~{}r_{\rm g}$ and located at $h\sim 3~{}r_{\rm g}$. A compact corona at $\sim 3$–$4~{}r_{\rm g}$ was also proposed by Dauser et al. (2012) and Caballero-García et al. (2018). Contrarily, Wilkins & Fabian (2012) studied the emissivity profiles due to various coronal geometries along with previous constraints from reverberation lags, and inferred a presence of the 1H 0707–495 coronal height at $\sim 2~{}r_{\rm g}$ that radially extends outwards to a radius of $\sim 30~{}r_{\rm g}$. The corona also seems to expand as the source luminosity increases (Wilkins et al., 2014). Dovčiak & Done (2016) suggested that the 1H 0707–495 corona must be large enough to intercept sufficient disc photons to produce the hard X-ray continuum and, during the maximum flux, it could be at $\sim 2-5~{}r_{\rm g}$ extending further outside to $\sim 20~{}r_{\rm g}$. Furthermore, Chainakun et al. (2019) developed a spherical corona model to explain the frequency-dependent time lags of AGN, and found that the 1H 0707–495 corona could be radially extended up to $\sim 10~{}r_{\rm g}$. Previous study using a dual lamp-post model also suggested the 1H 0707–495 corona could extend upwards, up to $\sim 20~{}r_{\rm g}$ (Hancock, Young, & Chainakun, 2023). Although this work assumes a lamp-post geometry, it suggests a dynamic corona of 1H 0707–495 that can move far beyond a few gravitational radii from the central black hole, and an increase of the coronal height with luminosity is also evidenced.

In conclusion, our study highlights the potential of machine learning, specifically regression techniques, for probing the coronal evolution in AGN. The constrained $h\sim 5-18~{}r_{\rm g}$ in both AGN here supports the height-changing corona and the light-bending scenario (Miniutti & Fabian, 2004), in which the coronal height is correlated to source luminosity. The fluctuations in the gravitational and magnetic fields caused by the black hole and the accretion disc can induce the corona to move closer or farther away from the black hole, leading to variations in the X-ray emission. The results are also in line with recent studies on other AGN such as NGC 5548 that, by probing the UV/optical PSD variation, revealed a dynamic corona chainging its position and probably also its spatial extension on various timescales (Panagiotou et al., 2022).

We thank the anonymous referee for comments that improved the paper. This research has received funding support from the NSRF via the Program Management Unit for Human Resources & Institutional Development, Research and Innovation (grant number B16F640076). N.M. thanks Suranaree University of Technology (SUT) and National Astronomical Research Institute of Thailand (NARIT) for the financial support.

The observational data were accessed from XMM-Newton Observatory (http://nxsa.esac.esa.int). The response function data analysed here were generated using the kynxilrev model available in https://projects.asu.cas.cz/stronggravity/kynreverb. The code for the developed RFR model to predict the coronal height is available in https://github.com/PChainakun/RanForPSD. Other derived data underlying this article will be shared on reasonable request to the corresponding author.