A Library of Synthetic X-Ray Spectra for Fitting Tidal Disruption Events

Tidal disruption events (TDEs) are a once-rare class of transient now discovered at a rapidly increasing pace. Current optical and soft X-ray wide-field surveys, such as ASAS-SN (Holoien et al. 2014, 2016), the Zwicky Transient Facility (van Velzen et al. 2021), and eROSITA (Sazonov et al. 2021), are discovering tens of new TDEs per year (see Saxton et al. 2020; van Velzen et al. 2020, for recent reviews), and almost one hundred candidate TDEs have been discovered to date. ⁵ Near-future surveys in the optical (Vera Rubin Observatory executing the Legacy Survey of Space and Time; Ivezić et al. 2019) and ultraviolet (ULTRASAT; Sagiv et al. 2014) will likely find hundreds to thousands more each year (Bricman & Gomboc 2020).

Such a large sample of TDEs raises the prospect of studying supermassive black hole (SMBH) mass and spin demographics. The canonical TDE forms when a star passes too close to an SMBH in a quiescent galactic nucleus, is torn apart, and then produces a luminous, multiwavelength flare (Hills 1975; Lidskii & Ozernoi 1979; Rees 1988; Evans & Kochanek 1989). The physical evolution and thermal emission from TDEs is governed by just a small number of free parameters (see Rossi et al. 2021 for a recent review). For the standard case of a main-sequence star disrupted on a nearly parabolic orbit, ⁶ the five controlling parameters are SMBH mass M_• and stellar mass (Rees 1988; Guillochon & Ramirez-Ruiz 2013; Ryu et al. 2020a, 2020b), SMBH spin a_• (Dai et al. 2015; Hayasaki et al. 2016; Wen et al. 2020), orbital pericenter (Guillochon & Ramirez-Ruiz 2013; Stone et al. 2013; Dai et al. 2015), and spin–orbit misalignment angle (Guillochon & Ramirez-Ruiz 2015; Hayasaki et al. 2016; Liptai et al. 2019). Other properties of the star, such as its metallicity, rotation rate, and age along the main sequence, can affect the evolution of debris, but are unlikely to do so at leading order (Gallegos-Garcia et al. 2018; Golightly et al. 2019; Kagaya et al. 2019; Law-Smith et al. 2019, 2020). The small number of free parameters associated with the flare suggests that M_• and a_• are usefully encoded in the resulting electromagnetic emission.

Measuring the distribution of M_• is a path toward quantifying the uncertain low-mass tail of the SMBH mass function, including elusive intermediate-mass black holes (IMBHs), and constraining SMBH seed formation scenarios (Volonteri 2010). Constraining the SMBH spin distribution would recover the growth history of SMBHs (Berti & Volonteri 2008). The a_• distribution is poorly understood at present, with most existing measurements arising from a single technique: iron Kα reflection spectroscopy (Reynolds 2014). The SMBH spin is imprinted on TDE observables in various ways, e.g., relativistic precessions (Stone & Loeb 2012; Guillochon & Ramirez-Ruiz 2015; Hayasaki et al. 2016; Liptai et al. 2019; Curd 2021; Andalman et al. 2022), the disk thermal X-ray spectrum (McClintock et al. 2006; Done et al. 2012), the TDE rate (Kesden 2012), and quasiperiodic oscillations in the X-ray spectrum (Dheeraj et al. 2019), raising the possibility of obtaining independent SMBH spin constraints, with different uncertainties and biases than those of the status quo.

Until recently, most attempts to measure M_• from TDE observations have relied on optical light curves (Mockler et al. 2019; Ryu et al. 2020c), the time evolution of which resembles theoretical mass fallback rates (Guillochon & Ramirez-Ruiz 2013). Unfortunately, first-principles models for TDE optical emission do not yet exist. More idealized models have not even converged on the underlying power source; at the time of writing, there is a vigorous debate about whether reprocessed accretion luminosity (Guillochon et al. 2014; Metzger & Stone 2016; Roth et al. 2016; Bonnerot & Lu 2020) or shock dissipation (Piran et al. 2015; Shiokawa et al. 2015) is the dominant power source for the optical photosphere.

The quasi-thermal soft X-ray emission observed in a fraction of TDEs is arguably better understood, as the only existing model for it invokes an origin in the innermost regions of the newly formed accretion flow (Ulmer 1999; Lodato & Rossi 2011). While uncertainties about the geometry of this inner accretion flow exist, particularly at early times, the X-ray-emitting material is generally matter that has dissipated orders of magnitude of energy. It is therefore reasonable to expect significant circularization. For example, the recent first-principles simulations of Andalman et al. (2022) find that even over the first six days of mass fallback alone, equatorial material accreting onto the SMBH passes through the innermost stable circular orbit with eccentricities e ≈ 0.4–0.8 (much smaller than typical post-disruption e ≈ 0.99, and likely to decrease further as the inner disk builds up and circularization progresses).

In recent work, we have combined Kerr metric ray-tracing (Psaltis & Johannsen 2012) with general relativistic slim-disk models (Sádowski 2009) to produce synthetic X-ray spectra for the inner regions of TDE disks. By fitting these models to X-ray observations, we have been able to constrain M_• and a_• for both SMBHs (Wen et al. 2020, hereafter W20) and IMBHs (Wen et al. 2021, hereafter W21). The IMBH case is particularly interesting as a novel probe of SMBH seed formation scenarios and unexplored mass ranges of ultralight dark matter candidates (W21). However, a practical weakness of our past work is the need to generate a large grid of computationally expensive synthetic X-ray spectra on the fly, for purposes of parameter estimation. The computational cost of generating custom spectral databases (∼10⁴ CPU hours for each TDE) would become prohibitive if applied at scale to much greater numbers of X-ray-bright TDEs, such as those expected in the eROSITA era (Khabibullin et al. 2014; Jonker et al. 2020; Sazonov et al. 2021).

