# RTFAST-Spectra: Emulation of X-ray reverberation mapping for active galactic nuclei

Benjamin Ricketts, Daniela Huppenkothen, Matteo Lucchini, Adam Ingram, Guglielmo Mastroserio, Matthew Ho, Benjamin Wandelt

► To cite this version:

Benjamin Ricketts, Daniela Huppenkothen, Matteo Lucchini, Adam Ingram, Guglielmo Mastroserio, et al.. RTFAST-Spectra: Emulation of X-ray reverberation mapping for active galactic nuclei. Monthly Notices of the Royal Astronomical Society, 2025, 538 (2), pp.1096-1117. 10.1093/mnras/staf337. hal-04903194

HAL Id: hal-04903194

https://hal.science/hal-04903194v1

Submitted on 18 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# RTFAST-Spectra: emulation of X-ray reverberation mapping for active galactic nuclei

Benjamin J. Ricketts<sup>1,1,2*</sup> Daniela Huppenkothen<sup>1,2</sup> Matteo Lucchini<sup>1</sup> Adam Ingram<sup>1,3</sup>
Guglielmo Mastroserio,<sup>4</sup> Matthew Ho<sup>1,5</sup> and Benjamin Wandelt<sup>5,6</sup>

Accepted 2025 February 24. Received 2025 February 21; in original form 2024 December 13

## ABSTRACT

Bayesian analysis has begun to be more widely adopted in X-ray spectroscopy, but it has largely been constrained to relatively simple physical models due to limitations in X-ray modelling software and computation time. As a result, Bayesian analysis of numerical models with high physics complexity have remained out of reach. This is a challenge, for example when modelling the X-ray emission of accreting black hole X-ray binaries, where the slow model computations severely limit explorations of parameter space and may bias the inference of astrophysical parameters. Here, we present RTFAST-Spectra: a neural network emulator that acts as a drop in replacement for the spectral portion of the black hole X-ray reverberation model RTDIST. This is the first emulator for the reltrans model suite and the first emulator for a state-of-the-art X-ray reflection model incorporating relativistic effects with 17 physically meaningful model parameters. We use principal component analysis to create a light-weight neural network that is able to preserve correlations between complex atomic lines and simple continuum, enabling consistent modelling of key parameters of scientific interest. We achieve a $O(10^2)$ times speed up over the original model in the most conservative conditions with $O(1$ per cent) precision over all 17 free parameters in the original numerical model, taking full posterior fits from months to hours. We employ Markov Chain Monte Carlo sampling to show how we can better explore the posteriors of model parameters in simulated data and discuss the complexities in interpreting the model when fitting real data.

**Key words:** black hole physics – methods: observational – methods: statistical – stars: black holes – radiative transfer – relativistic processes – software: machine learning.

## 1 INTRODUCTION

Compact objects remain one of the most intriguing phenomena in the known Universe. They are part of some of the most energetic systems in the universe and they present unique challenges in a regime where our understanding of fundamental physics may break down. In particular, black holes are both interesting laboratory for fundamental physics and unique objects in their own right.

One way of probing compact objects is via the study of their accretion flows. At its most extreme, accretion enables us to measure their fundamental physical properties such as the mass and spin of black holes (McClintock et al. 2011) and the equation of state of neutron stars (Özel & Freire 2016). Accreting black holes are also known to show energetic jets, and thus help us to better understand the physical conditions under which relativistic jets form (Nemmen et al. 2012). Characterizing the interaction between these accretion processes with the compact object, and the geometry of the system are crucial for our understanding of the underlying extreme physics.

Black hole accretion models have become significantly more complex as we come to better understand the phenomena that govern their behaviour and we gain access to increasingly varied data such as X-ray polarization (Weisskopf et al. 2016) and high-cadence time series (Jansen et al. 2001; Gendreau, Arzoumanian & Okajima 2012; Uttley & Casella 2014). Many of the foremost challenges in black hole physics arise when we try to combine the rich information in these data sets with sophisticated theoretical models (Remillard & McClintock 2006; Sobolewska & Papadakis 2009).

We are particularly interested in the inner-most regions of the accreting black hole, composed of an accretion disc (of which there are multiple interpretations of the thickness, truncation and magnetohydrodynamics) that emits thermally (Shakura & Sunyaev 1973), a region of energetic electrons around the black hole (known as the corona) that emits in a non-thermal power law (Thorne & Price 1975; Sunyaev & Truemper 1979), and, under certain conditions, a jet which also emits non-thermally (Blandford & Königl 1979; Levinson & Blandford 1996; Romero et al. 2017). A portion of the photons emitted by the corona illuminates the accretion disc, where it is reprocessed and generates what is called a reflection spectrum (Fabian et al. 1989). Modelling these phenomena requires

© 2025 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1"><sup>1</sup>Anton Pannekoek Institute, University of Amsterdam, Science Park 904, Amsterdam, NL-1098 XH, the Netherlands</footnote>

<footnote n="2"><sup>2</sup>SRON, Niel Bohrweg 4, Leiden, NL-2033 CA, the Netherlands</footnote>

<footnote n="3"><sup>3</sup>School of Mathematics, Statistics and Physics, Newcastle University, Herschel Building, Newcastle Upon Tyne, NE1 7RU, UK</footnote>

<footnote n="4"><sup>4</sup>Dipartimento di Fisica, University of Milan, Via Festa del Perdono, 7, I-20122 Milano MI, Italy</footnote>

<footnote n="5"><sup>5</sup>Institut d'Astrophysique de Paris, Sorbonne Université/CNRS, 98 bis Bd Arago, F-75014, Paris, France</footnote>

<footnote n="6"><sup>6</sup>Department of Physics and Astronomy, Johns Hopkins University, 3400 N Charles Street, Baltimore, MD 21218, USA</footnote>

<footnote>*E-mail: b.ricketts@sron.nl</footnote>
the computation of photon trajectories in a regime where relativistic effects dominate (Cunningham & Bardeen 1973; Dauser et al. 2013), which in turn strongly depend on the black hole spin. The different parts of the inner regions all interact with each other to contribute to a complex total spectrum that also varies with time. By observing the correlated variability of the emitted radiation as a function of time and energy, we are able to measure the light travel time difference of the reflected photons compared to the non-thermal coronal photons. Proper modelling of this emission allows us to measure the size of the inner region, thus the mass of the black hole.

A fundamental challenge in modelling accreting sources is the high degree of parameter degeneracy: different parameter combinations may fit a data set similarly well, and therefore choices such as starting parameters, random seeds, and optimization algorithm may determine the scientific conclusion. This is undesirable, as we require robust inferences of models and their parameters in order to advance our scientific understanding. One potential solution is to explore a wider range of parameter space using sampling algorithms like Markov Chain Monte Carlo (MCMC). This can be efficiently done in a Bayesian context, and allows us to report multidimensional probability distributions of the parameters given the data, rather than single-point best-fitting values and errors derived using approximations. Posterior distributions are particularly useful when the parameter space is complex and multiple combinations of parameters may generate models that fit the data similarly well (also known as multimodal probability distributions).

Some studies in X-ray spectral modelling have adopted Bayesian inference (Buchner et al. 2014; Dias et al. 2024) but have to contend with significant computational burdens. Models that calculate physical parameters self-consistently can be computationally expensive, and the high number of model evaluations required by sampling algorithms can make inferring posteriors for real data slow. As these models are also often high dimensional, the number of evaluations required to map the posterior space dramatically increases, creating an exponential increase in inference time. This leads to a trade-off between the complexity (and realism) of models and the computational feasibility of comparing these models to data.

The numerical model RTDIST (Ingram et al. 2022, henceforth I22) was developed to model the X-ray spectra and timing products from the corona and reflection spectrum of active galactic nuclei (AGNs) and X-ray binary black holes. RTDIST incorporates X-ray reflection and reverberation into one model that self-consistently calculates ionization within the disc according to the geometry of the system. We further expand on the features of RTDIST and assumptions in Section 2. RTDIST currently takes approximately 0.5 s to evaluate the time averaged spectrum once on a single core of a typical laptop. Timing products such as the lags-energy spectrum take slightly longer. When using MCMC sampling to obtain Bayesian posteriors, this means that modelling a single data set can take weeks to months of computation time.

An alternative to directly computing the numerical model often applied in X-ray astronomy are interpolation table models. Many models are computed on a grid of parameters which will then be collated into a table. The table is used in linear interpolation to obtain model estimates for parameters not present in the table itself (Wilms, Allen & McCray 2000; García et al. 2013). This approach is much faster computationally, but is no longer feasible for high-dimensional models like RTDIST as the required number of stored, pre-computed models expand to the power of the number of model parameters, leading to unfeasibly large interpolation tables to be read into memory.

