Towards more accurate synthetic reflection spectra: improving the calculations of returning radiation

We present a new model to calculate reflection spectra of thin accretion disks in Kerr spacetimes. Our model includes the effect of returning radiation, which is the radiation that is emitted by the disk and returns to the disk because of the strong light bending near a black hole. The major improvement with respect to the existing models is that it calculates the reflection spectrum at every point on the disk by using the actual spectrum of the incident radiation. Assuming a lamppost coronal geometry, we simulate simultaneous observations of NICER and NuSTAR of bright Galactic black holes and we fit the simulated data with the latest version of RELXILL (modified to read the table of REFLIONX, which is the non-relativistic reflection model used in our calculations). We find that RELXILL with returning radiation cannot fit well the simulated data when the black hole spin parameter is very high and the coronal height and disk's ionization parameter are low, and some parameters can be significantly overestimated or underestimated. We can find better fits and recover the correct input parameters as the value of the black hole spin parameter decreases and the value of the coronal height increases.


INTRODUCTION
Blurred reflection features are common in the X-ray spectra of stellar-mass black holes in X-ray binaries and supermassive black holes in active galactic nuclei (Fabian et al. 1989;Tanaka et al. 1995;Nandra et al. 2007;Miller et al. 2009;Walton et al. 2013).They are generated in the strong gravity regions of black holes by illumination of a cold accretion disk by a hot corona (Fabian et al. 1989;Reynolds & Nowak 2003;Bambi et al. 2021).X-ray reflection spectroscopy refers to the analysis of these relativistically blurred reflection features and can be a powerful technique to study the accretion process around black holes (De Rosa et al. 2019), measure black hole spins (Reynolds 2021;Bambi et al. 2021), and test Einstein's theory of general relativity in the strong field regime (Tripathi et al. 2019(Tripathi et al. , 2021a)).
The past decade has seen tremendous progress in Xray reflection spectroscopy, thanks both to a new gener-ation of reflection models and new observational facilities.Today the state-of-the-art in reflection modeling is represented by the relativistic models relxill (Dauser et al. 2013;García et al. 2014), reltrans (Ingram et al. 2019), reflkerr (Niedźwiecki & Życki 2008;Niedźwiecki et al. 2019), and kyn (Dovčiak et al. 2004), and by the non-relativistic models xillver (García et al. 2013) and reflionx (Ross & Fabian 2005).Despite the remarkable improvements with respect to the previous generation of reflection models, even these models present a number of simplifications, so caution is necessary when we attempt to infer very precise measurements of accreting black holes.Some simplifications are relatively well studied and they do not affect significantly the final measurements of a system if we select properly the spectra to analyze.For example, current models normally assume Keplerian and infinitesimally thin accretion disks.At least in the case of fast-rotating sources, such a simple model works quite well, without introducing appreciable systematic uncertainties in the final measurements, as it was shown by theoretical studies with GRMHD simulations (Shashank et al. 2022) as well by the analysis of specific sources with models admitting deviations from Keplerian motion and disks with finite thickness (Abdikamalov et al. 2020;Tripathi et al. 2020Tripathi et al. , 2021b;;Jiang et al. 2022).However, such a simple model breaks down and we can have unacceptably large systematic errors if we do not select properly the sources; for instance, if we use our reflection models to analyze sources with a mass accretion rate close to their Eddington limit (Riaz et al. 2020a,b).
The choice of the model for the emissivity profile is more subtle.If we knew the geometry and the properties of the corona, the emissivity profile could be calculated with ray-tracing techniques.For coronae of unknown geometry, it is common to employ a broken power law profile, or even a twice broken power law profile, where the values of the emissivity indices and the breaking radii can be determined by the fit.A twice broken power law profile should be able to approximate well most coronal geometries (Wilkins & Fabian 2012;Gonzalez et al. 2017), and in the analysis of high-quality data we should be able to infer the general behavior of the emissivity profile (see, e.g., Wilkins & Fabian 2011).
The returning radiation (or self-irradiation) is the radiation that is emitted by the accretion disk and returns to the accretion disk because of the strong light bending near the black hole.Cunningham (1976) was the first to present the calculations of the returning radiation of thermal spectra of thin disks.It turns out that a thermal spectrum in which the returning radiation is taken into account appears like a thermal spectrum without returning radiation but a slightly higher mass accretion rate (Li et al. 2005), and therefore the returning radiation is normally ignored in the spectral analysis of thermal components.On the other hand, the effect of the returning radiation is definitively important and non-negligible in the polarization of thermal spectra (Schnittman & Krolik 2009).
In the context of disk reflection, the returning radiation was first studied in Dabrowski et al. (1997), where the authors found that the effect is not important if the corona corotates with the disk.Niedźwiecki & Życki (2008) confirmed this conclusion but also pointed out that the returning radiation is instead important if the corona is static and close to the black hole event horizon.The importance of the returning radiation in the case of lamppost coronae was illustrated in Niedźwiecki et al. (2016) and Niedźwiecki & Zdziarski (2018), where the authors showed in which conditions the effect of the returning radiation can be significant.
The impact of the returning radiation on reflection spectra has been recently investigated in a few studies with different approximations (Wilkins et al. 2020;Dauser et al. 2022;Riaz et al. 2021Riaz et al. , 2023)).All these studies agree that the effect of the returning radiation is important only for fast-rotating black holes (a * > 0.9, assuming that the inner edge of the accretion disk is at the innermost stable circular orbit, or ISCO) and when the corona is compact and close to the black hole (i.e., when the corona illuminates mainly the very inner part of the accretion disk; for a lamppost corona h < 5 r g , where h is the coronal height and r g is the gravitational radius).These studies also agree on the fact that most measurements are not significantly affected if the theoretical model used to analyze the data does not include the effects of the returning radiation.In particular, the returning radiation does not seem to have a significant impact on the estimate of key-parameters like the black hole spin.However, all these studies employ some important simplifications in the calculation of the returning radiation and, at the same time, fast-rotating black holes with compact coronae close to them are the sources for which we can naturally get more precise measurements, as their relativistic reflection features can be stronger.Further studies are thus necessary to figure out the ac-tual impact of the returning radiation on X-ray reflection spectroscopy measurements.
In this manuscript, we present a new model to calculate reflection spectra of thin accretion disks in Kerr spacetimes produced by lamppost coronae.Our model does not exploit the existence of the Carter constant and therefore it can calculate reflection spectra of thin accretion disks in other stationary, axisymmetric, and asymptotically-flat spacetimes with minor modifications (Bambi 2022).If compared to the existing models, it presents two important improvements: i) the nonrelativistic reflection spectrum at every point on the disk is calculated by taking into account the actual spectrum of the incident radiation, which consists of a power law spectrum from the corona and a reflection spectrum from the returning radiation, and ii) the effect of the returning radiation is calculated up to the fourth order, when we find that further iterations are negligible.We note that point ii) is important to infer the actual spectrum of the incident radiation at point i), while it does not appreciably change the estimate of the emissivity profile.
We use the reflection spectra calculated by our new model to simulate simultaneous observations of bright Galactic black holes with NICER and NuSTAR.The simulated data are fit with the latest version of relxill, where the effect of the returning radiation is included only in the emissivity profile generated by lamppost coronae (Dauser et al. 2022).From the spectral analysis of the simulated data, we find that we cannot obtain good fits when the value of the black hole spin parameter is very high and the values of the coronal height and of the disk's ionization parameters are low, and in such a case the estimate of some parameters can be significantly biased.The quality of the fits improves and the systematic uncertainties decrease as we decrease the value of the black hole spin parameter and we increase the value of the coronal height.If we increase the value of the disk's ionization parameter, the quality of the fits improves, but mainly because the reflection features become weak, while we can have still problems to recover the correct input parameters.Our results do not contradict previous studies on the impact of the returning radiation (Dabrowski et al. 1997;Niedźwiecki & Życki 2008;Niedźwiecki et al. 2016;Niedźwiecki & Zdziarski 2018;Wilkins et al. 2020;Dauser et al. 2022;Riaz et al. 2021Riaz et al. , 2023)), but provide a more clear and solid picture of this issue on an important region of the parameter space.
The manuscript is organized as follows.In Section 2, we present our new model to calculate reflection spectra of thin disks in Kerr spacetimes and, in Section 3, we show its predictions and the impact of the returning radiation in reflection spectra.In Section 4, we simulate some observations of bright Galactic black holes with NICER and NuSTAR and we present the spectral analysis of these simulated data.Discussion and conclusions are reported in Section 5.In Appendix A and Appendix B, we report the key-equations to calculate the emissivity profile generated by a lamppost corona and how the returning radiation illuminates the disk, respectively.In Appendix C, we present an extension of our model in which the ionization parameter of the disk is calculated self-consistently as a function of the radial coordinate from the actual incident radiation and the disk's electron density.