In this paper, we address this problem, and thoroughly explore the parameter space of relativistic slim-disk models, by producing a large tabulated library of synthetic X-ray spectra. These spectra are derived from our original full ray-tracing model (W21). In Section 2, we review the underlying theoretical model and key assumptions. In Section 3, we present the tabulated library of synthetic X-ray spectra. We fit two previously studied TDEs by interpolating over this library in Section 4. In Section 5, we explore the different types of parameter degeneracy that arise in the multidimensional parameter space of X-ray-bright TDE disks, with the overall goal of determining which types of observations can constrain different underlying physical parameters. Finally, we conclude in Section 6. So that this library is widely usable as the TDE field transitions to large population studies, we have made it publicly accessible. ⁷

Our model is based on the work of W20 and W21, who used a general relativistic stationary slim-disk model (Abramowicz et al. 1996; Sádowski 2009) employing the ray-tracing code from Psaltis & Johannsen (2012) to calculate the synthetic spectra as seen by an observer. Here, we review the high-level details of the local X-ray emission and ray-tracing procedure. For the full details on relativistic slim-disk solutions and the ray-tracing techniques we use, we refer the reader to W20, W21, and the aforementioned citations in this section.

In order to calculate the local X-ray spectrum of the disk from the disk flux F(r), we follow Shimura & Takahara (1995), assuming a color-modified blackbody is emitted from each annulus of the disk. The specific intensity of the radiation, with frequency ν, can be expressed as

Here h is the Planck constant and k_B is the Boltzmann constant; f_c is commonly called the spectral hardening factor. T_e is the local effective temperature of an annular ring of the disk (assumed to be axisymmetric), which is determined by T_e = (F(r)/2σ)^1/4 (where F(r) is the flux from the disk at radius r and σ is the Stefan–Boltzmann constant). We use f_c values as estimated by the radiative transfer calculations of Davis & El-Abd (2019), and which are parameterized as

Here Σ(r) and Q(r) are the surface density and the strength of the vertical gravity at each annulus of the disk, respectively. Although this parameterization is fitted to the best theoretical models for f_c in the present literature, we note that they have only been tested/calibrated for accretion rate between 0.01 and 1 times the Eddington accretion rate, and only for nonspinning BHs. Throughout this paper, the Eddington accretion rate is given by

where η is the radiative efficiency of accretion, η₋₁ = η/0.1, and M₆ is the SMBH mass in units of a million solar masses, M₆ = M_•/(10⁶ M_⊙). Hereafter, we refer to accretion rate in dimensionless Eddington units, i.e., , and we assume η = 0.1 in all future references to Eddington ratios (even though the actual radiative efficiency of our slim-disk models can be significantly different).

We trace light rays as null geodesics from the image plane to the surface of the disk, which is truncated at an outer radius R_out and an inner radius R_in. We set R_out to two times the tidal radius R_t (Hills 1975),

Here, R_g = GM_•/c² is the gravitational radius and M_⋆ and R_⋆ are the mass and radius of the disrupted star, respectively. In calculating R_t (to compute R_out), we assume that the disrupted star has a solar mass and radius. When 2R_t > 600R_g, we reset it to be R_out = 600R_g. The error on the luminosity caused by this choice of R_out is always <1% (W21), because the outer disk is too cold to produce a significant X-ray luminosity. ⁸ On the other hand, when M_• > 4.2 × 10⁷ M_⊙, then R_out will usually be smaller than 15R_g. For such a case, the choice of R_out can affect the spectrum significantly. As the SMBH mass M_• increases, R_out/R_g will decrease quickly. If R_out < 15R_g, it is reset to R_out = 15R_g. Due to a singularity in the calculation of f_c for disk radii smaller than the innermost stable circular orbit (ISCO), we set R_in to the value of the ISCO. The error on the luminosity caused by different choices of inner edge is small. ⁹ A test in W21 shows that the error is <2% for different choices of spin and accretion rate for an intermediate-mass BH with inclination of 45°. W21 also finds that, for an accretion disk around a fast spinning BH, the error is less than 0.5% regardless of BH mass,.

If the peak accretion rate for TDEs is highly super-Eddington, the disk will be geometrically thick (H/R ∼ 0.5, where H is the half-disk thickness and R is the disk radius). When the disk inclination is close to edge-on, the edge of the disk will block most X-rays from the inner disk (W20). When the photon’s trajectory intersects the disk at any point during ray tracing, we stop further tracing and assume the photon is emitted from the first intersection point. In addition, we assume that the outer rim of the disk behaves as a single-temperature color-corrected blackbody, with an effective temperature equal to that of the standard annular effective temperature of the outermost annulus.

To summarize, in calculating the synthetic X-ray spectrum, we have assumed that:

• 1.  
  The inner disk is circular and equatorial.
• 2.  
  The disk structure is determined assuming a zero-torque inner boundary condition, large optical depth, and no self-irradiation. It is computed using the general relativistic slim-disk model of Sádowski (2009), which is similar to the classic Shakura–Sunyaev and Novikov–Thorne thin-disk models, but differs in its inclusion of (i) advective cooling and (ii) sub-Keplerian angular momentum when H/R is large (Abramowicz et al. 1988). These features allow slim disks to cover the super-Eddington accretion rates that may be present during the early stages of some TDEs.
• 3.  
  The local X-ray emissivity of the disk is consistent with modified multicolor blackbody radiation, e.g., Equation (1). The spectral hardening factor is calculated using Equation (2).
• 4.  
  The X-rays are generated at the photospheric surface of the disk at a finite height, which we assume to be equal to its scale height H. If a photon’s trajectory intersects the disk after emission, we assume that it will be absorbed completely. Photons follow null geodesics that can be strongly lensed by the Kerr metric they propagate through.
