Accretion Disk Size and Updated Time-delay Measurements in the Gravitationally Lensed Quasar SDSS J165043.44+425149.3

We analyze variability in 15-season optical lightcurves from the doubly imaged lensed quasar SDSS J165043.44+425149.3 (SDSS1650), comprising five seasons of monitoring data from the Maidanak Observatory (277 nights in total, including the two seasons of data previously presented in Vuissoz et al.), five seasons of overlapping data from the Mercator telescope (269 nights), and 12 seasons of monitoring data from the US Naval Observatory, Flagstaff Station at lower cadence (80 nights). We update the 2007 time-delay measurement for SDSS1650 with these new data, finding a time delay of ΔtAB=−55.1−3.7+4.0 days, with image A leading image B. We analyze the microlensing variability in these lightcurves using a Bayesian Monte Carlo technique to yield measurements of the size of the accretion disk at λ rest = 2420 Å, finding a half-light radius of log(r 1/2/cm) = 16.19−0.58+0.38 assuming a 60° inclination angle. This result is unchanged if we model 30% flux contamination from the broad-line region. We use the width of the Mg ii line in the existing Sloan Digital Sky Survey spectra to estimate the mass of this system’s supermassive black hole, finding M BH = 2.47 × 109 M ⊙. We confirm that the accretion disk size in this system, whose black hole mass is on the very high end of the M BH scale, is fully consistent with the existing quasar accretion disk size–black hole mass relation.


Introduction
At present there are only two observational techniques capable of measuring the optical/UV continuum source structure in quasars: the analysis of microlensing variability in gravitationally lensed quasar systems and reverberation mapping (e.g., Edelson et al. 2015;Cackett et al. 2018).Two somewhat different microlensing analysis techniques have been developed to make such measurements, broadly categorized as the single-epoch and multiepoch methods.The single-epoch technique requires simultaneous multiband images of a lensed quasar to determine whether deviations exist between the macroscopic lens models predictions of the flux ratios and the observed values and if the flux ratios are significantly dependent on wavelength (Pooley et al. 2007;Blackburne et al. 2011;Mosquera & Kochanek 2011).While the singleepoch technique is observationally efficient, its precision is limited by an unknown statistical prior, and its applicability is limited by the geometry of the lens.The (unknown) average mass of a lens galaxy star must be assumed, and the singleepoch technique is only appropriate for short time-delay systems because in systems with longer time delays, differentiating between the variability intrinsic to the quasar source itself and variability from extrinsic effects such as microlensing is not possible.The multiepoch technique requires the collection and analysis of many-season lightcurves, modeling the observed extrinsic variability resulting from the complicated magnification network produced by the lens galaxy stars' combined gravitational potential (see, e.g., Kochanek et al. 2006;Dai et al. 2010;Hainline et al. 2012Hainline et al. , 2013;;Cornachione et al. 2020, who utilize the methods of Kochanek 2004).This technique necessitates cosmological modeling of the effective transverse velocity as well as larger amounts of observing time; however, it does not necessitate the assumption of the mean mass of a lens galaxy star.
In this investigation, we present a multiepoch microlensing analysis of 15-season optical lightcurves from the doubly imaged gravitationally lensed quasar SDSS J165043.44 +425149.3(hereafter SDSS1650), which has an image separation of 1 2. The source redshift is z s = 1.547 and the lens redshift is z l = 0.58 (Morgan et al. 2003).The existing time-delay measurement, 49.5 ± 1.9 days (Vuissoz et al. 2007), was made using only two seasons of monitoring data; we use our significantly longer lightcurves to update this measurement.Unfortunately, this system has never been observed with the Hubble Space Telescope (HST), so we are forced to use ground-based imagery and the measured time delay to constrain the lens galaxy morphology.
In Section 2 we describe the new r-band observations of SDSS1650, and in Section 3 we detail our new time-delay measurement for SDSS1650 including our lightcurves, which have 13 yr of extended monitoring beyond the Vuissoz et al. (2007) data set.We describe our strong lensing model in Section 4, and in Section 5 we discuss our microlensing simulations.The results of the simulations are detailed in Section 6, and in Section 7 we discuss and summarize our findings.We adopt the following cosmological parameters for this work: Ω M = 0.3, Ω λ = 0.7, and H 0 = 70 km s −1 Mpc −1 (Hinshaw et al. 2009).