DESCRIPTION OF THE MODEL
This section describes how our model calculates reflection spectra of thin disks.In what follows, we assume that the spacetime metric is described by the Kerr solution, but our model does not exploit any specific property of the Kerr metric (e.g., the existence of the Carter constant) and we could implement another stationary, axisymmetric, and asymptotically-flat metric only by changing the metric coefficients.The accretion disk is supposed to be infinitesimally thin, on the plane perpendicular to the black hole spin.The inner edge of the disk is set at the ISCO radius, r in = r ISCO , and the outer edge is set at the radial coordinate r out = 500 r g .The motion of the material in the disk is Keplerian.We employ the lamppost coronal model, so the corona is a point-like source along the black hole spin axis and is only regulated by its height h.The emission of the corona is isotropic in its rest-frame and its spectrum is described by a power law with a high-energy cutoff, so there are two parameters: the photon index Γ and the high-energy cutoff E cut (which is measured in the rest-frame of the corona).With such a setup, we report below step by step the calculations of our model.

Illumination of the disk by the corona
To calculate how the corona illuminates the accretion disk, we modify the model blacklamp (Abdikamalov et al. 2019;Riaz et al. 2022), which is available on GitHub9 and Zenodo (Abdikamalov et al. 2024a).
The key-equations are reported in Appendix A. First, we identify 100 reference radii on the accretion disk.We employ the following algorithm where i = 0, 1, ..., i max − 1 and here i max = 100.If we want to change the number of reference radii and/or the radial coordinates of the inner/outer edges of the accretion disk, we have just to change the values of, respectively, i max , r in , or r out in Eq. (1).At this point, we consider the locally-Minkowskian reference frame of the corona and we fire photons to hit the disk at the 100 reference radii in Eq. (1).We require that the difference between the radial coordinate at which the photon hits the disk, say ri , and the desired reference radius, r i , is less than 10 −6 r g .When we meet this condition, we store the value of ri (which becomes our new reference radius i), we calculate the redshift g i , and we fire one more photon to hit the disk near ri and calculate the energy density illuminating the disk surface at ri (see Appendix A).
2.2.Non-relativistic reflection spectrum (0th order) For every reference radius ri , we calculate the nonrelativistic reflection spectrum in the rest-frame of the material of the disk.Here and in the next steps, we use the code of reflionx (Ross & Fabian 2005), which can calculate the non-relativistic reflection spectrum for an arbitrary input spectrum illuminating the disk.The spectrum of the radiation illuminating the reference radius i is described by the photon index Γ and the redshifted high-energy cutoff E cut,i = g i E cut , where g i is the redshift factor between the corona and ri .At the end, we have 100 non-relativistic reflection spectra, one for every reference radius.
A couple of comments are in order here.First, we use the code of reflionx, not the public table of reflionx, because the latter assumes that the incident radiation illuminating the disk has a spectrum described by a power law with a high-energy cutoff.While this would be enough at this step (non-relativistic reflection spectrum at the 0th order), in the next steps the illuminating spectrum includes a returning radiation component (see below).Second, the choice of the non-relativistic reflection model is not crucial for the purposes of this work because we do not analyze real data and we only compare theoretical predictions among different models that employ the same non-relativistic reflection calculations.We could have used the code of xillver (but, again, not the public table of xillver) or any other non-relativistic reflection model that can calculate the output spectrum for an arbitrary input spectrum.We do not discuss here the differences between the predictions of reflionx and xillver, but a work on the subject is in preparation.

