Introduction

Active galactic nuclei (AGN) are thought to be powered by the accretion of matter around supermassive black holes at the center of galaxies (E. E. Salpeter 1964;Y. B. Zel'dovich 1964). Quasars are bright AGN with unobscured accretion disks, and are some of the brightest objects in the Universe. They are observable at extreme cosmological distances, making them powerful probes of the early Universe (D. J. Mortlock et al. 2011;E. Bañados et al. 2018). Quasars are also thought to play an important role in galaxy evolution (A. Franceschini et al. 1999;G. Kauffmann & M. Haehnelt 2000;A. Hoshi et al. 2024). The variability of quasar brightnesses has been studied since their discovery (J. L. Greenstein 1963;C. Hazard et al. 1963;T. A. Matthews & A. R. Sandage 1963;J. B. Oke 1963;M. Schmidt 1963). These brightness variations are related to the physical properties of the black hole and accretion disk that power them. Furthermore, investigating the physics governing quasar light curves offers valuable insights into our understanding of cosmology (N. Khadka & B. Ratra 2020;B. Czerny et al. 2023).

The variability of quasars is most often modeled as an X-ray driving variability source corona above the black hole that illuminates a geometrically thin, optically thick accretion disk (N. I. Shakura & R. A. Sunyaev 1973;E. M. Cackett et al. 2007). The reprocessing of the driving variability to the UV/optical emitting regions of the accretion disk is represented by the so-called transfer functions, and introduces wavelength-dependent time lags (R. D. Blandford & C. F. McKee 1982), ranging from less than a day for small supermassive black holes (∼10 7 M ⊙ ) to several tens of days (∼10 10 M ⊙ ). These time lags can be measured through continuum reverberation mapping of UV/optical light curves to probe the relative size scales of the emitted regions, and are interconnected to properties of the accretion disk and black hole (E. M. Cackett et al. 2021;V. K. Jha et al. 2022; S. Wang et al. 2023).

Studies of quasar variability have found accretion disk sizes to be larger than the predictions of the standard thin-disk model (N. I. Shakura & R. A. Sunyaev 1973) by a factor of ∼2-4 (D. Mudd et al. 2018; W.-J. Guo et al. 2022; V. K. Jha et al. 2022). This is also consistent with accretion disk size measurements found using gravitational microlensing (S. Poindexter et al. 2008;S. Poindexter & C. S. Kochanek 2010;J. A. Blackburne et al. 2015;J. A. Muñoz et al. 2016; C. W. Morgan et al. 2018). New measurements and methods are needed to test accretion disk models and enhance our understanding of the physical processes governing quasar emissions.

Upcoming wide-field surveys such as the Rubin Observatory Legacy Survey of Space and Time (LSST) will observe an unprecedented quantity of data. The LSST main survey will cover 18,000 deg 2 and is projected to monitor tens of millions of quasars over a 10 yr period with six UV/ optical bandpass filters (ugrizy) at 55-185 samplings per band, or around 800 total visits across the 10 yr out to a redshift of z ∼ 7.5 (LSST Science Collaboration et al. 2009;A. Prša et al. 2023). A smaller sky area of 200 deg 2 known as the Deep Drilling Fields is expected to detect 40,000 additional ultrafaint AGN at higher cadences of about 1000 samplings per band (W. Brandt et al. 2018). Machine learning (ML) algorithms are well suited to analyze the vast amounts of data expected from LSST and other wide-field surveys. The quasar light curves from LSST pose challenges for traditional ML techniques due to being stochastic, multivariate, noisy, and irregularly sampled across bands with long seasonal gaps (see S. N. Shukla & B. M. Marlin 2021, for a review of ML methods on irregularly sampled time series).

UV/optical variability light curves are most commonly modeled using Gaussian process regression (GPR). In this framework, quasar light curves are fit with a specific kernel, e.g., the kernel associated with the damped random walk (DRW) process (Y. Zu et al. 2013). The optimized kernel parameters have been empirically shown to be related to properties of the quasar, such as the black hole mass (C. L. MacLeod et al. 2010; K. L. Suberlak et al. 2021).

The codebase JAVELIN (Y. Zu et al. 2011Zu et al. , 2016; M. M. Fausnaugh et al. 2016) uses a DRW kernel with tophat transfer functions to simultaneously model the variability and time delays. The bluest band is taken as the effective driving variability source and fit using the DRW, and the tophat transfer functions are used to measure the time delays between the other bands. This method relies on Markov Chain Monte Carlo (MCMC) to sample the DRW kernel and transfer function parameters. By default, JAVELIN measures the time delays but does not directly extract the physical properties of the quasar. Simple parametric models can be fit to the timedelay measurements as a secondary step to measure the accretion disk size and sometimes the temperature profile, although it is most frequently fixed to the thin-disk case. In D. Mudd et al. (2018), JAVELIN was modified to measure the accretion disk size directly by fixing the time delays to the thin-disk model.

The codebase CREAM (D. A. Starkey et al. 2015) is similar to JAVELIN but uses the thin-disk transfer functions directly instead of top hats. Instead of treating the bluest band as the effective driving variability, CREAM explicitly reconstructs the driving variability by modeling it as a Fourier series. MCMC is used to optimize the Fourier components of the driving variability and the accretion disk parameters, and both the driving variability and transfer function kernels are reconstructed. In theory, CREAM can recover the product of the black hole mass and accretion rate, MM , and the inclination angle, although the inclination cannot be recovered without very-high-fidelity data. For example, it was used in M. M. Fausnaugh et al. (2018) to model two Seyfert 1 galaxies to measure MM and to constrain the inclination for one galaxy.

Currently, JAVELIN and CREAM are the two most advanced methods of measuring time delays and accretion disk sizes in continuum reverberation mapping. However, both of them rely on computationally expensive MCMC sampling, which may be infeasible to apply to the tens of millions of quasar light curves expected from LSST. In addition, the DRW kernel in JAVELIN is only fit with respect to the bluest band; however, the information in each band is highly correlated to the other bands, so jointly modeling the light curves and time delays would improve the performance. This is especially the case for sparsely and irregularly sampled cadences like those expected from LSST. UV/optical light curves have also been shown to significantly deviate from a DRW (W. Yu et al. 2022), so a more flexible fitting method would be beneficial. CREAM does jointly model the driving variability and transfer functions; however, the driving variability is modeled as a Fourier series. This could be challenging for long light curves like what we expect from LSST, since the driving signal is expected to be stochastic across many timescales. Ideally, our methods would also be able to learn from features across an entire sample of light curves. For example, there is additional information in the mean brightnesses of each band that could be combined with time-delay measurements to better estimate accretion disk parameters and break some of the degeneracies. Our ML approach improves upon all these areas.

ML methods have been used to model quasar variability in several different ways. Y. Tachibana et al. (2020) used a recurrent auto-encoder to model quasar variability, and found it to perform better than a DRW model when applied to real data. P. Sánchez-Sáez et al. (2021) used a recurrent variational auto-encoder for anomaly detection to find changing-look AGN. J. W. Park et al. (2021) introduced a method to simultaneously reconstruct quasar light curves and predict accretion disk parameters using attentive neural processes. I. Čvorović Hajdinjak et al. (2022) introduced conditional neural processes to model quasar variability. X. Sheng et al. (2022) applied stochastic recurrent neural networks (RNNs) to reconstruct simulated LSST light curves. E. Danilov et al. (2022) developed a neural inference Gaussian processes method to fit quasar light curves. J. I.-H. Li et al. (2024) demonstrated how simulation-based inference could be used to predict accretion disk parameters. For microlensed quasar light curves, G. Vernardos & G. Tsagkatakis (2019) used a convolutional neural network (NN) to measure the accretion disk size and temperature profile in simulated light curves, and H. Best et al. (2024) predicted the black hole mass, inclination angle, and impact angle. J. Fagin et al. (2025) developed a method of predicting high-magnification microlensing events through real-time classification with simulated LSST light curves using an RNN.

J. Fagin et al. (2024) introduced latent stochastic differential equations (SDEs) as a method to reconstruct simulated LSSTlike quasar light curves and simultaneously predict the accretion disk and variability parameters. Latent SDEs are a type of generative neural network that can model continuoustime stochastic dynamics (X. Li et al. 2020). They are physically motivated by the fact that quasar light curves are generally well described by SDEs such as the DRW or higherorder continuous-time autoregressive moving-average (CARMA) processes (W. Yu et al. 2022). Latent SDEs can be viewed as infinite-dimensional variational autoencoders (D. P. Kingma & M. Welling 2013;D. J. Rezende et al. 2014) with an SDE-induced process as their latent state. J. Fagin et al. (2024) found the deep learning method to be superior to a multitask GPR baseline in reconstructing LSST light curves. In addition, their model simultaneously performed parameter inference based on the context vector of the encoder and the latent vector of the SDE. They were able to predict the black hole mass, temperature slope, and inclination angle, as well as the parameters of the DRW driving variability signal.

While the latent SDE method of J. Fagin et al. (2024) is able to simultaneously reconstruct the light curve and perform parameter inference, the reconstructed light curve and parameter predictions do not necessarily correspond to the same time delays. In this work, we are the first to combine the light-curve reconstruction and parameter inference into a selfconsistent, unified framework. This is achieved by developing the first auto-differentiable simulation of the accretion disk and including it into the architecture of our ML model. Within our NN, we use a latent SDE to generate the X-ray driving variability. We then predict the accretion disk parameters, which are converted to the corresponding transfer functions using our auto-differentiable simulation of the disk. The reconstructed driving variability is convolved with the transfer functions and then scaled to produce the mean best-fit reconstruction of each observed UV/optical band, along with their uncertainties. In addition, we predict the variability parameters of the driving signal and their uncertainties. The relative time delays between bands and their uncertainties are also predicted from the mean time delay of the reconstructed transfer functions.

A recurrent inference machine (RIM; P. Putzky & M. Welling 2017) is a technique that iteratively refines its predictions by feeding previous outputs back into the model. It has been used in several astrophysics applications (W. R. Morningstar et al. 2019;C. Modi et al. 2021;A. Adam et al. 2023;C. Rhea et al. 2023). Our goal in this work is to solve the blind deconvolution inverse problem of recovering the driving signal and transfer functions given the observations of our light curve. This is made particularly challenging given the stochastic nature of quasar variability and the fact that our observations are irregular and sparsely sampled. We use RIM in our NN to iteratively improve upon the light-curve reconstruction and accretion disk parameter estimation.

In Section 2, we describe how we build a realistic simulation of quasar light curves including an auto-differentiable version that is incorporated into our ML model. In Section 3, we present our NN architecture and training. In Section 4, we give our results on our model's performance in light-curve reconstruction and parameter inference. In Section 5, we discuss the results of our ML model and future prospects, and Section 6 gives our concluding remarks. Throughout this work, we assume a flat ΛCDM cosmology with H 0 = 70 km s -1 Mpc -1 , Ω m = 0.3, and Ω Λ = 0.7.

Light-curve Simulation

We train our ML model with realistic simulations of LSST 10 yr light curves. The reprocessing of the X-ray driving variability by the accretion disk to the UV/optical wavelength λ is modeled by

is the mean flux, ΔF λ (λ) is the amplitude of the variable flux, X(t) is the normalized driving variability (mean zero and variance one), and ψ(τ|λ) is the transfer function kernel (E. M. Cackett et al. 2007;D. A. Starkey et al. 2015;J. H.-H. Chan et al. 2025). Section 2.1 describes how we model the driving variability X (t). Section 2.2 introduces our auto-differentiable model of the accretion disk reprocessing and transfer functions ψ(τ|λ). In Section 2.3, we construct realistic spectra from our accretion disk model using templates for the spectral lines, host-galaxy flux, and extinction. We then integrate the spectrum across the filter response functions of each LSST band to get the mean flux and variability amplitude of each band, ¯( ) F and ΔF λ (λ). Section 2.4 gives the parameter ranges for building our training set, and Section 2.5 describes how the time series is degraded to mimic LSST observational cadences and noise.

Quasar Driving Variability Model

UV/optical quasar light curves are often modeled as a DRW, a type of Gaussian process also known as the Ornstein-Uhlenbeck process (C. E. Rasmussen & C. K. I. Williams 2006;Y. Zu et al. 2013). A DRW signal X(t) is governed by the SDE,

where ε(t) is a white-noise process with a mean of zero and variance of one, τ is the characteristic timescale, b is related to the mean of the process = X b , and σ is related to the standard deviation defined by the asymptotic structure function / = SF 2 (B. C. Kelly et al. 2009). The DRW can alternatively be characterized by its power spectral density (PSD), given by

for frequency ν. This model is useful because GPR can be used to measure the variability parameters τ and SF ∞ using the kernel of the Gaussian process:

SF exp , 4 2 where Δt is the time separation of two observations. The measured values of τ and SF ∞ have been empirically shown to relate to properties of the accretion disk such as the black hole mass (C. L. MacLeod et al. 2010; K. L. Suberlak et al. 2021).

More complex Gaussian processes have also been used, such as higher-order CARMA processes (W. Yu et al. 2022). These CARMA processes are the solution to higher-order SDEs than the DRW given in Equation (2). The X-ray driving variability has been empirically shown to be better modeled by a broken power-law (BPL) PSD than a DRW. We generate our simulated driving signal from a bended BPL PSD, given by ( ) ( )

where ν b is the break frequency between the lower power-law slope α L and the higher power-law slope α H . This PSD implies that ( ) P L at low frequencies when f ≪ ν b and ( ) P H when f ≫ ν b . A form of the bended BPL has been extensively used to model the X-ray variability (e.g., I. M. McHardy et al. 2004; P. M. O'Neill et al. 2005; P. Uttley & I. M. McHardy 2005; D. P. Summons 2007; A. Markowitz 2010; K. L. Smith et al. 2018; L. F. Sartori et al. 2019; H. Yang et al. 2022; B. Czerny et al. 2023;H. Yuk et al. 2025). Although the bended BPL is not a Gaussian process, its parameters can be measured by directly fitting the PSD. In the special case where α L = 0, α H = 2, and ν b = 1/(2πτ), the bended BPL recovers the PSD of the DRW in Equation (3). Typically, the BPL for X-ray variability has been measured with α L ∼ 1 and α H ∼ 3. Similarly to the DRW, the parameters of the BPL have been shown to relate to properties of the accretion disk and black hole (I. M. McHardy et al. 2005;P. Arevalo et al. 2024). For example, the break frequency decreases with an increased black hole mass. In this work, we do not correlate the variability parameters to the accretion disk parameters, since we are interested in measuring these correlations in data without bias.

We can generate light curves from any PSD using the method of J. Timmer & M. Koenig (1995) by taking the inverse Fourier transform with randomized complex phases. This method has been used to generate X-ray variability with the bended BPL (e.g., B. Czerny et al. 2023). We generate the driving variability at a length 8 times larger than our desired signal, and keep only the initial desired length. This is to avoid the boundary conditions of the Fourier transform, and so we can set the asymptotic mean and standard deviation to the longer time series without bias (described in more detail in Section 2.3). Generating the driving signal with 8 times the target length should be more than sufficient to achieve these goals. The driving variability is generated at daily intervals and in magnitude, so the flux is always positive.

The power spectrum in Equation ( 5) is in the rest frame of the quasar. We predict the break frequency in the observer frame to avoid degeneracies with the redshift. The posteriors of our break frequency and redshift can be combined after the fact to obtain the rest-frame break frequency:

Auto-differentiable Accretion Disk Model

Our goal is to use the accretion disk model to simulate our training set and as a direct component of our NN. We implement the accretion disk and transfer functions in PyTorch (A. Paszke et al. 2019), enabling automatic differentiation and allowing us to use the model with gradient-based optimization (i.e., we can use backpropagation and gradient descent to update the weights of our NN). The transfer functions ψ(τ|λ, η) are parameterized by a set of accretion disk and black hole parameters, given by the vector η.

We use a modified version of the Novikov-Thorne (NT) model (I. D. Novikov & K. S. Thorne 1973), the relativistic version of the Shakura-Sunyaev (SS) thin-disk model (N. I. Shakura & R. A. Sunyaev 1973). At large radii, the NT and thin-disk models predict viscous temperature slopes of T visc ∝ R -3/4 . Quasar microlensing studies have measured a general temperature profile of the form T ∝ R -β and favor shallower slopes β < 3/4 (M. A. Cornachione & C. W. Morgan 2020). Some accretion disk models predict shallower or steeper slopes than the thin-disk slope of β = 3/4. For example, the slim-disk model predicts a shallower slope of β ≈ 0.5 (M. A. Abramowicz et al. 1988), while the magnetorotational instability model of E. Agol & J. H. Krolik (2000) predicts a steeper slope of β = 7/8. The temperature slope can also be modified due to the presence of wind outflows, where the accretion rate becomes radially dependent (R. D. Blandford & M. C. Begelman 1999;B. You et al. 2016;Y.-P. Li et al. 2018;M. Sun et al. 2018;J. Huang et al. 2023;J. H.-H. Chan et al. 2025). We use the NT plus lamppost temperature profile:

where R is the radius away from the black hole on the disk, M is the black hole mass, ( ) M R is the radially dependent accretion rate, f NT (R, a) is the dimensionless flux factor from the NT model given in Appendix A, H is the corona height, η X is the X-ray radiative efficiency, and σ SB is the Stefan-Boltzmann constant. We could instead use the SS model with f SS , but we choose to use the NT model since it includes general relativistic (GR) corrections. To model the deviations of the viscous temperature due to wind inflows or other deviations from the NT model, we take the accretion rate to be a power law of the form

where s is the power-law slope of the accretion rate related to the asymptotic slope of the viscous temperature profile β by s = 3-4β, R in is the inner radius of the disk, equivalent to the innermost stable circular orbit (ISCO), and M in is the accretion rate at R in (J. H.-H. Chan et al. 2025). For reference, R in /R g = 9, 6, 1 for a = -1, 0, 1, respectively, where R g = GM/c 2 is the gravitational radius (I. D. Novikov & K. S. Thorne 1973). We want to define the accretion rate based on the Eddington ratio. To do this, we set the accretion rate for the thin-disk case to

( ) = = M s M 0 in 0 , where = M M 0 Edd Edd for Eddington ratio λ Edd and Eddington accretion rate ( / = M GMm c 4 Edd p T

2 ), where m p is the proton mass, σ T is the Thomson scattering cross section, and η is the overall radiative efficiency factor. We can then set M in such that the total viscous bolometric luminosity is fixed to the thin-disk case by

4 is the flux of the NT model at s = 0. The radiative efficiency η is set such that the total luminosity is constant with a and is given by

g in in the NT model. For reference, η = 0.0377, 0.0572, 0.4226 for spin a = -1, 0, 1, respectively. We note that for the SS model, the radiative efficiency should instead be set to η = R g /2R in . When β < 0.5, the bolometric luminosity diverges when integrating out to infinity (i.e., the denominator in Equation (8) becomes very large), but here we only consider the case when β ∈ [0.5, 1.0]. We include GR gravitational redshifting and Doppler shifting due to the Kerr black hole. In the rest frame of the disk, the effective temperature of each region of the disk is converted to disk flux through blackbody radiation:

where h is the Planck constant, k B is the Boltzmann constant, λ emit is the emitted photon wavelength, and g = 1/(1 + z Kerr ) is the redshift factor due to both gravitational redshifting and relativistic Doppler beaming effects near the Kerr black hole (C. T. Cunningham & J. M. Bardeen 1973). When simulating our transfer functions, we want to use the photons emitted at a redshift λ emit that will correspond to the observed wavelengths λ obs , in our case the effective wavelength of each LSST wave band. We must therefore account for the cosmic redshifting of the photon as well as the gravitational redshift and relativistic Doppler shifting from the Kerr black hole. An emitted photon will be observed at a wavelength

Kerr

where z is the cosmic redshift. The redshifting due to the Kerr black hole geometry is given by

where each g µν are the metric components of a Kerr spacetime, Ω K is the angular velocity, θ = θ inc is the inclination angle of the disk, and f is the azimuthal angle (J. P. Luminet 1979; S. Bhattacharyya et al. 2001;G. Mastroserio et al. 2018;M. Heydari-Fard et al. 2023). We approximate this effect without the need for GR ray tracing by ignoring the effects of light bending and assuming the accretion disk orbits the black hole at Keplerian angular velocity:

where the disk rotates in a circular orbit around the black hole (M. A. Abramowicz & P. C. Fragile 2013;G. Mastroserio et al. 2018). In this way, we can efficiently implement GR effects into our auto-differentiable accretion disk simulation. An example redshift map for different disk orientations is shown in the top panels of Figure 1. The spin has very little effect on the redshift map except in the inner region of the disk near the ISCO. The bottom panels show example transfer functions with and without all the GR effects. At typical rest- frame wavelengths and low inclination, these effects have minimal impact on the transfer functions. However, at high redshifts and large inclination, where the shorter wavelengths allow us to probe the inner regions of the disk, these effects can alter the shape of the transfer function, although the mean time remains nearly unchanged.

From the temperature profile defined in Equation ( 6), the blackbody flux in Equation (10), and standard time lags of the lamppost thin-disk model, we can generate the transfer functions (E. M. Cackett et al. 2007). See J. H.-H. Chan et al. (2025) for a full derivation of the transfer function. In order to be computationally efficient, we want to minimize the number of pixels we need in the grid to accurately calculate the transfer functions. We define an effective radius of the accretion disk by numerically solving for k B T eff (R c ) = hc/λ, where we use the wavelength of the reddest band (in our case the y band). To account for cases where there is no solution, we define a characteristic radius by

and calculate the transfer function using a 1000 × 1000 grid out to 100R c . This effectively calculates the disk out to an infinite outer radius, since the flux decays to an insignificant amount well before 100R c . The characteristic radius is the same as in J. H.-H. Chan et al. (2025) but using the full accretion disk model so cannot be solved for analytically.

The transfer functions introduce a wavelength-dependent time lag:

We normalize the transfer functions to represent a probability distribution such that ( ) = d 1 0 for all λ. The transfer functions are calculated out to 800 days, chosen to be sufficiently long such that going out further is insignificant even for the largest-mass black holes. We evaluate the transfer functions on daily intervals. The influence that each parameter has on the mean time delays and standard deviation of the transfer functions can be found in Appendix B.

Spectrum Model

In addition to the transfer functions, we model the quasar spectrum to obtain the mean flux in each broadband filter consistently with the accretion disk parameters. We do not need to model the spectrum in the auto-differentiable simulation, since the mean brightnesses of each wave band are independent free parameters fit to the observations, rather than being fixed by the accretion disk parameters.

To model the mean brightness in each band, we generate the quasar spectrum and then integrate it across the response function of each bandpass filter. Our accretion disk model that we defined in Section 2.2 can evaluate the mean flux of the disk at each wavelength for a given luminosity distance. We evaluate the mean fluxes corresponding to the observed wavelengths of the LSST filters (3000-11000 Å) to generate the continuum spectrum. To properly model the quasar spectrum, however, we need to also include the emission lines, flux from the host galaxy, spectral reddening due to extinction from dust along the line of sight, contamination in the Lyα forest, and the Lyman limit. These effects have been previously modeled by M. J. Temple et al. (2021) and calibrated for redshifts 0 < z < 5. We use the spectrum model of M. J. Temple et al. (2021), but replace the continuum from a BPL with the continuum generated from our accretion disk model. They use a set of templates to include the spectral bands, host-galaxy flux, and extinction. The Lyman limit is modeled by setting the flux at a rest-frame wavelength of λ < 912 Å to zero. Once the full spectrum is obtained, we use speclite (D. Kirkby et al. 2023) to integrate the flux across the LSST broadband filters and obtain the mean magnitude of each wave band. Figure 2 shows an example simulated quasar spectrum with the LSST response functions.

We first calculate the absolute magnitude of the i band M i from the continuum, which is used to set the scaling of the emission lines and host galaxy. We add an additional scatter to our spectrum by adding ( ) N 0, 0.1 mag to M i . The other main parameter is the extinction strength E(B -V ), which we randomly draw from ( ) N 0, 0.075 and take its absolute value. The extinction has a significant effect on the quasar brightness. We further add additional scatter by varying the emission lines by controlling the relative scale height of the Hα, Lyα, and near-line region, which we vary ( ) N 0, 0.05 . In real quasar data, individual emission lines can further vary in width and strength depending on the accretion disk parameters and the broad-line region of the disk. For example, if there is a stronger than average emission line in a wave band, then the quasar would appear brighter, and there would be longer time delays, so we may overestimate the black hole mass. We include an additional scatter in the mean brightness of each band drawn from ( ) N 0, 0.025 mag. We also include an additional constant offset to all bands drawn from ( ) N 0, 0.025 mag.

We limit the redshift to the range z ∈ [0.1, 5.0], which includes the vast majority of quasars that will be observed by LSST. Within our selected redshift range, due to the cutoff in the Lyman limit, we do not always observe across all bands. The bluest band of LSST starts at 3000 Å, so if the Lyman limit cutoff is redshifted past this range (starting at z = 2.29), parts of the LSST observations will be affected. At the highest redshifts, both the u and g bands will not be observed at all. To keep within the observational range of LSST, if the mean magnitude of the i band is greater than 27 mag or less than 13 mag, then we resimulate a new quasar with a random set of parameters. The 27 mag limit represents the faintest possible objects observable by LSST, and the 13 mag lower limit accounts for saturation with bright sources.

We obtain the mean brightness of each band both with and without the host galaxy, because the host-galaxy flux should not be variable. The driving signal is first produced in magnitude with mean zero and standard deviation σ, and added to the mean magnitude of each band without the host-galaxy contribution. Next, we convert from magnitude to flux and convolve with each transfer function kernel. The flux from the host galaxy is then added, and finally we convert back from flux to magnitude.

As mentioned in Section 2.1, the driving signal is initially generated 8 times larger than needed to asymptotically set the mean and standard deviation in the longer signal. We then take the first eighth of the signal, where the local mean and standard deviation will not be the same as their asymptotic values. When constructing the training set, the driving signal is set to the normalization of the reference i band's brightness before convolving with the transfer functions, since its mean brightness is arbitrary.

Parameter Ranges

The parameters to generate our light curves are given in Table 1. The light curves are simulated with 12 physical parameters, four for the driving variability and eight for the accretion disk. The driving variability depends on its standard deviation σ, the break frequency ν b , and the low-and highfrequency power laws α L and α H . The black hole geometry is determined by its mass M and dimensionless black hole spin a.

The quasar is at an inclination angle θ inc (with 0°being faceon) and a redshift z. We also vary the corona height H, corona X-ray strength f lamp , asymptotic viscous temperature slope β, and Eddington ratio λ Edd . When sampling the light curves, we draw each parameter uniformly from its minimum and maximum range. There are also additional parameters related to the quasar spectrum that influence the mean magnitudes, in particular the extinction strength E(B -V ), that are not predicted by our ML model.

We parameterize the ISCO height as (H -R in )/R g , since the corona must always be above the ISCO. We parameterize α H -α L instead of using α H directly, because we expect α H � α L . The X-ray radiative efficiency factor η X = (η/λ Edd ) (L X /L Edd ), where L X /L Edd ∼ 0.005 (F. Ursini et al. 2020). We parameterize the lamppost strength, which we define as

, since the albedo and X-ray luminosity ratio are degenerate. The albedo ranges from A = 0 for full absorption to A = 1 for full reflection, with typical values A ∼ 0.1-0.2 (F. Ursini et al. 2020). The driving variability parameters are chosen to be consistent with measurements from D. P. Summons (2007).

Mock LSST Observations

After the UV/optical light curves are simulated, we degrade them to mimic LSST-like observing cadences and noise in the same way as Section 2.4 of J. Fagin et al. (2024). The errors at each LSST observation are determined by

where σ sys is the systematic error, and σ rand is the photometric noise. The systematic error is set to 0.005 mag, the maximum value expected for LSST (Ž.

Ivezić et al. 2019; K. L. Suberlak et al. 2021). The photometric noise depends on the brightness of each observation, and is expected to follow ( ) ( ) ( ) ( )

where γ is a band-dependent factor, m is the magnitude of each observation, and m 5 is the 5σ depth of a point source observed at the zenith of the observation. We expect γ u = 0.038 and γ g,r, i,z,y = 0.039 (Ž. Ivezić et al. 2019;X. Sheng et al. 2022). We simulate LSST-like observations using rubin_simfoot_3 with the baseline_v2.1_10yrs rolling cadence, which gives the time and m 5 of each observation. We produce a random sample of 100,000 LSST-like observations by sampling anywhere in the sky with 750-1000 total observations across all bands. This filters the light-curve sample to include only the Eddington ratio -2 0

Note. The parameters are sampled uniformly between the minimum and maximum given in the table. The top four parameters are for the X-ray driving variability, while the bottom eight parameters relate to the black hole and accretion disk reprocessing.

Wide Fast Deep observations, the main LSST survey (see Figure 1 of A. Prša et al. 2023).

We randomly deviate each observation with variance LSST 2 . We combine observations to the nearest 1 day interval. We simulate 10.5 yr light curves with 10 yr of LSST observations, so our model is trained to reconstruct some time before and after the survey. The start time of the first LSST observation is randomly selected within the extra half-year. Any observation outside the range of 13-27 mag is excluded, which should be similar to the best-case magnitude limits of LSST. We discard observations in each band with fewer than 15 total observations across the 10 yr to account for bands near the magnitude limits.

Machine Learning Model

Overview

The overall goal of our ML model is to solve the inverse problem in Equation (1) related to reconstructing the driving variability X(t) and transfer function kernels ψ(τ|λ) given a set of noisy and sparsely sampled observations. In the context of quasar variability, this inverse problem involves deducing the underlying physical processes (e.g., the driving variability and accretion disk reprocessing physics) that produce the observed light curves. Unlike many inverse problems, where the forward process is well understood, this task is complicated by the fact that quasar variability itself is stochastic. This adds an additional layer of uncertainty to the problem, as the driving variability X(t) is not directly observable with UV/optical data and must be inferred by the ML model.

Since the reprocessing is formulated as a convolution of the driving variability and transfer function, we are essentially training a ML model to solve a blind deconvolution problem. Blind deconvolution is a classic problem in signal processing where both the input signal (in this case, X(t)) and the kernel (here, ψ(τ|λ)) must be inferred simultaneously from the observed data. This problem is particularly challenging because of the degeneracy between X(t) and ψ(τ|λ), i.e., different combinations of driving variability and transfer functions can produce the same observed flux. Moreover, the stochastic nature of X(t) exacerbates this challenge, as it introduces variability that cannot be precisely predicted.

In general, recovering the transfer function kernels is an intractable problem because of this degeneracy. We could instead treat the bluest band as an effective driving signal and define effective kernels with respect to it. This is what is done with JAVELIN for example, but then we sacrifice most of the physics of the reprocessing on the disk. To make the problem tractable, we parameterize the transfer functions through our auto-differentiable simulation ψ(τ|λ, η), where η is the vector of accretion disk parameters, given in Table 1. The driving signal is reconstructed using a latent SDE (X. Li et al. 2020;J. Fagin et al. 2024), parameterized by the context at each time and the latent vector ẑ. We then solve for the flux by embedding the physics of the reprocessing of the driving variability into our NN architecture by

, 18 flux predicted 0 latent SDE auto diff transfer function

where this convolution of the driving variability and transfer functions is evaluated numerically within our ML model using fixed-grid numerical integration. Our ML model also quantifies the uncertainty in the reconstructed F λ (t, λ), X(t), ˆ, the variability parameters, and the mean time delays ¯coming from each reconstructed transfer function ( ˆ) , . The input to the ML model is the brightness and error values (both in magnitude), for a total of 12 features at each time step (six bands and six errors). For training stability, each band is normalized to have mean zero and standard deviation one, but the mean and standard deviation are used to predict the latent space of the SDE and the parameter posteriors. The reconstructed driving signal and UV/optical variability are unnormalized after they are generated by the ML model. For each band that is not observed at a given time step, we set both the brightness and error to a dummy value to be masked by the ML model. The outputs of our ML model are the best-fit reconstruction of the UV/optical light curves, the reconstructed driving signal, the reconstructed transfer functions, the predicted accretion disk and driving variability parameters and uncertainty, and the predicted relative time delays between bands and uncertainties (where the i band is arbitrarily chosen as the reference band). The model is trained using supervised learning to best reconstruct the light curve, driving signal, parameters, and time delays.

Proof-of-concept Recurrent Inference Machine

We use a RIM to compare the reconstruction of our light curve with the observations, and iteratively adjust the accretion disk parameters by ˆˆ= +

and the latent space of the latent SDE by ˆˆ= + z z z i i i in Equation (18) with iteration i and initial values of zero. Due to constraints in GPU memory, only a small batch size can be used, and this process can be computationally costly. We therefore train the ML model first without this process and then test training using RIM with three iterations, although ideally we would use more iterations and train longer.

Model Architecture

A diagram of the model architecture is shown in Figure 3. The model has two main components to encode the light curve: the context network of the latent SDE, and the RIM network used only to adjust ˆand ẑ. Both contain bidirectional RNNs including a GRU-D layer (Z. Che et al. 2016), a type of gated recurrent unit (GRU; J. Chung et al. 2014) layer that is designed to handle the masking and irregular sampling. In addition, we combine the output of the RNN layers with a transformer encoder (A. Vaswani et al. 2023). The GRU-D layers work well at handling the irregular sampling in the observed light curves by masking the unobserved time steps, while the transformers provide improved ability to handle long-term dependencies. Both the context and RIM networks require GRU-D layers to encode the irregularly sampled observations. The use of transformers also allows us to efficiently scale the model to many more parameters compared to RNNs. We use bidirectional RNN layers since they process the time series both forwards and backwards.

After the light curves are encoded, the context and latent vector ẑ are used by the neural SDE solver and then projected with an RNN to produce the reconstructed driving signal. The context and RIM are also used to predict the estimated accretion disk parameters ˆthat are used by the autodifferentiable simulation to produce the reconstructed transfer functions. The driving signal and transfer functions are then convolved within the ML model, numerically evaluating Equation (18) to produce the mean reconstructed UV/optical light curves. The mean reconstruction of the driving signal and UV/optical light curve is scaled, and the uncertainties are quantified using another RNN. In addition, the mean time delays are found from the reconstructed transfer functions, and the uncertainty is quantified. We also predict the variability parameters of the driving signal, using additional information from analyzing the driving signal such as fitting its PSD with five linear fits. The inferred driving and accretion disk parameters are combined, and the uncertainty is quantified.

In summary, the ML model infers a posterior distribution on the accretion disk and variability parameters, a posterior of the time delay between wave bands, reconstructs transfer functions, and reconstructs the UV/optical light curve and the unobserved driving signal. Our model is built in PyTorch (A. Paszke et al. 2019) and has 38,551,054 parameters, the majority of which come from the two transformers. More specific details about the NN architecture are given in Appendix C.

Uncertainty Qualification and Loss

We train our model with supervised learning by minimizing a weighted sum of four loss components. The driving variability is modeled as a latent SDE. Unlike in J. Fagin et al. (2024), we do not sample the latent space of the SDE from a posterior distribution, but instead just directly predict ẑ. We find this to significantly improve the performance of the light-curve reconstruction, since otherwise the reconstructed driving signal is not as adaptable. Our light-curve reconstruction is found by convolving the driving variability with the predicted transfer function kernels and scaling it to best match the observed data. The first component of the loss is the negative Gaussian log-likelihood between the reconstructed and true light curve, including the unobserved driving signal and averaged across the observed bands:

where y is the true light curve, ŷ is the predicted mean, ˆ2 is the predicted variance, and N is the number of time steps. We predict ( ˆ) log 2 instead of ˆfor stability and to force the variance to be positive.

We include an additional component to the loss to ensure that our mean light-curve reconstruction closely matches the observations:

where σ LSST is the error of each observation, m t is a mask that is 1 if a band is observed at a time step t and 0 if it is not, and

is the number of observations. The additional term in the negative log-likelihood

does not depend on our predictions and is therefore not included.

For the parameter inference, we use a multivariate Gaussian mixture model. Each multivariate Gaussian has negative loglikelihood:

where y are the vectors of true values of our parameters, ŷi are the vectors of our predicted means, and Σ i are the covariance matrices of the prediction. To ensure that the covariance matrix is symmetric, positive semi-definite, and nonsingular, we predict the lower triangular matrices L i representing the Cholesky decomposition of the covariance matrices = L L i i i . The diagonal of L must be positive, so we take the softplus of each diagonal element. An m × m lower triangular matrix will have m(m -1)/2 free parameters, with m = 12 for our case. We include into the loss the negative loglikelihood of the Gaussian mixture model:

with mixture coefficients a i of our n = 5 multivariate Gaussians, normalized such that = = a 1 i n i 1 . We additionally predict the relative time delay of each band with respect to the reference i band. These time delays are based on our predicted transfer functions, but we also infer the lower triangular matrix to assign multidimensional uncertainty and include into the loss:

1 evaluated even for unobserved bands. The overall loss is the weighted sum of each term:

weighting L param twice to reflect our focus on parameter inference. When using a RIM, each term is additionally scaled by its iteration index to emphasize later iterations. We use uniform priors for the variability and accretion disk parameters when simulating our light curves, given in Table 1.

When training our ML model, we reparameterize the parameter labels from their physical values to between zero and one and then take the logit 12 to scale them from -∞ to ∞. We then evaluate the negative log-likelihood of our parameter posterior in this logit space. We take the sigmoid 13 after drawing samples from the posterior, scaling the predictions back to between zero and one before scaling back to the original physical range. This transformation prevents any posterior probability from being wasted in physically impossible parameter space or outside the range of the training set (e.g., we must restrict the spin to -1 < a < 1).

Training

We train our ML model with 100,000 light curves per epoch that are randomly regenerated on the fly. Regenerating the training set each epoch prevents overfitting by training with many more unique examples. This is especially important since our model has a large number of free parameters. We use fixed test sets of 10,000 light curves to evaluate the performance of our model after training.

Our ML model is trained for 30 epochs using the Adam optimizer (D. P. Kingma & J. Ba 2017). We train with four A100 GPUs (80 Gb) in multiple stages. We first train without RIM to train faster. We use an initial learning rate of 8 × 10 -4 , determined through empirical tuning to balance convergence speed and stability. The learning rate is exponentially decayed by 0.95 each epoch, and a batch size of 26 per GPU (the maximum that fits in GPU memory). In the final four epochs, we use three RIM iterations with a batch size of 14 to test the RIM technique. Throughout training, we use gradient clipping with a maximum gradient norm of 250 to help prevent exploding gradients. Training took around 3 weeks. We did not employ any systematic hyperparameter tuning due to computational constraints in training the NN.

Results

Light-curve Reconstruction Performance

An example reconstructed UV/optical light curve, X-ray driving signal, and transfer function from the nominal test set is shown in Figure 4. Our ML model is able to properly quantify the uncertainty in its reconstruction, and it reconstructs the driving signal and transfer functions using just the UV/optical observations. Only the transfer functions from the mean parameter posterior that are convolved with the driving signal are shown. Although we do not directly infer uncertainty in the reconstructed transfer functions, an ensemble of reconstructed transfer functions could be produced by sampling from the parameter posterior.

To evaluate the reconstruction performance of our model, we compare it to an exact multitask GPR baseline with DRW kernel on our test sets, the same baseline as in Section 4 of J. Fagin et al. (2024). We use a GPR baseline since it is the standard method of reconstructing quasar variability (e.g., Z. Stone et al. 2022), and has been previously used to compare performance to ML methods (e.g., Y. Tachibana et al. 2020;E. Danilov et al. 2022;J. Fagin et al. 2024).

To test the robustness of our ML model to irregular variability, we use additional test sets with different out-ofdistribution driving signals. Table 2 reports the median ± median absolute deviation of the negative Gaussian loglikelihood across each test set of 10,000 light curves to provide an outlier-robust summary. We find that our ML model outperforms the GPR baseline for each driving signal, despite being trained only using the bended BPL. The exact parameter space of each driving signal is given in Appendix D. Statistical significance is assessed by applying paired t-tests to the full set of 10,000 paired log-likelihood differences, yielding p < 10 -6 for each test set. We choose to test our model on driving signals that are variable (BPL, DRW), quasi-periodic (BPL +sine), and periodic (sine, sawtooth, square wave). In each case, we are still able to predict the accretion disk parameters despite using out-of-distribution driving signals. We demonstrate this in Appendix D, and show that we are able to predict the black hole mass, redshift, temperature slope, and Eddington ratio similar to the case for the bended BPL. We also show an example light-curve reconstruction using a sawtooth driving signal to demonstrate how our model fits an out-of-distribution light curve.

In the top panel of Figure 5, we show that the uncertainties we predict in our UV/optical light-curve reconstruction are well calibrated for the nominal test set, while the GPR baseline is misaligned. We may expect the GPR baseline to be misaligned, since it assumes a DRW kernel, while our driving signal is generated from a more general BPL PSD and convolved with the transfer function kernels.

Parameter Inference Performance

In Figure 6, we show the median prediction compared to the true value for each parameter across the test set, demonstrating the ability of our model to predict each parameter. In this 12

figure, we only show the median, but a full multivariate posterior is predicted for each light curve (see, e.g., Figure 12 of J. Fagin et al. 2024). In the bottom panel of Figure 5, we show that the uncertainties we predict in our parameter posteriors are well calibrated, although overall slightly underconfident.

For the accretion disk parameters, we can predict the mass and redshift with high accuracy and place meaningful constraints on the temperature slope and Eddington ratio. The ML model struggles to predict the inclination angle and black hole spin but can make some constraints, while the corona height and lamppost strength cannot be constrained at all. We show the influence that each parameter has on the transfer functions in Appendix B, which explains why some of the accretion disk and black hole parameters can be predicted more easily than others, although there is additional information in the mean brightnesses. The mass and temperature slope predictions are better than in J. Fagin et al. (2024) due to modeling the mean brightness of each band, and we can now also predict the redshift and constrain the Eddington ratio. Our ML model cannot predict the inclination angle well because the inclination angle has almost no effect on the mean time delays, except at very high inclinations. While there is a higher-order effect of the inclination on the standard deviation of the transfer function, the associated brightness suppression is less informative, as we properly model the brightness of each band. The difficulty in constraining the inclination angle is primarily why the light curve and mean time delays can be well reconstructed while the exact shape of the transfer functions can still be off, such as in Figure 4. The shape of the predicted redshift is not continuous due to the u and g bands being suppressed at high redshift (starting at z = 2.29), causing the predictions to be very accurate at specific redshifts where these bands are being cut off. In addition, we are able to predict the bended BPL PSD parameters of the driving signal, although it is more challenging than the DRW parameters (see J. Fagin et al. 2024). This is due to the bended BPL having two additional parameters corresponding to the low-and highfrequency limits of the power spectrum, and the parameters are approximately degenerate. We note that we restrict the posterior distribution to be within the uniform prior, as mentioned in Section 3.4. In Figure 6, we only show the median predictions, which by design will never reach the edges of the parameter space. In the worst case, the ML model will predict a median at the center of the prior, as is the case for the corona height and lamppost strength.

We demonstrate in Appendix D that our model can still infer the black hole mass with out-of-distribution driving signals, and the performance of all accretion disk parameters remains consistent with the results shown in Figure 6. In addition, we quantify the root mean squared error (RMSE) of each parameter for each test case. For instance, the nominal RMSE is 0.33 for ( ) / M M log 10 and 0.32 for the redshift.

Example Relative Time Delay Inference

The mean time delays between bands are found directly from the mean times of our predicted transfer functions so that our model is entirely self-consistent. The time-delay differences are defined with respect to the i band as a reference. Our model predicts time delays even in the case where we do not observe the bluest bands, since we still have the transfer functions corresponding to our predicted accretion disk parameters. Our ML model then produces uncertainties associated with the mean time delays by predicting the lower triangular matrix L of its covariance matrix Σ = LL ⊤ . An example posterior of the time-delay measurements is shown in Figure 7. The time-delay posterior is directly inferred by our ML model, instead of requiring MCMC sampling.

Recurrent Inference Machine Evaluation

We evaluate the performance of our model using the RIM for three iterations. The average loss across the nominal test set does decrease each iteration, given by 0.619, 0.580, and 0.579. Therefore, the RIM procedure has a minor improvement in the performance of the model compared to using just a single iteration.

Discussion

Comparison to Traditional MCMC-based Methods

We constructed the first ML approach to model the UV/ optical light curve, driving variability, accretion disk reprocessing transfer functions, and time delays in a single unified framework by embedding the accretion disk reprocessing model into our NN architecture. This makes our predictions more interoperable compared to modeling each process individually. Our method enables the fast inference of accretion disk and variability parameters, time delays between wave bands, and the reconstruction of the driving signal and UV/optical bands.

Unlike traditional methods of measuring time delays, our model can use information related to the mean brightness and variability to inform our ML model on the accretion disk parameters and time-delay measurements. We also model the driving variability in a free-form way (i.e., using a latent SDE) instead of using a DRW process like JAVELIN or a Fourier series like CREAM. In addition, we directly predict the variability parameters of the driving signal without requiring it to be a Gaussian process. Furthermore, we can predict more parameters of the accretion disk, including the mass, Eddington ratio, redshift, temperature slope, and variability parameters of the unobserved driving signal. Comparatively, CREAM typically only measures MM .

While our new method is more computationally demanding to train than J. Fagin et al. (2024) due to having many more parameters and the use of a RIM, during inference it is still quick enough to easily apply to the tens of millions of quasar light curves expected from LSST in several hours. Specifically, the average inference time on our GPUs using a batch size of 64 is about 27 minutes per million light curves using three iterations of the RIM, or 7 minutes using one iteration. Comparatively, using JAVELIN or CREAM on tens of millions of light curves would be infeasible. This is currently a major limitation of traditional MCMC-based methods; however, faster methods may be developed that take advantage of GPU acceleration to overcome this limitation.

Instead of being model dependent, we could also use analytic functions such as using simple top-hat transfer functions like JAVELIN defined with respect to the bluest band (i.e., only five transfer functions for six bands since the bluest band is treated as the effective driving variability). One may also reconstruct the effective transfer functions in a freeform manner, but this approach would require regularization or the use of a complete set of basis functions to make the inverse problem tractable. In such cases, the X-ray driving variability would not be reconstructed, which was a major goal of this work, but these methods could be explored in future studies.

Limitations of This Work

In this work, we did not account for time delays that may arise from the broad-line or narrow-line regions of quasars. Modeling the BLR is difficult since its geometry is not well known. Including the BLR into our simulation would introduce a large number of additional parameters (H. J. Best et al. 2025). The time delays caused by the BLR Table 1) across the nominal test set. The ideal case where the median prediction matches the truth is given by the black dashed line across the diagonal. can cause our method to overestimate the recovered black hole mass, since the BLR can introduce larger time lags. However, this is also the case for traditional sampling methods like JAVELIN and CREAM. Furthermore, there could be a radially dependent albedo or color-correction factors to the blackbody radiation. In addition, the thin-disk NT model may break down at high Eddington ratios, and different disk models such as the slim-disk model (M. A. Abramowicz et al. 1988) may be more accurate. Some of the deviations from the thin-disk model may already be partially accounted for by the inclusion of the accretion rate wind model (Equation ( 7)) that modifies the slope of the temperature profile. There may also be long negative time lags at the viscous timescale of the disk (P. Z. Yao et al. 2023;A. Secunda et al. 2023), which our model could be adapted to predict and could be used to study the vertical structure of the disk. We also assumed that each light curve is coming from type 1 quasars, but there may be contamination from type 2 quasars that are misidentified. We tested that our model is robust to different out-of-distribution driving signals, but it would be useful to test how our pretrained model behaves for different disk models. Moreover, we assumed that each driving signal was generated from a stationary process, but there could be nonstationary variability such as from flaring or tidal disruption events.

In our training set, we assumed we had the full 10 yr of LSST data. In future work, we plan to evaluate the model's performance using different stages of the survey, such as using only the first 1, 3, and 5 yr of LSST data. We may also compare the performance of our model between light curves from the Wide Fast Deep and Deep Drilling Fields. Furthermore, a fraction of the LSST quasar light curves will have observations from previous surveys that can be combined with LSST data to lengthen the light curve.

We find the driving signal parameters of the BPL to be more difficult to predict than the DRW parameters used in J. Fagin et al. (2024). This is because there are two additional parameters, and the additional degrees of freedom introduce degeneracies. Here we parameterized the driving variability as a bended BPL, since it has been found to fit X-ray variability data. Some authors such as M. Papoutsis et al. (2024) fix the lower-frequency slope to α L = 1, which would get rid of degeneracies, and our predictions for ν b and α H would improve. We choose to keep the bended BPL more general to account for the wide range of variability possible in LSST. We could, however, use any parameterization of the PSD to train our model.

Discussion of Machine Learning Architecture

As far as we are aware, this work represents the first application of transformers in modeling quasar variability. To process the irregularly sampled time series, we first use GRU-D and GRU layers and then transformer encoders concatenated with the output. The GRU-D handles the irregular sampling more effectively than transformers, which we found struggled to converge on their own. By combining both RNN and transformer architectures, we achieved stable convergence while leveraging the transformer's improved handling of longterm dependencies in lengthy time series. This approach also enabled us to scale our model to tens of millions of parameters, making it the largest deep learning model applied to quasar variability to date.

In future work, we may test improvements to the NN architecture. One major modification could be to reconstruct the light curve by repeatedly sampling driving signals and accretion disk parameters from the latent space. This would work like a variational auto-encoder to yield nonparametric posterior distributions. However, with our current training paradigm this would require repeated sampling during training to obtain approximate posterior distributions to use in the loss function, and this would be too computationally demanding to train. It could also be possible to replace all RNNs with just a single transformer encoder. Furthermore, we could explore other encoder architectures (e.g., M. Schirmer et al. 2022). In this work, we bin the time series and transfer functions to daily time intervals. Ideally, we would use smaller bins, but we are limited by GPU memory and training time. Using daily intervals could negatively impact the parameter inference of quasars with smaller-mass black holes (M ≲ 10 7 M ⊙ ), where the time delays can be on the order of 1 day or less.

We use our framework to iteratively improve the light-curve reconstruction and accretion disk parameter estimation with a RIM by analyzing the residuals between the predicted light curve and observations. We found the RIM to have only a minor advantage compared to using a single iteration. This is likely due to several factors. For example, we pretrained our model using a single iteration to save on training time. We also only use the RIM to adjust the mean accretion disk parameters and latent vector of the SDE, but perhaps all parameters should be iteratively adjusted. In future work, we plan to explore using a single encoder, instead of separate context and RIM modules, to iteratively adjust the entire context. Here we use a very large model with tens of millions of parameters, and it may be the case that a minimum is achieved with only a single iteration. In addition, our time series are stochastic, irregularly sampled, and noisy. Therefore, we cannot reconstruct the time series to the pixel level, as is the case for image reconstruction tasks where a RIM has traditionally been applied. The RIM technique may be more advantageous for well-sampled light curves, such as in the Deep Drilling Fields.

Outlook and Future Applications

Once LSST data becomes available, we can fine-tune our model using the data with self-supervised learning. This would help generalize our ML model trained with simulated data and supervised learning to the real LSST observations. This could be done simply by minimizing the loss with respect to how well our reconstructions match LSST observations (similar to Equation (20)). The light-curve reconstruction and parameter inference are linked through our auto-differentiable simulation of the accretion disk, so by fine-tuning the model on the observations, we can potentially improve all aspects of our predictions. Furthermore, semi-supervised learning could be used if there is a partial set of LSST light curves that have some parameter estimates from spectra. In addition, the parameters of the BPL driving signal may be related to properties of the accretion disk and black hole (P. Arevalo et al. 2024). In our training set, we choose to keep them independent so that these correlations could potentially be measured by our model without bias. However, our model could improve its predictions by finding these correlations when fine-tuning with the real data. Fine-tuning our model with data from the Deep Drilling Field quasar sample would be particularly beneficial. In future work, we could also explore adapting our method for fully self-supervised training.

The auto-differentiable simulation we develop enables new physics-based ML methods like ours to be developed to model quasar variability. In addition, more traditional methods that rely on gradient-based sampling, such as Hamiltonian MCMC (M. Betancourt 2018), can be developed using our simulation. Hamiltonian MCMC can greatly speed up the rate of convergence compared to traditional MCMC by using gradient-based guidance. Even without the use of autodifferentiation, our simulation can take advantage of GPU acceleration to speed up traditional MCMC methods. We benchmark the speed of our simulation with and without GPU acceleration. For a single transfer function, we find the GPU to be 1.5× faster than on the CPU (taking 0.54 s compared to 0.81 s). With batches of 50, we find this improves to a 67× speedup on a single GPU (taking 0.55 s per batch compared to 36.9 s). The batch size can be scaled to as many samples that can fit in GPU memory for maximum speedup. The driving variability can also be generated using PyTorch, and is included in the codebase although not directly used in this work.

Our accretion disk simulation includes a wind-based accretion rate and GR effects, and it can be readily adapted and expanded upon in future work. Our framework enables hypothesis testing by training multiple ML models, each incorporating a different disk model, and evaluating their inferred parameters and time delays against both observational data and general relativistic magnetohydrodynamics (GRMHD) simulations. In observational data, these quantities can be measured through independent methods, providing a complementary means of assessing model accuracy. In GRMHD simulations, the true underlying values are known, allowing for a direct comparison between model predictions and physically motivated disk structures (e.g., S. Koudmani et al. 2024). This approach provides a systematic way to evaluate competing accretion disk models and refine our understanding of disk variability. Our framework supports the integration of any auto-differentiable disk models, ensuring flexibility in exploring a wide range of accretion physics scenarios. Even complex simulations could be made autodifferentiable by training a NN to emulate them and embedding the fixed-weight pretrained model into our NN architecture. In addition, we model the quasar spectrum and integrate the mean brightness across the LSST response functions to obtain physically consistent mean brightnesses of each wave band. This is a major upgrade from J. Fagin et al. (2024), since this information can be used to break parameter degeneracies and infer photometric redshifts. This is especially important in cases with low variability.

The framework developed in this work offers significant potential for anomaly detection in various contexts. For instance, the latent space of the SDE can be used in an unsupervised way to identify out-of-distribution sources of variability by comparing it to the typical distribution of quasar light curves or employing methods like isolation forests. Additionally, by analyzing the reconstructed driving signal, we may detect quasi-periodic signals that could indicate the presence of a supermassive black hole binary. The model could also be extended through supervised learning to classify light curves as belonging to phenomena such as supermassive black hole binaries, changing-look AGN, flaring events, tidal disruption events, or gravitational microlensing.

LSST will observe tens of millions of quasars, and only a small fraction of them will get follow-up spectra to obtain precise redshift measurements. Therefore, most quasars will rely on photometric redshifts. In theory, the photometric redshift predictions of our ML model could outperform simple photometric redshift estimators based on the brightness differences between wave bands. This is because the photometric redshifts should depend most strongly on the asymptotic brightnesses of each band, and our model incorporates modeling the time-variable brightness. Furthermore, the redshift affects the time delays and driving variability timescale, so by analyzing all these processes in a unified framework, we can outperform methods relying on static data. Quasars that do have follow-up spectra can be better constrained by using the redshift as a prior in our ML model. The mass and other properties of the accretion disk may be constrained from measurements of the broad-line region (S. Panda et al. 2019), which could also be used to inform our ML model. Some of the accretion disk parameters are difficult to predict for individual quasars. Hierarchical inference could be used with the entire LSST sample to estimate the population-level distributions of the parameter space (S. Wagner-Carena et al. 2021). For example, the population-level distribution of the temperature slope β could be used to test accretion disk and wind models, despite being challenging to constrain for individual quasars with only six bands. Furthermore, we could measure the relationship between driving variability parameters ν b , σ, α L , and α H , and black hole properties like the mass.

Conclusions

Our model fits the UV/optical variability, reconstructs the driving variability signal, predicts accretion disk and variability parameters, and measures the relative time delay between bands, all self-consistently using a ML model that incorporates the physics of the accretion disk reprocessing into its architecture using an auto-differentiable simulation and latent SDEs. We incorporate transformers into our model to enhance its capacity for capturing long-term dependencies and to scale it to tens of millions of parameters. Our method is ready to be applied to the entire sample of tens of millions of monitored LSST quasars in a matter of hours. In comparison, using GPR frameworks such as Celerite (D. Foreman-Mackey et al. 2017) and traditional curve-shifting techniques like JAVELIN (Y. Zu et al. 2011Zu et al. , 2016) ) and CREAM (D. A. Starkey et al. 2015) will be computationally infeasible. Furthermore, we test the robustness of our ML model by comparing the performance of our pretrained model on light curves with out-of-distribution driving signals including variable, quasi-periodic, and periodic signals. We find our model outperforms a multitask GPR baseline in all cases and can still infer the accretion disk parameters.

