Extended corona models of X-ray reverberation in the AGN 1H 0707−495 and IRAS 13224−3809 Open Access

ABSTRACT

1 INTRODUCTION

The spectra of active galactic nuclei (AGNs) have a component of direct X-ray continuum emission, and that cold gas can reflect some of this X-ray continuum (Pounds et al. 1990). The cold, optically thick material seen through fluorescence and reflection is known to occur in the presence of an accretion disc where X-rays produced by inverse Compton scattering in a corona are reflected off the accretion disc producing a modified spectrum emission with fluorescent Fe K lines at 6.4 keV and other spectral features (George & Fabian 1991). The first hints of reflection or reprocessing time-delay due to the light travel time between the corona and disc were seen in XMM–Newton observations of Ark564 (McHardy et al. 2007), which led to the first robustly discovered delays in 1H0707−495 (Fabian et al. 2009) where the soft energy band (0.3–1 keV) lagged behind the hard band (1–4 keV) by 30 s. Many reverberation lags have since been discovered (e.g. Emmanoulopoulos, McHardy & Papadakis 2011; Kara et al. 2013b; Zoghbi et al. 2013).

The soft negative time lag is the signature of relativistic reflection that reverberates in response to continuum fluctuations and is interpreted as the light crossing time from the source to the reflecting region, correlating positively with the black hole mass (De Marco et al. 2013). Hard positive lags at lower frequencies are understood to originate from fluctuations of the accretion rate propagating from outer to inner radii, causing the hard X-rays produced at smaller radii to respond after soft X-rays produced at larger radii (e.g. Kotov, Churazov & Gilfanov 2001; Arévalo & Uttley 2006; Arévalo, McHardy & Summons 2008). Hard X-ray lags (i.e. hard photon variability lagging soft photon variability) were evident in X-ray binaries before they were discovered in AGNs (see e.g. Miyamoto et al. 1988; Nowak et al. 1999).

The reproduction of soft reverberation lags from a compact corona and the hard lags produced by propagating fluctuations through an extended region whilst maintaining the features of the energy dependence seen in the Fe K α line region is challenging and computationally intensive. Motivated by these phenomena, Chainakun & Young (2017), CY17 hereafter, developed an approximation of a vertically extended source using two X-ray point sources located on the rotation axis of the black hole. The model employs gravitational units where the gravitational radius r_g = GM/c² and gravitational time t_g = GM/c³ (where G is the gravitational constant, M is the black hole mass, and c is the speed of light). Instead of modelling the propagating fluctuations, a phenomenological function of expected source responses which react to these propagations is employed. The two X-ray sources are allowed to vary as they respond to primary intrinsic variations. The X-ray continuum variability depends on the response of the source and the X-ray reflection depends also on the disc response. Therefore, the model has to predict the time lags from both continuum X-ray sources and the associated disc responses. The model was fitted to a 120 ks timing XMM–Newton observation of the narrow-line Seyfert 1 galaxy PG 1244+026, examining the frequency and energy where the lags were found. The model revealed hard and soft X-ray sources at heights ∼6 r_g and ∼11 r_g, respectively, with the upper source producing small amounts of reflection which suggested a feasible geometry of a relativistic jet beaming away from the disc. CY17 suggested that the continuum flux from the upper source and the extra blackbody component that contribute significant flux at energies <1 keV are required to dilute the soft reverberation lags and to reproduce the absence of soft lag in the lag-energy spectrum of PG 1244+026.

The X-ray variability in 1H 0707−495 was found to be extreme and may have both intrinsic and environmental absorption origins (Parker et al. 2021). The geometry of 1H 0707−495 has been discussed by Szanecki et al. (2020) who developed a new extended lamp-post model which accounted for the spatial extent and rotation of the X-ray source. The investigation of the location and size of the corona indicated a compact corona at ∼2 r_g (highly centrally peaked rather than extended) regardless of the effects of ionized absorption from winds. They found no evidence that the size of the corona was correlated to the luminosity as reported by Wilkins et al. (2014). The time lags in IRAS 13224−3809 have been discussed by Alston et al. (2020) who studied short-time-scale variations and included all relativistic effects allowing for ultrafast outflows fitting multiple epochs where the source height changed. They tackled inherent degeneracies between the reverberation signal and black hole mass, inner disc radius, and height of the corona by tracking short-scale reverberation signatures where the source height changed. They found that the height of the corona increased with increasing luminosity and that black hole mass uncertainty estimates were comparable to the leading optical reverberation method by Peterson et al. (2004). In addition, Caballero-García et al. (2020) used all available data from the XMM–Netwon archive to simultaneously fit various flux states of the energy spectra and time lags using a new code which calculates reflection spectra from the accretion disc in response to an X-ray flare from a point source located above the black hole in accretion disc lamp-post geometry. The model strongly favoured a maximally spinning black hole and detected significant variations of the corona height, increasing from 3 − 5 r_g at lower flux states and extending to ∼10 − 20 r_g when the luminosity doubled. Recently, Chainakun et al. (2022b) constrained the reverberation signatures that appeared in the power spectral density of IRAS 13224−3809 and found that the lamp-post source height increased from ∼3 r_g to ∼25 r_g with the luminosity.