Returning radiation
We consider a point at the reference radius ri and we calculate how this point illuminates the disk.The key-equations are reported in Appendix B. We consider the locally-Minkowskian reference frame of this point and we fire photons isotropically (our photon grid is 4000×4000).The photon trajectories are calculated in Boyer-Lindquist coordinates, from our reference point to the accretion disk (returning radiation), the black hole event horizon or the plunging region (radiation lost into the black hole or the plunging region), or to some large radial coordinate faraway from the disk and the black hole (radiation escaping to infinity).The trajectories of the photons hitting the disk are also the trajectories that contribute to the returning radiation of our reference point (if we consider the time reversal process).To infer the total spectrum of the returning radiation, we sum up all redshifted non-relativistic reflection spectra from those trajectories, where the non-relativistic reflection spectra at the emission points are obtained by interpolation from the spectra at the 100 reference radii calculated in Subsection 2.2.We repeat these calculations for every reference radius.At the end, we have 100 spectra of the returning radiation, one for every reference radius ri .

Non-relativistic reflection spectrum (1st order)
We use the code of reflionx to calculate the nonrelativistic reflection spectra for the 100 reference radii generated by the combination of the radiation from the corona and the 1st order returning radiation.At the end, we have 100 non-relativistic reflection spectra, one for every reference radius.Note that at this step it would be incorrect to use the public table of reflionx because the spectrum illuminating the disk is not a power law with a high-energy cutoff as at the 0th order.
2.5.Non-relativistic reflection spectrum (higher orders) To calculate the returning radiation at higher order, we do not to need to repeat all calculations in Subsection 2.3, because the photon trajectories do not change.For every reference radius, we have to recalculate the total spectrum of the returning radiation, which is obtained by summing up all redshifted non-relativistic reflection spectra (1st order) from the emission points illuminating the points at the reference radii.At the end, we have 100 spectra of the returning radiation (2nd order), one for every reference point.
We consider two more iterations and for every reference point we obtain the spectra of the returning radiation at the 3rd and 4th orders.We do not go beyond the 4th order because we do not see any significant change in the spectra if we include higher order corrections.
2.6.Relativistic reflection spectrum At this point, we have a non-relativistic reflection spectrum at 4th order for every reference radius of the accretion disk.We use the ray-tracing code blackray10 (Abdikamalov et al. 2024b) to calculated the relativistic reflection spectrum of the whole disk as seen by a distant observer with viewing angle ι.The ray-tracing code has a grid on the image plane of the distant observer.From every point of this grid, we fire a photon and we calculate the photon trajectory backward in time from the image plane of the observer to the emission point on the accretion disk.We can thus calculate the redshift associated to that point of the grid and we can determine its reflection spectrum at 4th order.At the end, we have the image of the accretion disk and every pixel of this image has its own spectrum.The total spectrum of the source is obtained by summing up the spectra of all pixels.These calculations have been already extensively discussed in the literature and the reader can find all the details, for instance, in Bambi (2017).

Computational time
To calculate a full setup, our model takes 3-4 hours using 80 threads.The most time-consuming part is the calculation of the non-relativistic reflection spectra with the reflionx code, while the ray-tracing calculations with blacklamp and blackray are relatively fast.