X-ray astronomy is not the first field of research, both within and outside of astronomy, to have faced this problem. Interpolation models can be (relatively) easily replaced by creating an emulator composed of flexible, parametric or non-parametric functions that replicate the relationship of input parameters and model outputs. Emulators have been used in a wide variety of scenarios within astronomy including cosmological simulations (Heitmann et al. 2009; Bird et al. 2019; Kaushal et al. 2022), cosmological matter power spectra (Heitmann et al. 2013; Spurio Mancini et al. 2022), optical and UV energy spectra of galaxies and Type Ia supernovae (Alsing et al. 2020; Kerzendorf et al. 2021) and X-ray spectral studies of accretion disc winds in AGNs (Matzeu et al. 2022). Emulators have been built with a wide range of techniques, including Gaussian processes (Bastos & O'hagan 2009), artificial neural networks (ANNs; Himes et al. 2022), and kernel ridge regression (Vicent et al. 2018). In this paper, we utilize an ANN as an emulator.

ANNs have been shown to obey the universal approximation theorem (Baker & Patil 1998): according to which an arbitrary continuous function of a bounded set can be, within a given accuracy, approximated by an appropriate ANN. As such, they can straightforwardly model even very complex relationships between parameters and model outputs while being fast to compute and extremely memory-efficient in comparison with interpolation methods.

An ANN emulator is trained on a large training data set consisting of pairs of parameters and the corresponding outputs of the model to be emulated. Posterior sampling of one spectrum may take several million evaluations of the model. Using the same number of model evaluations to generate training data saves time in the long run, as the ANN model, once trained, is orders of magnitude more efficient than the original model. The ANN can also be easily loaded on to a GPU, and model evaluations can be parallelized to achieve even faster evaluation times.

A reduced model evaluation time allows for faster optimization, sampling and a more comprehensive exploration of parameter space. The latter is especially important to ensure robustness of the statistical inference. Reporting posteriors enables us to report degeneracies and correlations in parameter space. The differentiable nature of a neural network emulator also introduces the possibility of using modern, more efficient methods of optimization and sampling, developed outside of the field (Giolami & Calderhead 2011; Cranmer, Brehmer & Louppe 2020; Buchner 2023).

This work presents an emulator called RTFAST-spectra (subsequently RTFAST) as a drop-in replacement for the spectral portion of the RTDIST model (expanded upon in Section 2) which is considerably faster than the original model. We outline our methodology in Section 3, including dimensionality reduction via principal component analysis (PCA), the neural network architecture and a new, custom loss function to efficiently train the network on the PCA representation. In Section 4, we present a series of tests of RTFAST to demonstrate its reliability. In Section 5, we summarize our model and discuss its possible use cases and future work.

## 2 RTDIST

RTDIST is part of the RELTRANS package (Ingram et al. 2019) and that the new version of RELTRANS that includes RTDIST will soon be publicly available. A comprehensive introduction to RTDIST and its assumptions can be found in I22. RTDIST incorporates a considerable amount of theoretical development in X-ray reverberation and black hole modelling over the last decade. As such, we provide here a concise overview of RTDIST.
<table><caption>**Table 1.** List of parameters in RTDIST, the type of parameter space used in training and their units. All units indicated with an asterisk (*) are dimensionless. Note that $A_{norm}$ is a range calculated from a possible range of coronal fluxes specified in Section 3.1.2.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Spacing</th><th>Range</th><th>Units</th></tr></thead><tbody><tr><td>h</td><td>lampost corona height</td><td>Log</td><td>1.5–100</td><td>R<sub>g</sub></td></tr><tr><td>a</td><td>Spin</td><td>Linear</td><td>0–0.998</td><td>*</td></tr><tr><td>Inclination</td><td>Inclination</td><td>Log</td><td>1–80</td><td>deg</td></tr><tr><td>R<sub>inner</sub></td><td>Inner disc radius</td><td>Log</td><td>1–400</td><td>ISCOs</td></tr><tr><td>R<sub>outer</sub></td><td>Radius of the outer edge of the disc</td><td>Log</td><td>400–10<sup>5</sup></td><td>R<sub>g</sub></td></tr><tr><td>z</td><td>Redshift</td><td>Linear</td><td>0–0.1</td><td>*</td></tr><tr><td>Γ</td><td>Photon index</td><td>Linear</td><td>1.4–3.4</td><td>*</td></tr><tr><td>D<sub>kpc</sub></td><td>Angular diameter distance</td><td>Log</td><td>10<sup>2</sup>–10<sup>6</sup></td><td>kpc</td></tr><tr><td>A<sub>fe</sub></td><td>Iron abundance of disc in solar units</td><td>Log</td><td>0.5–10</td><td>*</td></tr><tr><td>logNe</td><td>Electron density of disc</td><td>Linear</td><td>15–20</td><td>N<sub>e</sub>: cm<sup>-3</sup></td></tr><tr><td>kTe</td><td>Observed electron temperature</td><td>Log</td><td>5–500</td><td>keV</td></tr><tr><td>1/B</td><td>Boost</td><td>Log</td><td>10<sup>-2</sup>–10<sup>-5</sup></td><td>*</td></tr><tr><td>Mass</td><td>Mass of BH</td><td>Log</td><td>10<sup>4</sup>–10<sup>11</sup></td><td>M<sub>☉</sub></td></tr><tr><td>Scale height</td><td>Aspect ratio of conical disc</td><td>Linear</td><td>0–0.176</td><td>*</td></tr><tr><td>b1</td><td>Angular emissivity parameter</td><td>Linear</td><td>0–2</td><td>*</td></tr><tr><td>b2</td><td>Angular emissivity parameter</td><td>Linear</td><td>-4–4</td><td>*</td></tr><tr><td>A<sub>norm</sub></td><td>Normalization of the total model</td><td>Log</td><td>10<sup>-4</sup>–4</td><td>*</td></tr></tbody></table>

RTDIST assumes a lamppost point source corona emitting a Comptonization spectrum located above a black hole, which has a thermally emitting disc with a fixed aspect ratio of height over radius. The spectrum emitted by the corona is calculated using the thermal comptonization model nthcomp (Zdziarski, Johnson & Magdziarz 1996). The power law index of this spectrum is assumed to vary with time (a feature inherited from other RELTRANS models), motivated by empirical constraints from observations (Kotov, Churazov & Gilfanov 2001). The spectral pivoting is implemented using a first-order Taylor expansion to linearize the model. This emission from the corona travels not only to the observer, but also illuminates the disc. When emission from the corona hits the disc, it is absorbed and reprocessed by material in the disc. This results in re-emission from the disc that is referred to as the reflection spectrum. This reflection spectrum has three key components that are apparent in the observed spectrum: the soft excess below 1 keV, the iron line complex around 6.4 keV and the Compton hump (Ross & Fabian 2005; García et al. 2013).

This reflection spectrum is also affected by relativistic effects due to the corona and reflecting material being situated near the black hole (Fabian 1989). These are modelled via ray-tracing that takes into consideration boosting from directed emission from the corona (rather than simply assuming isotropic emission), as well as blue- and red-shifting effects. This results in relativistic smearing of the entire reflection spectrum, which is most apparent in the narrow-line Fe K complex.

Distinct from other models (with the exception of the KYN package; Dovčiak, Karas & Yaqoob 2004), the ionization of the disc is calculated self-consistently from the distance to the source. The ionization of the disc affects the shape of the resultant reflection spectrum which, combined with the observed flux of the spectrum, allows for the constraining of distance to the source. In conjunction with an observed red-shift, one can theoretically constrain the Hubble constant from the spectra alone with a mass prior from observations in other wavelengths. Including the timing information into the model allows us to measure the distance even without a prior on the mass, since the mass becomes a parameter measured by the model.

Further extending this line of reasoning, considering multiple AGN simultaneously allows for the marginalization over peculiar velocity of individual sources. This allows for better constraints of

distance (and thus ionization) by sharing this common parameter within fits of data. This enables another independent method of constraining the Hubble constant with X-ray spectral information. This method also serves as a way to directly verify the geometric assumptions made by RTDIST if RTDIST agrees with cosmological measurements of the Hubble Constant. Theoretically speaking, this method can already be performed with observations made of nearby X-ray bright AGN with XMM-Newton (Mitchell et al. in preparation). Practically, the complex and long fitting process of an individual AGN makes this difficult to do while retaining the self-consistent calculations required to do this at scale. This particular scientific question directly motivated the development of RTFAST, and RTFAST has several design decisions related to this motivation, such as the modelled energy grid and restricted model parameters.

All parameters in RTDIST are listed in Table 1 as well as their ranges. RTDIST is a general model that can be applied to both AGN and X-ray binaries; however, we only train RTFAST for AGN so all parameter ranges are restricted to that use case.

## 3 RTFAST

In this section, we describe the framework for the training of RTFAST. We outline the choices we made when designing RTFAST as well as how we arrived at them. This section can be viewed as a way to approach similar modelling problems to those that we encountered.

### 3.1 Data generation

RTDIST can model a wide range of scenarios, ranging from AGN with low-temperature discs to super-Eddington accreting X-ray binaries. To model every scenario, no matter how unrealistic or unlikely, requires an enormous number of training examples. Thus, we constrain the modelling space to a subset of parameters physically meaningful in modelling AGN. We use sampling techniques, PCA and physical parameter constraints to achieve a manageable parameter space for neural network training. In the following, we will describe how we generate pairs of parameter sets $\theta_i$ and the associated model spectrum generated by RTDIST, $y_i$, $i \in N$. $\theta$ is a vector of parameters described in Table 1.
### 3.1.1 Parameter space

Many of the parameters of RTDIST span multiple orders of magnitude. We use the logarithm of those parameters in training data generation to ensure that all orders of magnitude are equally sampled and represented in the training data. In the case of both very small and very large numbers, this also ensures numerical stability within the neural network. In Table 1, we list each parameter, the type of parameter space used for that parameter, and its units. In actual usage of RTFAST, users input parameters in the conventional physical units of the original model.

Grid sampling is not very effective in higher dimensions due to the curse of dimensionality (Bellman 1966), and requires large amounts of data to effectively sample the entire parameter space. We used the implementation of Latin Hypercube sampling (LHS; McKay, Beckman & Conover 2000) from the SCIPY python package (Virtanen et al. 2020) to effectively sample the entire parameter space while preventing over/under sampling of areas of parameter space that can occur when uniformly sampling high-dimensional spaces. Latin hypercube sampling achieves this by considering the potential parameter space as components divided into strata of equal marginal probability. Each of these strata are then sampled once and are then randomly matched together. This ensures the entire space is covered effectively and all areas of the overall distribution are explored. This sampling results in all instances of parameter sets $\theta_i$.

### 3.1.2 Physically motivated parameter constraints

Simulating spectra with completely random parameter combinations within the very broad limits set in Table 1 resulted in the simulation of non-physical spectra, such as spectra with emission lines that exceeded the disc continuum emission by 10 orders of magnitude. This introduced additional complexity for the emulator to learn without any scientific gain. We imposed a series of physically motivated parameter constraints that limited our training space to an area that is physically realistic and useful to observers. We imposed six physical restrictions, which are the following:

(i) The redshift and angular diameter distance to the object must be consistent with a Hubble constant between 30 and $300 \text{ km s}^{-1} \text{ Mpc}^{-1}$. Current measurements indicate the Hubble constant to be between 67 and $73 \text{ km s}^{-1} \text{ Mpc}^{-1}$ at their most extreme (Freedman 2021). We allow for a wider range of possible inferred values to allow for flexibility in inference of the implied ionization of the disc in the final model.

(ii) The height of the corona, $h$, must be greater than the horizon radius as otherwise photons emitted from the corona do not escape the black hole.

(iii) The coronal flux of the object must be within $10^{-8}$ and $10^{-12}$ erg s$^{-1}$ cm$^{-2}$ at earth, the range currently observed from AGN.

(iv) The coronal flux of the object must be consistent with an accretion luminosity between $10^{-4}$ and $2.5 \ L_{edd}$ (Eddington luminosity). This aligns with X-ray fluxes of local AGN.

(v) All three of the following may not be true simultaneously:

(a) The photon index of the corona is greater than 2.75.

(b) The disc has an iron abundance of 4 solar units or greater.

(c) The electron density of the disc is $10^{19}$ or higher.

These particular restrictions incorporate knowledge from observations of X-ray binaries where these properties are at their most extreme (García et al. 2018; Dong et al. 2020; Ludlam et al. 2020). Disc densities have been shown to be anticorrelated with

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Corner plot of 4 of the training set parameters and their distributions. Possible training data to be generated is constrained by the physical constraints outlined in Section 3.1.2. The full distributions can be found in Fig. A1.</figcaption></figure>

black hole mass and higher disc densities have also been shown to decrease the inferred iron abundance of discs of AGN (Jiang et al. 2019), meaning these restrictions serve as a suitable upper limit for AGN.

(vi) The spin is forced to be pro-grade since no observational evidences of retrograde spin have been found (Draghis et al. 2024).

These restrictions are applied to all parameter sets $\theta_i$ generated in Section 3.1.1. This functionally necessitated the removal of instances of $\theta_i$ that were invalid, resampling the parameter space, and checking if these new parameter sets were valid. This was repeated until the desired number of parameter sets $\theta_i$ was reached.

These physical constraints are also visualized in Fig. 1. The full distributions can be found in Appendix A. By imposing these physical restrictions, we reduced the possible parameter space to 1 per cent of its original size (Table 1) and eliminated non-physical results. Fig. 1 shows clear correlations between parameters: for instance, $A_{norm}$ and mass show a dependence such that increasing $A_{norm}$ only allows for increasingly smaller range of masses: this is due to $A_{norm}$ being highly correlated with accretion luminosity and coronal flux. Lower values of $A_{norm}$ mean that AGN with smaller masses do not produce enough flux to be meaningfully observed at earth at the distances specified by our physical constraints. This effect is also seen in the correlations between $A_{norm}$ and redshift (and subsequently distance).

We note that, due to these restrictions, moving outside of these imposed physical bounds is likely to give incorrect results. Any usage of RTFAST should impose these physical constraints. We demonstrate this restriction in Section 4.5. If observations implied scenarios outside of these parameter constraints, it would be relatively trivial to loosen these restrictions by re-training the neural network on new data outside the original constraints.
### 3.1.3 Pre-processing of spectra

Neural networks training performs best when the values the network is supposed to reproduce are well-behaved and narrowly distributed. Therefore, we applied several reversible pre-processing steps. The neural network learns a relationship between the parameter and these pre-processed quantities, which can then be easily and efficiently expanded into a physically meaningful spectrum:

(i) All simulated fluxes below $10^{-11}$ photons cm$^{-2}$ s$^{-1}$ keV$^{-1}$ in the spectrum were fixed to $10^{-11}$ photons cm$^{-2}$ s$^{-1}$ keV$^{-1}$. This prevents the emulation of fluxes that we realistically will not observe. This threshold corresponds to an average of one photon per energy bin detected during a 1 Megasecond observation with the Athena X-IFU telescope (Kammoun et al. 2022). This is realistically extremely conservative and no current or planned telescope at the time of writing will detect photons below this limit.

(ii) Neural networks are known to struggle generating outputs that span multiple orders of magnitude. To mitigate this issue, we took the logarithm of the flux $y_i$ in each energy bin $E_i$ to obtain a new log-scaled quantity, $x_i = \log(y_i)$.

(iii) Each $x_i$ is then scaled such that the mean of $x_i$ for each bin $E_i$, taken over all N training examples, is 0 and the standard deviation is 1. This scaling enhances features such as changes in the iron absorption and reflection lines that deviate significantly from the average spectrum, both of which are particularly important in our science case.

(iv) We performed a PCA decomposition with 200 PCA components on the spectra. We go into more detail about this technique in Section 3.1.4.

(v) The 200 PCA components are scaled to a mean of 0 and standard deviation of 1 over the data set for identical reasons as in step 3, with the additional motivation that neural networks behave more stably when the target output is constrained around 0 with magnitudes of approximately 1.

This pre-processing was key in achieving and maintaining high precision over the entire spectrum by reducing the complexity of the model space and reducing the size of the neural network required to model the spectral behaviour.

### 3.1.4 PCA decomposition

The large number of energy bins present in X-ray spectra requires large neural networks by design. In turn, these networks struggle to capture small sharp features such as the iron line, as most spectral bins are dominated by the continuum emission, which thus dominated training performance. To reduce the number of outputs the neural network has to predict, we performed PCA decomposition to capture the relevant structure in the spectra in a lower dimensional space.<fnref n="1" /> This reduced the spectrum from 2017 energy bins (the XMM–Newton EPIC-PN energy grid above 0.1 keV) to K principal components. We found that by decomposing the spectra, correlations between energy bins were effectively preserved. This also led to considerably smaller neural networks as it reduced the problem to learning a relationship between 17 input parameters and K output features, which thus significantly reduced training time and complexity.

In order to determine K, we need to balance the trade-off between the reduction of the dimensionality of the output space to be predicted

by the neural network with the loss of accuracy implied by the use of a finite set of PCA components. To do this, we iteratively increased the number of PCA components K, and computed the average error in the reconstructed spectra across a set of test spectra that were not used in computing the PCA components. We stopped adding components when this average error fell below a defined threshold. We chose to set this threshold at 0.1 per cent over the entire parameter space, and defined the error as

$$x_i = (x_{i,NN} - x_{i,RTDIST})/x_{i,RTDIST}, \quad (1)$$

where $x_{i,NN}$ is the spectrum in bin i reconstructed from the PCA components, transformed as described in steps (i)–(iii) in Section 3.1.3, $x_{i,RTDIST}$ and denotes the corresponding transformed spectrum generated with RTDIST. This equation describes the fraction error $\chi_i$, which we will also use in later on to assess the reconstruction performance of the neural network combined with PCA.

At a threshold of 0.1 per cent reconstruction error, we effectively require that the PCA reconstruction must be an order of magnitude better than the performance we require from the neural network.

We found that $K = 200$ PCA components generally fulfilled our requirement, as shown in Fig. 2.

Fig. 3 shows the first 5 PCA components plotted as a function of energy. The original spectrum is recreated from a linear combination of these components, weighted by the PCA projection of the original data. It is these PCA projections that the neural network learns to emulate. Because PCA is a reversible operation, reconstructing the original spectra from the projections is trivial. The first component in Fig. 3 is largely flat and thus encapsulates the effect of a normalizing constant.<fnref n="2" /> The second PCA component is dominated by a power-law-like feature.

Notably, in Fig. 3, components 3 and onwards exhibit much more interesting features than the simple normalization and power-law component that the first two components can largely explain. Later components can be directly related to the soft excess below 1 keV and the iron line complex at 6.4 keV. These components show how the PCA decomposition directly captures physical features of interest. The remaining components that are not shown in this plot almost exclusively relate to better describing the soft excess and iron line.

Fig. 4 shows the cumulative explained variance ratio subtracted from 1. The explained variance ratio is a statistical measure of how much of the variation of a data set is captured by each principal component. If the explained variance ratio sums to 1, all variance in the data set is explained, and any input from the data set can be perfectly reconstructed. The first two components capture the majority of the explained variance, accounting for 99.7 per cent of the variance in the model. The following components capture increasingly smaller amounts of the overall explained variance despite being very important to recreating features in the spectrum such as emission lines.

### 3.1.5 Data set split

We generated a total of $1.1 \times 10^7$ spectra using the SHERPA python package (Freeman, Doe & Siemiginowska 2001) to interface with the original RTDIST model on the energy grid of the XMM–Newton PN CCDs from 0.1 to 20 keV. This data set amounts to a total of approximately 700 GB. We split the data into a training

<footnote n="2">²Ordinarily, we would divide out a normalizing constant from the spectrum and add it as a separate output to be modelled by the neural network. However, for RTDIST, the normalization is not independent from the shape of the spectrum, thus this route is not available to us.</footnote>

<footnote n="1">¹Alsing et al. (2020) is a historical example of using PCA decomposition to enable better emulation of spectra.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** A heat-map of the percentage residual difference between the PCA reconstruction and RTDIST sorted by the model parameter spin (a). Yellow to orange indicates 1–10 per cent error, green indicates 0.1–1 per cent error and blue indicates 0.01 per cent error and below. Black indicates an error above 10 per cent. White indicates that both the RTFAST and RTDIST output are below the detectable limit outlined in Section 3.1.3.</figcaption></figure>

set and validation set with a 90 per cent \ 10 per cent split. We found a reliable increase in model accuracy with increasing of data set size, but did not increase the training data set size further as we reached the $O(1$ per cent) accuracy requirement while also approaching storage limits on the computer we trained on. We also separately generated a test set of 1000 spectra for performance review using the same methodology as the training and validation set.

## 3.2 Neural network

A neural network consists of a large ensemble of neurons, arranged into several layers and nodes. The output of the $m^{th}$ node in the $n^{th}$

layer is $a_m^n = \sigma(\mathbf{a}^{n-1} \odot \mathbf{W}_m + \mathbf{b}_m)$, where $\mathbf{a}^{n-1}$ is the output of the $(n-1)$th layer, $\mathbf{W}_m$ and $\mathbf{b}_m$ are trainable arrays of hyper-parameters, and $\sigma()$ is the activation function.

Training of the neural network involves optimizing the parameters to minimize a loss function. This optimization is achieved using stochastic gradient descent, which allows for the stochastic approximation of gradient descent optimisation.

We explored a variety of strategies to train RTFAST including loss functions, different optimizers, learning rate schedulers, and architectures. We ultimately achieved best performance with a somewhat conventional approach with the exception of a custom loss function. We report how we reached this approach and our best-performing strategy in the following section.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The first 5 PCA components plotted as a function of energy. Each component is in arbitrary units and is scaled appropriately for their actual contribution in the final spectrum.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Cumulative summing of the explained variance by PCA component subtracted from 1.</figcaption></figure>

### 3.2.1 Architecture

The neural network in this paper is built using the PYTORCH package (Paszke et al. 2019). We utilized a feed-forward neural network. The architecture was deliberately kept methodologically very simple as previous work (Kwan et al. 2015; Rivera et al. 2015) has found high efficacy with feed-forward neural networks compared to those with more sophisticated architectures. Smaller networks also reduce overhead in computation and complexity in training.

After preliminary tests on a network that was hand-tuned, we initially performed grid sweeps of different architectural structures during training: changing the number of layers (4, 6, 8, 10, and 12), nodes (256, 512, and 1024) and type of activation function [Rectified error Linear Unit (ReLU) and Gaussian Error Linear Unit (GELU)] used in the network to find the best-performing model while keeping learning hyper-parameters fixed. ReLU activation functions can be defined as the non-negative part of its argument, also known as the ramp function. This is defined as

$$ReLU(x) = x^{+} = \max(0, x) = \frac{x + |x|}{2}. \qquad (2)$$

GELU activation functions can be considered a smoother form of the ReLU that are widespread in neural networks in the literature today. The GELU activation function is $x\Phi(x)$, where $\Phi(x)$ is the standard Gaussian cumulative distribution function. GELU weights inputs by their percentile rather than gating the inputs by their sign like in ReLU. The GELU($x$) is defined as

$$GELU(x) = x\Phi(x) = x \cdot \frac{1}{2}[1 + \erf(x/\sqrt{2})]. \qquad (3)$$

We found that a network of 8 fully connected linear layers with 256 nodes in each layer with GELU activation functions (Hendrycks & Gimpel 2016) between each hidden layer to perform the best for our use case. We show a visual representation of the best performing network in Fig. 5.

### 3.2.2 Hyperparameter grid sweeps

We then performed grid sweeps of hyperparameters, including batch size (256, 512, 1024), learning rate ($10^{-4}$, $5 \times 10^{-3}$, $10^{-3}$), optimizers [SGD (Sutskever et al. 2013), Adam (Kingma 2014), and AdamW (Loshchilov 2017)] and learning rate schedulers [no scheduling, Cosine annealing (Loshchilov & Hutter 2016), and ExponentialLR]. Grid sweeps in machine learning are the systematic training of a neural network with different neural networks and sets of pre-specified hyper-parameters (where hyperparameter is a parameter that defines the details of the learning process of a neural network) to then compare final performance and select the best performing combination. Hyperparameters include (but not exclusively): learning rates, the type of optimizer used during training, and batch size. This approach often results in faster development than trial and error development of a neural network at the cost of computational power.

We found that a learning rate of $10^{-4}$, a batch size of 1024, an Adam optimizer with otherwise default hyper-parameters and no learning rate scheduler yielded the best results. In general, we found larger

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Network architecture of RTFAST. RTFAST is composed of 12 hidden layers, each with 256 nodes, an input layer of 20 nodes and an output layer of 40 nodes. Between each layer, with the exception of the output layer, is a GELU activation function.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Training and validation loss as a function of training epoch. Due to the large training data size, initial validation loss is lower than training loss.</figcaption></figure>

batch sizes and lower learning rates to perform better. The addition of scheduling was negative for the Adam and AdamW optimizers.

## 3.3 Training

During training, we used a custom loss function that took into account the relative importance of each PCA component. We used a weighting factor to push the network to prioritize learning components that had the largest effect on the final reproduced spectrum. This also made the neural network training process more stable when compared to the non-weighted scenario.

We achieved this by weighting each component by a factor $w_i$ which was calculated as the following:

$$w_i = \log_{10}(\text{var}(\text{PC}_i)/\text{var}(\text{PC}_n)) + 1, \quad (4)$$

where var(PCᵢ) and var(PCₙ) denote the explained variance of the principal component to be weighted and the least informative principal component, respectively. The additional constant of 1 ensures that no component is totally ignored during training. The final loss function takes the form:

$$L = \sum_{i=0}^{n} w_i \cdot (PC_i - \hat{PC}_i)^2, \quad (5)$$

where PCᵢ is the scaled principal component and $\hat{PC}_i$ is the neural network output. This can be considered a weighted mean square error loss.

We trained all neural networks described in Section 3.2.1 for a maximum of 2000 epochs (where an epoch is one complete pass of the training set through the algorithm) with an early stopping criterion of requiring at least 1 per cent improvement over 100 epochs. In Fig. 6, we show the training and validation loss as a function of training epoch for the best performing network found by performing grid sweeps as described in Section 3.2.1. While both training and validation loss continued to improve over the last 100 epochs, we found the two were improving increasingly slowly, indicating that the network had learnt most of the information in the data set.

## 3.4 Ensemble

During training of the network, stochastic gradient descent optimization moves the network towards the global minimum of the loss landscape. While theoretically possible to reach the global

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparison of residual probability density functions on the test set between using a single neural network and the ensemble. The dashed red vertical lines indicate 1 per cent error and the dashed blue vertical line indicated 0 per cent error. The ensemble centres the emulator outputs directly around the truth, largely eliminating systematic over- or underestimation.</figcaption></figure>

minimum, this often takes a prohibitively long time and training is stopped before this global minimum is reached. A computationally less expensive strategy to improve performance of a neural network prediction is through the use of an 'ensemble' of neural networks (Ganaie et al. 2022). Each of the neural networks in an ensemble is trained separately either on the full data set or on a subset of it. Each of these neural networks will be trained to try and predict the correct solution but will give a slightly different solution to a given set of parameters. Averaging the outputs of each of the neural networks together has been shown to yield improvements in accuracy in regression problems like emulation (Mendes-Moreira et al. 2012; Lakshminarayanan, Pritzel & Blundell 2017).

We trained seven instances of the best performing neural network architecture from Section 3.2.1, following the same training pattern as in Section 3.3. The final validation losses of each of these networks were comparable to one another, indicating similar levels of convergence as for the initial network (Fig. 6).

We compared the improvement yielded from considering the ensemble of neural networks to a single instance of the neural network. To evaluate performance, we use the same fractional residual differences between the true spectrum and the prediction from RTFAST, parametrized in equation (1), as for the PCA reconstruction. In Section 3.1.4, we only compared how well PCA reconstructed the original spectra in order to understand the number of components required. Here, we compare the original spectra generated by RTDIST to emulated spectra, generated by letting the ANN predict PCA components from the same input parameters, which are then transformed as described in Section 3.1.3 into spectra comparable to RTDIST's output. For the test set, we obtain a distribution of residual differences per energy bin. Ideally, this distribution would be a delta function centred at 0 indicating that RTFAST perfectly reproduces the RTDIST output. In reality, this is not the case. In Fig. 7, the left-hand panel shows the distribution of residuals for a single neural network while the right-hand panel shows the distribution of residuals for the ensemble. What is apparent is that the single neural network overestimates the true output on average by about 3 per cent and exhibits a skewed distribution while the distribution of the ensemble instead is centred much closer to 0, implying that the ensemble on average estimates the correct output to better precision. The distribution of errors of the ensemble is tighter, with approximately 56 per cent more of the residuals sitting within the bounds of -1 to 1 per cent.

While evaluating seven neural networks does increase computation time, use of vectorization means that the total evaluation time of the
ensemble is only approximately 2.7 times longer than evaluating a
single neural network and the mitigation of error and systematic bias
is worth the increase in computation time. We adopt this ensemble
averaging method for the final version of RTFAST.

3.5 Computation speed

A key goal of RTFAST was to speed up model evaluation. RTFAST takes $1.62 \times 10^{-3}$ seconds for a single set of parameters on a Intel Xeon Gold 6342 2.80GHz processor. This is a 200 times speed up when directly compared to the 0.32 seconds it takes to call RTDIST on the same processor. This is the most conservative estimate of the speed up as the neural network also allows for vector computation, enabling parallel computation of models at much higher speed than in serial. For instance, the computation of 1000 different models on RTFAST takes $3.56 \times 10^{-2}$ seconds, while it took RTDIST 314 s called sequentially. This represents a speed up of 8000 times faster than RTDIST and 39 times faster than calling RTFAST sequentially.

The time it takes to perform a single evaluation of RTFAST in even the most conservative conditions is similar or smaller to the amount of time it takes to perform a single response convolution of XMM-Newton EPIC responses with the model output. This means that parameter inference is no longer computationally limited by the physical model but instead by the convolution of the instrument response and calculation of the likelihood.

4 ROBUSTNESS

In this section, we perform a suite of tests to assess the robustness
of the emulator. We consider the emulator to be robust if it reliably
reproduces the true output of RTDIST to within $O(3$ per cent) and
without any clear systematic biases. As we envision using RTFAST
with XMM-Newton data, we use the systematic error of the XMM-
Newton instrumentation as a threshold of how far off the truth
RTFAST can be. The most recent cross-calibrations estimate the
MOS and PN cameras to have an error in relative effective area to be
within ±3 per cent and ±2 per cent, respectively (Guainazzi et al.
2013). We subsequently adopt a criterion of 3 per cent error. We
explored multiple ways to show that RTFAST fulfils these criteria as
well as performing fits of simulated observations from RTDIST.

4.1 Direct comparison

In Fig. 8, we directly compare RTFAST to RTDIST for one example
from the test set. We observe that there are differences between the
models as a function of energy but note that all of these differences
are within the desired accuracy above 1keV. At the soft excess at low
energies, the difference between RTFAST and RTDIST becomes
larger and more jagged. While this is somewhat of a theme throughout
the overall parameter space for the emulator, the soft excess is
largely obscured by hydrogen column absorption and detector effects
in actual observations. These discrepancies at higher resolutions
become largely washed out under realistic observation conditions.
We discuss this further in a fit to simulated data in Section 4.4.

## 4.2 Expected model behaviour

In addition to a direct comparison from a randomly chosen parameter
set, we also tested that we observe the expected changes to the
spectrum from varying just one parameter. We vary Γ due to its easily
observable change in the continuum spectrum as well as absorption
and reflection lines becoming more apparent as Γ is increased. In

<figure><img src="image_9.png" /><figcaption>**Figure 8.** A direct comparison of RTFAST and RTDIST. The top panel shows the photon flux spectrum from RTFAST and RTDIST as a function of energy, with both lines on top of one another. The middle panel shows the same in $\text{keV}^2 \text{cm}^{-1} \text{s}^{-1} \text{keV}^{-1}$ units. The inset panel shows a zoom-in of the relativistically smeared iron line between 5.5 and 7keV. The bottom panel shows the fractional residual difference between RTFAST and RTDIST as a function of energy. The 3 per cent error region is shaded.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Model spectra generated by RTFAST are plotted on top of one another, coloured by increasing brightness with increasing Γ with all other parameters fixed to that of the parameters in Table 3.</figcaption></figure>

Fig. 9, we plot the resulting spectra from RTFAST [multiplied by TBABS to model the effects of a line-of-sight absorption column, using the abundances of Wilms et al. (2000)] as a function of Γ. We observe the expected behaviour of steepening of the spectrum as Γ is increased. We also observe the appearance of the Fe-K complex within the 6–8 keV range with the increasing steepness of the slope – a key component of reflection modelling within the model and a complex feature to model. Due to the relativistic smearing present in the original model, these are not sharp lines.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** A heat-map of the fractional residual difference between RTFAST and RTDIST (see equation 1). Yellow to orange indicates 1–10 per cent error, green indicates 0.01–3 per cent error and below. Black indicates an error above per cent. White indicates that the model output is below the detectable limit outlined in Section 3.1.3. An error of 3 per cent or less is comparable to instrumental error of current X-ray telescopes with the exception of XRISM. At the time of writing, RTDIST is not capable of XRISM resolution, so we did not attempt to achieve such resolution in RTFAST.</figcaption></figure>

Accurately reproducing key emission lines such as the Fe-K complex is of key importance to a robust data analysis with the emulator. These lines contain a wealth of physically relevant information that enable tight constraints on model parameters when applied to data. It is not an exaggeration to say that the failure to model these key areas of the spectrum would result in a practical inability to use RTFAST in modelling real observations.

## 4.3 Inaccuracy as a function of parameter space and energy

In Fig. 10, we show the fractional difference between RTFAST and RTDIST for the entire test set as a function of black hole spin (a).

We find that RTFAST generally performs within our requirements of 3 per cent or lower error in over 71 per cent of individual energy bins and extends to 15 per cent in the worst case scenarios. We do not find a consistent trend in increasing error as a function of any parameter, indicating that such errors are stochastically spread around the true model (as indicated previously in Figs 7 and 10). However, this means that small individual biases towards over- or underestimating the spectrum occur in small pockets of parameter space. This impacts the final posterior in the worst scenarios. We discuss this further with an example fit to a simulated spectrum in Section 4.4.

In Fig. 11, we show the distribution of percentage residual differences between RTFAST and RTDIST [see equation (1)]
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Distribution of the percentage residual difference between RT-FAST and RTDIST as a function of energy. The shade of region indicates the 25–75 per cent, 5–95 per cent, and 1–99 per cent region in increasing lightness, respectively.</figcaption></figure>

×100 per cent] as a function of energy for the same test set. This representation averages over the entire parameter space and instead shows general trends in errors for the emulator. The slight increasing of percentage error difference towards the highest and lowest energies is loosely related to the increasing steepness of the power law slope with higher Γ. As the most significant changes from the continuum (absorption edges and emission lines) do not deviate significantly in error from the errors in the rest of the continuum, RTFAST appears to capture much of the non-linear behaviour resultant from the general relativistic effects that affect the reflection spectrum. This partly shows the effectiveness of utilizing PCA decomposition to preserve correlation between energy bins by ensuring that the continuum does not dominate the spectrum in training and neglect the emulation of emission lines due to their small contribution to the overall spectrum.