This study builds on the timing and spectral analysis of reverberating AGNs reported by Hancock, Young & Chainakun (2022), HYC22 hereafter, using a variety of similar spectral flux levels as found in the spectra of each source. In addition we explore all suitable individual XMM–Newton observations of 1H 0707−495 and IRAS 13224−3809 between 2000 and 2016. Our aim is to develop the model initially created by CY17 to explore these time lag signatures in an extended corona scenario. These well-studied AGNs are selected due to their abundance of long observations, however, it should be noted that these are extreme narrow-line Seyfert 1 galaxies.

2 OBSERVATIONS AND DATA REDUCTION

The data for all observations outlined in Table 1 were downloaded from the XMM–Newton archive and processed using standard methods. The time lag estimates calculated between the soft (0.3–0.8 keV) and hard (1–4 keV) energy bands reported by HYC22 have been used throughout this study.

We initially assume 1H 0707−495 to have a black hole mass log M = 6.31 M_⊙ (Bian & Zhao 2003) and redshift z = 0.0411 (Leighly 1999). The luminosity distance from the NED data base is 187 Mpc. This AGN has been well documented to be dominated by relativistically blurred reflection at either a low or a moderate inclination angle (see e.g. Fabian et al. 2009; Zoghbi et al. 2010; Dauser et al. 2012; Kara et al. 2013a; Caballero-García et al. 2018). We employ the inclination angle i = 53° as derived from the emissivity profile by Wilkins & Fabian (2011). For IRAS 13224−3809, we assume a black hole mass of log M = 6.30 M_⊙ (Alston et al. 2018, 2020; Caballero-García et al. 2020). We also adopt appropriate values for redshift z = 0.0406 and inclination i = 64° as reported by Fabian et al. (2013). The luminosity distance from the NED data base is 310 Mpc.

3 THE EXTENDED CORONA MODEL (ECM)

The extended corona model (ECM) assumes a standard geometrically thin, optically thick accretion disc (Shakura & Sunyaev 1973) which extends from the innermost stable circular orbit (ISCO), or the radius of marginal stability r_ms, to 400 r_g around a central black hole. Note that accreting supermassive black holes were mostly found to be rapidly spinning (Reynolds 2021), therefore to limit the number of free parameters and to avoid the model degeneracy, we fix the black hole spin to be a = 0.998. The accretion disc is illuminated by two compact X-ray sources located on the symmetry axis at heights h₁ and h₂ which are the lower and upper source heights, respectively, whose amplitudes as a function of time are x₁(t) and x₂(t), respectively, as in CY17. We maintain the two source to investigate the extent of the corona and a basic sketch of this scenario is shown in Fig. 1. Physically, the lower X-ray source may represent the base of a jet-like structure or the lower region of a compact corona (Wilkins & Gallo 2015) and the upper source represents the farthest region from where a flare or response from the disc is detected hence it is interpreted as the upper extreme of the corona. A further plausible explanation of this region could be due to a periodic vertical collimation of the corona as a jet launching event subsides (Wilkins et al. 2015). These explanations also suggest that photon emission is beamed vertically away from the disc.

Figure 1.

A simple sketch of the extended corona scenario showing the two X-ray sources located on the rotation axis of an accreting black hole. The arrows show the trajectories of the continuum and accretion disc reflected photons.

Open in new tabDownload slide

We began by creating the disc model by integrating photon paths from the observer to the disc assuming a maximum black hole spin where a = 0.998 with the disc extending from ISCO out to 400r_g. The next task was to integrate photon paths from a single X-ray source to the disc, after which the spectra was calculated using reflionx for a given source height and disc inclination already computed in the previous steps. The final step was to calculate the reverberation signatures by looking at the already computed source-to-disc and observer-to-disc ray tracing runs. The disc response function, ψ_i(E_j, i), was computed and the variability of the disc reflection can be calculated via a convolution term: x_i(E_j, t)⊗ψ_i(E_j, i). The observed X-ray variability then can be written as the sum of the source and the disc variability. The lag-frequency spectra were calculated from the Fourier-phase differences between two energy bands, following the standard techniques (e.g. Cackett et al. 2014; Emmanoulopoulos et al. 2014; Chainakun & Young 2015).

