key: cord-0311396-5qk5zm0n
authors: Jiang, Jiachen; Dauser, Thomas; Fabian, Andrew C.; Alston, William N.; Gallo, Luigi C.; Parker, Michael L.; Reynolds, Christopher S.
title: XMM-Newton Observations of the Narrow-Line Seyfert 1 Galaxy IRAS 13224$-$3809 : X-ray Spectral Analysis II
date: 2022-04-21
journal: nan
DOI: 10.1093/mnras/stac1144
sha: 86312bd5987a0e5dffeaa72b1e10c303bb485f4b
doc_id: 311396
cord_uid: 5qk5zm0n

Previously, we modelled the X-ray spectra of the narrow-line Seyfert 1 galaxy IRAS 13224$-$3809 using a disc reflection model with a fixed electron density of $10^{15}$ cm$^{-3}$. An additional blackbody component was required to fit the soft X-ray excess below 2 keV. In this work, we analyse simultaneously five flux-resolved XMM-Newton spectra of this source comprising data collected over 2 Ms. A disc reflection model with an electron density of $n_{rm e}approx10^{20}$ cm$^{-3}$ and an iron abundance of $Z_{rm Fe}=3.2pm0.5Z_{odot}$ is used to fit the broad-band spectra of this source. No additional component is required to fit the soft excess. Our best-fit model provides consistent measurements of black hole spin and disc inclination angle as in previous models where a low disc density was assumed. In the end, we calculate the average illumination distance between the corona and the reflection region in the disc of IRAS 13224$-$3809 based on best-fit density and ionisation parameters, which changes from 0.43$sqrt{f_{rm AD}/f_{rm INF}}$ $r_{rm g}$ in the lowest flux state to 1.71$sqrt{f_{rm AD}/f_{rm INF}}$ $r_{rm g}$ in the highest flux state assuming a black hole mass of $2times10^{6}M_{odot}$. $f_{rm AD}/f_{rm INF}$ is the ratio between the flux of the coronal emission that reaches the accretion disc and infinity. This ratio depends on the geometry of the coronal region in IRAS 13224$-$3809. So we only discuss its value based on the simple `lamp-post' model, although detailed modelling of the disc emissivity profile of IRAS 13224$-$3809 is required in future to reveal the exact geometry of the corona.

IRAS 13224−3809 (IRAS 13224 hereafter) is a narrow-line Seyfert 1 galaxy (NLS1, Gallo 2018) hosting a supermassive black hole of BH = BH / ≈ 2 × 10 6 (Zhou & Wang 2005; Alston et al. 2020) . Detailed disc reflection modelling of the broad Fe K and Fe L emission lines in the X-ray band (Ponti et al. 2010; Fabian et al. 2013 ) suggests a maximally spinning black hole in the centre (Fabian et al. 2013; Chiang et al. 2015; Jiang et al. 2018 ). IRAS 13224 is well known for its complex and rapid X-ray variability that is dominated by varying continuum emission (Boller et al. 1997; Gallo et al. 2004; Fabian et al. 2013; Buisson et al. 2018; Jiang et al. 2018; Alston et al. 2019) . Besides, multiple blueshifted absorption features have been found in the X-ray spectra of IRAS 13224, which may correspond to a highly variable outflow with a line-of-sight velocity of ≈ 0.24 (Parker et al. 2017b; Parker et al. 2017a; Jiang et al. 2018 ; Pinto ★ E-mail: jcjiang@mail.tsinghua.edu.cn et al. 2018) or a thin layer of high-ionisation matter that is co-rotating with the innermost region of the accretion disc (Fabian et al. 2020 ).

In the first paper of this series (Jiang et al. 2018) , we analysed the spectra of IRAS 13224 that were extracted from the XMM-Newton observing campaign of this source in 2016. We followed the same approach as in Chiang et al. (2015) : a disc reflection model with a fixed disc density of e = 10 15 cm −3 was used on an earlier and smaller dataset. A super-solar iron abundance of Fe ≈ 24 was inferred by this model. An additional blackbody component was required to model the excess emission in the soft X-ray band. This blackbody component followed a flux-temperature relation of ∝ 4 , suggesting an emission region with a constant area, e.g. the inner accretion disc (Chiang et al. 2015) . Reverberation studies also support the disc origin of the soft excess emission in IRAS 13224 (Ponti et al. 2010) . Detailed reverberation lag analyses based on a simple lamppost model suggest an extended coronal region of ℎ ≈ 10 − 20 g in the high flux states and a more compact coronal region of ℎ ≤ 10 g in the lower flux states of IRAS 13224 (Alston et al. 2020), similar to measurements of micro-lensed quasars (Chartas et al. 2017 ) and X-ray occultation events (e.g. Risaliti et al. 1999; Gallo et al. 2021) .

In this work, we present a reflection-based model for the XMM-Newton spectra of IRAS 13224. In particular, a disc model with a variable disc density parameter ( e ) is considered (Ross & Fabian 2007; García et al. 2016) . The previous assumption in the disc reflection model of e = 10 15 cm −3 was only appropriate for the disc around a massive super-massive black hole both theoretically and observationally (e.g. BH ≥ 10 8 , Shakura & Sunyaev 1973; Jiang et al. 2019a,c) . At a higher density, free-free absorption become important up to X-ray energies, and the reflection spectrum of the disc produces a blackbody-like feature in the soft X-ray band (García et al. 2016 , or see Fig.1 ). Such a feature can explain the soft excess emission commonly seen in the X-ray spectra of unobscured Seyfert 1 galaxies (Mallick et al. 2018; Garcia et al. 2018; Jiang et al. 2019c Jiang et al. , 2020b . Moreover, the disc density parameter inferred by reflection spectroscopy also shows a correlation with the mass accretion rate of the accretion disc according to the studies of stellar-mass black holes (Jiang et al. 2019b (Jiang et al. , 2020a . Last but not least, a high density disc reflection model may be able to solve the problem of supersolar iron abundances obtained by previous models (Tomsick et al. 2018; Jiang et al. 2019c) . Lower inferred iron abundances are due to the model requiring a higher contribution of disc reflection to the total X-ray spectrum (Jiang et al. 2019c) .

In Section 2, we introduce our high density disc model for the XMM-Newton spectra of IRAS 13224. In Section 3, we estimate the size of the coronal region for different flux states based on the best-fit parameters given by our model. In Section 4, we discuss and conclude our results.

In this section, we introduce a high density disc-based model to explain the spectral variability of IRAS 13224. We first present an overview of the broad band spectra of IRAS 13224 in five different flux states. Then we focus on the spectral modelling using high density disc reflection model. 

We consider all the archival XMM-Newton observations of IRAS 13224 in this work, including the 1.5 Ms observing campaign in 2016 (PI: Fabian, A. C.) and previous observations with a length of around 600 ks. A full list of these observations can be found in Table 1 . We focus on the broad band spectral variability shown by the EPIC-pn (pn hereafter) observations. All the pn observations were operated in the Large Window mode. SAS v.1.3 was used to extract spectra from filtered clean event lists. We follow the standard data reduction method as already introduced in the first paper of this series (Jiang et al. 2018) . The same source regions as in the first paper of this series (Jiang et al. 2018) were used after using the EPATPLOT tool to check pile-up effects. IRAS 13224 shows rapid variability in the X-ray band (Boller et al. 2003; Fabian et al. 2013; Chiang et al. 2015) . In order to better understand the spectral variability of this source, we divided all the observations into five different flux intervals. Five intervals share a similar number of total photon counts. The TABGTI-GEN tool was used for this purpose. This method was also used in our previous analysis of the data (Parker et al. 2017b; Jiang et al. 2018) . The final spectra respectively have an average observed flux of log( /erg cm −2 s −1 ) = −11.948 ± 0.003 (F1, <1.1 cts s −1 ), −11.614±0.003 (F2, 1.1-1.7 cts s −1 ), −11.454±0.003 (F3, 1.7-2.4 cts s −1 ), −11.313 ± 0.002 (F4, 2.4-3.4 cts s −1 ) and −11.108 ± 0.002 (F5, >3.4 cts s −1 ) in the 0.5-10 keV band.

The spectra are grouped to have a minimum number of 20 counts per bin, and no group is narrower than 1/3 of the full width half maximum resolution of the pn instrument. XSPEC v.12.11.01 (Arnaud et al. 1985 ) is used to analyse the spectra. Fig. 2 shows these five flux-resolved spectra of IRAS 13224. They are unfolded using a constant model to remove the effect of instrumental response only for demonstration purposes. As shown in Fig. 2 , different energy bands show flux changes by a different factor between two states. For instance, F5 has a flux approximately 9 times higher than F1 in the 2-4 keV band. A lower flux increase factor is found in the other bands, e.g. a factor of 4 between 4-10 keV band a factor of 7 between 2-4 keV. Similar conclusions are found in Alston et al. (2020) where RMS spectra are analysed. Previous reflection modelling suggests that the 2-4 keV band of IRAS 13224 is dominated by highly variable coronal emission while the 4-10 keV band is dominated by less variable reflection from the innermost region of the disc (Fabian et al. 2009; Chiang et al. 2015; Jiang et al. 2018; Alston et al. 2020) . A higher fraction of disc reflection is found in low flux states than in high flux states (Jiang et al. 2018 ).

In this section, we introduce a high-density disc model for the reflection spectra in the five spectra of IRAS 13224 in different flux states 1 . The coronal emission of IRAS 13224 is modelled using the thermal Comptonisation model nthcomp (Życki et al. 1999) . We consider a disc-blackbody spectrum for the seed photons, which is characterised by the inner disc temperature db . This parameter determines the low energy rollover of a thermal Comptonisation spectrum and cannot be constrained by X-ray spectra alone due to Galactic absorption. We therefore fix this parameter at 10 eV, the typical value for an accretion disc around a supermassive black hole. The electron temperature of the corona is fixed at 100 keV 2 .

The rest-frame disc reflection spectrum is calculated using the reflionx code 3 (Ross & Fabian 2007 ). The nthcomp model is used as the illuminating spectrum (Jiang et al. 2020b ). db , e and photon index Γ of reflionx are all linked to the corresponding parameters in nthcomp. The other free parameters are the ionisation 1 We show a low-density disc model for the reflection spectrum of IRAS 13224 in Appendix A, where an electron density of 10 15 cm −3 is assumed for the disc. 2 In Appendix B, we show the model for the soft X-ray data is consistent when a low e = 10 keV is used, except for a slightly harder photon index. 3 Another popular model xillverd can be used to calculate rest-frame disc reflection spectra too (García & Kallman 2010) . In Appendix C, we present a fit for the reflection spectrum of IRAS 13224 using an xillverd-based model (relxilld). We find that the fit is unsatisfactory due to its limited range of e . Therefore, we choose reflionx over the public version of xillverd for a wider range of the disc density parameter in reflionx. state ( ), the iron abundance ( Fe ) and the electron density within the Thomson depth of the disc e . e is calculated using the solar abundances calculated in Morrison & McCammon (1983) .

The convolution model relconv is applied to the rest frame reflection model reflionx to account for relativistic effects (Dauser et al. 2010) . A broken power-law emissivity profile is considered. The spin of the central black hole, the iron abundance and the inclination angle of the disc are not expected to change on observable timescales. We therefore link these parameters for all five spectra.

Previous studies by Parker et al. (2017b) ; Jiang et al. (2018) ; Pinto et al. (2018) found a highly variable blueshifted ionised absorption line features, suggesting evidence for an ultra-fast outflow with a line-of-sight velocity of ≈ 0.24 . We model the absorption features using the photon ionisation plasma model xstar Kallman & Bautista (2001) . The xstar grid in this work is calculated using a power law-shaped illuminating spectrum with Γ = 2. Free parameters include the column density ( H ), the ionisation 4 ( ), the iron abundance ( Fe ) and the redshift parameter ( ). Note that solar abundance in the xstar model is calculated by Grevesse et al. (1996) , which is different from the values used by reflionx. Therefore, the iron abundance parameters of these components are not linked in our model.

In addition, the tbnew model is used to account for Galactic absorption (Wilms et al. 2000) . The cflux model is used to calculate the unabsorbed flux of each component in the full 0.5-10 keV band. The total model is tbnew*xstar* (cflux*relconv*reflionx + cflux*nthcomp) in the XSPEC format.

The model introduced above is able to describe the broad band spectra of IRAS 13224 very well with a total 2 / = 864.32/798. The best-fit models and corresponding data/model ratio plots are shown in Fig. 3 . The best-fit model parameters can be found in Table  2 . Note that the high density models do not require any additional component for the soft excess flux seen below 2 keV.

Our best-fit absorption model suggests an outflow with a lineof-sight velocity of 0.21 − 0.24 . The ionisation state of the outflow during F1 and F2 are significantly lower than that during F3 and F4. The highest flux state F5 shows no significant evidence for blueshifted absorption features. We therefore fix the ionisation parameter and the redshift parameter of xstar for F5 at the best-fit values for F4. Only an upper limit of the column density of the outflow is found (a 90% confidence range of H < 2 × 10 22 cm −2 ). Similar conclusions were found in previous analyses (Parker et al. 2017b; Jiang et al. 2018; Pinto et al. 2018) . Interested readers may refer to Parker et al. (2017b) for more discussion about the rapid variability of the absorption features.

The best-fit relconv model indicates a highly spinning black hole in IRAS 13224 (>0.98). The inclination of the disc is estimated to be within a 90% confidence range of 64-69 • . They are consistent with previous measurements when a fixed disc density of e = 10 15 cm −3 is assumed (Gallo et al. 2004; Fabian et al. 2013; Chiang et al. 2015; Jiang et al. 2018) .

We also find tentative evidence of a variable disc density parameter: during F1 and F3, the best-fit density parameter is lower than those for other flux states (F2, F4 and F5), although the density parameters for F2-F5 are consistent within their 90% confidence ranges. The uncertainty ranges of their density parameters and the iron abundance parameter are shown in Fig. 4 . When applying a high density disc model to the data of IRAS 13224, one can find significantly lower disc iron abundances of ≈ 3 than previous results (e.g. Chiang et al. 2015; Jiang et al. 2018 ).

Previously, we found the density of the inner disc in IRAS 13224 has to be higher than 10 19 cm −3 (Jiang et al. 2018) . In this paper, we are able to confirm the results and obtain a high disc density of around 10 20 cm −3 for IRAS 13224. Such a disc density is expected for a low-BH NLS1 like IRAS 13224 according to the thin disc-corona model (Svensson & Zdziarski 1994) . Using Eq. 8 in Svensson & Zdziarski (1994) , one can find a disc density of around 10 20 cm −3 at = 2 − 10 g for a 10 6 BH with / Edd = 1 and corona = 0.9. Similar comparison between density measurements for Sy1s and the prediction of the thin disc model can be found in Jiang et al. (2019c) . The parameter corona is the fraction of the disc energy that is radiated away in the form of non-thermal emission, including disc reflection and coronal emission. The value of corona is estimated to be close to 1 to explain the typical spectral energy distribution of Sy1s (e.g. Haardt & Maraschi 1991; Jiang et al. 2020b ). According to our best-fit models, the flux of the disc reflection component ( refl ) and the coronal emission ( pl ) increases with the X-ray luminosity as shown in the first panel of Fig. 5 . refl increases from 10 −12 erg cm −2 s −1 for F1 to 5 × 10 −12 erg cm −2 s −1 for F5 while pl increases from 3 × 10 −13 erg cm −2 s −1 for F1 to 4 × 10 −12 erg cm −2 s −1 for F5.

The second panel of Fig. 5 shows how the strength of the reflection component relative to the power-law continuum (reflection fraction) changes with X-ray flux. The reflection fraction in the figure is defined as the flux ratio between two components in the 0.5-10 keV band. A significantly higher fraction of reflection is shown in the lower flux states. A decreasing reflection fraction along with an increasing X-ray luminosity was also found in previous analyses of IRAS 13224 when a low density disc model was used and an additional blackbody model was required to model the soft excess (Jiang et al. 2018) .

In order to better demonstrate the spectral difference between the high and low flux states, we show a F1 data / F5 best-fit model ratio plot in Fig. 6 . The grey dashed horizontal line shows the 0.5-10 keV band flux ratio of these two states. The negative relation between the direct X-ray coronal emission and the reflection component can be explained by the strong light bending of the coronal emission when the corona is closer to the black hole: the flux of the continuum emission decreases as more coronal photons are lost to the event horizon and due to the large redshift closer to the black hole. The reflection spectrum of the disc is less affected by the light-bending effects (e.g. Miniutti et al. 2003; Dauser et al. 2013; Wilkins et al. 2021) . Reverberation studies of IRAS 13224 also support this conclusion (Alston et al. 2020 ). Other explanations may involve different geometries for the corona (e.g. extended sizes or outflowing geometries, Dauser et al. 2013; Steiner et al. 2017; Szanecki et al. 2020) . Returning radiation of the inner disc (Riaz et al. 2021) and small changes in the inner disc radius may introduce changes in the observed reflection fraction too. In the next section, we present calculations of the distance between the corona and the inner disc suggested by our best-fit spectral models. In particular, we consider the best-understood lamppost model because of the difficulty of calculating other complex geometries. Future efforts based on other theories are needed for the analysis of the X-ray data of IRAS 13224.

In this section, we calculate the average illumination distance ℓ between the corona and the illuminated disc by using the bestfit parameters obtained in Section 2. The estimation of is also determined by the geometry of the coronal region. In particular, we consider a lamppost geometry for simplicity, where a point-like isotropic source above the accretion disc is assumed (Matt et al. 1991) , to provide a consistent interpretation of our energy spectra as our reverberation analyses for the same set of data in Alston et al. (2020) .

We first calculate the average illumination distance between the corona and the disc by using the best-fit ionisation ( ) and density ( e ) parameters.

In the reflionx model, is defined as following by given the flux of the illuminating spectrum ( AD ) and the number density of the electrons ( e ),

where AD calculated between 0.01-100 keV (Ross et al. 1999 ).

AD can be calculated using

where AD is the fraction of the coronal emission that reaches the accretion disc divided by solid angle, pl is the luminosity of the corona, is the average illumination distance between the disc and the corona. Similarly, the observed coronal flux is

where INF is the fraction of the coronal emission that reaches infinity divided by solid angle, and is the distance of the source. By combining the equations above, we obtain the average illumination distance between the corona and the disc is

where = √︃ AD INF . In our work, the absorption-corrected coronal flux corona is calculated in the same energy band (0.01-100 keV) as the ionisation parameter.

The values of / are calculated by using Eq. 4 for all the flux states using their best-fit parameters, and the results are in Table 2 and the last panel of Fig. 5 . A luminosity distance of 288 Mpc reported on the NED website is used 5 (Jones et al. 2009 ).

The lowest value of ≈ 0.43 × 2×10 6 BH g = 8.6 × 10 5 BH g is found for the lowest X-ray flux state (F1). The highest flux state (F5) has the highest value of ≈ 1.71 × 2×10 6 BH g = 3.4 × 10 6 BH g .

So far, we have calculated the average illumination distance between the corona and the disc using the ratio = √︁ AD / INF . The exact values of this ratio depends on the geometry of the coronal region and the spacetime around the BH. In this work, we consider an isotropic point-like corona, namely the lamppost model (Matt et al. 1991) , for simplicity and consistency with our work in Alston et al. (2020) .

In order to calculate AD / INF , the flux incident on the accretion disc has to be calculated. Assuming an isotropic emission of the √︁ INF / AD in the fourth panel is estimated using the best-fit e , and the 0.01-100 keV band flux of the coronal emission. A black hole mass of 2 × 10 6 is assumed. corona, this ratio would be unity in flat space time. Due to the strong relativistic effects close to the black hole, however, several additional factors have to be taken into account.

First of all, light bending leads to an increasing fraction of photons bent towards the accretion disc, which is parametrized by the reflection fraction refl (see Dauser et al. 2014 ) and can easily reach values of 10 and larger for a compact corona 6 . Additionally, the energy shifts from the corona to the accretion disc AD and from 6 The reflection fraction here is different from obs refl shown in Table 2 Figure 7a shows the flux ratio AD / inf for the range of parameters allowed within the rage of our best fits. For very low heights of the corona, the flux incident on the accretion disc can be more than 100 times larger than the observed flux from the corona. In order to be able to compare our best fit parameters with these calculations, Fig. 7b plots lum (ℎ) = ℎ/ √︁ AD / INF .

Secondly, we estimate the size of the reflection region in the disc, based on which we argue that the average distance between the corona and the reflection region ( ) shares a similar value with the height of the corona (ℎ) when ℎ is sufficiently small. In Fig. 8 , we show the outer radius defining the area which sees most of the incident flux of the corona using the same codes as in Dauser et al. (2013) . For instance, using 50% of the enclosed flux as the averaged value of incident radiation, the corona can be seen that up to ℎ =5 g most of the radiation hits the disc between the ISCO and 2 g (see the dashed lines). The average distance between the corona and the reflection region is approximately = √ 5 2 + 2 2 = 5.38 g =1.08ℎ ≈ ℎ. In comparison, the difference between and ℎ becomes important when the corona is very far from the disc, e.g. = √ 2 + ℎ 2 > √ 2ℎ when ℎ > 20 g . Our best-fit spectral model suggest a compact region in the inner disc where reflection happens. Therefore, ℎ √︁ INF / AD can be approximated as √︁ INF / AD . In Fig. 7 , we show the values of √︁ INF / AD for the lowest and the highest flux states of IRAS 13224 that are calculated in Section 3.1. They correspond to a height of the corona that is varying between 3-6 g . Similar calculations have also been used to study the geometry of the innermost accretion region in Cyg X-1 (Zdziarski & De Marco 2020) . Further discussion of the lamppost model can be found in Niedźwiecki et al. (2019) where some minor model differences are discussed.

Previously, the X-ray spectra of IRAS 13224 were fitted by reflection models with super-solar iron abundances of Fe ≈ 20 and a low disc density of 10 15 cm −3 (e.g. Chiang Figure 8 . The outer radius, which encloses 90% (solid) and 50% (dashed) of the total flux incident on the accretion disc. To guide the eye, the solid line is plotted for = ℎ. An additional blackbody model and two relativistic disc reflection components were required for the stacked spectrum of IRAS 13224 (Fabian et al. 2013; Jiang et al. 2018 ). This double-reflection model was required because an averaged spectrum extracted from a wide range of flux states was considered 7 .

In this work, we consider all the archival XMM-Newton observations of IRAS 13224 and divide the data into 5 flux-resolved intervals. We find that only one high density disc reflection model with e ≈ 10 20 cm −3 and Fe ≈ 3 together with a thermal Comptonisation model is needed for each flux state. The decrease of inferred iron abundances in a high density model is due to an increase in the relative strength of the reflection component. Similar conclusions were found in other sources (e.g. Ton S180, Jiang et al. 2018) . A more detailed comparison between low-density and high-density models for the reflection spectrum of IRAS 13224 can be found in Appendix A.

Caballero- García et al. (2020) applied an intermediate density model to both the energy spectra and the lag spectra of IRAS 13224. The lamp-post geometry was used in their work. They concluded a low disc density of 10 18 cm −3 . This was because they did not include data below 3 keV in their spectral analysis. As shown in Ross & Fabian (2007) ; García et al. (2016) ; Jiang et al. (2019c) , the effects of the density parameter on the resulting reflection spectrum are, however, mostly confined to the soft X-ray band, e.g. <2 keV.

Despite a different density parameter, other spectral parameters are consistent with previous measurements, such as BH spin and disc inclination angle. The anti-correlation between the observed X-ray flux and the flux ratio of reflection and Comptonisation components suggests strong light-bending effects that take place in the innermost region of the accretion disc. Besides, we are also able to confirm the existence of variable UFO features that were identified in Parker et al. (2018) .

In this work, we also calculate the average illumination distance between the corona and the disc by using the best-fit density and ionisation parameters. When an isotropic point-like corona model is applied, our measurements suggest that the corona moves from 3× 2×10 6 BH g away from the disc during the lowest flux state (F1) to 6× 2×10 6 BH g during the highest flux state (F5). Alston et al. (2020) calculated the height of the corona by modelling the lag-frequency spectra extracted from the observations of IRAS 13224 in 2016 with a lamppost model: ℎ varies between 6 +4 −2 g and 15 ± 3 g for a similar BH ≈ 2 × 10 6 . These values of ℎ are 2-3 times higher than our estimations in Section 3.2.2. Similar measurements of ℎ were obtained by (Caballero- García et al. 2020) .

Although our spectral analysis assumes no particular coronal geometry by using a broken-power law emissivity profile, the disagreement between the lamppost interpretations for the reverberation lags and the observed e of IRAS 13224 suggests modifications to the simple lamp-post model may be needed (see Appendix D for more discussion concerning the geometry of the corona). For example, a coronal region that extends radially while moving along the rotation axis of the central black hole may exist in the NLS1 7 Appendix E shows an example where a double-reflection model similar to the one in Jiang et al. (2018) is used to fit the stacked spectrum of IRAS 13224.

1H 0707−495, a sister source of IRAS 13224 (Wilkins et al. 2014 ). If such a complex coronal geometry exists in IRAS 13224, the calculations of AD INF will have to change accordingly too. Besides, returning radiation might be another possible explanation for the difference of our measurements: some of the radiation of the disc, including thermal emission and reflection, may return to the inner disc due to strong light-bending effects. Returning radiation is particularly important around a rapidly spinning BH with * > 0.98 (Cunningham 1976 ), e.g. the supermassive BH in IRAS 13224. This process will consequently cause an underestimation of ℎ by increasing the ionisation state of the inner disc. Future analysis involving multiple reflection processes is needed for the observations of IRAS 13224 and other very soft NLS1s (e.g. Ross et al. 2002) . ACKNOWLEDGEMENTS This paper was written during the COVID-19 pandemic in 2020-2021. We acknowledge the hard work of all the health care workers around the world. J.J. acknowledges supports from the Tsinghua Astrophysics Outstanding Fellowship and Tsinghua Shui'Mu Scholarship Program.

The data underlying this article are available in the High Energy Astrophysics Science Archive Research Center (HEASARC), at https://heasarc.gsfc.nasa.gov. The reflionx model used in this work will be available upon reasonable request.

In this section, we compare the high-density disc model in this work with a low-density disc model for the reflection component in IRAS 13224.

As shown in Fig. 5 and Jiang et al. (2018) , the disc reflection component of IRAS 13224 makes more contribution to the total Xray luminosity in a low flux state than in a high flux state. Therefore, we focus on the F1 (lowest flux state) spectrum to demonstrate how the density parameter changes our model.

The model considering a variable disc density parameter as introduced in Section 2.2 is named 'Model 1' hereafter. Model 1 is able to provide a good fit for F1 spectrum with 2 / = 217.75/154, and we find a disc density of 9 × 10 19 cm −3 . The best-fit parameters are shown in the first column of Table A1 . They are consistent with the values obtained when we fit F1-F5 spectra simultaneously (see Table 2 ). F1 spectrum alone is not able to constrain the lineof-sight Galactic column density. We thus fix this parameter at H = 6.78 × 10 20 cm −2 (Willingale et al. 2013 ) in the following tests.

We then consider a low-density disc model assuming e = 10 15 cm −3 . An additional bbody model is required to model the soft excess emission of IRAS 13224 in combination with the relativistic disc reflection component (Model 2). The best-fit parameters are shown in the second column of Table A1 , and the best-fit model is shown in the left panel of Fig. A1 . Model 2 provides a slightly better fit than Model 1 with an insignificant improvement of Δ 2 ≈ 5 and one additional free parameter. We only obtain an upper limit for the Galactic column density at H < 9 × 10 20 cm −2 in Model 2. This parameter is fixed at H = 6.78 × 10 20 cm −2 calculated by Willingale et al. (2013) .

In comparison, Model 2 needs a harder power-law continuum (Γ ≈ 1.97) than Model 1 (Γ ≈ 2.07). Besides, Model 2 requires a significantly higher iron abundance ( Fe ). Because the inferred disc reflection component makes less contribution to the total X-ray luminosity ( obs refl = 1.27 +0.12 −0.11 ) when a low-density model is used. In comparison, obs refl for Model 1 is 3.5 +0.4 −0.3 . Fe > 20 is thus needed in Model 2 to fit the broad Fe K emission line in the data. Similar conclusions were found in other AGNs (e.g. Tons S180, Jiang et al. 2019c) .

Furthermore, Model 2 requires a disc ionisation parameter one order of magnitude higher than the one in Model 1. The resulting e is only 1.6 × 10 −4 of the value given by Model 1. This suggests a corona that is over 90 times larger than the estimated value in Section 3 (see eq. 4), e.g. over a hundred g , which is in tension with the observed reverberation lags of IRAS 13224 (e.g. Alston et al. 2020 ). is required to fit the soft excess emission. Right: a cool corona of e = 10 keV is assumed. The disc density parameter is allowed to vary in this model. The best-fit model for the data below 10 keV is consistent with the one based on the assumption of e = 100 keV. Table A1 . Best-fit parameters for F1 spectrum of IRAS 13224 using Model 1-3. See text for more details. Model 2 considers a grid of density parameters. An additional blackbody model is used to fit the soft excess together with the disc reflection model.

Lastly, similar to the tests in Section 4.3 in Jiang et al. (2018) , we investigate the possibility of fitting the spectrum of IRAS 13224 using a disc reflection model with an intermediate density parameter. Model 2, which includes an additional bbody model, is used for this purpose. We consider a grid of e by fixing this parameter at 10 16 , 10 17 , 10 18 and 10 19 cm −3 .

After adapting an intermediate disc density parameter, we find the fit is significantly worse than Model 1 although they provide acceptable fits justified by 2 / . See Table A1 for parameters and Fig. A2 for corresponding models. Model 2 with e = 10 15 cm −3 and Model 1 have 2 = 212.32 and 217.75. In comparison, 2 peaks at 265.71 when e = 10 17 cm −3 . The improvement of the fits in Model 1 and Model 2 with e = 10 15 cm −3 is in the soft X-ray band, e.g. <1 keV (see Fig. A3 ). At e = 10 15 cm −3 , the soft X-ray spectrum below 0.8 keV is modelled by bbody; at e > 10 19 cm −3 , this part of the spectrum is modelled by the disc reflection component. This is because the spectral shape of the soft excess emission requires the density parameter to be within a certain range of values, e.g. ≈ 9 × 10 19 cm −3 for IRAS 13224, to fit the data (see Fig. 1 ). Similar conclusion was found in our previous work (e.g. see Fig. 11 in Jiang et al. 2018 ). When we consider the best-fit density parameter of Model 1 by fixing e at 9 × 10 19 cm −3 in Model 2, the bbody component is not constrained. We fix of bbody at 46 eV, which is the best-fit value when assuming e = 10 19 cm −3 , we obtain only an upper limit of the normalisation of bbody (< 5 × 10 5 ).

We also note the differences in the parameters of Model 2 when different values of e are used. The most striking difference is the decrease of the inferred iron abundance in the reflection model along with the increase of observed reflection fraction (see Fig. A4 for a comparison of these parameters). A similar result was already achieved in Jiang et al. (2019c) , where a significantly lower iron abundance and a higher reflection fraction were found for the NLS1 Ton S180 when a high density model was used. The application of a high-density model increases the contribution of the reflection component indicated by the values of obs refl in the total emission. A lower iron abundance is thus required to fit the data.

During the 1.5 Ms XMM-Newton observing campaign in 2016, NuS-TAR also observed IRAS 13224 with a total exposure of around 500 ks. By fitting the 3-30 keV FPM spectra alone with a lowdensity model, Jiang et al. (2018) found a lower limit for the high energy cutoff cut at 15 keV 8 (2 ). We do not consider NuSTAR observations in this work, and a coronal temperature of e = 100 keV is used in Model 1. Because the NuSTAR observations only covered a fraction of the XMM-Newton orbits and we focus on the variability of the 0.5-10 keV emission of IRAS 13224.

To investigate how a low e may change our reflection model in the soft X-ray band, we fix e at 10 keV, the lowest value reflionx allows (Model 3). The best-fit parameters are shown in the third column of Table A1 , and the best-fit model in the 0.5-80 keV band is shown in the middle panel of Fig. A1 . When a low e is used, a low exponential cut-off is seen in the model above 20 keV. However, the fit for the soft X-ray emission is consistent except for the photon index of the continuum: a slightly lower value 8 cut is approximately 2-3 times the temperature of the corona ( e ). of Γ (1.982 +0.012 −0.018 ) is found in comparison with a highe model (2.072 +0.008 −0.002 ).

Another commonly used disc reflection model is relxilld (García et al. 2016) , which calculates reflection spectra from a photoionised slab illuminated by a power law-shaped spectrum and includes relativistic effects in the vicinity of a BH. The allowed range of e in relxilld is 10 15 -10 19 cm −3 . In Fig. C1 , we present a fit for the F1 spectrum of IRAS 13224 using relxilld, which has the same parameters as Model 1. We find that the fit is pegged at e = 10 19 cm −3 , the upper limit of the allowed parameter range in relxilld (see the right panel of Fig. C1 ).

Due to an inappropriate value of e , the relxilld model fails to reproduce the soft X-ray spectral shape in the data as shown in the ratio plot in Fig. C1 . Similar conclusions were found in Jiang et al. (2018) : when relxilld was used to model the stacked spectrum of IRAS 13224, the fit was pegged at e = 10 19 cm −3 and provided a significantly worse fit than other reflionx-based models. Therefore, we only use reflionx in this work.

In this work, we present a disc reflection model for IRAS 13224 based on a broken-power law disc emissivity profile. No particular geometry for the corona is assumed in our spectral model.

It has been calculated that a different coronal geometry, e.g. an aborted jet or a sphere enshrouding the inner disc, may lead to a very different disc emissivity profile and thus a different shape of the broad Fe K emission line in a disc reflection spectrum. A broken power law or a double-broken power law is a good approximation for the emissivity profiles of various geometries ( ∝ −q , e.g. Wilkins & Fabian 2012; Gonzalez et al. 2017 ): most of the coronal emission reaches the innermost region of the disc due to strong lightbending effects. Then the emissivity of the disc decreases quickly with radius. At even larger radius, the emissivity profile turns flatter and approaches −3 as in a flat, Euclidean spacetime.

In this section, we investigate whether a lamppost model can explain the spectrum of IRAS 13224. We consider F1 and F5 spectra in this section to demonstrate that a simple lamp-post model is not sufficient to explain the soft X-ray spectra of IRAS 13224.

We replace the relconv model in Model 1 with relconvlp (Dauser et al. 2014 ). The relconvlp model calculates the emissivity profile of a disc illuminated by an isotropic point-like source characterised by ℎ, the height of the corona on the rotation axis of the BH. By applying this model to the F1 spectrum, we find a significantly worse fit with Δ 2 = 86 and two fewer parameters in comparison with Model 1. The model is shown in the left panel of Fig. D1 . Significant residuals can be found below 1 keV and between 7-8 keV. We estimate the uncertainty of ℎ by using the steppar tool in XSPEC, and obtain ℎ < 2.5 g (Δ 2 < 2.706) for F1. Similar conclusions can be found for the highest flux state spectrum where residuals can be found below 1 keV.

In this paper, we adopt a phenomenological broken power-law emissivity as it provides a better fit to the data of IRAS 13224. 1 10 Figure C1 . A fit for the F1 spectrum of IRAS 13224 using the relxilld model. Due to the limited range of e , the fit is pegged at e = 10 19 cm −3 as shown in the right panel, while the reflionx model requires e ≈ 10 20 cm −3 (90% confidence range shown by the green shaded region). Due to an inappropriate value of e , we find it difficult to model the soft X-ray spectrum of F1 by using relxilld. A zoom-in on the data/model ratio plot in the soft X-ray band is shown in the left panel.

The stacked F1+F4 spectrum of IRAS 13224 is shown in the Panel A of Fig. E1 in comparison with F1 and F4 spectra. We follow the same analysing approach as in Jiang et al. (2018) by starting with a single-reflection, single-absorption model. An electron density of 10 15 cm −3 is used as in Chiang et al. (2015) and Jiang et al. (2018) . The full model is tbnew * xstar * ( relconv * reflionx + bbody + nthcomp) in XSPEC notation. The line-of-signt Galactic column density is fixed at H = 6.78 × 10 20 cm −2 (Willingale et al. 2013 ) in this test. We find that this single-reflection, singleabsorption model fails to reproduce the observed blue wing of the broad Fe K emission and the absorption features above 7 keV (see the Panel B of Fig. E1 ). An additional reflection component is able flux interval. We do not use F5 interval in this test because F5 does not show significant UFO absorptions as shown in Section 2.

to significantly improve the fit by Δ 2 = 30 and two more free parameters, which are the ionisation and the normalisation parameters of the second reflection component. The other parameters, including Fe , emissivity indices, * and are linked to the corresponding parameters of the first reflection component. However, a second absorption model is still required to fit the rest of the broad absorption feature above 8.5 keV. The width of the absorption line is due to the averaging effects of rapidly variable UFOs as seen in Parker et al. (2017a) .

In the end, we present the best-fit model for the stacked spectrum IRAS 13224 by considering the same model set-up as in Jiang et al. (2018) . A double-reflection, double-absorption model is required and able to provide a good fit to the stacked spectrum with 2 / = 188.18/156. The best-fit parameters are shown in Table E1 , and the best-fit model is shown in the Panel D of Fig. E1 Table D1 . Best-fit parameters for the low-flux (F1) and high-flux (F5) spectra of IRAS 13224 assuming a lamp-post emissivity profile. This model provides worse fits to the data than a model using broken power-law emissivity profiles. 1 Only upper limits are obtained for this parameter and shown in the brackets. Residuals are seen in the soft X-ray band, e.g. <1 keV. Solid lines: total models; dashed lines: reflection models; dash-dotted lines: Comptonisation models. Top panels show the best-fit models. The bottom two panels show corresponding data/model ratio and residual plots. responding data/model ratio plot can be found in the Panel E of Fig. E1 . So far, we have successfully reproduced the spectral fitting process in Jiang et al. (2018) by considering a stacked spectrum of IRAS 13224 extracted from low and high flux intervals combined. Due to the averaging effects, a double-reflection, double-absorption model is needed. We conclude that, in order to study the spectral variability of a variable AGN like IRAS 13224, one needs to avoid stacking spectra from a wide range of flux intervals. The spectra are unfolded using a model that is constant across the energy band for demonstration purposes only. Panel B: a data/model ratio plot for the F1+F4 spectrum using a single-reflection, single-absorption model. Significant positive residuals are seen in the 7-8 keV band and negative residuals are seen at 9 keV. Panel C: a data/model ratio plot for the same stacked spectrum using a double-reflection, single-absorption model. Negative residuals are still found at 9 keV. Panel D: the best-fit double-reflection, double-absorption model for the stacked spectrum of IRAS 13224. Corresponding data/model ratio plot is shown in the panel E. This model is similar to the ones in Jiang et al. (2018) .

Revisiting Narrow-Line Seyfert 1 Galaxies and their Place in the Universe. p

Astronomical Society of the Pacific Conference Series

simple lamppost model may be required for IRAS 13224. Similar conclusions were found in other AGNs. For example, 1H 0707−495, a sister source of IRAS 13224, shows similar rapid X-ray variability (e.g. Fabian et al. 2009 ) and steep double-broken power-law disc emissivity profiles (Wilkins & Fabian 2011) . Detailed modelling of its emissivity profiles indicates a coronal region that may extend in the radial direction while moving along the rotation axis of the central BH (Wilkins et al. 2014) . Similar modelling of emissivity profiles as in Wilkins et al. (2014) is required to test beyond-lamppost models for IRAS 13224.The reflection fraction is also an indicator of the coronal geometry as shown by both observations (e.g. Gallo et al. 2019 ) and theories (e.g. Dauser et al. 2014) . A consistent calculation of the reflection fraction parameter will need to be taken into account in future beyond-lamppost spectral models for IRAS 13224.Lastly, it is important to note that we estimate the values of √︁ INF / AD based on the measurements of the e and parameters.This parameter increases from 0.43× 2×10 6 BH g in the lowest flux state to 1.71 × 2×10 6 BH g in the highest flux state. These estimations are independent of the choice of the coronal geometry. In Section 3.2, we adapt the over-simplified 'lamp-post' model only to estimate the values of INF / AD . A different geometry may affect the values but will not change our spectral models for IRAS 13224.

IRAS 13224 shows rapid variability in the X-ray band on timescales of kiloseconds (e.g. Alston et al. 2020) . During the XMM-Newton observing campaign in 2016, the observed EPIC-pn count rate changed by up to 2 orders of magnitude on timescales of days (Jiang et al. 2018) , making it one of the most extreme known AGNs.To study its spectral variability, we divide the observations of IRAS 13224 into five flux intervals. In Section 2.2, we introduce a disc reflection model with an electron density of around 10 20 cm −3 to explain each flux-resolved spectrum of IRAS 13224. Alternatively, one may find an equally good fit using a low-density model as introduced in Section A. However, an additional bbody model is required to fit the soft excess emission when e = 10 15 cm −3 is assumed. Either a high-density or a low-density model has only one reflection component. Besides, one photoionisation model xstar is needed to model the blueshifted Fe / absorptions in the 8-9 keV band of flux-resolved spectra.In this section, we argue that the requirement for a doublereflection, double-absorption model as in Chiang et al. (2015) ; Jiang et al. (2018) is due to an averaged spectrum extracted from a wide range of flux states being used in the previous work. To test this idea, we consider a spectrum from F1 and F4 9 intervals combined. Table E1 . Best-fit parameters for F1+F4 spectrum of IRAS 13224 using the same model in Jiang et al. (2018) . The full model is tbnew * xstar1 * xstar2 * ( relconv * (reflionx1 + reflionx2) + bbody + nthcomp).