## 4.4 Recovering parameters

In this section, we show that RTFAST can retrieve posteriors for simulated data from RTDIST. As the original model contains multiple degenerate pairs of parameters (notably degeneracy between the distance and mass of the black hole in the inferred ionization of the disc in equation 10 of I22) and some of the parameters cannot be constrained from the spectra alone, we constrain the problem to a smaller suite of parameters in which we can retrieve reasonable posteriors.

We simulate an observation with 260ks XMM-Newton exposure with xspec's (Arnaud 1996) fakeit routine, basing our source parameters on Ark 564 (Kara et al. 2013, 2016; I22), listed in Table 3. We use the open source software NDSPEC [https://github.com/matteolucchini1/nDspec, Lucchini et al. (in preparation)] for implementation of instrument response convolution throughout the rest of this paper. We perform Bayesian parameter estimation using the EMCEE package (Foreman-Mackey et al. 2013).

We specify our priors in Table 2, but fix all parameters to their true values other than $h, a, i, \Gamma, D_{\text{kpc}}, A_{\text{fe}}, \log N_e, N_h$, and $A_{\text{norm}}$. We fix $z$ as redshift constraints for AGN can be found from observations in other wavelengths. For simplicity, we fix inner radius to the ISO as X-ray reflection models are generally sensitive to either spin or inner radius. We fix $h/r$ (scale height of the disc), $b_1$ and $b_2$ as these have been generally unconstrained in previous fits to data with RTDIST. We fix $kT_e$ as we are unable to observe the high-energy roll-over at approximately 50keV due to XMM-Newton only observing between 0.1 and 20keV.