• 5.  
  The inner edge of the disk is set to be the ISCO. The outer edge of the disk is twice the tidal disruption radius (2R_t). If the outer edge is >600R_g (<15R_g ), we reset it to be 600R_g (15R_g).

The tabulated synthetic spectra in the library are produced by the ray-tracing model of W21, which includes gravitational redshift, Doppler, and lensing effects self-consistently, as well as the additional effect of angular momentum loss by radiation (which was not considered in W20). In the context of relativistic slim disks, the disk X-ray spectrum can be affected by the following parameters: the black hole mass M_• and spin a_•, the disk accretion rate and inclination θ, viscosity α, and spectral hardening factor f_c. In W20, we showed that the impact of α on the spectrum is small compared with that of M_•, , and θ; as a result, for the model spectral library, we fix α = 0.1. We calculate f_c following Davis & El-Abd (2019); this prescription is calibrated for a range of sub-Eddington accretion disk properties, and extending it to highly super-Eddington accretion rates introduces a systematic uncertainty that should be quantified with future modeling. In general, for the tabulated model, there are only four free disk parameters, M_•, a_•, , and θ, while α is fixed at α = 0.1 and f_c is estimated using Equation (2).

We introduce two model functions, slimd and slimdz, for spectral fitting. slimd fits the spectra for a given luminosity distance and redshift, while slimdz fits the spectra by calculating the luminosity distance from the redshift with a flat cosmological model. We note that redshift can be a free parameter for slimdz. For more details about these two model functions, we refer the reader to the menu given with the synthetic spectral library.

The library contains model spectra calculated on a grid of parameters M_•, a_•, , and θ. The size of the grid is 85 × 33 × 100 × 119. One can find the specific grid values in the model files, and the ranges of values are listed in Table 1. We use linear interpolation to estimate the spectra for parameter values between the grid points. For the parameter θ, we use linear interpolation to estimate the actual flux at each energy bin, while, for the parameters M_•, a_•, and , we use linear interpolation to estimate the logarithm of the flux. ¹⁰ We then compare the flux predicted by interpolating from grid points on the table with that from direct, on-the-fly ray-tracing calculations. The grid is sufficiently fine that the error caused by interpolation is less than 10% for nearly all cases (more on this below). For each grid point, we calculate the spectral luminosity density in 500 photon energy bins spread over the 0.2–10 keV energy range.

Table 1. Parameter Ranges for Model Grid

M• [M⊙]	a•	[Edd]	θ [deg]
[103, 9 × 107]	[–0.998, 0.998]	[0.05, 600]	[2, 90]

Note. For calculating the Eddington accretion rate , we assume a constant radiation efficiency of η = 0.1 irrespective of the real, a_•-dependent efficiency.

Download table as:  ASCIITypeset image

Figure 1 shows the spectra and the error on the luminosity density caused by the linear interpolation between the grid points. For the upper panels, the spectra are calculated with the full ray-tracing code even if their parameters do not lie on the grid points of the library. We also produce spectra for disks with the same parameters, but by using linear interpolation on the library of pre-tabulated spectral models. We compare the full ray-traced spectra with those obtained through the linear interpolation method to evaluate the accuracy with which linear interpolation approximates the fully ray-traced spectra. The biggest error caused by the linear interpolation arises for spectra with a high-spin BH but low accretion rate. Those spectra are the most sensitive to spin and accretion rate, implying large changes in X-ray flux in response to small changes in parameters. As we can see from the lower panel of the left plot, the error in the interpolated luminosity density is <10% for all the cases considered over the 0.2–3 keV energy band. If we only consider the energy range over which the luminosity density decreases by a factor of <10⁴, the error is <5%.

Zoom In Zoom Out Reset image size

The right-hand side of Figure 1 shows the spectra and the fractional error caused by interpolation for disks with high accretion rates. As we can see from the spectra, the luminosity density decreases as the disk inclination becomes more edge-on. This is because the edge of the disk blocks X-rays emitted by the inner region of the disk, preventing them from reaching the observer (W20). The error can be up to 40% for an edge-on inclination. This effect was not seen in the left-hand panel of Figure 1, because of the smaller disk scale height at low Eddington ratio. Halving the distance between grid points for the inclination only reduces the error on the luminosity determination marginally (it is still around 40%). In the next section, we test how well high-inclination sources can be fit using the library of model spectra, taking ASASSN-15oi as a test case. We note that the discrepancy between the detailed model calculations and the linear interpolation calculations at high inclinations becomes smaller for larger M_• (>5.5 × 10⁶ M_⊙).

We have examined the flux errors caused by interpolating between library spectra and demonstrated that the procedure is generally accurate relative to fitting on the fly with our W21 ray-tracing model. In this section, we explore the effect on the M_• and a_• constraints from interpolating the library. We consider two cases: the TDE ASASSN-14li and the likely edge-on TDE ASASSN-15oi. We fit their X-ray spectra using both the library and the full W21 ray-tracing approach and compare the results for M_• and a_•.

In addition, we explore the constraints on parameters that were not investigated comprehensively in W20 (, θ, and N_H), because such a high-dimensional parameter sweep was computationally prohibitive using on-the-fly ray-tracing calculations. Now linear interpolation over the library of tabulated model spectra makes broader parameter surveys feasible. The fitting is done using XSPEC (Arnaud 1996).

For ASASSN-14li, we fit the 10 spectra with a slim-disk plus power-law model, where the power-law component, po, is only added to the fit function to describe the hard component detected during epoch 6. For the slim-disk component, we consider both the on-the-fly ray-tracing model of W21 and the library of pre-tabulated spectra from that model (slimd). We fit the spectra of the 10 epochs simultaneously, requiring the inclination θ, M_•, and a_• to be the same for all epochs (the priors are the same as listed in Table 1), but allowing the absorption (phabs) and accretion rate to vary between epochs.