Note that the light curve in each energy band always contains both continuum and reflection components. The contamination of the continuum flux in the reflection-dominated band, and the reflection flux in the continuum-dominated band cause dilution effects which can reduce the lag amplitude, meaning that the measured time lags are shorter than the intrinsic time lags. The dilution effects modify the shape of the lag-frequency spectrum without affecting the frequency at which the lags occur. For discussions on dilution, see e.g. Uttley et al. (2014), Kara et al. (2014), and Chainakun & Young (2015). In order to deal with dilution effects, the variability of the reflected photons of the disc is normalized using the reflected response fraction (defined as the reflected flux/continuum flux) of all energy bands E_j. We include a brightness parameter b (as measured in the frame of the observer) as a ratio of the brightness of the lower and upper source, b₁ and b₂, respectively, by defining b = b₂/b₁ and fix b₁ = 1 and allow b to vary between 1 and 3. The continuum flux of the lower source will be less than the upper source due to its closer proximity to the black hole, therefore photons will be subject to light bending effects towards the centre (Miniutti & Fabian 2004).

The model is capable of reproducing the prominent time lag features seen in AGNs when Γ₁ ≠ Γ₂ as seen in the top panel of Fig. 2. In addition Γ₁ > Γ₂ for positive low freqency lags and Γ₂ > Γ₁ for negaitive low frequency lags. As t_max increases, so does the amplitude of the hard lags as seen in the bottom panel of Fig. 2 and aliasing effects (phase wrapping) are moved to lower frequencies and positive lags will switch to negative lags. The low frequency hard lags increase with energy and the softer reverberation lags dominate at higher frequencies, consistent with the ‘two-mechanism’ features of propagation and reflection and was also consistent with the traditional spectral features that are widely explained by reflection from the inner disc.

Figure 2.

The frequency dependent lags varying with Γ and t_max. Both panels were set at i = 53°, h₁ = 2, h₂ = 3, Fe = 1.0, log ξ = 1.0, and q₂ = 0.5. The top panel shows the frequency dependent lags varying with Γ₁ and Γ₂. The positive low frequency fluctuation lags and the negative soft reverberation lags are produced when Γ₁ ≠ Γ₂. The bottom panel shows lag behaviour with different values of t_max, where Γ₁ = 2.0 and Γ₂ = 1.5, b = 2.0, and t_shift = 10.

Open in new tabDownload slide

We generated the initial model by running the code for a course range of parameter values to fit the lag-frequency spectra of the fully combined observations 1H0707−495 using an inclination angle i = 53°. We set the density profile to constant (p = 0) and the response decay function q₂ = 0.5 and iron abundance Z_Fe = 1.0. The lower source height was fixed at 2 r_g and upper source height was allowed to vary between 3 and 20 r_g thus obtaining the extent of the corona. While fixing the lower source at 2 r_g is a pragmatic choice that simplifies the fitting and limits the number of free parameters, having one source this close to the horizon is comparable to the limits by previous spectral or timing modelling, (e.g. Kara et al. 2013a; Caballero-García et al. 2018, 2020). In addition, Γ₁ and Γ₂ were allowed to vary between 1.4 and 3.3 in line with reflionx, the ionization parameter specified at the ISCO could vary from 0 to 3 for neutral and highly ionized exploration respectively. The brightness parameter was also allowed to vary from 1 to 3. The remaining free parameters were t_max and t_shift and we initially set these ranging from 50 to 600 t_g and 0 to 120 t_g, respectively. The black hole mass for each source was initially tested using sensible values from literature (Bian & Zhao 2003; Alston et al. 2020) and allowed to vary during the fitting procedure.

4 RESULTS

4.1 Model testing on 1H 0707−495

The ECM fitting was developed using the Interactive Spectral Interpretation System (isis) Version 1.6.2-27. Generally, much better fits were obtained using single observations than combining the observations into groupings. Although spectral combinations are useful for obtaining snapshots of the various epochs, the true variability is best revealed by examining each observation in turn, however this is time-consuming and computationally expensive. First, we independently examined the posterior distribution for the mass using combined observation best fit. For this inspection we use the isis_emcee module outlined at Remeis-Wiki (2018) by adopting the methods of the mcmc hammer (Goodman & Weare 2010; Foreman-Mackey et al. 2013). We initiate 500 walkers and run the chain with 10 000 iterations. The walkers converged tightly after 1000 steps and settled into an acceptance rate of ∼0.7 after ∼3000 iterations to obtain a mass of log (M/M_⊙) = 6.22 ± 0.01, where the errors are calculated at the 90 per cent confidence limit. The posterior density and contour plot is shown in the top row of Fig. 3. This value is comparable to that reported in the literature (e.g. Zhou & Wang 2005; Zoghbi et al. 2010) although uncertainties in the estimations have been acknowledged. This initial model fit was reasonable where the χ² was 2.32 and the upper source was located at |$20.00^{+0.00}_{-0.11}\ r_g$|⁠.

Figure 3.

Independent (emcee) tests of the posterior density for the combined observations of 1H 0707−495 (top row) with mean mass μ = 6.22 (solid blue line) and σ = 0.01 (black dashed lines). IRAS 13224−3809 is also shown (bottom row) with mean mass μ = 6.34 and σ  = 0.01. The contour plots of the mass and normalization are shown in the right column. All x-axes are identical for comparison.

Open in new tabDownload slide