PREDICTIONS
In this section, we present the predictions of our new model.The ray-tracing calculations are compared with those in Dauser et al. (2022) and we show that the predictions of the two models are in good agreement.The two models solve different equations (our model does not exploit the properties of the Kerr solution and solves equations valid for any stationary, axisymmetric, and asymptotically-flat spacetime), but the theoretical framework is almost the same.The key-feature of our new model is that it calculates the non-relativistic reflection spectrum at every point on the disk by using the actual spectrum of the incident radiation, which is the combination of a power law spectrum from the lamppost corona and of the sum of reflection spectra from different parts of the disk returning to the disk.-Fraction of photons returning to the accretion disk (returning radiation, red curves), falling onto the black hole or the plunging region (black curves), and escaping to infinity (blue curves) as a function of the radial coordinate of the emission point on the disk.The solid curves are for the predictions of the model presented in this manuscript and the dashed curves are for the predictions of the model presented in Dauser et al. (2022).The photon trajectories depend on the spacetime metric and the results shown here are for a Kerr spacetime with spin parameter a * = 0.998.
First, from the ray-tracing calculations described in Subsection 2.3, for every point on the disk emitting reflection photons, we can evaluate the fraction of photons that returns to the disk (returning radiation), falls onto the black hole or the plunging region, and escapes to infinity.Fig. 1 compares our results with those in Dauser et al. (2022).The dashed curves for Dauser et al. (2022) are difficult to see in the plot because of the very good agreement between the two codes.The figure shows the calculations for a Kerr spacetime with spin parameter a * = 0.998, where r in = 1.237 r g .Most of the photons emitted from the very inner part of the accretion disk (r < 2 r g ) return to the disk or fall onto the black hole and the plunging region.Already for r > 2 r g , most of the photons can escape to infinity and their fraction approaches 1 as the radial coordinate increases.
Figs. 2 and 3 show some non-relativistic reflection spectra as calculated in Subsections 2.2 and 2.5 assuming that the spectrum of the corona is described by a power law with photon index Γ = 1.7 and high-energy cutoff E cut = 300 keV, the ionization of the disk is ξ = 100 erg cm s −1 (in Fig. 2) and 1000 erg cm s −1 (in Fig. 3), and the black hole spin parameter is a * = 0.998.The figures show the results for three different coronal heights (h = 2, 5, and 10 r g ) at three different radial coordinates on the accretion disk (R = 1.3, 3, and 10 r g ).The blue solid curves are the non-relativistic reflection spectra without returning radiation (Subsection 2.2) and the red dashed curves are for the non-relativistic reflection spectra with the returning radiation calculated at the 4th order (Subsection 2.5).In general, the iron line region is not particularly affected by the returning radiation, while larger differences between the spectra without and with returning radiation appear at lower energies (below 2-3 keV) and around the peak of the Compton hump.
Figs. 2 and 3 show even the spectra of the incident radiation: the black solid curves are for the direct incident radiation from the corona, while the dashed black curves are for the total incident radiation illuminating the disk (i.e., direct radiation from the corona + returning radiation from the disk).As we can see from both figures, the total incident spectrum can be significantly different from a power law spectrum, so it should not be a surprise that the non-relativistic reflection spectrum can be appreciably affected by the returning radiation.In some cases, at some energies the flux of the incident radiation with returning radiation can be more than an order of magnitude higher than that without returning radiation.For ξ = 100 erg cm s −1 , the contribution of the returning radiation is mainly at very low energies and around the peak of the Compton hump.For ξ = 1000 erg cm s −1 , the returning radiation seems to contribute appreciably over the whole X-ray spectrum.
Fig. 4 shows the direct incident X-ray flux from the corona (blue solid curves) and the incident flux of the returning radiation calculated at the 4th order (red dashed curves) for a coronal height h = 2, 5, and 10 r g in Kerr spacetime with a * = 0.998.The spectrum of the corona is described by a power law with photon index Γ = 1.7 and high-energy cutoff E cut = 300 keV.Fig. 4 can be compared with Fig. 5 in Dauser et al. (2022) and, even in this case, the two models seem to agree.If the corona is very close to the black hole, the flux of the returning radiation is lower than the direct flux from the corona in the inner part of the accretion disk and higher in the outer part (left panel in Fig. 4).For higher values of the corona height, we can see the opposite effect, namely the incident flux of the returning radiation can exceed that from the corona in the inner part, even if only in a very small region of the accretion disk, and is lower in the outer part (right panel in Fig. 4).
Figs. 2, 3, and 4 can be qualitatively explained as follows; see also Dauser et al. (2022), where such a behavior was found and discussed for the first time.If the corona is close to the black hole (h is low), the effect of light bending is strong and most of the radiation illuminates the inner part of the accretion disk, while at larger radii the disk is not illuminated well.If the corona is not close to the black hole (h is high), the effect of light bending is weak and there are no dramatic differences between the illumination of the disk at small and large radii.The returning radiation is mainly produced by the inner part of the accretion disk (as shown even in Fig. 1), but such a radiation does not necessarily return to the inner part of the disk.The combination of these effects leads to what we see in Figs. 2, 3, and 4.
• For h = 2 r g , the corona strongly illuminates the inner part of the accretion disk, while the contribution of the returning radiation is subdominant (in Figs. 2 and 3, the difference between the spectra without and with returning radiation is small for R = 1.3 r g and in Fig. 4 the red dashed curve is below the blue solid curve at small radii).Since the corona directly illuminates well the inner part of the accretion disk, there are many photons that can return to the disk at large radii, where the illumination from the corona is weak and therefore the returning radiation is dominant (in Figs. 2 and 3, the difference between the spectra without and with returning radiation is large for R = 10 r g and in Fig. 4 the red dashed curve is above the blue solid curve at large radii).
• For h = 10 r g , the corona illuminates the disk more uniformly.At small radii, the contribution of the returning radiation is relatively high because such a region is not strongly illuminated by the corona (in Figs. 2 and 3, the difference between the spectra without and with returning radiation is relatively large for R = 1.3 r g and in Fig. 4 the red dashed curve is above the blue solid curve at small radii).At large radii, we have the opposite effect because the inner disk is not strongly illuminated and therefore there are not many photons returning to the disk and, at the same time, the outer part of the disk is already illuminated well by the corona (in Figs. 2 and 3, the difference between the spectra without and with returning radiation is small for R = 10 r g and in Fig. 4 the red dashed curve is below the blue solid curve at large radii).
Relativistic reflection spectra without returning radiation and with returning radiation calculated at the 1st, 2nd, 3rd, and 4th orders are shown in Fig. 5, Fig. 6, Fig. 7, and Fig. 8. Every figure shows the results for three different coronal heights (h = 2, 5, 10 r g ) and three different disk's ionization parameters (ξ = 100, 1000, and 10000 erg cm s −1 ).The inclination angle of the disk with respect to the line of sight of the observer is ι = 60 • in Fig. 5 and Fig. 6 and ι = 30 • in Fig. 7 and Fig. 8.The photon index of the coronal spectrum is Γ = 1.7 (Fig. 5 and Fig. 7) and 2.5 (Fig. 6 and Fig. 8).The high-energy cutoff is always E cut = 300 keV.We note that all spectra are normalized to have the peak of the iron Kα line at νF ν = 1 in order to compare better the shapes of the spectra (without such a normalization, the flux of the spectra without returning radiation would be lower than the flux of the spectra with returning radiation, as the latter increases the incident flux illuminating the disk).The height of the corona is the most important parameter to regulate the difference between the spectra without and with returning radiation.This is understandable, because photons returning to the disk are mainly emitted from the very inner part of the accretion disk, and the emissivity of the inner part of the accretion disk increases as the coronal height decreases.