Figure 2 compares the constraints on M_• and a_• obtained through fitting the observed spectra with our original on-the-fly ray-tracing model (W21) and with the model spectral library. In both analyses, we fit the same observed spectra assuming the same parameter priors. The 1σ and 2σ contours overlap closely, indicating that the results from the library and ray-tracing fits agree. We have achieved this strong agreement despite a dramatic reduction in computational cost: generating the confidence contours in Figure 2 takes ∼10⁴ CPU hours for the full, on-the-fly ray-tracing calculation, but only ∼10 CPU hours when interpolating over our pre-tabulated model library.

Zoom In Zoom Out Reset image size

For completeness, we also compare the M_• and a_• constraints from the ray-tracing models of W21 and W20 (which differ only by inclusion of angular momentum losses through radiation in W21). The two models give consistent results for the 2σ confidence levels (CLs), and W21 produces tighter constraints for the 1σ CL (see the Appendix).

Figure 3 shows the constraints on additional slim-disk fit parameters for ASASSN-14li obtained from using the library of model spectra. Due to the resulting computational efficiency, we can derive, for the first time, constraints on , N_H,i (in units of 10²² cm⁻²), and θ for each epoch while allowing M_• and a_• to float. The upper left panel shows the best-fit values and error bars for . We obtain a poor constraint (at 2σ CL) on the dimensionless mass accretion rate for most epochs, because the individual constraints on M_• and a_• are weak. As a result, spectra with different mass accretion rates and {a_•, M_•} combinations can fit the data equally well.

Zoom In Zoom Out Reset image size

The upper right panel shows the constraints on N_H for all epochs. Here, we constrain N_H well at most epochs, and we find a decline in N_H from epoch 1 to epoch 10. The fitted values of N_H in epochs 1 through 4 are larger than the value for the host galaxy plus the Galactic absorption combined (4 × 10²⁰ cm⁻²; Miller et al. 2015) at a 2σ CL, while the fitted values of N_H in epochs 5 through 10 are consistent with it at 1σ CL. We also plot the contours of N_H,1 versus N_H,10 in the lower right panel. From these contours, we find a decrease in N_H from epoch 1 to epoch 10 with a 2σ CL. This indicates that the additional gas liberated in the TDE might be responsible for the enhanced absorption at early times, and that it disappears at late times.

The lower left panel shows the joint constraint on a_• and θ. The 1σ contour is small, because we constrain the spin to be a_• > 0.95. θ is poorly constrained at 2σ CL due to the large range in allowed {a_•, M_•} values.

For the case of ASSASN-15oi, there are only two epochs of observations. We fit those two spectra with the slim disk plus a broken power-law function (bknpow in XSPEC). For the slim-disk component, as for ASSASN-14li, we use both the full on-the-fly ray-tracing model (W21) and the library of model spectra (slimd). The bknpow fit is used to describe the spectrum of the background photons; the bknpow parameter values are determined through separate fitting of the background spectra. During the fitting of the two source spectra, the parameters of the broken power law are kept fixed at their best-fit values. We note that the normalization of the broken power law is allowed to float (this is consistent with the procedure employed in W20). We use the floating normalization of the broken power law to account for the difference in area over which the background and source photons are extracted.

In W20, we considered a variety of scenarios to explain the unusual time evolution (Gezari et al. 2017) of ASASSN-15oi. Here we focus only on the “slimming disk” scenario of W20, so as to investigate the accuracy of parameter estimation with the library in the special case of nearly edge-on TDEs. In W20, we used the mass fallback law in Guillochon & Ramirez-Ruiz (2013) as a prior for this scenario. However, for our reanalysis using the library of model spectra, it is computationally feasible to keep as a free parameter. To remain close to the other priors used in W20, we allow to float, but set Edd and 70° < θ < 90°. With these assumptions, we fit the spectra with the ray-tracing model of W21 and the library of model spectra.

Figure 4 compares the constraints on M_• and a_• obtained through fits using the library of model spectra versus on-the-fly ray-tracing (W21). The 1σ and 2σ constraints from the two approaches overlap, demonstrating the benefit of fitting faster with the pre-tabulated library.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the constraints on additional slim-disk model parameters for ASASSN-15oi that we obtain from fitting the library model spectra to the two epochs of XMM-Newton data. As we did for ASASSN-14li, we constrain these parameters for the first time while allowing M_• and a_• to float. The best-fit values of the accretion rate are only loosely constrained, and a_• is unconstrained. The constraint on θ, which is well within the assumed prior range for θ, indicates that the disk is in fact nearly edge-on. The poor constraint on spin may be due to the X-ray photons from the inner part of the edge-on disk being blocked during the early epoch. For both epochs, there are strong constraints on N_H. We find, at >3σ confidence, the decrease in N_H over time also seen by W20; this behavior may be related to the source being viewed close to edge-on (Lu & Bonnerot 2020).

Here we have refit the events ASASSN-14li and ASASSN-15oi with the library of model spectra. We successfully reproduce the constraints on M_• and a_• obtained from the more computationally costly method of on-the-fly ray-tracing. In addition, we obtain new constraints on the , θ, and N_H parameters. This indicates that fits obtained by linearly interpolating over the library of model spectra are reliable representations of the results that would have been obtained had the models been calculated using the full ray-tracing code at all parameter values.

In this section, we investigate how the accuracy of estimated parameters—M_•, a_•, mass accretion rate, and disk inclination—is affected by the number of observing epochs and the number of X-ray photons detected. We also explore the impact of photoelectric absorption by including model phabs. In addition, we survey the parameter space to determine which types of disks yield stronger or weaker constraints on M_• and a_•. For instance, holding all else equal, does a highly super-Eddington mass accretion rate allow for a stronger or weaker constraint on M_• and a_•?