For this mass and geometry, the photon indices Γ₁ and Γ₂ were well constrained at |$2.6^{+0.01}_{-0.07}$| and |$2.1^{+0.72}_{-0.01}$|⁠, respectively, where the ionization log ξ ∼ 3.0. The upper source was exactly twice as bright as the lowers source given by parameter b ∼ 2. The source response time from disc fluctuations |$t_{\rm{max}} = 250.00^{+67.82}_{-0.04}\ t_g$| and the time taken for the fluctuations to propagate from the lower to upper source |$t_{\rm{shift}} = 20.00^{+11.05}_{-11.92}\ t_g$|⁠. Although this initial fit to the combined data is promising and can provide a good description of the data, we acknowledge the reasonably large errors returned for t_max, and the loosely constrained errors found for t_shift. In general, this is typical of the fits found for all data.

Further complexity was revealed when attempting to standardize the mass value for the remainder of the data to the value obtained from model fitting via the University of Bristol high performance computer BlueCrystal. A fixed mass of log (M/M_⊙) = 6.22 generally returned much poorer fits and the fitting procedure remained computationally intensive.

To develop this model further, the reverberation signatures were read into a table model that can be used to fit the data in isis and xspec. This step required the use of astropy and heasp, the latter of which is a component of the standard heasoft installation (NASA 2020). The table model was initially generated using python and the parameter space was given finer step sizes between parameter values to obtain a finer grid. The output model was huge (>40 GB) and the model loading time was too long for general usage, therefore the parameter space was constrained further. The final table model was 19 GB and the parameter space is outlined in Table 2. Although this is still very cumbersome, further constraints at this stage were halted due to dynamic benefits of the model parameter space.

The first table model fit was conducted on the combined data for 1H 0707−495, achieving a very good model fit to the data where χ² = 0.50.¹ This fitting was applied to all of the combination flux groups for low, medium, high etc and the results are presented in Table A1. The combined data and model fit along with the lag-frequency spectrum file predicted by the model is presented in the top panel of Fig. 4. Of course, a full table model listing every possible integer value of parameter space would be extremely difficult to achieve given the intensive computational power required, so only the closest lag-frequency, that is the predicted reverberation signatures are shown unbinned at full resolution by the solid grey wavy lines.

Figure 4.

The combined lag-frequency Table model fits for 1H 0707−495 (top panel) and IRAS 13224−3809 (bottom panel) showing the best model fit (red) where the mass is variable. The closest reverberation signature is shown by the grey wavy line. Note that the lags are estimated between the soft, 0.3–0.8 keV, and hard, 1–4 keV, bands. Furthermore, there is significant cancellation of positive and negative lags due to phase wrapping and binning.

Open in new tabDownload slide

4.2 Model testing on IRAS 13224−3809

Once again we explore the initial model fit using isis_emcee with 500 walkers and 10 000 iterations and walkers converged tightly after only ∼200 steps, settling into an acceptance rate of ∼0.7 after ∼300 iterations. We obtained a mass of log (M/M_⊙) = 6.34 ± 0.01 for IRAS 13224−3809, where the errors are calculated at the 90 per cent confidence limit. Again we note the model capability of returning a mass value closely comparable to log (M/M_⊙) ∼ 6.3 as reported in, e.g. Alston et al. (2020). The mass posterior density and contour plot, along with the model fit is shown in bottom row of Fig. 3. The ECM fit was statistically good where χ²  = 1.47 with the upper source located at |$5.00^{+5.52}_{-0.02}\ r_g$|⁠. For this mass and geometry, the photon indices Γ₁ and Γ₂ were well constrained at |$2.60^{+0.02}_{-0.59}$| and 1.90 ± 0.01, respectively, where the ionization log ξ ∼ 1.0. The upper source was about 3 times as bright as the lowers source (b ∼ 3). The source response time from disc fluctuations |$t_{\rm{max}} = 1300.00^{+0.00}_{-304.36}\ t_g$|⁠. Also note that the ranges of the parameters t_max and t_shift required adjusting up to 1500 t_g and 150 t_g, respectively, for this source as summarized in Table 2. The model was able to provide a good statistical description of all data when fitted via the table model. The best-fitting result for the combined observations is presented in the bottom panel of Fig. 4. Note that there is significant cancellation of positive and negative lags dues to phase wrapping and binning.

4.3 Constraining the black hole mass

Before exploring the simultaneously fitted ECM model, we inspected the variable mass model fits as a function of the negative (reverberation) time lags. For each source this relationship is shown in Fig. 5, suggesting that whilst the black hole mass and time lag relationship is not evident for 1H 0707−495 (p > 0.05), it is moderately anticorrelated in IRAS 13224−3809. Since the black hole mass for each source should be constant during these observations, the variable mass found in IRAS 13224−3809 could be induced by other geometric effects in an extended corona environment that are not related to the central mass. It could also be due to the fact that the height and the mass are degenerate (Caballero-García et al. 2018, 2020), since the mass and the source height affect the lags in a similar way. In any case, our results suggest that it might be better to fix the black hole mass when performing reverberation analysis.