Due to the exceptionally high counts and to speed up computations, we use a Gaussian likelihood. Errors are generated by xspec's fakeit

<table><caption>**Table 2.** Priors for a more realistic fitting scenario. $\mathcal{U}(a, b)$ and $\log \mathcal{U}(a, b)$ indicate uniform and log uniform priors bounded between lower bound $a$ and upper bound $b$.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Prior</th></tr></thead><tbody><tr><td>$h$ ($R_g$)</td><td>Height of the corona in black hole gravitational radii</td><td>$\log \mathcal{U}(1.5, 10)$</td></tr><tr><td>a</td><td>Spin</td><td>$\mathcal{U}(0.5, 0.998)$</td></tr><tr><td>$i$ (°)</td><td>Inclination (degrees)</td><td>$\log \mathcal{U}(30, 80)$</td></tr><tr><td>$r_{\text{in}}$ (ISCOs)</td><td>Inner radius of disc</td><td>$\log \mathcal{U}(1, 20)$</td></tr><tr><td>Γ</td><td>Photon index</td><td>$\mathcal{U}(2, 2.75)$</td></tr><tr><td>$D$ ($10^5$kpc)</td><td>Distance in kpc</td><td>$\log \mathcal{U}(3.5, 500)$</td></tr><tr><td>$A_{\text{fe}}$</td><td>Iron abundance in solar units</td><td>$\log \mathcal{U}(0.5, 3)$</td></tr><tr><td>$\log(n_e/\text{cm}^{-3})$</td><td>Electron density of the disc</td><td>$\mathcal{U}(15, 20)$</td></tr><tr><td>$N_h$ ($10^{20}\text{cm}^{-2}$)</td><td>Hydrogen column density</td><td>$\log \mathcal{U}(10^{-3}, 1)$</td></tr><tr><td>$A_{\text{norm}}$ ($10^{-4}$)</td><td>RTDIST normalization</td><td>$\log \mathcal{U}(10^{-4}, 0.1)$</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Posteriors of fitting for a more realistic scenario where $2\sigma$ contours are plotted. The plot limits indicate the $3\sigma$ sigma boundary. We show five of the eight parameters varied. The full posteriors can be found in Appendix B. The blue lines show the true parameters used to generate the simulated observation.</figcaption></figure>