We aim to incorporate as much physics into our ML model as possible. The latent SDE captures the physics of the stochastic driving variability, while the auto-differentiable simulation models the reprocessing of the driving signal on the accretion disk. In previous works, there has been no mechanism to ensure that the inferred accretion disk parameters correspond to the time delays in the light-curve reconstructions (J. W. Park et al. 2021;J. Fagin et al. 2024). In addition, our model fits the power spectrum of the reconstructed driving signal with several linear segments to better estimate the driving signal parameters. By embedding these physical processes into the NN, we achieve a model that is more robust and interpretable compared to traditional blackbox parameter estimators. Furthermore, by linking these processes the model can refine its parameter estimates when using self-supervised learning on real observations, such as those from LSST. On real data, we can also compare the predicted time delays of our model to those obtained using traditional curve-shifting techniques, providing a way to validate the time-delay measurements and search for anomalies. Additionally, this method can test reprocessing models by comparing the reconstructed driving signals to X-ray data. The ML approach we present is highly general and can be adapted to other multivariate time series with irregular sampling, particularly for blind deconvolution or inverse problems. Our ML model and auto-differentiable and GPU-accelerated accretion disk simulation are open-sourced and available on GitHub,foot_4 with a copy preserved in Zenodo (doi:10.5281/ zenodo.15446208).