We use the model slimd and the XMM-Newton pn response file to generate mock observational spectra. These mock spectra have no background component, but they include statistical noise to better simulate the observations. The parameters used to generate the mock spectra are given in Table 2. We rebin the data such that a bin contains at least one X-ray photon (i.e., count). We use the library of model spectra to fit these mock spectra, employing Poisson statistics (Cash 1979).

Table 2. Input Parameters and Photon Number Counts for Mock Observational Spectra

Spectrum	Counts	NH [1022 cm−2]	[Edd]	θ [deg]	M• [M⊙]	a•	Exp. Time [ks]
1	46,393	0.03	3	45	106	0.9	3
2	46,676	0.03	3	...	106	0	16
3	46,022	0.03	3	...	104	0.9	55
4	47,697	0.03	0.2	...	...	...	200
5	49,621	0.03	1	45	106	0.9	9
6	47,981	0.03	0.9	...	...	...	10
7	48,328	0.03	0.8	...	...	...	12
8	45,978	0.03	0.7	...	...	...	14
9	48,466	0.03	0.6	...	...	...	19
10	46,089	0.03	0.5	...	...	...	25
11	47,314	0.03	0.4	...	...	...	40
12	47,827	0.03	0.3	...	...	...	75
13	192,620	0.03	1	...	...	...	35
14	18,987	0.03	1	45	104	0.9	50
15	17,242	0.03	0.9	...	...	...	50
16	15,500	0.03	0.8	...	...	...	50
17	13,533	0.03	0.7	...	...	...	50
18	3032	0.03	10	80	106	0.9	2
19	20,769	0.03	2	...	...	...	4.2
20	21,342	0.03	1	...	...	...	9
21	20,684	0.03	0.9	...	...	...	10
22	21,085	0.03	0.8	...	...	...	12
23	20,134	0.03	0.7	...	...	...	14

Note. We assume an absorbed slim-disk model, phabs*slimd, to generate the mock spectra using the XMM-Newton pn response file. We adopt D_L = 100 Mpc (corresponding to z = 0.023) for all spectra, and we use grppha to rebin the data so that there is at least one count per spectral energy bin.

Download table as:  ASCIITypeset image

We first investigate whether there are parameter degeneracies using a single-epoch mock spectrum. Then, we explore the effect of the total number of counts. Moreover, we show that multiepoch spectra help to break degeneracies between parameters. We also demonstrate that the early-time epochs of edge-on TDEs have important constraining power.

5.1. Parameter Degeneracies

Figure 6 shows the best-fit confidence contours for different parameters derived from library fits to different mock observed spectra. Here, we consider four mock spectra with different combinations of , a_• and M_•, sharing a common source inclination of θ = 45°. ¹¹ The input parameters for these mock spectra can be found in the first four rows of Table 2 and are shown below each of the figure panels. For each mock spectrum, we adjust the exposure time so that all have a similar number of counts in the 0.3–10 keV band.

Zoom In Zoom Out Reset image size

The four top left panels show the fits to a super-Eddington, high-a_• spin, and relatively high-M_• source spectrum. The M_• and a_• contours show that M_• and a_• are well constrained, indicating that the spectral shape is sensitive to these two parameters. From the other contours, we find that and θ are constrained poorly, and that both are highly degenerate with other disk parameters. The 2σ range of extends to values larger than 10 (in Eddington units), due to the fact that the X-ray luminosity is weakly sensitive to accretion rate for highly super-Eddington slim disks (Abramowicz et al. 1988; W20). The confidence contours of θ extend up to ∼80°, but not higher, because an even higher inclination would cause the accretion disk to block most X-rays coming from the inner disk, reducing the observable flux.

The four top right panels show results for the same input parameters as the four top left panels except for accretion rate, which has now been set to the sub-Eddington value of . From the M_• and a_• contours, we see that the model fits still constrain M_•, but no longer a_•. The low accretion rate makes the disk cooler, moving the spectrum to lower X-ray energies, which weakens the constraint on a_•. The constraints on and θ become even weaker when compared with the case of the super-Eddington spectrum. In order to constrain a_• using a single epoch of data when the mass accretion rate is sub-Eddington, one needs to know the value of M_• so that it can be held fixed in the fit.

The four bottom right panels show the results for a mildly super-Eddington disk with a_• = 0. The constraints on all parameters are weaker than in the high-spin case shown in the upper left panels. Similar to the sub-Eddington mass accretion case, a lower spin makes the disk cooler, resulting in a softer spectrum and poorer parameter estimates.

The four bottom left panels show the fit results for an IMBH: M_• is reduced to 10⁴ M_⊙, with and a_• = 0.9 as in the top left panels. The confidence contours on each parameter are the tightest for all of the four cases shown in this figure. This is because a smaller M_• results in a hotter disk (roughly, ), meaning that a larger portion of the source spectrum falls in the X-ray band. We note that for such a low M_•, the Eddington luminosity is 100 times lower than for the other three panels (where M_• = 10⁶ M_⊙), implying that the system has to be ≈10 times ¹² closer to yield the same number of X-ray counts for fixed observing time, as we have assumed in making these comparisons.

In all of these efforts to fit mock spectra with the library of model spectra, the input values of all parameters (M_•, a_•, , θ) are recovered within the 1σ confidence contours. In all cases, we find that: (1) M_• can be constrained relatively well; (2) both and θ are poorly constrained, as they are highly degenerate with other parameters. The spin a_• can be well constrained for cases with a high BH spin, super-Eddington accretion rates, or a BH in the intermediate-mass regime (although, in the IMBH case, it is more observationally challenging to obtain a large number of counts). As a rough heuristic summarizing these results, we find stronger constraints on all parameters when the information-rich peak of the quasi-thermal spectrum lies in (or at least closer to) the observing band, and weaker constraints when the observing band only contains the information-poor Wien tail.

5.2. Effect of X-Ray Photon Counts