SIMULATIONS AND FITS
In the latest version of the relxill package, the effect of the returning radiation is taken into account only in the emissivity profile, while the non-relativistic reflection spectra are still those calculated assuming that the incident radiation has a power law spectrum with a highenergy cutoff (Dauser et al. 2022).Such an approximation is motivated by the fact that the actual incident spectrum, which is the result of the combination of the direct spectrum from the corona and the returning radiation, has a shape rather close to a power law; see Fig. 7 in Dauser et al. (2022).The key-parameters are the black hole spin, the coronal height, and the disk's ionization parameter.The black hole spin and the coronal height regulate the intensity of the returning radiation.The -Synthetic reflection spectra in the rest-frame of the material in the disk at different radial coordinates (R = 1.3, 3, and 10 rg) without returning radiation (blue solid curves, label "0") and with returning radiation calculated with four iterations (dashed red curves, label "4").In every panel, we also show the spectrum of the incident radiation from the corona (black solid curves, label "Direct") and the spectrum of the total incident radiation (direct radiation from the corona + returning radiation from the disk, black dashed curves, label "Total").The calculations assume that the spacetime is described by the Kerr metric with spin parameter a * = 0.998, the spectrum of the corona at the emission point is described by a power law with photon index Γ = 1.7 and high-energy cutoff Ecut = 300 keV, and the disk's ionization parameter is ξ = 100 erg cm s −1 .The y-axis shows the quantity νFν in arbitrary units, where ν is the photon frequency and Fν is the observed flux at the frequency ν.The coronal height is h = 2, 5, and 10 rg for the upper, central, and lower panels, respectively.See the text for more details.
ionization parameter regulates the shape of the reflection spectrum: if the ionization is low, there are many features in the reflection spectrum, but these features are washed out when the value of the ionization parameter becomes very high (and there are no features when the disk is fully ionized).
To quantify the capability of the relxill package to analyze real sources, we simulate some observations of bright Galactic black holes.We assume that the energy flux of the source is 10 −8 erg cm −2 s −1 in the 1-10 keV energy range.In XSPEC language, the model of the simulations is tbabs×(cutoffpl + reflection) , where tbabs takes the Galactic absorption into account (Wilms et al. 2000), cutoffpl describes the Comptonized spectrum of the corona, and reflection is the reflection spectrum calculated by our model.The input parameters are summarized in Tab. 1.We consider two coronal heights (h = 2 and 5 r g ) and three possible disk's ionization parameters (ξ = 100, 1000, and 10000 erg cm s −1 ).The other input parameters are the same for all simulations.The black hole spin parameter is set to a * = 0.998 to maximize the impact of the returning radiation.
For the same input parameters, we run a first simulation without returning radiation and a second simulation with returning radiation calculated at the 4th order.We assume simultaneous observations of NICER and NuS- TAR and we consider an exposure time of 5 ks for NICER and of 30 ks for NuSTAR (for simplicity, we run the simulations assuming 60 ks with one detector rather than 30 ks with two detectors).
We fit the simulated data with the model where relxilllp reflionx is the lamppost flavor of the latest version of the relxill package (Dauser et al. 2022), in which the returning radiation can be turned off and on, and has been modified to read the reflionx table instead of the xillver one.In this way, we use the same model as in the simulations to calculate nonrelativistic reflection spectra and we include the flux correction factor discussed in Appendix B in Dauser et al. (2022).If we used relconvlp×reflionx, the fitting model would not include such a flux correction factor, while such an effect is automatically included by construction in the calculations of our new model.
In relxilllp reflionx we turn the returning radia-tion off to fit the simulations without returning radiation and we turn it on to fit the simulations with returning radiation calculated at the 4th order.For the simulations without returning radiation, we recover the correct input parameters and the fits are good, in the sense that we do not see unresolved features in the plots of the ratios between the simulated data and the best-fit models.This is just a healthy check for our new model and we do not report here the results because the latter are perfectly consistent with our expectations of good fits and correct estimate of the model parameters.
For the simulations with returning radiation, the quality of the fits is not always good and some parameters are significantly underestimated or overestimated.The results of these fits with returning radiation are summarized in Tab. 2 and the residuals are shown in Fig. 9.In the data to the best-fit model ratios, there are large residuals for simulations A1 and A2, there are still some residuals for simulations B1 and B2, while there are no clear residuals for simulations A3 and B3 with a high value of the ionization parameter ξ.The discussion of -Direct incident X-ray flux from the corona (blue solid curves) and incident X-ray flux generated by the returning radiation as calculated at the 4th order (red dashed curves) for a coronal height h = 2 (left panel), 5 (central panel), and 10 rg (right panel).These results are obtained assuming that the spacetime geometry is described by the Kerr solution with spin parameter a * = 0.998 and that the spectrum of the corona is described by a power law with photon index Γ = 1.7 and high-energy cutoff Ecut = 300 keV.
these fits is postponed to the next section.