with the hidden state of the GRU cell being updated with each RIM iteration. The context is also used to predict the parameters, the latent space of the SDE, and the normalization of the light curve.

There is an additional RNN without the transformer to produce uncertainty in the reconstructed driving and UV/ optical variability with the same architecture but no transformer. We use one more RNN to project the output of the latent SDE to the mean and standard deviation on the observation space (see J. Fagin et al. 2024). This RNN has a linear skip connection between the output of the SDE and a two-layer bidirectional GRU-based RNN, so the RNN layers can be easily bypassed if they are unnecessary.

Multi-layer perceptrons (MLPs) produce the posterior parameters of the accretion disk and variability, the latent space of the SDE, the normalization of the light curve, and the uncertainty in the time-delay estimates. While η and ẑ are adjusted iteratively, the rest of the inferred parameters are produced directly each iteration. Each MLP takes as input the context, the mean and standard deviation of the input light curve, as well as the other previously predicted parameters of the network, while η and ẑ also use the output of the RIM network. To adjust η and ẑ, we use an additional two-layer network, depending on the RIM network output and iteration number, that scales the updates by a factor dependent on the iteration count, allowing the network to gradually suppress changes as the RIM approaches convergence.

There are two MLPs in the latent SDE: the posterior drift function that decodes the context and latent vector, and the diffusion network that is applied element-wise to satisfy the diagonal noise (X. Li et al. 2020). The latent SDE solves the equation ( ) ( ( ) ) ( ( ) ) ( ) ( ) µ = + µ d z t t z t d t t z t d W t , ; , ; , C1 where µ(t, z(t);θ µ ) is the drift term that determines the deterministic component of the latent dynamics with learnable parameters of the NN θ µ , σ(t, z(t);θ σ ) is the diffusion term that represents the stochastic component of the latent dynamic with learnable NN parameters θ σ , W(t) is a Wiener process (Brownian motion) that captures the random fluctuations over time, and ( ) = z z 0