routine that incorporate the XMM-Newton PN background spectrum as well as Poissonian uncertainties. We use the spectrum between 0.3 and 10 keV. We employ 100 walkers with 70 000 steps each for our MCMC, and otherwise use the default parameters.

Fig. 12 shows the posteriors for five parameters: height, spin, $\Gamma$, $\log N_e$, and $A_{\text{norm}}$ (a corner plot showing the posteriors of all the free parameters is presented in Appendix B). Table 3 shows the full list of constrained parameters. The one-dimensional (1D) posteriors for the majority of parameters include the truth within the $3\sigma$ bounds, and sometimes fall within the $2\sigma$ bounds of the posterior. The two-dimensional (2D) posteriors show a clear bias in the some of the inferred parameters and correlations between parameters. Notably, the relationship between $A_{\text{norm}}$, a proxy for the accretion rate, and $h$, the height of the corona, shows a clear ‘banana’ degeneracy. This is realistic since a higher accretion rate is required to create the necessary reflection observed in the spectrum. There is also a bias towards lower values of spin and $\log N_e$ for this fit, leading the true value to lie outside of the $3\sigma$ bounds of the posterior.
<table><caption>**Table 3.** Input parameters and posterior fitted parameters for the simulated observations with $3\sigma$ uncertainties as recovered from fitting with RTFAST.</caption><thead><tr><th>Parameter</th><th>Input</th><th>Fitted parameters</th></tr></thead><tbody><tr><td><i>h</i> (<i>R</i><sub>g</sub>)</td><td>6</td><td>4.81<sup>+0.60</sup><sub>-0.48</sub></td></tr><tr><td><i>a</i></td><td>0.9</td><td>0.64<sup>+0.05</sup><sub>-0.06</sub></td></tr><tr><td><i>i</i> (°)</td><td>57</td><td>57.78<sup>+0.40</sup><sub>-0.42</sub></td></tr><tr><td><i>r</i><sub>in</sub> (ISCOs)</td><td>1</td><td>1</td></tr><tr><td><i>r</i><sub>out</sub> (<i>R</i><sub>g</sub>)</td><td>2 × 10<sup>4</sup></td><td>2 × 10<sup>4</sup></td></tr><tr><td><i>z</i></td><td>0.024917</td><td>0.024917</td></tr><tr><td>Γ</td><td>2.45</td><td>2.44<sup>+0.01</sup><sub>-0.01</sub></td></tr><tr><td><i>D</i> (10<sup>5</sup> kpc)</td><td>1</td><td>0.77<sup>+0.087</sup><sub>-0.073</sub></td></tr><tr><td><i>A</i><sub>fe</sub></td><td>1</td><td>0.98<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>log(<i>n</i><sub>e</sub>/cm<sup>-3</sup>)</td><td>17</td><td>16.41<sup>+0.08</sup><sub>-0.07</sub></td></tr><tr><td><i>kT</i><sub>e</sub></td><td>50</td><td>50</td></tr><tr><td><i>1</i>/<i>b</i></td><td>1</td><td>1</td></tr><tr><td><i>M</i><sub>BH</sub></td><td>3 × 10<sup>6</sup></td><td>3 × 10<sup>6</sup></td></tr><tr><td><i>h</i>/<i>r</i></td><td>0.02</td><td>0.02</td></tr><tr><td><i>b</i><sub>1</sub></td><td>0</td><td>0</td></tr><tr><td><i>b</i><sub>2</sub></td><td>0</td><td>0</td></tr><tr><td><i>N</i><sub>h</sub> (10<sup>20</sup>cm<sup>-2</sup>)</td><td>5</td><td>4.82<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td><i>A</i> (10<sup>-4</sup>)</td><td>2.2</td><td>3.1<sup>+0.6</sup><sub>-0.3</sub></td></tr></tbody></table>