DISCUSSION AND CONCLUSIONS
For a qualitative discussion of the impact of the returning radiation on observed reflection spectra, we can refer to Fig. 5, Fig. 6, Fig. 7, and Fig. 8.After fixing the black hole spin parameter to a * = 0.998, the two key-parameters are the coronal height h and the disk's ionization parameter ξ, while the impact of the values of the photon index Γ and of the inclination angle of the disk ι is weaker.Within the lamppost setup, the coronal height h is the parameter regulating the intensity of the returning radiation simply because it determines the emissivity profile of the direct radiation from the corona: as shown in Fig. 1 (where the three curves depend only weakly on the black hole spin parameter except for the fact a * determines the inner edge of the disk, r in ) only a significant fraction of photons emitted at small radii returns to the disk, and therefore the effect of the returning radiation increases/decreases when the emission at small radii increases/decreases.The ionization parameter ξ determines scattering probabilities and energy lines, but its impact on the spectra is more complicated.For ξ = 100 erg cm s −1 , the returning radiation decreases the relative intensity of the reflection features in the soft X-ray band and increases the intensity of the Compton hump, while we see the opposite trend for ξ = 10000 erg cm s −1 .On the other hand, for ξ = 1000 erg cm s −1 , it seems like the intensity of the reflection features in the soft X-ray band and the Compton hump are both increased with respect to the iron line.
For a quantitative discussion of the impact of the returning radiation, we have the simulations and fits presented in Section 4. We remind the reader that our simulations are for a Kerr spacetime with spin parameter a * = 0.998, which represents the situation in which the observed spectra are more strongly affected by the returning radiation.Moreover, we are assuming a static corona: in the case of a rotating corona, the impact of the returning radiation would be weaker (Dabrowski et al. 1997;Niedźwiecki & Życki 2008).First, we note that the quality of the fits with h = 2 r g and a low or moderate value of the ionization parameter ξ is bad: the reduced χ 2 is unacceptably high and there are large residuals in the ratio plots.For simulation A1 with ξ = 100 erg cm s −1 , the coronal height is clearly overestimated and the ionization parameter ξ, the iron abundance A Fe , and the black hole spin parameter a * are underestimated.For simulation A2 with ξ = 1000 erg cm s −1 , the fit recovers relatively well the coronal height and the black hole spin parameter, even if the quality of the fit is still bad.If we increase the coronal height to h = 5 r g , the quality of the fits improves: in simulations B1 and B2, the reduced χ 2 is around 1.5 and in the ratio plots we do not see the large residuals of the fits of simulations A1 and A2.The value of the coronal height is recovered quite well, while the estimate of the spin is not consistent with the spin parameter of the simulations (see Tab. 2).In the two cases with a high ionization parameter (A3 and B3), the quality of the fits is good, but mainly because the reflection features are weak.In the fit of simulation A3, the coronal height is clearly overestimated, the spin parameter is unconstrained, and the inclination angle of the disk is underestimated.In the case of the fit of B3, the spin parameter is poorly unconstrained and we do not recover the correct value at the 90% confidence level.So for high values of the ionization parameter we may not recover the correct input parameter without a clear signal in the residuals.
We would like to stress that the main feature of our model is that we calculate the reflection spectrum at ev- -Synthetic reflection spectra without returning radiation (blue solid curves, label "0") and with returning radiation calculated with one (red solid curves, label "1"), two (red dashed curves, label "2"), three (green dashed-dotted curves, label "3"), and four (black dotted curves, label "4") iterations.We consider three different values of coronal height (h = 2, 5, 10 rg) and three different values of ionization parameter of the disk (ξ = 100, 1000, and 10000 erg cm s −1 ).The calculations assume that the spacetime is described by the Kerr metric with spin parameter a * = 0.998, the spectrum of the corona at the emission point is described by a power law with photon index Γ = 1.7 and high-energy cutoff Ecut = 300 keV, and the inclination angle of the disk with respect to the line of sight of the observer is ι = 60 • .The y-axis shows the quantity νFν in arbitrary units (ν is the photon frequency and Fν is the observed flux at the frequency ν) and normalized to have the peak of the iron Kα line at 1 in order to facilitate the comparison of the shapes of the spectra.See the text for more details.
ery radial coordinate of the accretion disk taking into account the actual spectrum illuminating the disk, which is a combination of a power law spectrum from the corona and a reflection spectrum returning to the disk.The four iterations described in Subsection 2.5 are necessary to find the final reflection spectrum, while higher order effects have a negligible impact on the emissivity profile of the disk.To prove such a statement, we rerun our model assuming that the incident spectrum illuminating the disk is always a power law with a high-energy cutoff, and we take the returning radiation into account only in the intensity of the radiation, which is the approximation adopted in Dauser et al. (2022).The results for Γ = 1.7 and ι = 60 • are shown in Fig. 10, where we clearly see that the impact of the returning radiation is much weaker than in Fig. 5 and the 1st order calculations are enough, in agreement with the conclusion reported in Dauser et al. (2022).We note that we have checked our emissivity profile with such an approximation with the emissivity profile of relxilllp reflionx with returning radiation and the two profiles agree.
In the previous section, all simulations assumed a * = 0.998 to maximize the impact of the returning radiation.If we decrease the value of the black hole spin parameter, the inner edge of the accretion disk (which is supposed to be at the ISCO radius in this work) moves to larger radii.Since it is mainly the radiation emitted at very small radii to return to the disk (see Fig. 1), we automatically decrease the impact of the returning radiation on the observed spectrum.To be more quantitative, we run another cycle of simulations with a * = 0.9.The input parameters are summarized in Tab. 3 and they are exactly the same as in Tab. 1 except for the value of the black hole spin parameter.We assume the same properties for the sources and for the observations of NICER and NuS-TAR.We fit the simulated data with tbabs×(cutoffpl + relxilllp reflionx) and the ratio plots are shown in Fig. 11 (left panels for h = 2 r g , right panels for h = 5 r g ).
Tab. 4 shows the best-fit values of the fits with returning radiation.The quality of the fit is definitively better than the fits for the simulations with a * = 0.998: the reduced χ 2 is relatively close to 1 even in simulations C1 and C2, in Fig. 11 we do not see the large residuals of Fig. 9, and we substantially recover all the input parameters (with some exceptions).
In Riaz et al. (2021), it was shown that, when we fit synthetic spectra that include the returning radiation with a reflection model that does not include the returning radiation, the results of the fits can depend significantly on the considered energy band.That study was limited to neutral disks.With the model presented in the current paper, we can check if such a conclusion remains true for a more realistic ionized disk.We thus repeat some fits assuming that we have only NICER observations.We refit simulations A1 (h = 2 r g , a * = 0.998, and ξ = 100 erg cm s −1 ), A3 (h = 2 r g , a * = 0.998, and ξ = 10000 erg cm s −1 ), C1 (h = 2 r g , a * = 0.9, and ξ = 100 erg cm s −1 ), and C3 (h = 2 r g , a * = 0.9, and ξ = 10000 erg cm s −1 ).The results of these new fits are reported in Tab. 5 and Fig. 12.The results are quite different from the fits of the simulated observations with NICER+NuSTAR.We often recover the correct input parameters and there are no very large residuals in the ratio plots.In particular, we recover the correct black hole spin parameter in A1, A3, and C1, while a * cannot be constrained in C3.The values of the coronal height h and of the inclination angle i are recovered correctly (but in the fit of simulation A3 the value of h is poorly constrained).We thus recover the conclusions of Riaz et al. ( 2021) for ionized disks: if we have only data in the soft X-ray band, the impact of the returning radiation is weak.
We conclude this section with two considerations.First, in the model presented in this paper, the returning radiation is only represented by reflection photons.This is a good approximation for accretion disks of supermassive black holes (where the thermal spectrum of the disk is in the UV band) and of Galactic black holes in the hard state and at low luminosities (where the thermal component is negligible with respect to the continuum from the corona and the reflection spectrum).In general, even the thermal radiation could contribute to the returning radiation, and it is also possible that the returning radiation of the thermal spectrum acts as a corona in sources in the soft state (see, for instance, Connors et al. 2021).In a future paper, we plan to present an extension of our model that includes the calculations of the thermal spectrum and where the returning radiation is generated by both the reflection and thermal components.
Second, our new model can calculate reflection spectra with returning radiation calculated with the correct spectrum illuminating the disk at every radial coordinate, but it is too slow to be used to analyze real X-ray spectra.With the current architecture of relxill and the other reflection models, there is no straightforward way to include the returning radiation as calculated in the present work.This is because in all these models reflection spectra are obtained by calculating the integral (see, for instance, Bambi 2017) where g is the redshift factor, g * is the relative redshift factor, f is the transfer function, and I e is the specific intensity of the radiation in the rest-frame of the material in the disk.The transfer functions (encoding all de- tails about photon trajectories and particle motion in the disk) are tabulated in a FITS file for a grid of black hole spins and viewing angles.Non-relativistic reflection spectra and lamppost emissivity profiles, which determine I e , are tabulated in other two FITS files.These quantities are tabulated in FITS files to be able to calculate many spectra quickly during the data analysis process, when we need to scan the parameter space and find the bestfit values of all free parameters.Especially the calculation of the non-relativistic reflection spectra, which requires to solve radiative transfer equations, is very timeconsuming.In order to include the calculations of the returning radiation with the actual incident spectrum, we would need to add to the model a very heavy FITS file, beyond the capability of the RAM of current normal computers.Our plan for the near future is to develop a reflection model based on a library of reflection spectra.
While the construction of such a library is certainly timeconsuming, it can be done on a large computer cluster.
The model can then use this library to calculate quickly many spectra and analyze real data.The calculation of the emissivity profile produced by a lamppost corona is already discussed in the literature.Here we simply report the main steps and the key-equations.-As in Fig. 5 but calculating the reflection spectrum assuming that the spectrum of the radiation illuminating the disk is described by a power law with a high-energy cutoff as in Dauser et al. (2022).
In the locally-Minkowski reference frame of the corona, the initial conditions of the photon 4-momentum can be written as where δ is the angle between the black hole spin axis and the direction of the photon 4-momentum.In Boyer-Lindquist coordinates, the initial conditions of the photon 4-momentum become where E is the conserved photon energy and can be inferred from the condition k µ k µ = 0.The redshift factor of the photon emitted from the corona and hitting the disk at the radial coordinate r is where k d µ is the photon 4-momentum on the disk, u µ d is the 4-velocity of the material in the accretion disk, k c µ is the photon 4-momentum at the corona, and u µ c is the 4-velocity of the corona.The intensity of the radiation illuminating the disk is TABLE 3 Summary of the input parameters of the simulations with a * = 0.9.In simulations C, the coronal height is h = 2 rg, while in simulations D we have h = 5 rg.The disk ionization parameter is ξ = 100, 1000, and 10000 erg cm s −1 for, respectively, simulations 1, 2, and 3   where dΩ o is the solid angle in the source rest-frame, dA o is the surface area in the reference frame of the disk patch, γ is the Lorentz factor of the material in the accretion disk as measured in the locally-non-rotating reference frames    and Ω K is the angular velocity of the material in the accretion disk.Last, the energy density illuminating the disk surface (= emissivity) at the reference radius ri is given by The calculation of the returning radiation is described in Subsection 2.3.For every reference radius, we consider a point on the disk and we fire photons isotropically.In the locally-Minkowski reference frame of the emission point, the initial conditions of the photon 4-momentum can be written as k (α) = (−E 0 , E 0 sin θ e cos ϕ e , E 0 sin θ e sin ϕ e , E 0 cos θ e ) , where θ e and ϕ e are the emitting angles.In the locally-non-rotating reference frames, the photon 4-momentum k where γ is still the Lorentz factor of the material in the accretion disk given in Eq. (A5) and u (ϕ) is the ϕ component of the 4-velocity of the material in the disk as measured in the locally-non-rotating reference frames.The photon 4-momentum in Boyer-Lindquist coordinates are where E is the conserved photon energy and can be inferred from the condition k µ k µ = 0.The total X-ray flux of the returning radiation at the k-th order illuminating the reference point at the radial coordinate r is F k ret (r, E) = j g 3 j I k e (r j , E g j ) ∆Ω j cos θ j , (B4) where the summation is over all the photons that return to the disk.The photon j returns to the disk from the radial coordinate r j and with redshift factor g j .θ j is the incident angle of the photon j at the radial coordinate r.The total X-ray flux illuminating the reference point at the radial coordinate r at the k-th order is the sum of the direct X-ray flux from the corona and X-ray flux from the returning radiation at the k-th order In the calculations presented so far in this manuscript, the disk electron density is always assumed to be n e = n 0 = 10 15 cm −3 at any radial coordinate.The value of the ionization parameter ξ changes (we have presented most results for three cases, ξ = 100, 1000, and 10000 erg cm s −1 ), but is also constant over the whole disk.However, for a real accretion disk we can expect that the electron density has a non-trivial radial profile, so n e = n e (r).At the same time, the ionization parameter is defined as where F X is the total X-ray flux illuminating the disk (direct X-ray flux from the corona + returning radiation).In general, both F X and n e are functions of the radial coordinate r and therefore even ξ should be a function of r.For the lamppost model, we can calculate F X (r).If we assume some radial profile for the disk electron density, we can determine ξ = ξ(r) and use the correct ionization parameter at every point of the disk in the calculation of reflection spectra 11 .We note that, depending on the radial profiles of L X and n e , ξ(r) may not be a monotonic function of the radial coordinate r, and the maximum value of ξ, say ξ max , may not be at r in .
It is straightforward to implement a non-trivial disk electron density profile in our new model and calculate the ionization parameter at every radial coordinate of the disk with Eq. (C1).Fig. 13 shows synthetic reflection spectra without and with returning radiation for two electron density profiles.In the upper panels, we still assume a trivial profile, where the electron density is constant over the disk.In the lower panels, we consider the density profile of the Shakura-Sunyaev model (Shakura & Sunyaev 1973) where the minimum of n e (r) is set to n 0 = 10 15 cm −3 .The height of the corona determines the profile of F X , but not its normalization, which depends on the coronal luminosity.In our model, we regulate the coronal luminosity by changing the maximum value of the disk's ionization parameter, ξ max .In Fig. 13, we assume ξ max = 10000 erg cm s −1 .The calculations are done assuming that the spectrum of the corona is described by a power law with photon index Γ = 1.7 and high-energy cutoff E cut = 300 keV, and that the inclination angle of the disk with respect to the line of sight of the observer is ι = 60 • .Fig. 14 shows the ionization parameter profiles for the two choices of the electron density profile of Fig. 13 and the coronal heights h = 2, 5, and 10 r g .The exact choice of the electron density profile has a strong impact on the profile of ξ, while the impact of the returning radiation is much weaker.The returning radiation slightly increases the value of the ionization parameter resulting from the illumination of the corona.5 with the ionization parameter ξ calculated with Eq. (C1).In the upper panels, we still assume the same electron density over the whole disk, ne = n 0 = 10 15 cm −3 .In the lower panels, we assume the electron density profile of the Shakura-Sunyaev model, ne(r) ∝ r 3/2 [1 − (r in /r) 1/2 ] −2 .See the text for more details.13.We have two models of the disk electron density profile: ne = n 0 = 10 15 cm −3 (dashed curves) and the Shakura-Sunyaev profile ne(r) ∝ r 3/2 [1 − (r in /r) 1/2 ] −2 (solid curves).The blue curves are the ionization parameter profiles without returning radiation and the red curves the profiles calculated taking the returning radiation into account.
Fig.1.-Fraction of photons returning to the accretion disk (returning radiation, red curves), falling onto the black hole or the plunging region (black curves), and escaping to infinity (blue curves) as a function of the radial coordinate of the emission point on the disk.The solid curves are for the predictions of the model presented in this manuscript and the dashed curves are for the predictions of the model presented inDauser et al. (2022).The photon trajectories depend on the spacetime metric and the results shown here are for a Kerr spacetime with spin parameter a * = 0.998.
Fig.2.-Synthetic reflection spectra in the rest-frame of the material in the disk at different radial coordinates (R = 1.3, 3, and 10 rg) without returning radiation (blue solid curves, label "0") and with returning radiation calculated with four iterations (dashed red curves, label "4").In every panel, we also show the spectrum of the incident radiation from the corona (black solid curves, label "Direct") and the spectrum of the total incident radiation (direct radiation from the corona + returning radiation from the disk, black dashed curves, label "Total").The calculations assume that the spacetime is described by the Kerr metric with spin parameter a * = 0.998, the spectrum of the corona at the emission point is described by a power law with photon index Γ = 1.7 and high-energy cutoff Ecut = 300 keV, and the disk's ionization parameter is ξ = 100 erg cm s −1 .The y-axis shows the quantity νFν in arbitrary units, where ν is the photon frequency and Fν is the observed flux at the frequency ν.The coronal height is h = 2, 5, and 10 rg for the upper, central, and lower panels, respectively.See the text for more details.
Fig.4.-Direct incident X-ray flux from the corona (blue solid curves) and incident X-ray flux generated by the returning radiation as calculated at the 4th order (red dashed curves) for a coronal height h = 2 (left panel), 5 (central panel), and 10 rg (right panel).These results are obtained assuming that the spacetime geometry is described by the Kerr solution with spin parameter a * = 0.998 and that the spectrum of the corona is described by a power law with photon index Γ = 1.7 and high-energy cutoff Ecut = 300 keV.
Fig.5.-Synthetic reflection spectra without returning radiation (blue solid curves, label "0") and with returning radiation calculated with one (red solid curves, label "1"), two (red dashed curves, label "2"), three (green dashed-dotted curves, label "3"), and four (black dotted curves, label "4") iterations.We consider three different values of coronal height (h = 2, 5, 10 rg) and three different values of ionization parameter of the disk (ξ = 100, 1000, and 10000 erg cm s −1 ).The calculations assume that the spacetime is described by the Kerr metric with spin parameter a * = 0.998, the spectrum of the corona at the emission point is described by a power law with photon index Γ = 1.7 and high-energy cutoff Ecut = 300 keV, and the inclination angle of the disk with respect to the line of sight of the observer is ι = 60 • .The y-axis shows the quantity νFν in arbitrary units (ν is the photon frequency and Fν is the observed flux at the frequency ν) and normalized to have the peak of the iron Kα line at 1 in order to facilitate the comparison of the shapes of the spectra.See the text for more details.
Fig. 9.-Data to best-fit model ratios for simulations A1-A3 (h = 2 rg and a * = 0.998) and B1-B3 (h = 5 rg and a * = 0.998) with returning radiation.The red and blue colors represent, respectively, the NICER and NuSTAR data.ξ in units erg cm s −1 .APPENDIX APPENDIX A: CALCULATION OF THE EMISSIVITY PROFILE PRODUCED BY A LAMPPOST CORONA Fig.10.-As in Fig.5but calculating the reflection spectrum assuming that the spectrum of the radiation illuminating the disk is described by a power law with a high-energy cutoff as inDauser et al. (2022).
g , with returning radiation

)
APPENDIX B: CALCULATION OF THE RETURNING RADIATION