is the latent vector or initial condition of the SDE. The latent SDE consists of an Itô SDE solver using the Euler-Maruyama numerical approximation scheme broken up into 2000 time intervals of 2.3 days each.

We use a Gaussian mixture model to parameterize the posterior of the accretion disk and variability parameters with five multivariate Gaussians (discussed further in Section 3.4). For the mixture coefficients, we use a softmax activation function to normalize them to probabilities. The variability parameters take as additional input the results from analyzing the reconstructed driving signal by conducting five linear fits to different pieces of the power spectrum and using the slope and y-intercept. We also use the mean, standard deviation, mean absolute deviation, total variation, and total square variation, as well as the results of the two latent SDE MLPs at ẑ. To obtain a vector of accretion disk parameters ˆfrom the posterior distribution, we sample from the posterior and take its mean. We use the mean times in the reconstructed transfer functions to obtain the relative time delays with respect to the reference i band, and we assign uncertainties to the prediction. To normalize the light curve, we include three terms: a mean and standard deviation for each band's brightness in magnitude, as well as an additional bias term in the flux. The bias term is required since the host galaxy can contribute nonvariable flux. We note that when a band is not observed, its mean brightness is set to the maximum limit of 27 mag and standard deviation to 1 mag to inform the ML model.

Each MLP consists of four fully connected layers. The RNNs use tanh activation, transformers use GELU (D. Hendrycks & K. Gimpel 2023), and fully connected layers use LeakyReLU (A. L. Maas et al. 2013). Whenever appropriate, we include residual skip connections (K. He et al. 2015) and layer normalization (J. L. Ba et al. 2016) to enhance gradient flow, stabilize training, and improve model convergence. We use a hidden size of 256, transformer encoder size of 512, context size of 128, and latent size of 16. The transformer encoders have five layers, eight heads, and use a sinusoidal time embedding.