Observations and Data Reduction
We compiled a 626-epoch Sloan Digital Sky Survey (SDSS) r-band lightcurve using observations from several sources.We include 277 nights of data from Maidanak (see Ehgamberdiev 2018 for a description of the observatory and its suitability for observing gravitational lenses May and 2016 July, using three 600 s exposures in both cases.The pixel scale for OneChip is 0 18 and that of Tek2k is 0 33.A summary of all observations is provided in Table 1. We stack the three subexposures to create a high signal-tonoise ratio frame at each epoch, discarding any subexposures contaminated by satellite trails or telescope jitter, etc., leaving 80 nights of usable observations.Four reference stars were used to identify the point-spread function (PSF) and flux normalization for every image.The lens galaxy is too faint to be detected in the r band (Vuissoz et al. 2007;Mosquera et al. 2011) so we do not attempt to subtract any lens galaxy flux from the PSF fits of the images.Because we have data from different detectors (for which the same reference stars were not used), we normalized the magnitude scale to that from the Maidanak lightcurve.The USNO magnitude offsets are −13.312and −11.519 for images A and B, respectively, and for Mercator, the offsets were −0.013 and 0.010.An example image for OneChip with the reference stars labeled is shown in Figure 1.There were 23 nights of data obtained with OneChip, and 62 nights with Tek2k.We provide the reduced lightcurves in Table 2, and in Figure 2 we show the total combined lightcurves, with 626 epochs of data extending over 15 seasons.(This table is available in its entirety in machine-readable form.)

Time-delay Measurement
Before performing a microlensing analysis it is imperative to identify and eliminate variability intrinsic to the quasar source itself, which requires the shifting of the lightcurves by the time delay.The previous measurement of the time delay for SDSS1650 was performed in Vuissoz et al. (2007) with polynomial fitting and minimum dispersion methods with only two seasons of data.Both methods generated small error bars, as the polynomial fit utilized models for the intrinsic and extrinsic variability, which had a low number of degrees of freedom, and the minimum dispersion technique simply does not account for microlensing.
We update the time-delay measurement for SDSS1650 with our new data (15 seasons), finding a time delay of 4.0 days.This estimate is obtained by running the free-knot spline method implemented in PyCS3 (Millon et al. 2020b) and described in detail in Millon et al. (2020a).We applied this technique to the entirety of our data from the Mercator and Maidanak observatories, but we chose to only include USNO observations up to HJD 2457922 in the time delay analysis.The infrequent USNO observations following that season require excessive interpolation and add no meaningful statistical weight to the measurement.This low cadence is not a problem for the subsequent microlensing analysis, however, so we analyze the full lightcurves in Section 5. We used the following mean spacing between the knots: η ä [35,45,55,65] for the intrinsic variability and η ml ä [300, 600] for the microlensing variability.The estimates resulting from these different choices of model hyperparameters are combined in the same procedure as the one described in Millon et al. (2020a).As a consistency test, we run our time-delay estimation pipeline on the individual Maidanak, Mercator, and Maidanak+Mercator data sets.We obtain estimates consistent with our baseline estimate: Δt AB,Mercator = - 3.8 .The precision of our time-delay result appears to be somewhat worse than, though consistent with, the Vuissoz et al. (2007) measurement.However, it is likely that Vuissoz et al. (2007) underestimated the uncertainties in their two-season measurement.The minimum dispersion techniques and polynomial fitting methods used in this previous work either did not account for the additional uncertainties introduced by the microlensing variability or modeled it with a simple linear model.We tested this hypothesis by applying our free-knot spline method to the lightcurves from Vuissoz et al. (2007).Analyzing only the first two seasons using the free-knot spline technique increases the uncertainties in the time-delay measurement by a factor of ∼5.7× relative to those reported by Vuissoz et al. (2007).
Here, we implement flexible models for both the quasar intrinsic variability and the microlensing variability, yielding more robust but slightly less precise time-delay measurements (Bonvin et al. 2016).An illustration of the free-knot spline fits to the data is shown in Figure 3.Both the extrinsic microlensing variability and the intrinsic quasar variability are modeled with free-knot splines; we simultaneously optimize for the time delay.

Models for Strong Lensing
We utilized the GRAVLENS software (Keeton 2001) to model the lens galaxy mass profile.We modeled SDSS1650 with two components: a Navarro, Frenk, and White (NFW; Navarro et al. 1997) profile for the dark matter halo and a deVaucouleurs profile for the baryonic component.Because the fraction of baryonic mass to dark matter is not known a priori, we made a ten-model sequence for which 0.1-1.0 of the total mass is produced by a baryonic component, with steps of 0.1 (the parameter f M/L ).With these models we cover the plausible range of stellar mass fractions (κ * /κ) and shear; both variables are demonstrated by Schechter & Wambsganss (2002) and Vernardos & Fluke (2014) to have a strong influence on microlensing statistics.The two components were forced to be concentric; however, we permitted the position of the galaxy to float within the range of the observational uncertainty.
Since this system has never been observed with HST, there are no observational constraints on the lens galaxy morphology, forcing us to model the DeVaucouleurs and NFW components as spheres.This limitation affected the precision of our final results, but it did not preclude our analysis from converging.Additionally, our f M/L series allows us to cover a wide range of likely lens galaxy parameters.Following Vuissoz et al. (2007), we used H 0 , the positions of images A and B, and the measured time delay to constrain the effective radius of the lens galaxy to 0 8.The convergence (κ), stellar convergence fraction (κ * /κ), and shear (γ) at the location of each image are provided in Table 3.

Magnification Patterns
Following Kochanek (2004) we made use of the inverse rayshooting technique of Wambsganss et al. (1992) to construct our magnification patterns.For each magnification pattern, the stellar mass ranges were consistent with the Gould (2000) Galactic initial mass function extending with a - dN dM M 1.3 over a maximum to minimum ratio of 50.The average mass of a star within the lens galaxy (the mean microlens mass) 〈M * /M e 〉 is inferred during the Monte Carlo simulations.Each magnification Figure 3.In the top panel, we display the results of our free-knot spline fitting technique using PyCS as applied to the lightcurves of SDSS1650.The black curve is fit to the intrinsic variability of the source, and the black ticks on the curve are the final positions of the spline knots after running the PyCS analysis.The spline knots are given a uniform initial spacing, but they are free to move during the analysis.In the gaps in the lightcurve where the data do not provide any constraints on the knot positions, the knots maintain their initial spacing and we enforce linear interpolation between them.The orange (image A) and blue (image B) points are the data shifted by the best-fit time delay.The orange and blue curves are the models for the microlensing variability in images A and B, respectively.Fit residuals are shown in the bottom two panels.
where θ E is proportional to 〈M * /M e 〉 1/2 .In this system, the Einstein radius R E of a 1M e star is 4.54 × 10 16 cm.The pixel scale (2.22 × 10 14 〈M * /M e 〉 1/2 cm) was selected to allow for resolution of the minimum expected accretion disk size and the map dimensions were selected to permit constraint of our larger source sizes over the long time duration of our lightcurves.In total, 40 sets of magnification patterns per each image per macro model (800 magnification patterns in total) were generated, which removes any possible influence due to systematics in the magnification patterns.

Microlensing
Using linear interpolation, we shifted the lightcurve for image A to align with observation epochs in curve B. The two lightcurves are then divided to produce a "difference lightcurve" in which the remaining variability is only extrinsic (viz.microlensing).This is shown by the blue points in Figure 4. We use the Monte Carlo method described in Kochanek (2004) to retrieve the values of our physical parameters, which were the most likely to be capable of reproducing the observed microlensing variability in our lightcurves.We first convolve the magnification patterns with a Gaussian kernel over a source size range with 17 equal intervals in log space over log r cm s ( ˆ) = [14.5,18.5] 〈M * /M e 〉 1/2 cm.It has been determined by Mortonson et al. (2005) and Vernardos & Tsagkatakis (2019) that the half-light radius r 1/2 is what primarily sets the strength of microlensing fluctuations, not the particular surface brightness profile used; therefore, for computational efficiency, we selected a Gaussian profile for the convolutions.The accretion disk scale radius r s is measured at the center of the rest-frame monitoring band (Equation ( 2) from Morgan et al. 2010): where L/L E is the Eddington fraction, η is the accretion efficiency  h = L Mc 2 , and M BH is the black hole mass.A calculation of the black hole mass for SDSS1650 is presented in Section 6.The scale radius in Einstein units r s ˆis related to the scale radius in physical units by ˆ☉ .In essence, each Monte Carlo trial attempts to recreate the observed lightcurves by moving a model quasar source across a magnification pattern.Equipped with the convolved magnification patterns we randomly select an effective velocity at a unique initial position (x i , y i , with direction θ i ) for every trajectory, with the assumption of independent and uniformly distributed values.A log-uniform prior on the effective transverse velocity in Einstein units v e ˆwas used, ranging from [10, 10 6 ] 〈M * /M e 〉 1/2 km s −1 .
We made use of the Advanced Research Computer9 highperformance cluster at the United States Naval Academy for our simulations.We identify a reduced χ 2 cut (χ 2 < 2) that produces a sufficiently large number of solutions for the Bayesian analysis, while still producing a preponderance of good fits.A total of 10,000,000 trials per magnification pattern were attempted; 2.03 × 10 8 solutions were retained following the reduced χ 2 cut.In Figure 4 we display the five best simulated difference lightcurves for SDSS1650.
We use Bayesian statistics to construct likelihood functions for our variables of interest, including the physical extent of the quasar and the transverse velocity v e ˆof the system.We construct probability densities by marginalizing over the remaining model parameters, using priors covering every parameter's expected range.For example, the probability density for the scale radius of the accretion disk in Einstein units r s ˆis found as follows, where ξ represents the additional variables, including the effective velocity v e ˆ, f M/L , and our trajectory variables (x i , y i , θ i ).The statistical priors we used for the Bayesian integration are encapsulated in the expression π(ξ).
We follow Equation (5) of Mosquera & Kochanek (2011) and the method of Kochanek (2004) to model the true effective transverse velocity of the observer, source, and lens v e , along with the observer's velocity projection v CMB across the line of sight to the quasar, the bulk motion of the lens σ pec(zl) and the source σ pec (zs) , and the velocity dispersion of the lens galaxy, σ * .The velocity dispersion, σ * , was estimated from the mass of a singular isothermal sphere lens model.The east and north components of the observer's velocity (−365.1 km s −1 and −56.9 km s −1 , respectively) were found using cosmic microwave background (CMB) dipole projection along the line of sight to SDSS1650.We estimated the components of the system's peculiar velocity σ pec (z s ) = 91.73km s −1 and σ pec (z l ) = 135 km s −1 using the models of Tinker et al. (2012).The probability density for this velocity model is shown in Figure 5. Since ˆ, we can convolve the source velocity probability density and our velocity model to derive a probability density for 〈M * /M e 〉: We display this result in Figure 6.We were unable to eliminate extremely small microlens masses from our probability distribution; however, that is not unexpected given the limited observational constraints (e.g., no HST imagery) available for the generation of our lens galaxy model.To confirm the robustness of our transverse velocity model, we also estimated the accretion disk size with the assumption of a prior on the mean microlens mass extending uniformly from 0.1 £ M * /M e £ 1.0.The results that follow are consistent with those derived using this (more physically plausible) prior on the microlens masses.

Microlensing Sizes
By convolving the probability density for the mean microlens mass with that for the accretion disk scale radius in Einstein units r s ˆ, we produce a probability density for the scale radius r s in unscaled, physical units.We display this result in Figure 7.The dashed curve in Figure 7 is the resulting probability distribution when the uniform prior on the mean microlens mass is implemented.Integrating the probability density in Figure 7 yields an expectation value for the accretion disk size at λ rest = 2420 Å of log(r 1/2 /cm) = -+ 16.19 0.58 0.38 when assuming a 60°inclination angle, where the half-light radius is related to the scale radius as r 1/2 = 2.44 r s .
ˆ.The effective velocity v e model in unscaled physical units (km s −1 ) is plotted by the thin black curve; this model is the statistical prior on the effective source velocity.One concern regarding the accuracy of this result is the possibility of contaminating (unmicrolensed) flux from the broad-line region (BLR) included in the lightcurves affecting our measurement of the accretion disk size (e.g., Fian et al. 2018Fian et al. , 2021;;Paic et al. 2022).When we tested for this possibility by rerunning the simulations including 30% flux contamination from the BLR, our results remained consistent.See Rivera et al. (2023) for further discussion on the likelihood and limitations of BLR contamination biasing multiepoch microlensing results.

Black Hole Mass
We used the PyQSOFit code (Guo et al. 2018) to measure the FWHM of the Mg II line using SDSS spectra taken on 2001 June 19 (Richards et al. 2002).The spectral fits are shown in Figure 8.The ground-based Apache Point Observatory 2.5 m lacks the angular resolution to measure each spectrum separately, so the SDSS spectra include the magnified flux from both quasar images.To find the unmagnified flux, we divided the measured flux by the sum of the magnifications of the two images yielded by the macroscopic mass model, which is consistent with the measured time delay.We measured the FWHM of Mg II and used the measurement of the continuum luminosity at λ rest = 3000 Å (5953.65 km s −1 and 8.73 × 10 45 erg s −1 Å −1 , respectively) to calculate the black hole mass using the fitting form:

Discussion and Summary
In Figure 9 we display the quasar accretion disk size-black hole mass relation (Morgan et al. 2010(Morgan et al. , 2018) ) with the present SDSS1650 measurement shown in blue and updated measurements included for the quasars SBS0909+532, Q0957+561, and FBQ0951+2635 using the values from Cornachione & Morgan (2020) and Rivera et al. (2023).Not only is our new SDSS1650 measurement fully consistent with the relation, but it represents the most massive black hole studied with the multiepoch microlensing technique.In Figure 9 we also display a size estimate for the accretion disk based on the observed Figure 8. Line fits of the SDSS spectra for SDSS1650 produced using PyQSOFit (Guo et al. 2018).The data are shown in black, the overall fit is shown in blue, the continuum is in orange, and the red lines show the individual components for each line.The gray bars shown in the top panel show the regions used to calculate the continuum fit.The bottom three panels show the fits to the C IV, C III, and Mg II lines (from left to right, respectively).9.The accretion disk size-black hole mass relation.Microlensing (black crosses) and luminosity (gray x's) size measurements vs. black hole mass for a subset of gravitationally lensed quasars.The values for SDSS1650 are in blue and have been adjusted to 2500 Å for comparison with the sample from Morgan et al. (2010).Both the microlensing and luminosity size measurements for this quasar are consistent with the expected relation (the solid black line gives the microlensing relationship, and the black dashed line gives that for the luminosity sizes) within the expected errors.
) taken from 2004 May to 2008 October.The data from 2004 May to 2005 September were previously published in Vuissoz et al. (2007).For most of the epochs, six dithered images were taken, each with an exposure time of 300 s.The pixel scale for the CCD was 0 266.An additional 269 observations from Mercator were taken between 2004 August and 2008 October.Most of the data consist of five exposures per night, each with an exposure time of 360 s.The pixel scale of the CCD was 0 193.Our new observations were taken using the 1.55 m Kaj Strand Astrometric Reflector at the US Naval Observatory Flagstaff Station (NOFS) with both the 2048 × 4096 EEV (OneChip) and 2048 × 2048 Tek2k CCDs.We took observations with Tek2k from 2008 April to 2022 June and OneChip observations between 2010

Figure 1 .
Figure 1.Example of a night of OneChip data taken with the 1.55 m Kaj Strand Astrometric Reflector at NOFS on 2010 May 8.The field of view is ¢ ´¢ 2.6 3.6.The reference stars used in the analysis are shown.

Figure 2 .
Figure 2. r-band lightcurves for SDSS1650.The lightcurve for image A is shown in the top panel; the lightcurve for image B is shown in the bottom panel.The yellow (OneChip) and purple (Tek2k) data were observed with NOFS, the orange data with Maidanak, and the pink data with Mercator.

Figure 4 .
Figure 4.The best five model lightcurves from our Monte Carlo simulation (black solid lines) are overlaid against the observed r-band difference lightcurves (blue error bars) for SDSS1650.

Figure 5 .
Figure 5. Probability density for the effective transverse velocity.The thick curve is the posterior probability density for the effective velocity in Einstein units of 

Figure 6 .
Figure 6.Posterior probability density for the mean mass of a microlens 〈M * 〉 in the lens galaxy.

Figure 7 .
Figure7.Probability density for the source size in centimeters at λ rest = 2420 Å.The dashed curve is the probability density yielded by the application of a uniform mass prior extending over [0.1, 1]M e , and the solid curve is the probability density arising instead from the prior on the effective velocity.The vertical lines show the Schwarzchild radius and the innermost stable circular orbit (ISCO) for our calculated black hole mass.

Table 1
Summary of SDSS r-band Monitoring Data Taken with the Maidanak Observatory 1.5 m Telescope, Uzbekistan, the Mercator Observatory 1.2 m Telescope, La Palma, and the USNO 1.55 m Kaj Strand Telescope,