Figure 5.

The variable mass from the Table model fits as a function of the reverberation lag for 1H 0707−495 (top) and IRAS 13224−3809 (bottom). Other physical processes related to the coronal properties, rather than geometry, may involve in manifesting the lag-mass scaling relation under the extended corona environment especially in IRAS 13224−3809 where the anticorrelation between the lags and the mass is seen. Of course, there is only one true value of the mass for each AGN.

Open in new tabDownload slide

To obtain the best mass value for all data, we use similar procedures to Alston et al. (2020) by simultaneously fitting the data. We allowed the mass to vary between values of log (M/M_⊙) = 5.5–8.0 and loaded the first data set and ran the table model fit, then added the second data set where the identical parameters were tied for each set. Essentially we leave the first mass as a free parameter and subsequently tie the remaining data masses to the first data set. This method of building one large parameter file was computationally intensive. The results confirmed our expectations that not all data sets would obtain individually excellent fits and, for 1H 0707−495, the χ² statistics ranged from 0.19 to 2.98. This method, however, helped to reduced the large errors for t_max and t_shift for the majority of the data. The single mass value obtained for 1H 0707−495 was log (M/M_⊙) = 6.04 ± 0.01. The 1H 0707−495 mass constrained here is slightly smaller than the mass of log (M/M_⊙) = 6.37 as reported in Zhou & Wang (2005). In fact, Caballero-García et al. (2018) also adopted the mass from Zhou & Wang (2005) and showed that if the mass is larger than this, the lamp-post model will not be able to fit the reverberation time-lag data. Our obtained mass for 1H 0707−495 is then still in the acceptable regime found by Caballero-García et al. (2018). We repeated these methods to obtain the single mass value for IRAS 13224−3809. All 17 data sets were reasonably well constrained statistically where χ² ranged 0.17–2.75, achieving a single mass value of log(M/M_⊙) = 6.30 ± 0.01, which is closely consistent with what obtained by Alston et al. (2020). These results are presented in Table A2. Fig. 6 shows an overview of the behaviour of the time lags as a function of frequency for different black hole mass values and upper X-ray source locations. While the positive lags in previous literature were usually modelled independently using a power-law function, the ECM model can simultaneously produce both negative and positive lags.

Figure 6.

Top panel: The variable mass behaviour for IRAS 13224−3809 combined table model where h₂ = 13.0 r_g and logM  = 6.0 (black), 6.3 (red) and 6.7 (blue). Specific model parameters of the reverberation signature are also provided. Bottom panel: The variable height behaviour is shown for the same signature parameters and mass values as seen in the top panel, where h₂ = 5 r_g (black), 10 r_g (red), and 20 r_g (blue). The other model parameters were Γ₁ = 2.6, Γ₂ = 1.7, log ξ = 3.0, b = 2.0, t_max = 300, t_shift = 20 (for both panels). Note that in the ECM model, the mass scales not only the negative reverberation lags but also the positive lags.

Open in new tabDownload slide

4.4 ECM correlations

We explore the model for correlations using the Spearman’s rank method. For this we use all individual observation table model fits returned from the single mass value phase for each AGN to maximize the available data. We found moderate correlations between Γ₁ and t_max and between Γ₂ and t_max and t_shift. On the other hand, large errors are still a feature of t_max therefore we acknowledge that these correlations are not well constrained. Further investigation reveals these correlations are often much stronger in IRAS 13224−3809 than they are in 1H0707−495. The brightness parameters b was strongly anticorrelated with t_max in IRAS 132224−3809, where no correlation was evident in 1H 0707−495. Note that Caballero-García et al. (2018) also suggested that the 1H 0707−495 mass should be fixed at values below ∼(2–3) × 10⁶ M_⊙, otherwise the observed reverberation lags cannot be explained. The time lag amplitude was moderately correlated with the upper photon index Γ₂ in all cases. The results are presented in Table 3, which shows the parameters of interest and the overview of the Spearman’s rank correlations r_s and their p values for each source in columns 3 and 4 and for any global correlations in the final column. Note that the covering fraction is also tested from the results reported in HYC22. We found no correlation between the upper source height and the luminosity, however, the time lags do correlate moderately with the upper X-ray source height in 1H 0707−495, and a stronger relationship is evident in IRAS 13224−3809. Further consideration of the lags seen ≲250 s in the latter source provides a much stronger correlation where r_s = 0.68 with p-value = 0.005. The six strongest correlations found for these sources are shown in Figs 7 and 8, respectively. The panels within these figures also show where the limits of each parameter were reached as denoted by the blue arrows.

Figure 7.

The 1H 0707−495 ECM results showing moderate to strong correlations where the p-value <0.05. For clarity of the h₂–ξ_ECM correlations (upper middle panel) have been heavily averaged into four groups and shown by the yellow points for data points located in quadrants identified by the grey dashed lines.

Open in new tabDownload slide

Figure 8.