Appendix D Robustness Test

To test the robustness of our trained ML model to out-ofdistribution variability, we compare its performance on test sets of light curves with BPL, DRW, BPL+sine, sine, sawtooth, and square wave driving signals (see Section 4.1). The BPL is our nominal test set and uses the same parameter range as our training set, given in Table 1. The DRW is the same as the BPL but with α = 0 and α H = 2, which is outside the parameter range of our training set. We use a sine function with period and relative width of each rectangular pulse chosen between [0.1, 0.9]. For each periodic signal, we select a random phase to shift the start time. We show an example reconstruction of a light curve with sawtooth driving variability in Figure 9. Despite the power spectrum of a sawtooth wave being very different than a BPL, our model still fits the observations and the signal is contained within the predicted 2σ credible interval. Since our model is trained using only stochastic signals, it will not recognize that this light curve is actually periodic but is still able to reconstruct the accretion disk parameters and time delays between wave bands.

To evaluate our model's ability to predict the accretion disk parameters for out-of-distribution light curves, we compare the median black hole mass prediction to the true values for each test set in Figure 10. We also give the normalized RMSE for each test set in Table 3. Our model can estimate each accretion disk parameter similarly to Figure 6, although there are more outliers when predicting the redshift, temperature slope, and Eddington ratio. Furthermore, it can still infer the time delays between wave bands, as demonstrated in Figure 9. Note. We note that predicting the mean of our uniform prior would yield a normalized RMSE of / 1 12 28.87%.