This bias in some recovered parameters is physically explainable. Fig. 13, top panel, shows an observation simulated with RTDIST along with the spectrum predicted by the emulator for the exact same set of input parameters. We also show example spectra generated from posterior draws, and the mean of those posterior draws. The bottom panel shows the error-weighted residuals of the mean of the posterior draws alongside the same error-weighted residuals of the ‘correct’ spectrum given the true parameters predicted by the emulator. The RTFAST spectrum overpredicts a considerable soft excess. As a result, the model will underestimate $\log N_e$ to compensate for this bias during inference.

The inference of lower spin is better understood when freezing all parameters to the true values and then evaluating the likelihood on a 2-d grid of mass and spin. In Fig. 14, we took a grid of 10 points in mass and spin and evaluated the log-likelihood of each point on the grid. We took the normalized logarithm of the negative log-likelihood and plotted the contours. For a specified mass of $3 \times 10^6 M_\odot$, the spin with the highest likelihood is around 0.7. In contrast, Fig. 15 plots the likelihood as a function of spin when all other parameters are fixed to the true value. The point of maximum likelihood is clearly located extremely close to the truth, identifying a high, but not maximal spin as the most likely outcome. This implies a bias in the emulator such that it produces a spectrum that would be identified at lower spin at a slightly smaller mass for the original numerical model. When mass is fixed, the effects of spin on the spectrum are well-modelled and results in the correct inferred spin.

We plot the 2D likelihoods in the same manner in Appendix C for a number of other parameters combinations which show this bias in the emulator predominately drives a bias in the inferred mass spin likelihood for this particular simulated data. We also plot similar figures as Fig. 15 for all parameters in Appendix C.

This significant bias in the modelling is partly caused by the emulator’s unusually bad performance for this particular set of parameters of Ark 564. This particular spectrum has deviations of more than 10 per cent in the soft X-ray which places it within the 5 per cent worst performing spectra for the emulator. This amount

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Plots of the fitted spectrum, the emulated spectrum when the true parameters are inputted and the true spectrum generated by RTDIST. Posterior draws are plotted in orange, while the mean of these draws is plotted in red. The true simulated spectrum is plotted in black, with errors from Poisson noise and systematics. The emulator output for the true parameters is plotted in blue.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** The two dimensional log-likelihoods of fitting the spectrum, varying mass and spin, plotted in contours through parameter space. Brighter contours indicate higher likelihood. The red lines indicate the true parameters and what should be the peak of likelihood is represented by the red dot where they intersect.</figcaption></figure>

of bias in inferred parameters is thus an upper limit to the emulator’s worst performance. The bias in the model could be mitigated in this instance by fitting the spectrum between 1 and 10 keV, as the model is otherwise accurate for the rest of the spectrum. Notably, modelling the soft excess in the X-ray below 1–2 keV has been a well-known
<figure><img src="image_16.png" /><figcaption>**Figure 15.** The log likelihood as a function of spin when all other parameters are kept fixed. The true spin is plotted as a vertical dashed line at 0.9.</figcaption></figure>