The IRAS 13224−3809 ECM results showing moderate to strong correlations where the p-value <0.05. Each correlation has been divided into bin sizes as shown by the vertical grey dashed lines and the resultant mean plot (and associated errors) has been shown by the large yellow data points for clarity of each correlation.

Open in new tabDownload slide

Both AGNs have a moderate relationship for Γ₂ and t_shift which can be seen as the strongest global correlation where r_s = −0.594 and p-value = 0.00034. The lower source photon index Γ₁ also correlates strongly with t_max in IRAS 13224−3809 where r_s  = 0.716 and p-value = 0.001 with no counterpart correlation seen in 1H 0707−495. For the upper source Γ₂ a moderate correlation is seen with the brightness parameter b in 1H 0707−495 and another strong inverse correlation emerged with t_max in IRAS 13224−3809 where r_s = −0.725 and p-value = 0.001. Finally, a strong correlation found in IRAS 13224−3809 and not in 1H 0707−495 was ΔΓ, (Γ₁ − Γ₂), with t_max where r_s = 0.82 and p-value = 5.86 × 10⁻⁵. It is interesting to note that unlike 1H 0707−495, IRAS 13224−3809 simultaneous fitting does not always lead to Γ₁ > Γ₂. These results suggest that the spectral properties of the corona may be specific to the source. The Γ–L/L_Edd and Γ–h₂ relationships for these AGNs are also presented in Fig. 9. It can be seen that lower and higher limits of Γ and h₂ often reached their limits and these parameters may be wider as indicated by the blue arrows in each panel.

Figure 9.

The evolution of Γ as a function of the Eddington fraction and the source height for 1H 0707−495 (top panels) and also for IRAS 13224−3809 (bottom panels). The Γ₁ and Γ₂ values of IRAS 13224−3809 are more separated and occupied at more extreme values than those of 1H 0707−495, especially when 5 r_g ≲ h₂ ≲ 10 r_g.

Open in new tabDownload slide

In fact, the photon index Γ₁ tends to increase sharply with increasing luminosity above ∼2.5 × 10⁴² erg cm s⁻¹ for both sources. A similar trend in AGN data against the luminosity ratio L_x/L_Edd has been reported by Yang et al. (2015) and interpreted as a two-phase accretion flow model, although more data would be required to investigate this scenario. For 1H 0707−495, the correlation is flat at lower luminosity then a mild positive correlation kicks in when the luminosity is 10^−1.7 ≲ L/L_Edd ≲ 10^−1.5. Contrarily, for IRAS 13224−3809, we see the progressively flat profile of Γ instead, and their Γ₁ and Γ₂ values are extreme and more separated. The behaviour of Γ may suggest that contrasting mechanisms are driving the variability in these AGNs (within the model constraints).

5 DISCUSSION

The initial fitting of the ECM to the various flux groups achieved reasonable fits and model descriptions, however, these were not always well constrained and error ranges often tended to the model maximum and minimum allowed values. The model can provide a good statistical fit to the observed data when the mass is employed as a free parameter and is capable of finding the upper X-ray source heights between 3 and 20 r_g. Fixing the mass at the single best value obtained from the simultaneous fits also achieved similar results and whilst the black hole mass of 1H 0707−495 dropped from 6.22 (obtained by the independent emcee fit) to 6.04 ± 0.01log (M/M_⊙), both model fits remain statistically reasonable. The black hole mass for IRAS 13224−3809 was estimated at 6.30 ± 0.01log (M/M_⊙) and remained reasonably consistent with the independent measurement. These obtained mass values contained very small errors and they fall within the limits of the variable mass values from the model fits presented in Table A1.

Caballero-García et al. (2020) fix the mass of IRAS 13224−3809 at log (M/M_⊙) = 6.30 and their model looks at the reverberation from a single point source and predicts only the negative reverberation lags, whilst the positive lags at low frequency are modelled separately using the KYNXILREV phenomenological power law. Our model is mass dependent and considers both soft and hard lag mechanisms simultaneously. Nevertheless, we note that a limiting factor may be that we have fixed the lower source height at 2 r_g; an assumption required to ease the heavy computations required to create the predicted time lags and hence we are estimating the extent of the corona rather than its true size that may be obtained by employing h₁ as a free parameter. This would be extremely intensive with the current model and further development should maintain its complexity whilst enhancing model computational performance.

Despite this, some meaningful relations among the ECM parameters can still be inferred. For 1H 0707−495, we find significant correlation (p < 0.05) between h₂ and Γ₁, and significant anticorrelation between h₂ and ξ. This means that when the 1H 0707−495 corona extends upwards (larger h₂), the coronal emission seems to be softer (larger Γ) while producing less overall ionization on the disc (smaller ξ). This is expected since the vertically extended corona should produce less intense illumination pattern on the inner disc (e.g. Wilkins & Fabian 2012). The tendency of softer corona with increasing its vertical extent is also found in IRAS 13224−3809. This is in line with Alston et al. (2020) and Chainakun et al. (2022b) where the lamp-post geometry was used and the photon index of the continuum was found to increase with the source height. In fact, Chainakun et al. (2022b) also found that during the source height increases, the disc itself generates more high-frequency variability, suggesting that the inner disc becomes more active. Here, we find the significant anticorrelation between Γ₂ and both t_max and t_shift (r_s = −0.724 and −0.648, respectively). Therefore, when the extending corona becomes softer, the source response as well as the signal propagation occurs faster. This may agree with Chainakun et al. (2022b) that as the source height increases the innermost region is more active so that the corona requires shorter time in response.