In this section, we explore the dependence of our joint {M_•, a_•} constraints on the total number of X-ray counts. In order to quantify this dependence, we use a figure of merit (FoM). Following the Dark Energy Task Force team (Andreas et al. 2006), we define the FoM as the reciprocal of the area of the error ellipse in the plane. The FoM can be calculated as

where is the covariance matrix of and a_•. A larger (smaller) FoM indicates a better (weaker) constraint on and a_•.

In order to explore the relationship between the FoM and the number of counts in the spectra, we consider four different sets of spectra: super-Eddington IMBH (M_• = 10⁴ M_⊙, ), super-Eddington SMBH (M_• = 10⁶ M_⊙, ), sub-Eddington SMBH (M_• = 10⁶ M_⊙, ), and low-spin SMBH (M_• = 10⁶ M_⊙, ). We use the slimd model to generate 11 mock spectra for each set, each with the same a_• = 0.9 (except for the low-spin case) and θ = 45°. Each set undergoes mock observations with a range of different exposure times ¹³ that are chosen for equal numbers of counts among the four sets. We fit these spectra using the slimd model, allowing all the parameters to vary freely. We do not get the covariance matrix directly from the fitting result, as it is not accurately determined in XSPEC, which only estimates the Fisher matrix (Arnaud 1996). However, we run the steppar command for the parameters M_• and a_•, and then calculate the covariance matrix with the likelihood function (derived from Cstat) and the corresponding values of M_• and a_•.

Figure 7 shows the FoM versus the number of X-ray counts for the four sets of spectra that we simulated. In all cases, a power-law function fits the data well. The power-law index varies depending on the underlying parameters: it is highest for the IMBH case and lowest for the sub-Eddington SMBH case. Interestingly, the normalization of the IMBH case’s power-law is lower than that of the super-Eddington SMBH case, so the FoM of the IMBH only bests the super-Eddington SMBH once the number of counts exceeds several thousand. In summary, the FoM power-law indices fitted here quantify the trend noted heuristically in Section 5.1: hotter disks generally provide more constraining power per X-ray photon.

Zoom In Zoom Out Reset image size

5.3. Constraints from Fits to Multiepoch Data

Unlike in most active galactic nuclei (AGNs), the accretion rate varies strongly with time in TDE disks. Therefore, multiepoch X-ray spectroscopy has the potential to obtain more stringent constraints on the best-fit parameters for IMBH/SMBH accretion than it does for typical AGNs. In this section, we explore how multiepoch spectroscopy influences parameter estimation. We first examine the constraints obtained for a source with a low or high accretion rate. In order to do so, we generate nine epochs of mock spectroscopic data for SMBH TDEs, using phabs*slimd as an input model (see rows 5–13 in Table 2). Next, we study the case where multiepoch spectra are obtained for IMBH TDEs. We create four mock spectra with an exposure time of 50 ks (rows 14–17 in Table 2). In this analysis, when fitting the simulated spectra, we fix the value of N_H as a prior and do not allow it to float or vary in time.

Figure 8 shows the constraints on M_• and a_• obtained when jointly fitting multiepoch spectra together. For three different sets of spectra, we simultaneously fit all epochs with the model phabs*slimd. We require the same M_•, a_•, and θ within each set of spectra, because those parameters are not expected to change between observations. The figure shows that a fit to multiepoch spectra obtained when the mass accretion rate is nearly super-Eddington (upper left panel) yields tighter constraints ¹⁴ on both M_• and a_• than do multiepoch spectral fits to only sub-Eddington epochs (upper right panel). The explanation for these tighter constraints is twofold. First, a higher accretion rate yields a harder spectrum, increasing the fraction of X-ray emission that falls in the observed band and thus also boosting its constraining power. Second, a higher accretion rate pushes the effective temperature maximum of the disk closer to the ISCO, where the influence of the BH spin has a stronger effect on the emergent X-ray emission.

Zoom In Zoom Out Reset image size

In our past work (W20, W21), fits to multiepoch spectra constrained both M_• and a_• better than fits to individual, single-epoch spectra did. However, as multiepoch spectra usually have more counts in total, we need to test whether the multiepoch nature of multiple observations is intrinsically valuable for breaking parameter degeneracies. To do this, we have generated a single mock spectrum with and set the exposure time such that the spectrum has the same number of counts as the four spectra combined in each of Figure 8's two upper panels. The lower panel of Figure 8 shows the M_• and a_• constraints from this single-epoch spectrum. It is clear both by eye and by FoM analysis that this single-epoch, trans-Eddington spectrum provides weaker constraints than multiepoch spectra including at least one trans-Eddington epoch (top left panel), but stronger constraints than multiepoch spectra that only span sub-Eddington epochs.

The lower right panel shows the constraints on M_• and a_• from the fitting of the four epochs of simulated IMBH TDE spectra. As we showed in Section 5.1, the model fits to a single epoch of a super-Eddington IMBH spectrum can already constrain both M_• and a_• reasonably well. We generate four spectra with a 50 ks exposure time each (see rows 14–17 in Table 2). The number of counts for each epoch is about 1/2 to 1/3 that in row 3 of Table 2 (our previous IMBH analysis, where the Eddington ratio was higher) and also in rows 5–12 (our previous multiepoch analysis for SMBHs, which are usually brighter at fixed ). We retrieve the input M_• and a_• values to within 1σ CL, and the CL contours are small relative to the other panels. The constraints would be even better for slim-disk fits to more highly super-Eddington accretion spectra, a situation that arises frequently for IMBH TDEs.

5.4. Effect of Edge-on Disk

In this section, we explore how an edge-on inclination (θ) of the inner disk with respect to our line of sight affects the constraints from model fits to multiepoch spectra. When the accretion rate decreases, the disk height is reduced and more X-rays from the inner disk have a free line of sight toward the observer. We generate six mock SMBH spectra with θ = 80° and ranging from 10 to 0.7. (See Table 2, rows 18 to 23, for the values of other parameters.) The 0.3–10 keV X-ray flux peaks for disks with . For , the disk aspect ratio increases, and the X-ray emission of the inner disk is increasingly blocked by the disk edge. For , the sky-integrated luminosity in the 0.3–10 keV band decreases quickly with decreasing mass accretion rate.