problem for X-ray reflection and often requires additional model components (Lewin et al. 2022; Masterson et al. 2022; Xu et al. 2022). Using RTFAST between 1 and 10 keV more generally may be a better approach to mitigate biases caused by the emulator (particularly as the emulator performs the worst between 0.1–1 keV and 10–20 keV (which can be seen in Fig. 11) as well as avoid issues in modelling of the soft excess.

Modelling actual data is much more difficult. Maximum likelihood fits to the data may be unstable, and one may achieve different best fits for different starting parameters, indicating that degeneracies are present in the data. Thus, it often requires making assumptions about the geometry of the system to eliminate degeneracies by fixing certain parameters, similar to other analyses with RTDIST when fitting to real data (Nathan et al. 2024). Including timing products alongside the spectrum in the inference assists greatly with reducing degeneracy (as seen in I22) but cannot solve all of these challenges. RTFAST will allow for greater exploration of these degeneracies by enabling fast posterior calculations, even for complex physical models.

## 4.5 Poor extrapolation

In this section, we show the pitfalls of calling RTFAST outside of the training data bounds. Neural networks and other machine learning techniques have been repeatedly shown to extrapolate extremely poorly outside of their training bounds (Pastore & Carnini 2021) and RTFAST is not an exception to this. In Fig. 16, we generate a sample with a coronal height, *h*, of 400 $R_g$ from both RTFAST and RTDIST. We can clearly see that RTFAST gives a completely nonsense result, overestimating the spectrum by many orders of magnitude. Note that for other parameters than *h*, the effect can be more subtle and is more likely to be misleading. Users should never, under any circumstance, utilize parameter bounds outside those imposed in Section 3.1.2.

# 5 DISCUSSION AND CONCLUSIONS

## 5.1 Emulator

We have designed and implemented RTFAST, an emulator that acts as a drop-in replacement for the model RTDIST for AGN. RTFAST is able to recreate theoretical model spectra to within 3 per cent precision of the original model with an $O(10^2)$ speed up in conservative conditions. This is the first emulator for the reltrans

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Spectrum generated by RTFAST outside of the training bounds compared to RTDIST. We set coronal height (*h*) to 400 $R_g$ and we observe that RTFAST produces extremely unrealistic photon fluxes. This is a clear demonstration of RTFAST being unable to extrapolate behaviour.</figcaption></figure>

model suite and the first emulator to achieve this level of precision for this many physically meaningful model parameters within X-ray spectroscopy. The wide range of magnitude in model spectra has previously made this difficult to achieve when training neural networks to this level of precision.

By enabling computations on the GPU, the emulator allows us to harness parallelization to yield an 8000 times speed up over the original model when called sequentially. Because neural networks are differentiable models, where we can compute gradients of the model output as a function of input parameters, the emulator also allows the use of modern optimization and sampling algorithms to greatly speed up analysis. This in turn enables full Bayesian analysis on AGN and, in future, X-ray binaries. We will use this in future work to perform more robust measurements of black hole properties and perform simultaneous modelling of multiple sources to retrieve information about the black hole population as a whole.

Upon extensive testing in Section 4, we have found that there are small biases to inferred parameters introduced by RTFAST. These small biases are physically explainable and the emulator seems to capture important physical relationships between parameters. The biases depends on the section of parameter space being investigated, as well as what selection of parameters are being inferred, and any users of RTFAST should be careful about where these biases can creep into their interpretation of an observation. We suggest that users should always check the fitted model and plot posterior draws from the original RTDIST to assess whether they should trust the final results coming from RTFAST.

It should also be noted that the computational gains of RTFAST are relatively independent of the inclusion of additional physics in the case of RTDIST. Introducing complexity into the physical model is likely to introduce only small differences to the final model output, thus the complexity of RTFAST will likely have to increase only marginally to be able to model these scenarios. Even if future iterations of RTDIST were to take 100 s to evaluate instead of ≈ 0.5 seconds, RTFAST will still take on order of $10^{-3}$ seconds to evaluate.
We suggest that emulators such as RTFAST will allow us to probe far more complex physical scenarios while still being able to perform the millions of evaluations required for proper inference of the posteriors.

For those in the community who wish to implement this type of work on their own models, we recommend reducing the complexity and size of the model's output (here: the spectrum) using algorithms that preserve important relationships in the model (e.g. PCA). For instance, for models that feature large suites of absorption and emission lines, preserving their appearance as a function of input parameters in a way that the neural network does not need to learn will allow developers to keep the network size low. Domain knowledge of both the model and astronomical sources modelled was key in constraining the problem as highlighted in Section 3.1.2. When developing emulators, we recommend considering both observational and theoretical constraints to the model parameter space and the model itself. We discuss this in more detail in Appendix D.

## 5.2 Bayesian analysis

In Section 3.5, we showed that RTFAST is, at its most conservative estimation, 200 times faster than the original RTDIST model when called sequentially. This becomes even faster once considering the ability to evaluate multiple models in vector form – calculating spectra up to $O(10^4)$ times faster than the original model. Current MCMC inference using RTDIST can take up to a month or more, particularly if multiple data sets are fit jointly. Using RTFAST reduces those fits to a few hours at most. This allows for further in-depth exploration of the posteriors of the degenerate parameter space of RTDIST, illustrated in Section 4.4.

RTFAST was developed to make it possible to properly explore the large parameter space of RTDIST. RTFAST enables us to do this via multiple avenues: Monte Carlo Markov chain (which already has widespread use in the field), nested sampling (which has already been somewhat adopted by some of the community; Buchner et al. 2014), Hamilton Monte Carlo (a Monte Carlo variant that uses gradients in respect to parameters to achieve faster optimization; Duane et al. 1987) and simulation-based inference (a machine learning technique that approximates the posterior of an observation via simulations of a theoretical model; Cranmer et al. 2020). All of these techniques (with the exception of Hamilton Monte Carlo) are already possible to use with RTDIST but are prohibitively computationally expensive to use practically. We have implemented all of these techniques in an accompanying notebook to RTFAST in the RTFAST repository https://github.com/SRON-API-DataMagic/RTFAST.

It is possible to use RTFAST to produce MCMC chains to a considerable length towards convergence and then switch to RTDIST for the final steps in the chain towards convergence to avoid the biases potentially present in sampling posteriors with RTFAST. This allows for a dramatic reduction in computation time while still retrieving the complex posteriors.

## 5.3 Future work

RTFAST approximates the spectrum generated by the RTDIST model, but RTDIST comprises a number of different physical processes (i.e. reflection from the disc (García et al. 2013), relativistic effects on photons emitted by the corona (Dauser, García & Wilms 2016), etc). Our next-generation emulator will focus on reproducing

these individual effects using the design philosophy of RTFAST as a template. This would allow for a modular design that enables easy enabling and re-enabling of the true underlying physical model. An example of this would be to emulate the disc illumination profile resulting from the geometry of the system rather than computing it directly with ray-tracing. Retraining a single component may also result in less required training data when updating the emulator compared to retraining the emulator on the total output.

RTFAST currently contains systematic errors of order of 3 per cent which is just accurate enough to be comparable with the systematic errors of current generation X-ray instruments such as XMM-Newton. We would like to achieve better precision in recreating the spectrum in the future, especially to reduce or even eliminate bias in the inferred parameters in the future. We expect to be able to reduce the error by at least one order of magnitude by changing our strategy slightly.

Currently, we emulate all components in the spectrum at once. The RTDIST spectrum is mainly made up of 2 components: a comptonized power-law approximating the corona (nthcomp), and a reflection spectrum resulting from reprocessed coronal emission from the disc. The reflection spectrum itself is by far the most computationally expensive component of RTDIST as it requires the calculation of relativistic effects. This component typically composes approximately 10 per cent of the total flux in observed spectra from both AGN and X-ray binary sources. Emulating the reflection spectrum alone would mean that the resultant error in the total spectrum would be an order of magnitude smaller and prevent strong biases in parameters inferred from the continuum. This would likely require returning to the reltransDCp flavour of reltrans, due to the entanglement of normalization and the spectrum by $A_{\text{norm}}$ in RTDIST.

Another potential extension would be to emulate the model after convolution with the instrument response. This would allow for the bottleneck of convolution of the instrument response to be avoided entirely, a topic that is particularly important for instruments such as XRISM which feature extremely large, high-resolution instrument responses. We opt to not adopt this strategy in this paper to keep the emulator somewhat instrument independent.

RTDIST contains functionality to compute timing products caused by reverberation between the corona and disc. In the future, we will design an instrument-independent emulator that will allow for the computation of RTDIST timing products to further constrain black hole parameters. This is outside of the scope of this paper due to the complexity of emulating the cross-spectrum, caused by its higher dimensionality and greater range in model outputs that could not be constrained in the same way as in Section 3.1.2.

While RTFAST is an ensemble where the final outputs of its constituent neural networks are simply averaged together, there are potentially better approaches to this ensemble technique such as the weighting of models by their relative accuracy to the final spectrum or explicitly training a considerably larger ensemble and pruning networks to find a more optimal combination of networks to better approximate RTDIST.

During development, we attempted use of an active learning technique called Query by Dropout Committee (Ducoffe & Precioso 2015). This technique attempts to select the most effective data to generate during training to reduce the total amount of data required to train RTFAST. This ultimately was not used in RTFAST as generating data during training became unwieldy when performing architecture and learning sweeps in Sections 3.2.1 and 3.3. This strategy will likely be more appropriate for the cross-spectrum as it both takes more time to compute and will be larger in data-size than the spectrum alone. As such, we would like to re-attempt this
technique when emulating the cross-spectrum. Another potential technique for reducing the amount of required training data is active subspace sampling (Constantine, Dow & Wang 2014). Active subspace sampling uses evaluations of the gradient to find the areas of strongest variability, and subsequently uses these directions to construct a sampling surface on a lower dimensional subspace.

## ACKNOWLEDGEMENTS

We would like to thank the referee for their constructive comments. This work used the Dutch national e-infrastructure with the support of the SURF Cooperative using grants EINF-9197 and EINF-5505. We'd like to express our thanks to D.Burke in particular and the CIAO help desk for their assistance in the use of the SHERPA package within this paper. AI acknowledges support from the Royal Society.

## DATA AVAILABILITY

The training data generated for this project is too large to be open access; however, the parameters for all generated spectra in this project can be found at https://doi.org/10.5281/zenodo.14264681. The RTDIST model used to generate the training data is currently not available publicly, but is available on request. The training program is publicly available and can be found at https://github.com/SRON-API-DataMagic/rtdist-emulator. If a user wishes to use the trained emulator for fitting of observations, they can find the trained model and the required program to install and use it at https://github.com/SRON-API-DataMagic/RTFAST.

## REFERENCES

* Alsing J. et al., 2020, ApJS, 249, 5

* Arnaud K., 1996, Astronomical Society of the Pacific (ASP), MS 179 Utah Valley University, 800 W. University Parkway, Orem, UT 84058-5999 USA, in Jacoby G. H., Barnes J., eds., Astronomical Data Analysis Software and Systems V, ASP Conference Series. Vol. 101, p. 17

* Baker M. R., Patil R. B., 1998, Reliable Comput., 4, 235

* Bastos L. S., O'hagan A., 2009, Technometrics, 51, 425

* Bellman R., 1966, Science, 153, 34

* Bird S., Rogers K. K., Peiris H. V., Verde L., Font-Ribera A., Pontzen A., 2019, J. Cosmol. Astropart. Phys., 2019, 050

* Blandford R., Königl A., 1979, ApJ, 232, 34

* Buchner J., 2023, Stat. Surv., 17, 169

* Buchner J. et al., 2014, A&A, 564, A125

* Constantine P. G., Dow E., Wang Q., 2014, SIAM J. Sci. Comput., 36, A1500

* Cranmer K., Brehmer J., Louppe G., 2020, Proc. Natl. Acad. Sci., 117, 30055

* Cunningham C. T., Bardeen J. M., 1973, ApJ, 183, 237

* Dauser T., Garcia J., Wilms J., Böck M., Brenneman L. W., Falanga M., Fukumura K., Reynolds C. S., 2013, MNRAS, 430, 1694

* Dauser T., García J., Wilms J., 2016, Astron. Nachr., 337, 362

* Dias S. D., Vaughan S., Lefkir M., Wynn G., 2024, MNRAS, 529, 1752

* Dong Y., García J. A., Steiner J. F., Gou L., 2020, MNRAS, 493, 4409

* Dovčiak M., Karas V., Yaqoob T., 2004, ApJS, 153, 205

* Draghis P. A., Miller J. M., Costantini E., Gallo L. C., Reynolds M., Tomsick J. A., Zoghbi A., 2024, ApJ, 969, 40

* Duane S., Kennedy A., Pendleton B. J., Roweth D., 1987, Phys. Lett. B, 195, 216

* Ducoffe M., Precioso F., 2015, preprint (arXiv:1511.06412)

* Fabian A., 1989, MNRAS, 238, 41P

* Fabian A. C., Rees M. J., Stella L., White N. E., 1989, MNRAS, 238, 729

* Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

* Freedman W. L., 2021, ApJ, 919, 16

* Freeman P., Doe S., Siemiginowska A., 2001, in Astronomical Data Analysis. p. International Symposium on Optical Science and Technology, San Diego, CA, United States, 76

* Ganaie M. A., Hu M., Malik A. K., Tanveer M., Sughanthan P. N., 2022, Eng. Appl. Artif. Intell., 115, 105151

* García J., Dauser T., Reynolds C., Kallman T., McClintock J., Wilms J., Eikmann W., 2013, ApJ, 768, 146

* García J. A. et al., 2018, ApJ, 864, 25

* Gendreau K. C., Arzoumanian Z., Okajima T., 2012, in Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray. p. 322, SPIE Astronomical Telescopes + Instrumentation, Amsterdam, Netherlands

* Girolami M., Calderhead B., 2011, J. R. Stat. Soc. Ser. B: Stat. Methodol., 73, 123

* Guainazzi M., Kirsch M., Haberl F., de Juan Ovelar M., 2013, Technical Report, XMM-Newton Calibration Technical Note. XMM-SOCCAL-TN-0018, Calibration status document, ESA-ESAC, Villafranca del.

* Heitmann K., Higdon D., White M., Habib S., Williams B. J., Lawrence E., Wagner C., 2009, ApJ, 705, 156

* Heitmann K., Lawrence E., Kwan J., Habib S., Higdon D., 2013, ApJ, 780, 111

* Hendrycks D., Gimpel K., 2016, preprint (arXiv:1606.08415)

* Himes M. D. et al., 2022, Planet. Sci. J., 3, 91

* Ingram A., Mastroserio G., Dauser T., Hovenkamp P., van der Klis M., García J. A., 2019, MNRAS, 488, 324

* Ingram A. et al., 2022, MNRAS, 509, 619(122)

* Jansen F. et al., 2001, A&A, 365, L1

* Jiang J. et al., 2019, MNRAS, 489, 3436

* Kammoun E. S., Barret D., Peille P., Willingale R., Dauser T., Wilms J., Guainazzi M., Miller J. M., 2022, A&A, 664, A29

* Kara E., Fabian A., Cackett E., Uttley P., Wilkins D., Zoghbi A., 2013, MNRAS, 434, 1129

* Kara E., Alston W., Fabian A., Cackett E., Uttley P., Reynolds C., Zoghbi A., 2016, MNRAS, 462, 511

* Kaushal N., Villaescusa-Navarro F., Giusarma E., Li Y., Hawry C., Reyes M., 2022, ApJ, 930, 115

* Kerzendorf W. E., Vogl C., Buchner J., Contardo G., Williamson M., Van der Smagt P., 2021, ApJ, 910, L23

* Kingma D. P., 2014, preprint (arXiv:1412.6980)

* Kotov O., Churazov E., Gilfanov M., 2001, MNRAS, 327, 799

* Kwan J., Heitmann K., Habib S., Padmanabhan N., Lawrence E., Finkel H., Frontiere N., Pope A., 2015, ApJ, 810, 35

* Lakshminarayanan B., Pritzel A., Blundell C., 2017, Adv. Neural Inform. Process. Sys., 30

* Levinson A., Blandford R., 1996, ApJ, 456, L29

* Lewin C. et al., 2022, ApJ, 939, 109

* Loshchilov I., 2017, preprint (arXiv:1711.05101)

* Loshchilov I., Hutter F., 2016, preprint (arXiv:1608.03983)

* Ludlam R. et al., 2020, ApJ, 895, 45

* Masterson M. et al., 2022, ApJ, 934, 35

* Matzeu G. et al., 2022, MNRAS, 515, 6172

* McClintock J. E. et al., 2011, Class. Quantum Gravity, 28, 114009

* McKay M. D., Beckman R. J., Conover W. J., 2000, Technometrics, 42, 55

* Mendes-Moreira J., Soares C., Jorge A. M., Sousa J. F. D., 2012, Acm Comput. Surveys (csur), 45, 1

* Nathan E. et al., 2024, MNRAS, 533, 2441

* Nemmen R. S., Georganopoulos M., Guiriec S., Meyer E. T., Gehrels N., Sambruna R. M., 2012, Science, 338, 1445

* Özel F., Freire P., 2016, ARA&A, 54, 401

* Pastore A., Carnini M., 2021, J. Phys. G: Nucl. Part. Phys., 48, 084001

* Paszke A. et al., 2019, Curran Associates, Inc, p. 8024

* Remillard R. A., McClintock J. E., 2006, Annu. Rev. Astron. Astrophys., 44, 49

* Rivera J. P., Verrelst J., Gómez-Dans J., Muñoz-Marí J., Moreno J., Camps-Valls G., 2015, Remote Sensing, 7, 9347

* Romero G. E., Boettcher M., Markoff S., Tavecchio F., 2017, Space Sci. Rev., 207, 5

* Ross R., Fabian A., 2005, MNRAS, 358, 211
Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Sobolewska M., Papadakis I., 2009, MNRAS, 399, 1597

Spurio Mancini A., Piras D., Alsing J., Joachimi B., Hobson M. P., 2022,
MNRAS, 511, 1771

Sunyaev R., Truemper J., 1979, Nature, 279, 506

Sutskever I., Martens J., Dahl G., Hinton G., 2013, In International conference on machine learning. p. 1139

Thorne K. S., Price R. H., 1975, ApJ, 195, L101

Uttley P., Casella P., 2014, Space Sci. Rev., 183, 453

Vicent J., Verrelst J., Rivera-Caicedo J. P., Sabater N., Muñoz-Marí J., Camps-Valls G., Moreno J., 2018, IEEE J. Selected Topics Applied Earth Observ. Remote Sensing, 11, 4918

Virtanen P. et al., 2020, Nature Methods, 17, 261

Weisskopf M. C. et al., 2016, Results Phys., 6, 1179

Wilms J., Allen A., McCray R., 2000, ApJ, 542, 914

Xu Y. et al., 2022, MNRAS, 513, 1910

Zdziarski A. A., Johnson W. N., Magdziarz P., 1996, MNRAS, 283, 193

APPENDIX A: PHYSICAL PRIOR DISTRIBUTIONS

In Fig. A1, we plot the distributions of all training parameters as specified in Section 3.1.2.

<figure><img src="image_18.png" /><figcaption>**Figure A1.** Corner plot of all training set parameters and their distributions. Possible training data to be generated is constrained by the physical constraints outlined in Section 3.1.2.</figcaption></figure>
## APPENDIX B: FULL POSTERIORS

In Fig. B1, we report the full posteriors of the fit in Section 4.4. We also plot the MCMC chains in their entirety in Fig. B2.

<figure><img src="image_19.png" /><figcaption>**Figure B1.** Posteriors of fitting. The $2\sigma$ contours are explicitly plotted. The blue lines show the true parameters used to generate the simulated observation.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure B2.** MCMC chains as a function of steps of each walker. The blue dotted line indicates the true parameters.</figcaption></figure>

## APPENDIX C: LIKELIHOOD SPACE

In this section, we show and further discuss the 2D biases for the inference in Section 4.4. In Figs C1 and C2, the maximum likelihood is very close to the true parameters when all other parameters are fixed. This supports that the emulator correctly captures the relationship between spin and $\Gamma$ and $N_h$.

<figure><img src="image_21.png" /><figcaption>**Figure C1.** 2D grid of the likelihood for spin versus $\Gamma$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure C2.** 2D grid of the likelihood for spin versus $n_H$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure C3.** 2D grid of the likelihood for coronal height versus spin.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure C4.** 2D grid of the likelihood for spin versus inclination.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure C5.** 2D grid of the likelihood for spin versus $A_{fe}$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure C6.** 2D grid of the likelihood for spin versus $A_{norm}$.</figcaption></figure>

Figs C3 and C4 show that the likelihood space is somewhat more complicated for inclination and height, which could potentially bias the inference to values of incorrect spin (with not as much an effect as Fig. 14). Fig. C3 shows that, while the point of highest likelihood is near the truth, there are areas of higher likelihood at higher values of height where spin is lower. This is somewhat

<figure><img src="image_27.png" /><figcaption>**Figure C7.** 2D grid of the likelihood for spin versus z.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Figure C8.** 2D grid of the likelihood for spin versus b1.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure C9.** 2D grid of the likelihood for spin versus b2.</figcaption></figure>

realistic (as increasing height allows for more photons to escape the influence of the black hole but decreasing spin also allows for more photons to escape the black hole due to less extreme general relativistic effects), but this likelihood space should be more clearly correlated. Fig. C4 shows a small split in the possible inferred inclination and spin of the system. The truth once again fits in the
<figure><img src="image_30.png" /><figcaption>**Figure C10.** 2D grid of the likelihood for spin versus boost.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Figure C11.** 2D grid of the likelihood for spin versus distance.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Figure C12.** 2D grid of the likelihood for spin versus scale height of the disc.</figcaption></figure>

space of highest likelihood, and the correct degenerate relationship of high spin and low coronal height and vice versa resulting in similar reflection spectra due to relativistic effects is recovered. Figs C5–C15 show the remaining parameters correlated with spin, excluding

<figure><img src="image_33.png" /><figcaption>**Figure C13.** 2D grid of the likelihood for spin versus kTe.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>**Figure C14.** 2D grid of the likelihood for spin versus log N<sub>e</sub>.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>**Figure C15.** 2D grid of the likelihood for spin versus inner radius.</figcaption></figure>

outer radius due to it being very uninformative to fitting of the spectrum.

Fig. C16 shows the one dimensional likelihoods for each parameter. We observe that, when every other parameter is fixed, the emulator recovers the correct parameter value as the point of maximum likelihood.
<figure><img src="image_36.png" /><figcaption>**Figure C16.** Plots of the log likelihoods as a function of each parameter while all parameters are fixed at the true value. The true parameter is plotted as a vertical dashed line.</figcaption></figure>

## APPENDIX D: BUILDING EMULATORS FOR PHYSICAL MODELS

This section contains the insights we developed building RTFAST-spectra. These insights are intended to be at least somewhat spectral model neutral, though every problem in the machine learning context has its own quirks that prevent one-to-one translation between problems.

Feed-forward neural networks were overwhelmingly the best strategy in terms of emulation quickly and accurately. Other types of layers such as convolutions and long short-term memory were generally more inaccurate, or in the cases where they were of similar

accuracy, slower than the methodologically simple linear network. Convolutional layers performed particularly poorly. We posit that, for many spectra, different emission lines in the spectrum are often too far apart within the spectrum to be encoded within the space of a single convolutional kernel. This means that important correlations are not captured by this type of neural network layer. A sufficiently wide enough kernel to capture these correlations quickly approaches that of a fully connected neural network, and loses any unique benefits offered by convolutional kernels.

Most notably, improvement in performance was overwhelmingly dominated by the addition of more data in training. We attempted several strategies in improving choice of data generation for training,
but found that all more intelligent strategies involving attempted
probing of the neural network's inaccuracies tended towards produc-
ing almost uniformly distributed points across the parameter space
chosen once enough data had been generated.

After a very large amount of data generated, we generally found
that the neural networks improved in training logarithmically, with
validation loss tightly linked to the training loss – not exhibiting
evidence of over-fitting. While it is likely that we could have trained
for much longer until true over-fitting had occurred, we found it more
effective to build an ensemble to yield the improvement we desired.
This also removed some systematic biases we found neural networks
frequently exhibited when emulating spectra (this is particularly
obvious in Fig. 7). The ensemble generally mitigated systematic
bias in the individual neural networks, and reduced the extremes of
outlying bad performance. Adding more models into the ensemble
seems to lead to a saturation in performance gains, but we could not
test this extensively due to computational resource constraints. This
generally implies that ensembles for these types of problems act much
closer to the theoretical best case for ensembles. Ensembles work
most effectively when errors in the neural network are from under-
fitting rather than due to systematic bias in the training data. If the
errors are approximately Gaussian due to the random initialization

of weights at the beginning, this means that simple averaging of
the outputs of multiple neural network outputs will more closely
approximate the truth.

Notably, the choice of loss functions also led to considerable
improvements compared to the choice of model architectures.
Choosing the standard mean-squared error loss to train the network
to emulate PCA components, we found the neural network to be
almost unable to learn any relation between the parameters and
PCA vectors. Weighting vectors allowed the neural network to first
fit simple correlations and then learn the higher order correlations
encoded in the PCA vectors as they became more relevant in the
weighting during training. Some experiments we performed on
non-decomposed forms of the original model were improved when
we encoded the distribution of the entire spectrum into the loss
via Kullback-Liebler divergence as an additional component to
the mean-squared error. As Kullback-Liebler divergence calculates
the difference between the distributions, this encodes correlations
between energy bins mathematically and punishes the neural network
for not taking these correlations into account.

This paper has been typeset from a TeX/LATEX file prepared by the author.