Fig
Fig.14.-Radial ionization profiles for the models shown in Fig.13.We have two models of the disk electron density profile: ne = n 0 = 10 15 cm −3 (dashed curves) and the Shakura-Sunyaev profile ne(r) ∝ r 3/2 [1 − (r in /r) 1/2 ] −2 (solid curves).The blue curves are the ionization parameter profiles without returning radiation and the red curves the profiles calculated taking the returning radiation into account.
Summary of the input parameters of the simulations considered in this work with a * = 0.998.In simulations A, the coronal height is h = 2 rg, while in simulations B we have h = 5 rg.The disk ionization parameter is ξ = 100, 1000, and 10000 erg cm s −1 for, respectively, simulations 1, 2, and 3.The values of the other parameters do not change among different simulations.
Best-fit values of simulations A1-A3 and B1-B3 with returning radiation.The reported uncertainties correspond to the 90% confidence level for one relevant parameter (∆χ 2 = 2.71).(B) means that the 90% confidence level reaches the boundary of the parameter or the best-fit value is stuck at the boundary of the parameter.The values of the input parameters are reported in Tab. 1.
. The values of the other parameters do not change among different simulations.

TABLE 4
Best-fit values of simulations C1-C3 and D1-D3 with returning radiation.The reported uncertainties correspond to the 90% confidence level for one relevant parameter (∆χ 2 = 2.71).(B) means that the 90% confidence level reaches the boundary of the parameter or the best-fit value is stuck at the boundary of the parameter.The values of the input parameters are reported in Tab. 3.

TABLE 5
Best-fit values for simulations A1, A3, C1, and C3 with returning radiation, assuming that we have only NICER data.The reported uncertainties correspond to the 90% confidence level for one relevant parameter (∆χ 2 = 2.71).(B) means that the 90% confidence level reaches the boundary of the parameter or the best-fit value is stuck at the boundary of the parameter.The values of the input parameters are reported in Tab. 1 and Tab. 3.