Upcoming wide field surveys such as the Rubin Observatorys Legacy Survey of Space and Time LSST will monitor thousands of strongly lensed quasars over a 10year period. Many of these monitored quasars will undergo high magnification events HMEs through microlensing as the accretion disk crosses a caustic places of infinite magnification. Microlensing allows us to map the inner regions of the accretion disk as it crosses a caustic even at large cosmological distances. The observational cadences of LSST are not ideal for probing the inner regions of the accretion disk so there is a need to predict HMEs as early as possible to trigger highcadence multiband or spectroscopic followup observations. Here we simulate a diverse and realistic sample of 10year quasar microlensing light curves to train a recurrent neural network RNN to predict HMEs before they occur by classifying the location of the peaks at each time step. This is the first deep learning approach to predict HMEs. We give estimates at how well we expect to predict HME peaks during LSST and benchmark how our metrics change with different cadence strategies. With LSSTlike observations we can predict approximately 55 of HME peaks corresponding to tens to hundreds per year and a false positive rate of around 20 compared to the number of HMEs. Our network can be continuously applied throughout the LSST survey providing crucial alerts to optimize followup resources.
['Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Instrumentation and Methods for Astrophysics']
Predicting High Magnification Events in Microlensed Quasars in the Era of LSST using Recurrent Neural Networks
{'Predicting High Magnification Events in Microlensed Quasars in the Era of LSST using Recurrent Neural Networks': "Joshua Fagin, 1, 2, 3 Eric Paic, 4 Favio Neira, 4 Henry Best, 1, 2, 3, 5 Timo Anguita, 6, 7 Martin Millon, 8 Matthew O'Dowd, 1, 2, 3 Dominique Sluse, 9 and Georgios Vernardos 1, 2, 3  \n1 The Graduate Center of the City University of New York, 365 Fifth Avenue, New York, NY 10016, USA  \nDepartment of Astrophysics, American Museum of Natural History, Central Park West and 79th Street, NY 10024-5192, USA  \n2  \n3 Department of Physics and Astronomy, Lehman College of the CUNY, Bronx, NY 10468, USA  \n- 4 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland\n- 5 Department of Theoretical Physics and Astrophysics, Faculty of Science, Masaryk University, Kotlářská 2, CZ-611 37 Brno, Czech Republic\n- 6 Instituto de Astrofísica, Departamento de Ciencias Fisicas, Facultad de Ciencias Exactas, Universidad Andres Bello, Av. Fernandez Concha 700, Las Condes, Santiago, Chile  \n7 Millennium Institute of Astrophysics, Nuncio Monseñor Sótero Sanz 100, Providencia, Santiago, Chile  \n8  \nKavli Institute for Particle Astrophysics and Cosmology and Department of Physics, Stanford University, Stanford, CA 94305, USA 9 STAR Institute, Université de Liège, Quartier Agora Allée du six Août, 19c B-4000 Liège, Belgium", 'ABSTRACT': "Upcoming wide field surveys such as the Rubin Observatory's Legacy Survey of Space and Time (LSST) will monitor thousands of strongly lensed quasars over a 10-year period. Many of these monitored quasars will undergo high magnification events (HMEs) through microlensing as the accretion disk crosses a caustic, places of infinite magnification. Microlensing allows us to map the inner regions of the accretion disk as it crosses a caustic, even at large cosmological distances. The observational cadences of LSST are not ideal for probing the inner regions of the accretion disk, so there is a need to predict HMEs as early as possible to trigger high-cadence multi-band or spectroscopic follow-up observations. Here we simulate a diverse and realistic sample of 10-year quasar microlensing light curves to train a recurrent neural network (RNN) to predict HMEs before they occur by classifying the location of the peaks at each time step. This is the first deep learning approach to predict HMEs. We give estimates at how well we expect to predict HME peaks during LSST and benchmark how our metrics change with different cadence strategies. With LSST-like observations, we can predict approximately 55% of HME peaks corresponding to tens to hundreds per year and a false positive rate of around 20% compared to the number of HMEs. Our network can be continuously applied throughout the LSST survey, providing crucial alerts to optimize follow-up resources.  \nKeywords: Quasars(1319) - Active galactic nuclei(16) - Quasar microlensing(1318) - Gravitational lensing(670) - Neural networks(1933) - Time series analysis(1916)", '1. INTRODUCTION': "The most widely accepted unified model of Active Galactic Nuclei (AGN) postulates that their structure revolves around a central Super Massive Black Hole (SMBH) at the center of galaxies (Antonucci 1993; Urry & Padovani 1995). The tidal stretching and viscous  \nCorresponding author: Joshua Fagin  \njfagin@gradcenter.cuny.edu  \nfriction experienced by matter in the immediate vicinity of the SMBH creates hot plasma, forming an accretion disk that emits a bright continuum of light across a wide range of wavelengths with stochastic variability. Quasars are bright AGN with unobscured accretion disks. Their radial energy profile is believed to follow the thin-disk model (Shakura & Sunyaev 1973; Page & Thorne 1974), but recent observations from microlensing (Poindexter et al. 2008; Poindexter & Kochanek 2010; Blackburne et al. 2015; Muñoz et al. 2016; Mor-  \nan et al. 2018) and reverberation mapping (Mudd et al. 2018; Jha et al. 2022) have found larger accretion disk sizes by a factor of ∼ 2-4. Further out, clouds of ionised gas form the broad and narrow line regions (BLR and NLR), reverberating accretion disk light and inducing broad and narrow line emissions in the quasar spectra depending on the gas velocity (Blandford & McKee 1982; Sluse et al. 2007; O'Dowd et al. 2015; Grier et al. 2017). Quasars also influence the evolution of their host galaxy due to the energy outflow they create (e.g. Di Matteo et al. 2008; Fabian 2012; Kormendy & Ho 2013). Quasars are among the most luminous objects in the Universe, making them valuable cosmological probes for measuring the expansion rate of the Universe (Lusso et al. 2019; Wong et al. 2020) and understanding the re-ionization mechanism (Grissom et al. 2014). These, however, require a precise understanding of quasar structure and the underlying physical processes.  \nStrongly lensed quasars occur when light from a quasar source is bent by a galaxy along the line of sight, forming multiple images of the quasar. Each image of the strongly lensed quasar can be independently magnified through microlensing due to stars or other compact objects in the lensing galaxy. This effect leads to a micromagnification of the light profile of the accretion disk. It is possible to measure the size of microlensed accretion disks using single epoch, multi-waveband observations using the different apparent size of the disk in different wavebands (e.g. Mediavilla et al. 2015). Alternatively, long-term monitoring of a strongly lensed quasar allows us to detect microlensing by measuring the micromagnification variation as the alignment between the microlenses and the background quasars varies over time. This long-term monitoring requires suppressing the intrinsic variability of the quasar by taking the pairwise difference of two strongly lensed images, corrected for the time delay between the images. If the intrinsic variability is perfectly subtracted, the resulting difference light curve, that we will refer to as the microlensing light curve , only contains the microlensing signal. These microlensing light curves provides constraints not only on the size of the disk (Kochanek 2004; Morgan et al. 2008; Cornachione & Morgan 2020; Rivera et al. 2024) but also on the geometry of the Broad Line Region (Sluse et al. 2011; Sluse & Tewes 2014; Paic et al. 2022; Savić et al. 2024) and can reveal the existence of substructures in the accretion disk, such as secondary supermassive black holes (Yan et al. 2014; Millon et al. 2022).  \nA high magnification event (HME) occurs when a quasar approaches or crosses a microcaustic, a region of  \ninfinite magnification. As a result, a specific area of the quasar's accretion disk becomes highly magnified. As the quasar moves relative to the microcaustic, the differential magnification provides a way of scanning though these different emission regions. In order to constrain the structure of the AGN engine though microlensing, however, we require high-cadence and multi-wavelength observations. These microlensing events can last from as little as a few weeks to a few years. It has been shown that X-ray monitoring of the HME can place constraints on the spin, inclination, and size of the innermost stable circular orbit (ISCO) of the black hole (e.g. Chartas et al. 2017). Furthermore, optical monitoring has the opportunity to measure the temperature profile of the inner regions of the disk (e.g. Eigenbrod et al. 2008). High-cadence and multi-wavelength observations near caustic-crossing events are crucial to extract the detailed accretion disk structure encoded in microlensing light curves (e.g. Krawczynski et al. 2019; Vernardos & Tsagkatakis 2019; Best et al. 2024). In the best case scenario, the entire caustic-crossing event will be observed by expensive and high-cadence follow-up to place the most stringent constraints on the inner workings of the AGN. Therefore, being able to predict these events ahead of time is a critical step in order to trigger appropriate follow-up observations.  \nThe most furnished microlensing light curve data set to date was released by COSMOGRAIL (Millon et al. 2020b) with 20-year long microlensing light curves. Out of the 23 observed lensed quasars (17 doubly lensed and 6 quadruply lensed) and a total of 203 seasons of monitoring, only three microlensing light curves display microlensing events higher than 1 magnitude. These light curves only contain observations from the r -band, severely limiting the constraints that can be placed on the AGN structure and the ability to predict HMEs.  \nThe advent of wide-field monitoring surveys such as the Rubin Observatory's Legacy Survey of Space and Time (LSST) is expected to increase the sample of monitored lensed quasars by at least tenfold. LSST is expected to simultaneously monitor hundreds to thousands of strongly lensed quasars (Oguri & Marshall 2010) in six UV/optical wavebands ( ugrizy ) at around 55-185 samplings per band or around 800 observations across the ten years (Abell et al. 2009; Prša et al. 2023). Tens to hundreds of quasar images should undergo HMEs of more than 1 magnitude per year (Neira et al. 2020). We cannot accurately constrain the innerstructure of the accretion disk using microlensing light curves at LSST cadences alone due to the sparse and irregular sampling. However, LSST microlensing light curves will still allow us to detect the characteristic  \nsteady rise of a microlensing event occurring up to years before the peak. Furthermore, by using different wavebands, we can potentially determine the precise timing of the HME peak as it affects various regions of the disk causing the wavebands to cross.  \nAn HME alert system was first developed for the Optical Gravitational Lensing Experiment (OGLE; Woźniak et al. 2000), and a triggering mechanism for multi-band follow-up observations was proposed by Wyithe et al. (2000). These efforts focused on predicting HMEs in Q 2237 + 0305 due to its brightness and rapid microlensing time scale but have not been used extensively. Thus far, there has been no prediction and follow-up of an HME as it unfolds. This is likely due the small sample of lensed quasars and a lack of long-term, multi-band monitoring that could distinguish microlensing from the intrinsic variability of the quasar. With LSST and other wide field surveys, however, this will no longer be the case. New and more sophisticated methods must be developed to quickly and autonomously be applied to the thousands of monitored lensed quasars expected in the near future.  \nIn this work, we aim to create a neural network able to predict the occurrence of microlensing HME peaks using the multi-band microlensing light curves of LSST. This tool will be crucial to optimize follow-up observations of such HMEs and caustic crossings. Machine learning was first applied to microlensed quasar light curves in Vernardos & Tsagkatakis (2019) to measure the accretion disk size and temperature profile using a convolutional neural network for simulated data. Best et al. (2024) also used a convolutional neural network to measure the black hole mass, inclination angle, and impact angle of simulated caustic crossing events. As a pilot application, they measured the mass and inclination angle from archival caustic crossings of Q 2237 + 0305 and found them to be consistent with previous estimates. Machine learning has also been extensively used to model the intrinsic quasar variability (e.g. with simulated LSST quasar lightcurves, see Fagin et al. 2024). Here we train a recurrent neural network (RNN) using binary classification to identify the peaks of HMEs. Our RNN is able to account for the irregular cadences and seasonal gaps expected from LSST, and could be continually applied to LSST microlensing light curves in real time as new data is obtained.  \nIn Section 2, we describe how we build our simulated data set of microlensing light curves. Section 3 presents our machine learning model architecture, training, and an example application to an LSST microlensing light curve. In Section 4, we benchmark the performance of our RNNs using different observational strategies includ-  \ning the baseline case of simulated LSST-like observations. In Section 5, we discuss the results of our network and its applicability to LSST. Section 6 gives our concluding remarks. In this work, we assume a flat Λ CDM cosmology with H 0 = 72 . 0 kms -1 Mpc -1 , Ω m = 0 . 26 , and Ω Λ = 0 . 74 .", '2. MICROLENSING LIGHT CURVE SIMULATION': 'In order to train a machine learning model to predict HME peaks, we require a library of microlensing light curves that is representative of the expected LSST sample. Since the goal is to predict the peaks of the HMEs identified on the light curves, we focus on constructing a library where the microlensing durations and amplitudes cover the range of what could be expected in current and future systems. Each light curve in the training set has two main components: the microlensing signal and the associated correlated noise. The noise in the microlensing lensing light curves is typically highly temporally correlated, often referred to as red noise as opposed to uncorrelated Gaussian white noise (Millon et al. 2020c; Paic et al. 2022). The red noise is mainly due to imperfect subtraction of the intrinsic variability of the quasar or actual physical processes such as the intrinsic variability being echoed by the BLR and differentially magnified by microlensing (Sluse & Tewes 2014; Paic et al. 2022). Since these physical processes are not well understood and difficult to model, we rely on an empirical approach, where we treat all these effects as red noise added to our simulated microlensing light curves.', '2.1. Microlensing simulation': 'Microlensing brightness variations are uncorrelated between the multiple images of a strongly lensed quasar, unlike the intrinsic variability which is identical in each image. Since the latter can be subtracted, albeit imperfectly leaving the so-called red noise (see next section), we focus on simulating the microlensing variations. To model the microlensing variability, we use the tool presented in Neira et al. (2020), which is able to generate simulated microlensing light curves for each image of the lensed quasars.  \nThe tool takes as an input an accretion disk, velocity, and macro lens model and outputs a library of light curves. The light curves are generated as follows:  \n- · The accretion disk model is convolved with magnification maps selected from the GERLUMPH project (Vernardos et al. 2014).\n- · In the convolved map, multiple tracks are defined to simulate the net transverse movement of the disk, accounting for the movement of the lens, mi-  \nolenses, and observer, based on the velocity parameters and probing time.  \n- · The pixel values along these tracks correspond to the microlensing light curves (see section 2 in Neira et al. (2020) for a description of the tool, models, and parameters).  \nWe use the simulated microlensing light curves from Neira et al. 2024 (in prep.), which makes use of the catalogue of lensed quasars presented in Oguri & Marshall (2010). A total of 2,821 mock systems are simulated. While the detail of the models and parameters used to generate the microlensing light curves from these systems can be found in Neira et al. 2024 (in prep.), here we briefly describe them.  \nThe mass distribution of each lensing galaxy is modeled as a singular isothermal ellipsoid (SIE). For an SIE, the convergence and shear at the position of the quasar images is defined as:  \nκ ( x, y ) = γ ( x, y ) = 1 2 θ E √ q ω ( x, y ) , (1)  \nwhere θ E is the Einstein radius, q is the axis ratio of the lensing galaxy, and ω ( x, y ) is the elliptical radius where x, y define the coordinates of the lens plane with the x-axis aligned with the major axis of the lens. An additional term is added to the shear that corresponds to the lens environment, resulting in a small scatter in the final κ and γ distribution, shown in Figure 1. The smooth matter fraction, s , is computed by assuming a light distribution and mass-to-light ratio of the lens (see Foxley-Marrable et al. 2018; Vernardos 2019, and Neira et al. 2024 (in prep.) for details). This allows us to select a magnification map from the GERLUMPH library (Vernardos et al. 2014) corresponding to the selected κ , γ , and s . Each map has 10,000 pixels per side with a resolution of 0 . 025 R E /pixel, where R E is the Einstein radius of the microlensing objects and is defined as:  \nR E = √ D S D LS D L 4 GM c 2 , (2)  \nwhere D L , D S , and D LS are the angular diameter distances from the observer-to-lens, observer-to-source, and lens-to-source respectively, G is the gravitational constant, c is the speed of light, and M = 0 . 3 M ⊙ is the mass of the microlenses (the approximate mean of a Salpeter initial mass function (Salpeter 1955)).  \nThe effective transverse velocity that defines the track along the magnification map is defined as:  \nυ e = υ o 1 + z l D LS D L + υ ⋆ 1 + z l D S D L + υ g , (3)  \nFigure 1. Distribution of the lensing macromodel parameters γ and κ for our quasar sample.  \n<!-- image -->  \nFigure 2. Distribution of transverse velocity υ e in our simulated light curves.  \n<!-- image -->  \nwhere υ o is the transverse velocity of the observer as measured with respect to the cosmic microwave background (CMB) velocity dipole (Kogut et al. 1993), υ ⋆ is the bulk velocity of the microlenses, and υ g is the combined peculiar velocities of the lens and source. The distribution of υ e across all the light curves is shown in Figure 2.  \nFigure 3. Distribution of accretion disk radius r disk in our simulated lensed quasar sample. The sizes are in units of the Einstein radius of their respective system.  \n<!-- image -->  \nThe size of each accretion disk is assumed to follow a standard thin-disk model (Shakura & Sunyaev 1973):  \nR λ = 9 . 7 × 10 15 ( λ rest µm ) 4 / 3 ( M BH 10 9 M ⊙ ) 2 / 3 ( f E η ) 1 / 3 [ cm ] , (4)  \nwhere R λ is the size of the disk at λ rest , M BH is the mass of the black hole, f E is the Eddington ratio and η is the accretion efficiency. We adopted fixed typical values of f E = 0 . 25 and η = 0 . 15 (e.g. Blackburne et al. 2011). For the black hole mass we have adopted an empirical model from MacLeod et al. (2010):  \nlog 10 ( M BH ) = 2 -0 . 27 · M i , (5)  \nwhere M i is the absolute magnitude of the source in the i -band. Given a black hole mass estimate, the size of the accretion disk can be computed at any wavelength through Equation (4). The influence of the specific shape of the brightness profile has only minor influence in the microlensing signal (Mortonson et al. 2005; Vernardos & Tsagkatakis 2019). Instead, the signal is mostly influenced by the source size. Thus, we adopt a Gaussian-shaped two-dimensional brightness profile, where we match the size from Equation (4) to the half light radii of the profile. In Figure 3 we show the distribution of the disk sizes across all the systems with respect to their microlensing Einstein radius (see Equation (2)).  \nWith the above, we are able to generate microlensing light curves by subtracted the light curves between each pair of quasar images. For each pair of images for all simulated systems, we generate 100 light curves in the six  \nLSST bands. We label all peaks that have a minimum amplitude of 0 . 5 magnitude in the u -band, which is the bluest and thus most prone to microlensing variations. This value is chosen to be roughly above the maximum amplitude of the added red noise. After the difference light curves are generated and the red noise is included from Section 2.2, we apply an additional threshold for event labels to avoid including event peaks with little microlensing in the difference light curves that could lead to our network predicting more false positives. We apply a threshold of 1 magnitude from the maximum to minimum of any time step in the light curve across all bands. If the light curve is below the threshold, the labels are set to zero. This threshold is applied in order to remove the initially labeled peaks in the cases where the difference light curve shows little brightness variability. We justify this because the typical length of the events are of the order of years, and the length of the light curves are set to 10 years. This makes it a common occurrence that two events from different light curves overlap, thus yielding a difference light curve that has less variability. We produce a data set comprising of 282,100 light curves. We randomly select 80% of the light curves for training and reserve 10% for validation and 10% for testing.', '2.2. Red noise generation': 'The reverberation of the continuum by the BLR with a time lag induces an echo of the intrinsic variability within the observed light curve. If the broad emission lines fall in the monitored band and the microlensing is not identical in both images, this echo will appear in the microlensing light curve (e.g. Sluse & Tewes 2014). There may also be extended continuum emission coming from other regions of the disk (Sluse 2024). The frequency and amplitude of this imprint is characterized by the size of the BLR and the variability of the continuum, which is well described by a damped random walk (MacLeod et al. 2010). Furthermore, the variability of the light curve (hence the microlensing light curve) can be affected by observational effects such as the seasonal change of airmass and contamination by the lensed arc or the lensing galaxy when measuring the photometry of the images (see section 3.3.3 of, Sluse et al. 2006). Since these effects all require precise knowledge about the physical properties of the lens system, we choose to assimilate them into the red noise that is added to the immaculate microlensing light curve.  \nFor this work, we use data from the COSMOGRAIL program (Courbin et al. 2005), which released the longest quasar microlensing light curves (Millon et al. 2020b), observed in the r -band with the Swiss 1.2-metre  \nFigure 4. The top panel shows a COSMOGRAIL r -band microlensing light curve of RXJ 1131 -1231 over a period of 18 years computed by subtracting image B from image A after shifting it by the time delay of 2.8 days (Millon et al. 2020b). A spline fit with η = 300 days is applied to the microlensing light curve. The bottom panel displays the residuals of the spline fit, where the correlated noise is clearly visible.  \n<!-- image -->  \nLeonhard Euler Telescope. In particular, the microlensing light curve of RXJ 1131 -1231 displays the events with the highest amplitude and the strongest correlated noise. Hence, we use this data to quantify a conservative red noise that can be added to LSST-like light curves. As shown by Figure 4, the microlensing light curve is first fit with a free-knot spline using the PyCS3 package (see Millon et al. 2020a, for details of the implementation). These piece-wise polynomials allow for a smooth fit of a targeted time scale of variation by constraining the distance η between two consecutive knots. Microlensing events are believed to occur on times scales longer than several years (e.g. Mosquera & Kochanek 2011), we therefore choose η = 300 days to prevent the spline from fitting intra-season features while recovering the microlensing variations, the fit obtained is shown on the top panel of Figure 4, and the residuals are shown on the lower panel. As shown in Figure 5, the power spectrum of the data and of the spline fit are identical for frequencies below 1 / 200 days -1 , which sets the boundary between high and low frequencies. At higher frequencies, the power spectrum of the data is identical to the residual one. Our goal is then to generate time  \nseries with the same power spectrum as the residuals to add it to the simulated microlensing light curves. We use the red noise generator implemented in PyCS3 to fit the slope β and amplitude σ of the residuals power spectrum. To account for the change in slope of the power spectrum of the residual above 1 / 10 days -1 , we fit β and σ separately in the frequency windows [ 1 / 200 , 1 / 10 ] days -1 and [ 1 / 200 , 1 / 10 ] days -1 . The addition of a generated red noise to the spline fit is shown in Figure 5 and has a power spectrum compatible with original data. The same parameters are then used to add red noise to each simulated microlensing light curve.', '2.3. Baseline observation strategy': 'We first simulate microlensing light curves with weekly cadence. We then want to degrade the cadences to mimic different observation strategies. Our baseline is LSST-like observations which are simulated using the observation times produced by rubin\\_sim 1 and the baseline\\_v2.1\\_10yrs cadences. We sample 100,000  \n<!-- image -->  \nFigure 5. The left panel shows the power spectrum fitting of the residuals shown in the bottom panel of Figure 4. The vertical dashed lines at 1 / 200 days -1 and 1 / 10 days -1 delimit the frequency range on which the red noise is added. Adding residual-like noise to the spline fit (top panel of Figure 4) creates a new realization of the data with the same power spectrum. The right panel shows an example microlensed light curve with one realization of noise added with the same power spectrum as the left panel.  \n<!-- image -->  \ndifferent sky positions that have between 750 and 1,000 total observations across the ten years to include only light curves from the main LSST survey (see section 2.4 of Fagin et al. (2024)). The time of each observation is combined to the nearest weekly interval to maintain a fixed number of time steps in our light curves. Times that are not observed in a band are masked by the RNN. Microlensing events happen over time scales of months to years, so combining the light curves to weekly intervals should not affect our ability to predict HMEs. Separate sky positions are used for the training, validation, and test sets to avoid bias.', '3.1. Network architecture': 'Our neural network 2 is trained to classify HME peaks of quasar microlensing light curves. The goal is to train our model so that it can be applied every week throughout the lifetime of LSST, including during seasonal gaps where there are no new observations. This is done by including only observations made at or before the current time of the prediction, and then masking any further time steps. The network is trained to classify if there is an event peak within a 21 week time window before or after each time, chosen to be a bit more than half a season. We then make a prediction for each time step  \nstarting when there are 85 weeks of data so the network has at least two seasons or so of data to inform its predictions. We normalize the light curve by subtracting each magnitude by the first observation in the r -band (or the g -band for some observation strategies we test with no r -band observations). Data augmentation is used during training to avoid overfitting the training set. This comprises of selecting a random cadence strategy for the light curve (i.e. for the nominal LSST light curves, we select the 10-year cadence from a specific sky position) and a random time step of the observation. The validation and test sets remain fixed for comparison.  \nThe input into the RNN is the relative brightness at each time step for each band. There are 522 time steps corresponding to the brightness at each week across the ten years. Nominally we include all six LSST bands ugrizy , but we also explore situations where less bands are observed. Unobserved points are set to a dummy value of zero to be masked by our network. We include an additional feature that is set to one up to and including the current time step and then zero at later time steps. This feature is necessary to inform the network of the current time step, including when there are no new observations.  \nOur RNN processes the time series in two separate paths, both forwards and backwards, since bidirectional RNNs can have improved performance. Each RNN path first contains a GRU-D layer (Che et al. 2016) a modified version of the gated recurrent unit (GRU; Chung  \nFigure 6. Diagram of our RNN model. The numerical value in each layer represents the output size. We have two RNN paths, one for the forward light curve (LC) and one for the reverse. The last time step of the final RNN layers are input into the FC layers. Each RNN layer is followed by layer normalization, and the first two FC layers are followed by LeakyRelu and layer normalization. The final FC layer uses sigmoid activation function to normalize the output of our network to a probability for classification, i.e. restrict it between 0 and 1. We also include residual skip connections between the output of each layer and the next layer.  \n<!-- image -->  \net al. 2014) that is designed to handle masking and irregular time intervals. This is followed by two GRU layers. Each GRU layer has a tanh activation function and a hidden size of 128. Each RNN layer is followed by layer normalization (Ba et al. 2016). We also employ residual skip connections between the output of each RNN layer and the next layer which has been shown to improve training stability and performance (He et al. 2015). The output of each RNN path is combined and followed by two fully-connected (FC) layers. Each FC layer has a hidden size of 256 and are proceeded by a LeakyReLU activation function (Maas et al. 2013) and layer normalization. The FC layers also include skip connections to the previous layer. We then have one final fully-connected layer with a hidden size of 256 and output size of 1 that is followed by a sigmoid activation function. The sigmoid activation function forces the final output of the network to be a value between 0 and 1 representing the probability of the HME peak being within 21 weeks of the time step. A simple diagram of our RNN is shown in Figure 6. Each network has a to-  \ntal of 642,673 trainable parameters and is built using PyTorch (Paszke et al. 2019).', '3.2. Network training': 'The network is trained by minimizing a binary crossentropy (BCE) loss function given by:  \nL ( y, ˆ y ) = -1 N N ∑ i =1 y i log(ˆ y ) + (1 -y i ) log(1 -ˆ y i ) , (6)  \nwhere y i is the training label, ˆ y i is the prediction of the network, and N is the number of training examples. We minimize the loss function using an Adam optimizer (Kingma & Ba 2014) with a batch size of 512 and a learning rate of 0.002 that is exponentially decreased to 0.0002 over the course of training. For training stability, we employ gradient clipping with a maximum gradient norm of 1. We train the network for multiple passes of our training set known as epoch. We train for a total of 300 epochs. As mentioned previously, we use data augmentation to regenerate the training set each epoch by choosing a random cadence strategy and time of the observation. We train separate networks for a variety of different observational strategies to see how the number of observations, seasonal gaps, and different bands affects the performance of our model. The BCE loss, accuracy, precision, and recall compared to the epoch are shown in Figure 7 for the LSST baseline observations. We find further training beyond the 300 epochs to have very little improvement on the performance of our networks. The other networks, trained with the different observational strategies, show similar convergence after 300 epochs. Training each network took around four days with an NVIDIA A10 GPU (24 GB).', '3.3. Application to microlensing light curves': 'We propose a three stage triggering system for rapid follow-up of HMEs. At each time step, the network makes a prediction to classify if an HME peak is within a 21 week window of the latest data. We quantify a positive prediction if the network predicts a probability above 0.5. The three triggers are colored from lowest to highest alert as green, yellow, and red. Green alert is triggered after any positive prediction. Yellow alert is triggered if there is a green alert continuously for 7 weeks in a row. Red alert is triggered if there is a yellow alert for 7 weeks in a row and continues until a negative prediction. As soon as there is a negative prediction, the count is reset to avoid false positives. Figure 8 shows an example prediction of our network for an LSST-like light curve. In this case, our network is able to predict the HME peaks ahead of time. The network is still able to make predictions during the seasonal gaps even  \nFigure 7. BCE loss (shown in the top left panel and given in Equation (6)), accuracy (top right panel), precision (bottom left panel), and recall (bottom right panel) as a function of the training epoch for the training (blue) and validation (orange) sets.  \n<!-- image -->  \nthough no additional observations take place. Ideally the alert system should match up with the green label windows and we would predict the HME peaks 21 weeks beforehand, identifying their temporal locations. Another example is given in Appendix A that includes a false positive, correct, and missed prediction of our network.', '4. RESULTS': 'We explore a variety of observational strategies to assess how important high-cadence and multi-band data is in predicting and analyzing HMEs. For each observational strategy, we train a separate network and then evaluate its performance on the test set. We use the same test set but with different observational cadences to fairly compare the performance of each model. Our baseline observing strategy represents LSST-like observations and is described in Section 2.3. We also evaluate our model on a cadence strategy where there are no season gaps in the data and an observational mask is randomly sampled for each band at every time step with mean number of visits [57,72,186,194,169,174] across the ten years in the ugrizy bands respectively. We further consider a hybrid case where we have LSST observations but a subset of bands are observed with the random ca-  \ndence mask to test how our network would perform with partial observations during the seasonal gaps. In addition, we evaluate the performance in the idealistic case where we have regularly sampled data at every weekly time step. Moreover, we evaluate the performance when we have regular sampling but including season gaps every year that range from 140-180 days, chosen randomly for each light curve. For each cadence strategy, we additionally evaluate the performance using only a subset of the LSST bands. All the different strategies are summarized in Table 1 with various metrics to evaluate the performance of each model.  \nBased on the different observational cadences we tested, we find the color information to be the most important factor for predicting HMEs. We expect this to be the case since the different wavelengths probe different radii of the accretion disk, and therefore, are microlensed by different regions of the caustic map. With regular weekly observations in just the r -band, we are unable to make high fidelity predictions since we correctly identify just 24% of the HME peaks. When we have regular sampling in both the g and i bands, the network performs much better, correctly identify 58% of the peaks. Additional bands further increase the per-  \nFigure 8. Example prediction of our network with LSST-like observations. The bottom panel represents the microlensing light curve observations for the six LSST bands, normalized such that the first observation in the r -band is zero. The middle panel shows the output of our neural network by classifying if there is an event peak within a 21 week window per side at each time step. The first 85 weeks and last 21 weeks are grayed out since we do not make predictions for those times. The orange label represents the peak of the event while the wider, green label is the window that the network is trained to predict if there is an HME peak. The top panel represents a mock alert system where the first 7 weeks of a continuous positive prediction (i.e. the probability is greater than 0.5 at the time step) are set to green alert followed by 7 weeks of yellow alert and then greater than 14 weeks is red alert.  \n<!-- image -->  \nformance of our model, but having at least two bands is most important.  \nWe also find that the season gaps play a significant role in the performance of our models. In the example in Appendix A, the missed HME peak at around 3,200 days could be predicted if the characteristic rise and band crossing did not fall within the season. In addition, the exact location of the HME peak at around 2,000 days is not identified because it falls directly in the season gap. The network clearly knows there should be an HME peak somewhere in the season gap but is unsure about its exact location. With all LSST bands but no season gaps, we identified 69% of the event peaks instead of the LSST baseline case of 55%. With regular sampling and no season gaps, this increases to 72%. Thus the seasonal gaps play a larger role than the irregular sampling. We do find that with season gaps but regular sampling, we identify 66% of the peaks, so the performance loss from the season gaps can be made up with enough observations. Furthermore, we can make  \nup some of the performance gap if there are some observations without seasonal gaps, such in one or two bands.', '5. DISCUSSION': "Our trained network is ready to be applied to the entire sample of lensed quasars of LSST. In the initial stages of survey before making predictions, we can finetune the network to calibrate it to the exact observational cadences and red noise properties of the lensed quasar sample. It will also be necessary to measure the time delays between lensed images in order to produce the microlensing light curves. We begin making predictions after an initial time of 85 weeks so the network has at least two seasons or so of data and to reduce the potential uncertainty in the time delay measurements. A shorter initial time could be used if the time delays can be measured, although we would expect worse performance in the early time of the survey, since there would fewer observations to infer the microlensing time scales of the system.  \nOur baseline LSST observations come from selecting sky positions within the main LSST survey. A limited  \nTable 1. Summary of metrics comparing the performance of each separately trained network on their test set using different observation strategies. The baseline for LSST-like observations is given at the top. The accuracy, precision, recall, and F1 score come from evaluating our network on random time steps throughout the survey across the test set. We note that the positive labels represent only about 3.6% of the data set, so an accuracy of 0.964 could be achieved by just predicting the negative class. The correct peak, red in label, and false positive metrics come from evaluating each time step of the entire 10-year light curve for 1,000 examples in the test set, i.e. making predictions like in Figure 8, and weighted by the total number of HME peaks. Correct peak means the peak was correctly identified (i.e. the orange label in the middle panel of Figure 8 has predicted probability at least 0.5). Red in label means that there is a red alert within the 21 week event peak window (i.e. the green label in the middle panel of Figure 8 overlaps with the red alert in the top panel at some point). A false positive represents a red alert outside any event peak window (i.e. if the red alert in the top panel of Figure 8 did not overlap at all with the green label in the middle panel).  \nnumber of strongly lensed quasars may be monitored in the Deep Drilling Fields which will have much higher cadences. Our predictions would significantly improve from the increased number of observations. We expect the performance metrics for lensed quasars in the Deep Drilling Fields to be similar to the regular sampling with seasons in Table 1.  \nOur goal is to train a general HME predictor for LSST. Thus, we use a diverse training set that should be representative of the entire expected LSST quasar sample. For specific systems that have already been well studied prior to LSST, we could train individual models to make more tailored predictions. For example, an HME predictor could be developed specifically for Q 2237 + 0305 like in Woźniak et al. (2000). We would expect these models specially trained for individual systems to have improved performances, since it can specialize within a much narrower parameter space. We could also give the network priors on the κ , γ , s , and velocity distributions of systems if they have been estimated. Knowing these values could constrain the time scales of microlensing events and help the network identify where the peak of the HMEs are with more clarity. The alert system itself could also be tailored to individual systems or populations with similar microlensing timescales. We choose a 21 week window around the HME peaks, so the training labels extent to around a season long. Some systems will be easier to trigger ahead of time, and we may then  \ndefine a longer window, while others will be more difficult and may benefit from a narrower window. Fine tuning the length of the training labels can be explored in future work.  \nIn our microlensing simulation, we use a simple Gaussian brightness profile scaled to the disk size of each waveband. We expect that at the observational cadences of LSST, the shape of the inner region of the brightness profile has little effect on the microlensing light curve. Furthermore, we are particularly interested in predicting the HMEs ahead of time. At these time scales and cadences, only the size of the disk chosen as its half-light radius is important, as the exact shape of the brightness profile has little impact on the light curve (Mortonson et al. 2005). This is no longer the case during caustic crossing events, as the inner region of the accretion disk can be highly magnified. The analysis of the highcadence follow-up, where the inner detail of the accretion disk becomes encoded in the light curve and may be able to constrain system parameters (Best et al. 2024), is beyond the scope of this work.  \nWe build our training set by adding red noise from the residual power spectrum measured in the r -band of RXJ 1131 -1231 because this is currently one of the best sampled, highest signal-to-noise ratio light curves with well-measured time delays. Since the photometric Gaussian noise is fairly low, this system is ideal to isolate and quantify the correlated red noise. Since the  \nCOSMOGRAIL data only has r -band observations, we assumed the red noise in each wavelength to be equivalent. In reality, this might not necessarily be the case, depending on the physical origin of the red noise. If it originates from a flux contamination of the multiple image, we can expect a similar amplitude in all photometric bands. On the other hand, if it originates from the quasar variability echoed by the BLR, we may expect more correlated noise in the photometric bands overlapping with a broad emission line of the quasar. This would therefore depend on the redshift of the quasar. We choose a conservative approach for this work by using RXJ 1131 -1231, since it has one of the largest amplitudes of red noise. Once LSST data becomes available, the residual power spectrum of each LSST band can be measured from all monitored systems, and we can fine-tune our model based on the real population of red noise. The photometric errors are included in the red noise generated from RXJ 1131 -1231 from the Leonhard Euler Telescope, which should be similar to LSST's photometric errors for objects this bright. Nevertheless, by using the entire LSST population, we will have access to more realistic photometric errors for each system. For bright quasars of (<19 mag), we expect the photometric errors to be only ∼ 0.01 magnitude and play a minor role. For the fainter quasars (>20 mag), the photometric errors can be significant and our network should be adapted to better account for them. There may also be additional uncertainty captured within the red noise related to uncertainty in measuring the time delay between lensed images. RXJ 1131 -1231 is a well studied system and there are only small uncertainties on its time delay measurements. This will not be the case for all systems monitored in LSST, and this additional uncertainty will be incorporated into the population of red noise. Recent methods of image deconvolution could largely reduce these effects, hence reducing the amount of correlated noise in the extracted microlensing light curves (Millon et al. 2024).  \nIn our training set, we apply two thresholds to the HME labels. The first is that an event needs to be at least 0.5 magnitude to be considered a HME. This value is chosen to be just above the maximum amplitude of the added red noise. Furthermore, in the subtracted microlensing light curves, we apply an additional threshold of 1 magnitude between the maximum and minimum of the light curve across all bands. This second threshold helps exclude labels where long-term microlensing events cancel each other out, leaving little indication of an HME in the light curve. These thresholds are designed to minimize true false positives and avoid triggering expensive follow-up observations unnecessarily. In-  \ndeed, the choice of these thresholds will affect how the network performs. In Appendix A, we show an example where if the threshold would have been lower, the first false positive would have been considered a correct prediction, as this peak would have been labeled. Thus our reported false positive rate is a conservative estimate, since edge cases where the HME is just below the threshold will often be misclassified by the network. Our focus is on predicting the HMEs that have the most significant impact on the microlensing light curves. Future studies could explore adjusting these thresholds, potentially lowering them to capture more subtle events. This approach could be particularly useful in scenarios where the amplitude of the red noise is reduced, allowing for a broader range of HMEs to be detected.  \nIdeally we would simulate the intrinsic quasar variability and use a variable disk when convolving with the microlensing maps instead of using a data-driven approach. Then we would measure the time delays between each image using the initial simulated data of LSST and subtract the images to produce our simulated microlensing light curves. This would then include the effects from the uncertainty in measuring the time delays, the photometric errors, and red noise in a more self-consistent way. We would, however, need to model the BLR reverberation which can have very complicated geometry. There may also be extended continuum emission coming from other regions of the disk that we may not be able to simulate (Sluse 2024). Since accounting for all the aspects of the intrinsic variability may be infeasible, we prefer a data-driven approach and defer exploring the inclusion of these effects directly into the simulation for future work.  \nHere we fix the threshold for a positive prediction of our network to a probability of 0.5, as is typically done in binary classification. By raising this threshold, however, we would predict fewer HME peaks but lower the false positive rate. We could also adjust the balance of classes in the training set. During training, we mimic observing at a random time throughout the 10-year survey, and we find a rate of positive labels to be around 3.6%. We could instead choose to bias the training set by including more or less mock observations near the HME peaks. Including more time steps near the HME peaks would enable the network to make more positive predictions, increasing the number of correctly predicted HME peaks but also raising the false positive rate. Within our metrics, we include predictions that fall within seasonal gaps. By focusing low-cost photometric follow-up near the predicted HME, we could improve our prediction and reduce the number of false positives before trigger-  \ning more expensive follow-up observations, such as X-ray or spectroscopic monitoring.  \nWe may also improve upon our machine learning model in future work. We use RNNs that included masking of missing values through GRU-D layers (Che et al. 2016). We may also explore other RNN architectures such as stochastic RNNs, which have been used to model UV/optical quasar variability (Sheng et al. 2022). Continuous recurrent units (Schirmer et al. 2022), which are a probabilistic recurrent architecture where the hidden state evolves as a linear stochastic differential equation, are a promising way of dealing with the irregular sampling. We may also consider using latent stochastic differential equations which has worked well at modeling the intrinsic quasar variability for simulated LSST light curves (Fagin et al. 2024). Transformers (Vaswani et al. 2023) such as multi-time attention networks (mTANDs; Shukla & Marlin 2021), where the time embedding is fed to the attention mechanism to incorporate irregularly sampled time series, may also show improved performance compared to RNNs.  \nThe largest current catalog of monitored lensed quasars is COSMOGRAIL which has only r -band observations (Millon et al. 2020b). We show in Table 1 that for regular observations including seasonal gaps in the r -band, only around 15.7% of event peaks are correctly identified with a false positive rate of 12.2%. These metrics indicate that even with the most ideal case of regular weekly sampling, without multi-band data the network fails to make meaningful predictions. There are currently no catalog of long-term monitored lensed quasars with multi-band data, so we defer testing our model on real data for future work.", '6. CONCLUSIONS': 'We trained an RNN to predict HME peaks in LSSTlike microlensed quasar light curves. We make our model and training procedure open source and available at https://github.com/JFagin/Quasar\\_Event\\_Predictor. A pretrained neural network like ours could be quickly and autonomously applied to the entire sample of monitored LSST strongly lensed quasars each week as new data is obtained in less than a minute on a single GPU.  \nThis is the first deep-learning method of predicting HMEs.  \nFor the baseline LSST cadences, we can expect to correctly identified around ∼ 55% of HMEs with a false positive rate of ∼ 20%. We therefore expect to correctly identify tens to a couple hundred HMEs per year (Neira et al. 2020). These metrics can be improved as highcadence follow-up is triggered ahead of the peak of the event.  \nWe envision our network to be continuously applied throughout the entire LSST survey to give weekly notifications to experts on potential incoming HME peaks. These alerts will be crucial for optimizing follow-up observational resources. High-cadence multi-band or spectroscopic follow-up observations of these HMEs will offer unique insight into the inner structure of AGN as imprinted in the light curves during caustic crossing events.  \nSoftware: PyTorch (Paszke et al. 2019), Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), Astropy (Astropy Collaboration et al. 2018),', 'ACKNOWLEDGMENTS': 'This work originated in the Lensing Odyssey 2021 workshop, and so we would like to acknowledge the organizers and attendees for the fruitful discussions. This work used resources available through the National Research Platform (NRP) at the University of California, San Diego. NRP has been developed, and is supported in part, by funding from National Science Foundation, from awards 1730158, 1540112, 1541349, 1826967, 2112167, 2100237, and 2120019, as well as additional funding from community partners. Support was provided by Schmidt Sciences, LLC. for JF, HB, MO, and GV. MM acknowledges support by the SNSF (Swiss National Science Foundation) through mobility grant P500PT\\_203114 and return CH grant P5R5PT\\_225598. TA acknowledges support from ANID-FONDECYT Regular Project 1240105, ANID Millennium Science Initiative AIM23-0001 and the ANID BASAL project FB210003. HB acknowledges the GAČR Junior Star grant No. GM24-10599M for support.', 'A. ANOTHER EXAMPLE PREDICTION': 'In Figure 9 we show another example of the predictions of our network. In this case there is a false positive, correct prediction, and missed prediction from left to right.  \nFigure 9. Example of the prediction of our network on LSST-like observations like Figure 8.  \n<!-- image -->', 'REFERENCES': "Abell, P. A., Burke, D. L., Hamuy, M., et al. 2009, Lsst science book, version 2.0, Tech. rep., LSST Antonucci, R. 1993, ARA&A, 31, 473,  \ndoi: 10.1146/annurev.aa.31.090193.002353 Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f Ba, J. L., Kiros, J. R., & Hinton, G. E. 2016, Layer Normalization, arXiv, doi: 10.48550/ARXIV.1607.06450 Best, H., Fagin, J., Vernardos, G., & O'Dowd, M. 2024, Monthly Notices of the Royal Astronomical Society, 531, 1095, doi: 10.1093/mnras/stae1182 Blackburne, J. A., Kochanek, C. S., Chen, B., Dai, X., & Chartas, G. 2015, ApJ, 798, 95, doi: 10.1088/0004-637X/798/2/95 Blackburne, J. A., Pooley, D., Rappaport, S., & Schechter, P. L. 2011, ApJ, 729, 34, doi: 10.1088/0004-637X/729/1/34 Blandford, R. D., & McKee, C. F. 1982, ApJ, 255, 419, doi: 10.1086/159843 Chartas, G., Krawczynski, H., Zalesky, L., et al. 2017, ApJ, 837, 26, doi: 10.3847/1538-4357/aa5d50  \nChe, Z., Purushotham, S., Cho, K., Sontag, D., & Liu, Y. 2016, Recurrent Neural Networks for Multivariate Time Series with Missing Values, arXiv, doi: 10.48550/ARXIV.1606.01865  \nChung, J., Gulcehre, C., Cho, K., & Bengio, Y. 2014, Empirical Evaluation of Gated Recurrent Neural Networks on Sequence Modeling, arXiv,  \ndoi: 10.48550/ARXIV.1412.3555  \nCornachione, M. A., & Morgan, C. W. 2020, ApJ, 895, 93, doi: 10.3847/1538-4357/ab8aed  \nCourbin, F., Eigenbrod, A., Vuissoz, C., Meylan, G., & Magain, P. 2005, in Gravitational Lensing Impact on Cosmology, ed. Y. Mellier & G. Meylan, Vol. 225, 297-303, doi: 10.1017/S1743921305002097  \nDi Matteo, T., Colberg, J., Springel, V., Hernquist, L., & Sijacki, D. 2008, The Astrophysical Journal, 676, 33, doi: 10.1086/524921  \nEigenbrod, A., Courbin, F., Meylan, G., et al. 2008, A&A, 490, 933, doi: 10.1051/0004-6361:200810729  \nFabian, A. C. 2012, ARA&A, 50, 455, doi: 10.1146/annurev-astro-081811-125521  \nFagin, J., Park, J. W., Best, H., et al. 2024, The Astrophysical Journal, 965, 104,  \ndoi: 10.3847/1538-4357/ad2988  \nFoxley-Marrable, M., Collett, T. E., Vernardos, G., Goldstein, D. A., & Bacon, D. 2018, MNRAS, 478, 5081, doi: 10.1093/mnras/sty1346  \nGrier, C. J., Trump, J. R., Shen, Y., et al. 2017, ApJ, 851, 21, doi: 10.3847/1538-4357/aa98dc  \nGrissom, R. L., Ballantyne, D. R., & Wise, J. H. 2014, A&A, 561, A90, doi: 10.1051/0004-6361/201322637  \nSluse, D., & Tewes, M. 2014, Astronomy & Astrophysics, 571, A60  \nSluse, D., Claeskens, J.-F., Altieri, B., et al. 2006, A&A, 449, 539-550, doi: 10.1051/0004-6361:20053148  \nSluse, D., Schmidt, R., Courbin, F., et al. 2011, A&A, 528, A100, doi: 10.1051/0004-6361/201016110  \nUrry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630  \nVaswani, A., Shazeer, N., Parmar, N., et al. 2023, Attention Is All You Need. https://arxiv.org/abs/1706.03762  \nVernardos, G. 2019, MNRAS, 483, 5583,  \ndoi: 10.1093/mnras/sty3486  \nVernardos, G., Fluke, C. J., Bate, N. F., & Croton, D. 2014, The Astrophysical Journal Supplement Series, 211, 16, doi: 10.1088/0067-0049/211/1/16 Vernardos, G., & Tsagkatakis, G. 2019, Monthly Notices of the Royal Astronomical Society, 486, 1944, doi: 10.1093/mnras/stz868 Wong, K. C., Suyu, S. H., Chen, G. C. F., et al. 2020, MNRAS, 498, 1420, doi: 10.1093/mnras/stz3094 Woźniak, P. R., Udalski, A., Szymański, M., et al. 2000, The Astrophysical Journal, 540, L65, doi: 10.1086/312874 Wyithe, J. S. B., Webster, R. L., Turner, E. L., & Agol, E. 2000, Monthly Notices of the Royal Astronomical Society, 318, 1105, doi: 10.1046/j.1365-8711.2000.03746.x Yan, C.-S., Lu, Y., Yu, Q., Mao, S., & Wambsganss, J. 2014, ApJ, 784, 100, doi: 10.1088/0004-637X/784/2/100"}