Nevertheless, in the ECM environment, the time lags in IRAS 13224−3809 correlate stronger with h₂ than those found in 1H 0707−495. Consideration of the lags seen under 250 s provides a much stronger correlation where r_s = 0.68 with p-value  = 0.005, hence there are hints that much longer time lags will still follow this correlation albeit with a shallower slope. Furthermore, Chainakun et al. (2019) pointed out that in order reproduce the time lags, their spherical corona model required higher coronal temperatures for a lower optical depth, τ, supporting the τ and coronal temperature anticorrelation argument discussed by Tortosa et al. (2018). Recently, Chainakun et al. (2022a) presented a new method of predicting the black hole mass of X-ray reverberating AGNs using artificial neural networks and pointed out that the inconsistency of the lag–mass relationship may be due to the lag amplitudes being more strongly affected by other geometric effects that are not related to the mass of the black hole. Also, there is a significantly larger number of parameters in the ECM than in the lamp-post model. Perhaps, this might explain the observed anticorrelation between the lags and the mass seen in IRAS 13224−3809 under the extended corona environment when using individual observations (Fig. 5), that the lag–mass relation is also modulated by other coronal parameters in the model.

In the lamp-post geometry, the X-ray source heights in IRAS 13224−3809 have been reported to correlate positively with the luminosity in the 2–10 keV energy range (Alston et al. 2020; Caballero-García et al. 2020; Chainakun et al. 2022b). However, there was no evidence of this relationship in this study. Whilst this is consistent with the findings of Szanecki et al. (2020), it is in contrast to their findings that the source is compact within a few gravitational radii (when modelled using the relativistically smeared reflection from the accretion disc). Variations of the inner regions and accretion flow have been explained via transitions between the jet emitting disc (JED) and standard accretion disc (SAD) framework (Marcel et al. 2022) and the framework of coupled hot accretion and jet model for GBHBs and AGNs (Yang et al. 2015), suggesting that a two-phase accretion flow could be driving the observed variability.

For 1H 0707−495, the profile of Γ is flat at lower luminosity before showing a mild positive correlation at 10^−1.7 ≲ L/L_Edd ≲ 10^−1.5 (Fig. 9). Although this is a small fraction of the Bolometric luminosity, similar findings have been reported by the two-phase accretion Type II luminous hot accretion flow (LHAF) and disc–corona regime, after which the evolutionary turnover point is reached when the accretion rate changes. Furthermore, the evolutionary behaviour of Γ₁ and Γ₂ is very similar when tracked against the Eddington fraction. Although there is not a multitude of individual AGNs categorized in the two-phase accretion scenario due to the requirement for multiple long observations, a theoretical explanation of this behaviour was found in AGNs and reported by Zdziarski et al. (2003) who suggested that the emission of a cold medium irradiating hot plasma could be due to the cooling via reprocessed emission of the hot plasma causing the observed X-ray spectrum to soften, leading to a higher value of Γ.

Contrarily, the Γ profile for IRAS 13224−3809 is relatively flat, which is comparable to the Γ − L_Bol/L_Edd behaviour of the jet-phase regime of the JED–SAD framework where the Γ is determined by the energy distribution of power law electrons in the jet and roughly constant for sources with different luminosity’s, neither source follows the expected Γ ≈ 2.1. Furthermore, for IRAS 13224−3809, the current model values of Γ₁ are extreme although not uncommon (see e.g. Wilkins et al. 2014) and Γ₂ may drop well below 1.4 approaching lower extreme values. There appears to be no relationship with either source Γ and the accretion rate in IRAS 13224−3809. Similar Γ behaviour was reported from RXTE observations of MCG-6-30-15 where the spectral index increased with luminosity, reaching a final upper value after which it remained roughly static as the luminosity continued to increase (McHardy, Papadakis & Uttley 1999), possibly due to a harder spectral component that does not change within the observation duration as the continuum component steepens the spectra.

A steeper spectrum is expected for increasing luminosity and this is evident in the spectral model for IRAS 13224−3809 as reported in HYC22. The ECM scenario, however, shows very limited spectral index variations for its flux variations. This has been previously reported for NGC 5548 (Sobolewska & Papadakis 2009) and PG 0804+761 (Papadakis, Reig & Nandra 2003). Furthermore, Caballero-Garcia et al. (2012) found no significant spectral variability in NGC 4151 and NGC 2110, but reported the significant flux variations that indicated the intrinsic variability of the central source in NGC 4388, NGC 4945, and IC 4329. This raises the question of different spectral states in AGNs. However, 1H 0707−495 and IRAS 13224−3809 are both classified narrow-line Seyfert 1 galaxies where the 2–10 keV flux variations are almost always associated with spectral variations (Papadakis et al. 2002). The modelling applied thus far (spectral and ECM models) suggests that different mechanisms may be contributing to the variability that may be explained by the differences in the geometry of each source.