Figure 9 shows the constraints on M_• and a_• obtained by fitting our model spectral library to the mock observational spectra. The upper left panel shows the constraints on M_• and a_• from fitting a single spectrum with . M_• is constrained to within a factor of ≈10 only, and there is no constraint on a_•. The immediate reason for this is the very low number of X-ray counts (see Table 2), but the deeper explanation is that the X-ray photons from the inner disk are almost entirely blocked by the geometrically thick disk and are unable to reach the observer.

Zoom In Zoom Out Reset image size

The upper right panel shows the constraints on M_• and a_• obtained by fitting multiepoch spectra. In this case, the multiepoch fit to five spectra includes two with high accretion rates (), but not so high as to heavily obscure the inner, X-ray-emitting disk annuli. As a result, both M_• and a_• are strongly constrained.

The lower panel shows the result combining all six spectra. Interestingly, the addition of just one more spectrum improves the joint constraint on M_• and a_• by a factor of ≈2.5 compared to the top right panel. This improvement arises because the additional, spectrum provides a strong constraint on the disk inclination similar to that shown in the lower right panel of Figure 5.

5.5. Effect of N_H

In this section, we quantify the importance of absorption on the X-ray spectrum. As the X-ray emission may be subject to absorption, the total amount of which is not known a priori, one needs to fit the absorption together with the other TDE parameters. The absorbing column N_H may be portions of the interstellar medium located in the Milky Way or in the TDE host galaxy, in which case N_H should be the same from epoch to epoch. Alternatively, a large fraction of N_H could arise from local gas produced in the TDE itself, in which case N_H would be expected to vary in time, likely experiencing an eventual decline. This locally produced gas could include bound but poorly circularized debris (Loeb & Ulmer 1997; Guillochon et al. 2014; Roth et al. 2016) or outflowing ¹⁵ material (Metzger & Stone 2016; Dai et al. 2018; Roth & Kasen 2018; Lu & Bonnerot 2020; Piro & Lu 2020).

Here we consider the effect of absorption on the best-fit parameters using fits both to a single-epoch spectrum and to multiepoch spectra. As before, we define the mock spectrum with an absorbing column of N_H = 0.03 × 10²² cm⁻², but we now allow N_H to float in our fits, rather than specifying it as a prior. Figure 10 shows the effect of absorption on the parameter constraints from fits to single-epoch spectra (in one case super-Eddington, in the other sub-Eddington). From the contours of , we see that there is a degeneracy between N_H and M_•, a_•, and . As a result, we find poorer constraints on M_•, a_• and compared with the fit where N_H is known a priori (see plots with in Figure 6). However, the best-fit value of N_H is constrained to 0.035 ± 0.01 (10²² cm⁻²), and thus the input value is regained at a 1σ CL. From the CL contours on , we see that this degeneracy (of N_H with M_•, a_•, and ) remains for sub-Eddington spectra, weakening constraints on parameters of interest. While the uncertainty in N_H can affect constraints on other disk parameters, N_H itself can be well fit for both super- and sub-Eddington cases, at least for the values we have explored here.

Zoom In Zoom Out Reset image size

Figure 11 shows the effect of allowing the absorption parameter to float freely when fitting multiepoch data. The left panel shows the constraints on M_• and a_• with fixed N_H,i , while the right panel shows the M_• and a_• constraints from fitting the same spectra but allowing N_H,i to float freely. We see that the additional parameters N_H,i make the constraints on M_• and a_• significantly worse at both the 1σ and 2σ levels. The substantial degeneracies of absorption with parameters of (generally) greater interest may motivate other approaches for constraining this “nuisance parameter,” such as fitting X-ray absorption lines (Miller et al. 2015), although this goes beyond the scope of the present work.

Zoom In Zoom Out Reset image size

X-ray observations of TDEs have a rich history. They are not only the first means by which TDEs were detected (e.g., Bade et al. 1996; Komossa & Bade 1999), but by probing regions of spacetime close to the event horizon of the black hole, they have the potential to measure fundamental parameters such as black hole mass M_• and spin a_•. Our recent work (W20, W21) aimed to systematize the process of parameter estimation in TDE X-ray spectra, but one key problem was computational expense. Applying full, on-the-fly general relativistic ray-tracing calculations to estimate confidence contours for the {M_•, a_•} plane can take over ∼10⁴ CPU hours. In this paper, we have addressed this problem by constructing a pre-tabulated library of slim-disk model spectra that can be used by the community to fit soft X-ray observations of TDEs over only ∼10 CPU hours.

We use the model of W21 to construct this tabulated model. First, we solve the general relativistic stationary slim-disk model, including the effect of angular momentum loss by radiation. Second, we use the relativistic ray-tracing code to calculate synthetic spectra. In these two steps, we have included gravitational redshift, Doppler, and lensing effects self-consistently. When calculating the spectrum, we determine the spectral hardening factor f_c using Equation (2), and set the inner and outer disk edges to the ISCO and 2R_t, respectively (except for very high-mass SMBHs). There are four free parameters for the model: M_•, a_•, , and the inclination θ. Their ranges are listed in Table 1. We employ linear interpolation to calculate the spectra with parameters between the grid points. We evaluate the error on the calculated fluxes caused by the interpolation.

We refit the X-ray spectra of ASASSN-14li and ASASSN-15oi with the library of model spectra and compare the results with those presented in W20. Using simulated mock X-ray spectra, we explore the full range of degeneracies between model parameters. We further evaluate how multiepoch X-ray spectra, an edge-on disk configuration, and absorption can affect the constraints on parameters. Our conclusions are:

• 1.  
  The error on the flux in the 0.2–10 keV band caused by linear interpolation is typically <5%. The error increases for higher-energy photons and can rise to 10% far out on the Wien tail, when the flux is a factor of 10⁶ lower than at the peak of the X-ray spectrum. The error can be even larger, up to 40%, for an edge-on, high disk configuration (see Figure 1). The constraints on M_• and a_• for ASASSN-14li and ASASSN-15oi obtained using the library of model spectra to fit the data are consistent with those from the full ray-tracing model (W20). We also confirm the decrease found by W20 in N_H with time at the >2σ CL for ASASSN-14li and at >3σ for ASASSN-15oi.
• 2.  
  Both super- and sub-Eddington spectra can constrain M_•. High spin, low BH mass, and/or highly super-Eddington spectra allow a stronger constraining power on {M_•, a_•} than low spin, high BH mass, and/or sub-Eddington spectra (see Figure 6). Typically, a single-epoch spectrum at a highly super-Eddington mass accretion rate can constrain {M_•, a_•} well.
• 3.  
  The constraining power of M_• and a_• increases as a power law with the number of spectral counts. The index of this power law is higher for high , higher a_•, and low M_•.
• 4.  
  Multiepoch spectra can help to break degeneracies in parameter estimation (see Figure 8). Nevertheless, multiple epochs of observations of sub-Eddington mass accretion rate provide weaker constraints on {M_•, a_•} than do super-Eddington epochs with similar total counts.
• 5.  
  Fitting a single spectrum of a system with edge-on inclination constrains {M_•, a_•} poorly, but simultaneous, multiepoch spectral fitting can place much tighter constraints (see Figure 9).
• 6.  
  The level of absorption N_H can be recovered from fitting spectra with either super-Eddington or sub-Eddington mass accretion rate. However, the uncertainty in N_H adversely affects the constraints obtained on M_• and a_•.

Our relativistic slim-disk library provides a systematic, computationally efficient framework for fitting and interpreting the growing wealth of observational data. As noted in our previous papers (W20, W21), our model disks are inherently idealized, because complete, first-principles models for the hydrodynamic evolution of TDEs do not yet exist (see, e.g., Bonnerot & Stone 2021 for a recent review). Most notably, at early times, TDE disks may possess substantial eccentricity or tilt that we currently ignore. Likewise, the slim-disk solutions could break down for highly super-Eddington epochs.

These potential problems are mitigated for sufficiently late-time observations, once the inner disk has fully circularized, internal and external torques have aligned it with the black hole equatorial plane, and the accretion rate has become sub-Eddington. These issues may also be reduced in certain regions of parameter space: higher spin values and more massive SMBHs likely favor fast alignment of initially tilted disks (Franchini et al. 2016), and deeply plunging TDEs may experience more rapid circularization (Andalman et al. 2022). If rapid circularization is a prerequisite for high X-ray luminosity, as has been suggested in past work (Dai et al. 2015), then residual disk eccentricity may not be a problem for most X-ray-bright TDEs. As the theory of TDEs progresses, the library framework we have constructed can be extended to incorporate more complex aspects of TDE disk physics and thereby eliminate some of these systematic uncertainties.

The pre-tabulated model library presented here largely reduces the computational cost of fitting single- and multiepoch TDE X-ray spectra with the slim-disk model. We anticipate that these publicly available tools will be more and more useful to the community as it transitions to an era of population studies on X-ray-bright TDEs. The eROSITA telescope has already found dozens of new X-ray TDEs (Sazonov et al. 2021) and is likely to find hundreds if not thousands more before its mission is over (Khabibullin et al. 2014; Jonker et al. 2020). The planned Einstein Probe mission has similar potential that may be realized this decade (Yuan et al. 2015). This near-future wealth of X-ray TDEs will carry unprecedented information on the demographics of {M_•, a_•} for massive black holes in quiescent galaxies, which can be explored efficiently with the framework we have introduced here.

We thank Chi-Kwan Chan for his help in accessing computational resources and in constructing the model. We also thank Keith Arnaud for his guidance for adjusting our model to the XSPEC table model. S.W. and A.I.Z. thank Steward Observatory and the UA Department of Astronomy for postdoctoral support for S.W. S.W. also thanks RU Department of Astrophysics/IMAPP for postdoctoral support. A.I.Z. acknowledges additional support from NASA ADAP grant #80NSSC21K0988. N.C.S. received support from the Israel Science Foundation (Individual Research grant 2565/19). Our disk calculations were carried out at UA HPC systems, which are supported by the National Science Foundation under grant No. 1228509. Our calculations of spectral energy distribution were carried out at Open Science Grid (Pordes et al. 2007; Sfiligoi et al. 2009), which is supported by the National Science Foundation award #2030508. Discussions during the 2020 Yukawa Institute for Theoretical Physics workshop on “Tidal Disruption Events: General Relativistic Transients” at Kyoto University were important to the completion of this work.

In W21, we updated the slim-disk solution by including the effect of angular momentum loss by radiation (Abramowicz et al. 1996). We found that this effect becomes strong for cases of high spin and low accretion rates, and can be ignored for low spin or high accretion rates. Here, we compare the constraints on M_• and a_• between the slim-disk models of W21 and W20, i.e., with and without the inclusion of angular momentum loss through radiation, respectively. We use the same spectra and parameter priors as in W20.

Figure 12 shows the constraints on M_• and a_• with the slim-disk models of W21 and W20. Using the W21 model, we recover the best-fit and the 2σ contour on the uncertainties for M_• and a_• from W20. However, the 1σ contour obtained from W21 is much smaller than that from W20. This difference arises because the χ² value of the best fit is smaller by about 1.0 for W21. The decrease in χ² can be understood as due to the high best-fit value of the BH spin making the inclusion of angular momentum loss through radiation important.

Zoom In Zoom Out Reset image size