Note that variations in the optical depth τ of the corona cause spectral changes in the Comptonized emission leading to a steeper spectrum and the temperature decreases with increasing τ (e.g. Haardt, Maraschi & Ghisellini 1997), therefore the geometry of 1H 0707−495 should consist of a smooth corona where the intrinsic flux varies in unison with Γ₁ and Γ₂. IRAS 13224−3809 however, may contain a slightly different geometry since the Γ₁ and Γ₂ values are more clearly separated, generally remaining relatively constant and occupied at the extreme values of the model especially when 5 r_g ≲ h₂ ≲ 10 r_g (Fig. 9). This suggests that IRAS 13224−3809 corona is more patchy (i.e. having clear, distinct spectral zones of spectral Γ₁ and Γ₂ that may be related to two distinct zones of corresponding τ and coronal temperature), especially at locations where the corona ≲ 10 r_g.

Global coronal correlations across AGN samples have remained elusive and/or confusing (see e.g. Sarma et al. 2015; Hinkle & Mushotzky 2021; Kamraj et al. 2022) the latter authors cautioning the use of coronal parameters and the Eddington fraction relations to infer properties of black hole systems. In contrast, this work finds that strong relations are evident in the grouped data and when drilling down to individual observations. These differences are possibly due to the choice of spectral models and parameters explored.

Many of the χ² values fell below unity, suggesting that the model may be overfitted due to too many free parameters. This suggests that deeper model assumptions leading to the freezing of at least one more parameter may also be appropriate for future modelling. Some parameter ranges may need to be slightly wider to accommodate better statistics (whilst remaining feasible), especially for h₂, Γ₁ and Γ₂. On the other hand, the large errors seen in t_max and t_shift will be difficult to constrain with current observations and the large errors associated with the low frequency fluctuation lags. More AGNs need to be explored by this model to test its capabilities against a wider sample with higher quality time lags versus frequency allowing a better statistical comparison of the data and models.

6 CONCLUSIONS

This work has explored the extended corona scenario using two X-ray point sources. The lower source may represent the base of a jet-like structure or the lower regions of a compact corona with the upper source representing the extended region where the flares are detected due to the periodic vertical collimation of a jet structure or outer region of the corona. The model is capable of fitting the lag-frequency of 1H 0707−495 and IRAS 13224−3809, where their black hole mass can be constrained to log (M/M_⊙) = 6.04 ± 0.01 and 6.30 ± 0.01, respectively.

This work also supports the advantages of exploring corona correlations using individual observations and also suggests that the physics and geometry of the corona may diverge between different sources. For the first time, we have gained an insight to the behaviour of Γ for each of the X-ray sources located above the black hole axis through the use of reverberation mapping. We find the tendency of softer corona with increasing its vertical extent in both AGNs. This produces a much less-ionized disc which is evident in 1H 0707−495. On the other hand, for IRAS 13224−3809 we also find the hint of increasing Γ₂ with decreasing t_max and t_shift, suggesting that shorter propagating fluctuations and faster coronal response may be evident when the corona extends vertically upwards. The intrinsic flux of IRAS 13224−3809 show less variability with Γ and its corona may be more patchy in a sense that contains clearer separate spectral zones of distinct Γ₁ and Γ₂ that may also link to two clearer distinct zones of coronal temperature than that of 1H 0707−495.

ACKNOWLEDGEMENTS

The calculations in this work were carried out using the high performance computer BlueCrystal of the Advanced Computing Research Centre, University of Bristol, UK. SH thanks the Science and Technology Facilities Council (STFC) for funding and the Bristol, Cardiff & Swansea Centre for Doctoral Training (CDT) Team for support. PC thanks funding support from (i) Suranaree University of Technology (SUT), (ii) Thailand Science Research and Innovation (TSRI), and (iii) National Science Research and Innovation Fund (NSRF), project no. 160355. This research has made use of isis functions (ISISscripts) provided by the Erlangen Center for Astroparticle Pysics (ECAP)/Remeis observatory and the Massachusetts Institute of Technology (MIT) (http://www.sternwarte.uni-erlangen.de/isis/). We wish to thank the anonymous reviewer for comments and suggestions which has improved the quality of this manuscript.

DATA AVAILABILITY

The original data underlying this article were obtained from the XMM–Newton Observatory (http://nxsa.esac.esa.int). The derived data generated in this research can be downloaded via http://www.star.bris.ac.uk/steff/hancock.html.

Footnotes

REFERENCES

APPENDIX A: TABLE MODEL ECM RESULTS