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

HAL
open science

# Dust Reverberation Mapping in Distant Quasars from Optical and Mid-infrared Imaging Surveys

Qian Yang, Yue Shen, Xin Liu, Michel Aguena, James Annis, Santiago Avila, Manda Banerji, Emmanuel Bertin, David Brooks, David Burke, et al.

► To cite this version:

Qian Yang, Yue Shen, Xin Liu, Michel Aguena, James Annis, et al. Dust Reverberation Mapping in Distant Quasars from Optical and Mid-infrared Imaging Surveys. The Astrophysical Journal, 2020, 900, 10.3847/1538-4357/aba59b. insu-03748046

HAL Id: insu-03748046

https://insu.hal.science/insu-03748046v1

Submitted on 10 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Dust Reverberation Mapping in Distant Quasars from Optical and Mid-infrared Imaging Surveys

Qian Yang<fnref n="1" />, Yue Shen<fnref n="1" />,<fnref n="2" />,<fnref n="40" /> ID, Xin Liu<fnref n="1" />,<fnref n="2" /> ID, Michel Aguena<fnref n="3" />,<fnref n="4" />, James Annis<fnref n="5" /> ID, Santiago Avila<fnref n="6" /> ID, Manda Banerji<fnref n="7" />,<fnref n="8" />, Emmanuel Bertin<fnref n="9" />,<fnref n="10" /> ID, David Brooks<fnref n="11" /> ID, David Burke<fnref n="12" />,<fnref n="13" /> ID, Aurelio Carnero Rosell<fnref n="14" />,<fnref n="15" /> ID, Matias Carrasco Kind<fnref n="1" />,<fnref n="2" /> ID, Luiz da Costa<fnref n="4" />,<fnref n="16" /> ID, Juan De Vicente<fnref n="17" /> ID, Shantanu Desai<fnref n="18" /> ID, H. Thomas Diehl<fnref n="5" /> ID, Peter Doel<fnref n="11" />, Brenna Flaugher<fnref n="5" /> ID, Pablo Fosalba<fnref n="19" /> ID, Josh Frieman<fnref n="5" />,<fnref n="20" /> ID, Juan Garcia-Bellido<fnref n="6" /> ID, David Gerdes<fnref n="21" />,<fnref n="22" /> ID, Daniel Gruen<fnref n="12" />,<fnref n="13" />,<fnref n="23" /> ID, Robert Gruendl<fnref n="1" />,<fnref n="2" /> ID, Julia Gschwend<fnref n="4" />,<fnref n="16" /> ID, Gaston Gutierrez<fnref n="5" /> ID, Samuel Hinton<fnref n="24" /> ID, Devon L. Hollowood<fnref n="25" /> ID, Klaus Honscheid<fnref n="26" />,<fnref n="27" /> ID, Nikolay Kuropatkin<fnref n="5" /> ID, Marcio Maia<fnref n="4" />,<fnref n="16" /> ID, Marisa March<fnref n="28" />, Jennifer Marshall<fnref n="29" /> ID, Paul Martini<fnref n="26" />,<fnref n="27" /> ID, Peter Melchior<fnref n="30" /> ID, Felipe Menanteau<fnref n="1" />,<fnref n="2" /> ID, Ramon Miquel<fnref n="31" />,<fnref n="32" /> ID, Francisco Paz-Chinchon<fnref n="2" />,<fnref n="7" /> ID, Andrés Plazas Malagón<fnref n="30" /> ID, Kathy Romer<fnref n="33" /> ID, Eusebio Sanchez<fnref n="17" /> ID, Vic Scarpine<fnref n="5" /> ID, Michael Schubnell<fnref n="22" /> ID, Santiago Serrano<fnref n="34" />,<fnref n="35" /> ID, Ignacio Sevilla<fnref n="17" /> ID, Mathew Smith<fnref n="36" /> ID, Eric Suchyta<fnref n="37" /> ID, Gregory Tarle<fnref n="22" /> ID, Tamas Norbert Varga<fnref n="38" />,<fnref n="39" /> ID, and Reese Wilkinson<fnref n="33" />

## Abstract

The size of the dust torus in active galactic nuclei (AGNs) and their high-luminosity counterparts, quasars, can be inferred from the time delay between UV/optical accretion disk continuum variability and the response in the mid-infrared (MIR) torus emission. This dust reverberation mapping (RM) technique has been successfully applied to ~70 z < 0.3 AGNs and quasars. Here we present first results of our dust RM program for distant quasars covered in the Sloan Digital Sky Survey Stripe 82 region combining ~20 yr ground-based optical light curves with 10 yr MIR light curves from the WISE satellite. We measure a high-fidelity lag between W1 band (3.4 µm) and g band for 587 quasars over 0.3 < z < 2 (⟨z⟩ ~ 0.8) and two orders of magnitude in quasar luminosity. They tightly follow (intrinsic scatter ~0.17 dex in lag) the IR lag–luminosity relation observed for z < 0.3 AGNs, revealing a remarkable size–luminosity relation for the dust torus over more than four decades in AGN luminosity, with little

<footnote n="1">¹ Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA; qany@illinois.edu, shenyue@illinois.edu</footnote>

<footnote n="2">² National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA</footnote>

<footnote n="3">³ Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, CP 66318, São Paulo, SP, 05314-970, Brazil</footnote>

<footnote n="4">⁴ Laboratório Interinstucional de e-Astronomia—LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil</footnote>

<footnote n="5">⁵ Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA</footnote>

<footnote n="6">⁶ Instituto de Física Teorica UAM/CSIC, Universidad Autonoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="7">⁷ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="8">⁸ Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="9">⁹ CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France</footnote>

<footnote n="10">¹⁰ Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France</footnote>

<footnote n="11">¹¹ Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="12">¹² Kavli Institute for Particle Astrophysics & Cosmology, P.O. Box 2450, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="13">¹³ SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA</footnote>

<footnote n="14">¹⁴ Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="15">¹⁵ Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="16">¹⁶ Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ—20921-400, Brazil</footnote>

<footnote n="17">¹⁷ Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain</footnote>

<footnote n="18">¹⁸ Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India</footnote>

<footnote n="19">¹⁹ Institut d'Estudis Espacials de Catalunya (IEEC), E-08193 Barcelona, Spain</footnote>

<footnote n="20">²⁰ Kayli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="21">²¹ Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="22">²² Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="23">²³ Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA</footnote>

<footnote n="24">²⁴ School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia</footnote>

<footnote n="25">²⁵ Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA</footnote>

<footnote n="26">²⁶ Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="27">²⁷ Department of Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="28">²⁸ Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA</footnote>

<footnote n="29">²⁹ George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA</footnote>

<footnote n="30">³⁰ Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA</footnote>

<footnote n="31">³¹ Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain</footnote>

<footnote n="32">³² Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona), Spain</footnote>

<footnote n="33">³³ Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton BN1 9QH, UK</footnote>

<footnote n="34">³⁴ Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain</footnote>

<footnote n="35">³⁵ Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain</footnote>

<footnote n="36">³⁶ School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, UK</footnote>

<footnote n="37">³⁷ Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831</footnote>

<footnote n="38">³⁸ Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, D-85748 Garching, Germany</footnote>

<footnote n="39">³⁹ Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians Universität München, Scheinerstr. 1, D-81679 München, Germany</footnote>

<footnote>Received 2020 June 18; revised 2020 June 26; accepted 2020 July 5; published 2020 September 1</footnote>

<footnote n="40">⁴⁰ Alfred P. Sloan Research Fellow.</footnote>
dependence on additional quasar properties such as Eddington ratio and variability amplitude. This study motivates further investigations in the utility of dust RM for cosmology and strongly endorses a compelling science case for the combined 10 yr Vera C. Rubin Observatory Legacy Survey of Space and Time (optical) and 5 yr Nancy Grace Roman Space Telescope 2 µm light curves in a deep survey for low-redshift AGN dust RM with much lower luminosities and shorter, measurable IR lags. The compiled optical and MIR light curves for 7384 quasars in our parent sample are made public with this work.

*Unified Astronomy Thesaurus concepts: Active galactic nuclei (16); Quasars (1319); Dust continuum emission (412); Reverberation mapping (2019)*

*Supporting material:* figure set, machine-readable tables

# 1. Introduction

In the widely accepted unified model of active galactic nuclei (AGNs),<fnref n="41" /> the inner regions surrounding the accreting super-massive black hole (SMBH) include (roughly from inside out but with potential spatial overlaps) a hot X-ray-emitting corona, an accretion disk, a broad emission line region, a dusty toroidal structure, and a narrow emission line region. The toroidal dusty region, commonly referred to as the dust torus, plays the central role in the AGN unification scheme (Antonucci 1993; Urry & Padovani 1995) that unifies type 1 (unobscured, broad-line) and type 2 (obscured, narrow-line) AGNs.

The spectral energy distributions (SEDs) of AGNs show a pronounced peak in the mid-infrared (MIR; e.g., Sanders et al. 1989; Elvis et al. 1994), which is interpreted as thermal emission from hot dust in the torus region. The dust grains in the torus absorb the incident UV/optical continuum emission from the accretion disk and reradiate in the infrared (Rieke 1978), with a time lag corresponding to the average light-travel time from the accretion disk to the dust torus. The torus extends from the dust sublimation radius outward (Barvainis 1987), with the near-infrared (NIR) radiation arising at the inner edge near the sublimation radius and longer wavelength radiation arising from the outer regions with lower dust temperatures (e.g., Nenkova et al. 2008; Netzer 2015, and references therein). Figure 1 shows a schematic for the broadband SED of unobscured broad-line AGNs.

The compact size (~subparsec to parsec for typical quasars) of the torus is difficult to spatially resolve directly. There are only a handful of nearby bright AGNs for which we can marginally resolve the torus structure using IR interferometric techniques (e.g., Swain et al. 2003; Jaffe et al. 2004; Tristram et al. 2007; Kishimoto et al. 2009, 2011; Weigelt et al. 2012; Burtscher et al. 2013; GRAVITY Collaboration et al. 2019). Beyond the nearby universe it becomes difficult for IR interferometry to directly resolve the torus owing to the reduced angular size and brightness of the distant source.

Dust reverberation mapping (RM) offers an alternative route to infer the size of the AGN torus without the need for spatial resolution. The echo of reprocessed dust emission to UV/optical continuum variations measures the average light-crossing time (hence a typical size) of the dust torus to the central engine. The dust RM lag has been measured between the UV/optical and the NIR (e.g., Penston et al. 1971; Clavel et al. 1989; Glass 1992, 2004; Oknyanskii 1993; Oknyanskij et al. 1999; Nelson 1996; Minezaki et al. 2004, 2019; Suganuma et al. 2006; Koshida et al. 2009, 2014), or UV/optical versus MIR (Vazquez et al. 2015; Lyu et al. 2019).

These previous dust RM studies have confirmed correlated variability between the optical and NIR (K-band) emission, justifying the RM technique. The dust torus radius, *R*, inferred from the lag between the optical and NIR variability, was found to scale with the AGN luminosity, *L*, as $R \propto L^{1/2}$ (e.g., Suganuma et al. 2006; Koshida et al. 2014). Comparisons between the dust RM lags and direct interferometric measurements found general agreement between the size measurements from the two methods (e.g., Kishimoto et al. 2011).

To date, most of the dust RM works on torus size measurements were limited to relatively low redshift ($z < 0.3$). It is important to extend dust RM to high-redshift quasars to better sample the high-luminosity end of the $R-L$ relation and to investigate any redshift evolution in the physical properties (such as size and grain physics) of the dust torus. At higher redshifts, we can only observe luminous quasars, whose observed-frame torus lags are longer than their low-$z$ counterparts owing to both the higher luminosity and the cosmic time dilation. Therefore, a sufficiently long time baseline is required to measure the dust RM lag for these distant and luminous quasars.

In this paper, we present dust RM measurements for a large sample of quasars at a median redshift of $\langle z \rangle \sim 0.8$. We make use of the multiepoch MIR imaging from the all-sky Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010; Mainzer et al. 2014), combined with optical light curves from all available ground-based facilities, including multiepoch data from the Dark Energy Survey (DES; Abbott et al. 2018), the Sloan Digital Sky Survey (SDSS; York et al. 2000), Pan-STARRS (PS1; Chambers et al. 2016), the Catalina Real-time Transient Survey (CRTS; Drake et al. 2009), the Palomar Transient Factory (PTF; Law et al. 2009), the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014), and the Zwicky Transient Facility (ZTF; Bellm et al. 2019). The WISE MIR data cover a baseline of $\sim 10$ yr with a cadence of $\sim 6$ months, and the optical data cover a combined baseline of $\sim 20$ yr. These multiyear baselines provide a unique opportunity to measure the dust echoes in luminous and distant quasars. Lyu et al. (2019) already demonstrated the power of combining WISE light curves with ground-based optical light curves in measuring dust echoes in luminous Palomar-Green (PG) quasars at $z < 0.5$. Here we extend this exercise to even higher redshifts with SDSS quasars. At the median redshift of our sample $\langle z \rangle \sim 0.8$, WISE $W1$ (3.4 µm) and $W2$ (4.6 µm) data mainly probe the dust emission at rest frame 2 µm, allowing us to directly compare the MIR lags of these quasars with the NIR (mostly K-band) lags in nearby AGNs (e.g., Suganuma et al. 2006; Koshida et al. 2014).

The structure of this paper is as follows. In Section 2, we describe our quasar sample and the photometric data (in MIR and optical). We describe the dust RM measurements in

<footnote n="41">⁴¹ In this paper AGNs refer to SMBHs that are accreting efficiently, with a standard optically thick, geometrically thin accretion disk (Shakura & Sunyaev 1973).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic representation of the broadband SED of an unobscured broad-line quasar, including contributions from the accretion disk (Shakura & Sunyaev 1973), dusty torus (e.g., Hönig et al. 2006; Nenkova et al. 2008), stellar emission (Bruzual & Charlot 2003), hot corona (Haardt & Maraschi 1991), and soft excess (e.g., Done et al. 2012 and references therein). The upper right inset shows the rest-frame wavelengths in the W1 and W2 bands as a function of redshift. This schematic is for demonstration purposes, and the relative contributions are not exact and vary from object to object.</figcaption></figure>

Section 3. We present the results of MIR lags in Section 4, and we discuss the relation between the dust radius and quasar luminosity in Section 5. We conclude in Section 6 with an outlook for future work. Throughout this paper, we adopt a flat $\Lambda$CDM cosmology with parameters $\Omega_{\Lambda} = 0.7$, $\Omega_{m} = 0.3$, and $H_{0} = 70$ km s$^{-1}$ Mpc$^{-1}$. All quoted uncertainties are $1\sigma$.

## 2. Sample and Data

### 2.1. The S82 Quasar Sample

Our parent quasar sample includes the 9258 spectroscopically confirmed broad-line quasars in the Stripe 82 (S82) region (MacLeod et al. 2012) that are included in the SDSS DR7 quasar catalog (Shen et al. 2011). There are several advantages of this sample for our dust RM study for high-z quasars: (1) these quasars cover a broad range in redshift and quasar luminosity and are representative of the luminous quasar population at high redshift; (2) they are bright enough for reliable photometric measurements from ground-based imaging surveys with small-aperture telescopes and from WISE; (3) they have ~20 yr of optical light curves combining all available photometry from various surveys described in Section 1; and (4) they have well measured spectroscopic properties, such as BH mass and Eddington ratios ($L/L_{Edd}$) from the Shen et al. (2011) catalog.

Our targets well sample the high-luminosity regime in the torus size–luminosity relation, compared with earlier dust RM measurements in low-redshift AGNs and PG quasars (e.g., Suganuma et al. 2006; Koshida et al. 2014; Lyu et al. 2019). Many of these S82 quasars have observed-frame MIR lags of more than a few years; therefore, it is necessary to have decade-long light curves to meaningfully measure the lag. There are still selection biases for the highest-redshift/luminosity quasars owing to the duration of our light curves, which will be discussed in Section 5.

### 2.2. WISE Light Curves

WISE scanned the full sky from 2010 January to July in four bands centered at wavelengths of 3.4, 4.6, 12, and 22 µm (W1, W2, W3, and W4). The secondary cryogen survey and Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE; Mainzer et al. 2011) Post-Cryogenic Mission mapped the sky from 2010 August to 2011 February. The NEOWISE Reactivation Mission (NEOWISE-R; Mainzer et al. 2014) surveys the sky in the W1 and W2 bands from 2013 twice a year. WISE obtains ~10–20 observations within a 36 hr window in each visit. We calculate the median magnitude and magnitude error, specifically the semiamplitude of the range enclosing the 16th and 84th percentiles of all flux measurements within a 6-month window. We limit to good-quality single-epoch data points with the best frame image quality score ($qi\_fact = 1$), observed far away from the South Atlantic Anomaly ($saa\_sep \ge 5$), with no contamination from the Moon ($moon\_masked = 0$), and excluding spurious detection ($cc\_flags = 0$). The WISE magnitudes are profile-fitting magnitudes and are converted from Vega to AB magnitude as $m_{AB} = m_{Vega} + \Delta m$, where $\Delta m$ is 2.699, 3.339, 5.174, and 6.620 in the W1, W2, W3, and W4 bands, respectively (Jarrett et al. 2011).

We extract the WISE light curves from the latest observations up to 2019 December 13 (released on 2020 March 26) for 9258 S82 quasars, using a matching radius of 2". Since W1 is the most sensitive band with the highest fraction of WISE detection of S82 quasars, we focus on the W1-band MIR light curves in this work. The results of multiband (W1 and W2) MIR lags and implications on dust torus properties (such as the radial temperature profile) will be presented in a future paper.

### 2.3. Optical Light Curves

We compile all available optical photometric data from various ground-based imaging surveys that cover the S82 region, including SDSS, PS1, DES, CRTS, ASAS-SN, PTF, and ZTF (see Table 1). These ground-based optical surveys cover different epochs and have different bandpasses and depths. We homogenize these optical light curves as described in detail in Section 2.4. Table 1 summarizes the characteristics of these optical surveys. The default optical magnitude type is the PSF magnitude.

Since WISE spans from 2010 to 2019, DES data (covering ~2013–2018) are crucial for measuring the MIR time lag. We therefore restrict our analysis to the subset of 7582 quasars from the S82 sample located at R.A. < 46 deg or R.A. > 316 deg within the footprint of the DES wide survey. A total of 7384 of these quasars have available WISE light curves in W1. We summarize the sample statistics in Table 2.

### 2.4. Optical Photometric Calibration

To calibrate the optical data from different surveys onto the same flux scale, we apply additive corrections to the optical magnitudes taking into account differences in filter curves and reported photometric magnitudes. To correct for different filter curves, we convolve the SDSS spectrum with the PS1/DES/ZTF filter curves to obtain synthetic magnitudes and compare to those derived with the SDSS filters to derive the corrections. CRTS data are observed through an unfiltered wide band, so we apply a constant offset to the CRTS magnitudes to match
<table><caption>**Table 1**
Survey Information</caption><thead><tr><th>Survey</th><th>Filter</th><th>Time</th><th>Cadence</th><th>N<sub>epoch</sub></th><th>Coverage</th><th>Depth</th></tr></thead><tbody><tr><td>SDSS(S82)</td><td>g</td><td>1998–2007</td><td>5 days</td><td>~60</td><td>300 deg<sup>2</sup></td><td>22.2</td></tr><tr><td>PS1</td><td>g (r, i)</td><td>2011–2014</td><td>2/season</td><td>~10</td><td>3π</td><td>22.0</td></tr><tr><td>DES</td><td>g (r, i)</td><td>2013–2018</td><td>1–4/season (wide-field)</td><td>~10</td><td>5100 deg<sup>2</sup></td><td>23.57</td></tr><tr><td>PTF</td><td>g</td><td>2009–2014</td><td>5 days</td><td>~4</td><td>11,233 deg<sup>2</sup></td><td>19</td></tr><tr><td>ZTF</td><td>g</td><td>2018</td><td>3 days</td><td>~30</td><td>2.5–3π</td><td>20.5</td></tr><tr><td>CRTS</td><td>unfiltered</td><td>2005–2013</td><td>3 weeks</td><td>~40</td><td>33,000 deg<sup>2</sup></td><td>19–21</td></tr><tr><td>ASAS-SN</td><td>V</td><td>2012–2019</td><td>3–4 days</td><td>~200</td><td>all sky</td><td>17</td></tr><tr><td>WISE</td><td>W1 (W2)</td><td>2010–2019</td><td>6 months</td><td>12–15</td><td>all sky</td><td>17.6 (Vega)</td></tr></tbody></table>

Note. The cadence of WISE light curves is 6 months per visit. Optical data were obtained in annual “seasons.” For example, SDSS S82 data were obtained annually within a 2 to 3 month window, and the cadence effectively samples timescales from days to years. For a single season, the median SDSS cadence is ~5 days. Since most quasars in our sample are fainter than 19 mag in *r* band, we bin the CRTS and PTF data annually (inverse-variance-weighted mean) to increase the signal-to-noise ratio. The exact cadences in most of these optical surveys are more complicated than the quoted approximate values, given various observing constraints.

<table><caption>**Table 2**
Sample Statistics</caption><thead><tr><th>Cut</th><th>Number</th><th>Section</th></tr></thead><tbody><tr><td>All</td><td>9258</td><td>Section 2.1</td></tr><tr><td>Matched in WISE</td><td>8990</td><td>Section 2.2</td></tr><tr><td>Covered in DES</td><td>7384</td><td>Section 2.3</td></tr><tr><td>Lag quality cuts</td><td></td><td></td></tr><tr><td>S/N<sub>MIR</sub> > 4</td><td>1315</td><td>Section 3.1</td></tr><tr><td>autocut</td><td>1064</td><td>Section 3.5</td></tr><tr><td>finalcut</td><td>587</td><td>Section 3.5</td></tr></tbody></table>

the median CRTS magnitude to the contemporary calibrated PS1 magnitudes.

All optical data are then cross-calibrated to DES g band and converted to physical fluxes for our lag measurements. For surveys with multiband coverage, we also include *r*-band data (converted to g band) to increase the cadence; we only include PS1 and DES *i*-band data (converted to g band) when there are no g-band or *r*-band data within ±1 yr; no other bands were used given uncertainties in color transformations. We can safely ignore the small delays (≤a few days) across optical continuum bands owing to accretion disk RM for our quasars (e.g., Jiang et al. 2017; Mudd et al. 2018; Homayouni et al. 2019; Yu et al. 2020).

In this work we do not correct for host galaxy contamination in the photometry since the vast majority of our quasars are at $z > 0.5$ and are luminous enough to dominate the emission in both optical and MIR. Constant host stellar emission also does not affect lag measurements. The intrinsic MIR variability required for successful lag detection (see Section 3.5) is significantly higher than any systematic magnitude uncertainties due to host galaxy photometry.<fnref n="42" /> Even in the worst-case scenario where the epoch-by-epoch host galaxy photometry introduces significant, random variability in the light curves, it will only make the lag more difficult to detect rather than bias the lag in any particular direction. We also do not consider the

<table><caption>**Table 3**
Table Format for Compiled Light Curves</caption><thead><tr><th>Column</th><th>Units</th><th>Description</th></tr></thead><tbody><tr><td>DBID</td><td></td><td>Object ID of SDSS S82 quasars</td></tr><tr><td>MJD</td><td>days</td><td>Modified Julian Date</td></tr><tr><td>SURVEY</td><td></td><td>Name of the imaging survey</td></tr><tr><td>BAND</td><td></td><td>Photometric filter</td></tr><tr><td>MAG</td><td>mag</td><td>Optical magnitude in AB; WISE magnitude in Vega</td></tr><tr><td>MAG_ERR</td><td>mag</td><td>Uncertainty in magnitude</td></tr><tr><td>FLUX</td><td>mJy</td><td>Flux density</td></tr><tr><td>FLUX_ERR</td><td>mJy</td><td>Uncertainty in flux density</td></tr></tbody></table>

Note. This table compiles light curves for all 7348 S82 quasars with both DES and WISE coverage. Each row corresponds to a single epoch in a given survey and filter for a single object. All optical magnitudes and fluxes are converted to DES g band as described in Section 2.4.

(This table is available in its entirety in machine-readable form.)

contamination of broad emission line flux in the broad photometric bandpasses. The delayed response in broad-line flux to optical continuum variations occurs on much shorter light-crossing times than the torus lag, and therefore neglecting this complication will not affect the MIR lag measurements.

The final merged optical light curves and WISE W1 light curves for all 7384 quasars are available through the electronic version of this paper as a machine-readable table. The format of the table is described in Table 3. While these light curves have the longest duration compared to data sets used in earlier studies, their cadences are typically insufficient for BLR RM. Fortunately for dust RM, the much-extended torus entails much longer time delays and broader transfer functions, which makes MIR lag measurements possible even with the sparse sampling of the WISE light curves.

### 3. Lag Measurements

Robustly measuring the lag between two sets of light curves is a nontrivial task, which depends on the quality of the light curves (e.g., duration, cadence, signal-to-noise ratio), as well as the intrinsic variability of the light curves. The success of the lag measurement critically depends on the prominent features in the variable light curve; thus, only for quasars with significant variability during the monitoring period can we measure a reliable time lag. Furthermore, quasar variability is stochastic in nature, and given insufficient time baselines or

<footnote n="42">⁴² To investigate the stability of WISE epoch-by-epoch photometry for galaxies, we use ~15,000 $z < 0.3$ SDSS star-forming galaxies within the S82 region. The measurement-uncertainty-subtracted, median rms for these galaxies is 0.02 mag (0.0006 mag) in W1 (W2). This systematic uncertainty is far too small to have any impact on our lag measurements, where the quasar light largely dominates over host galaxy light in W1. A small fraction (0.8%) of these galaxies do show significant intrinsic variability, many of which are due to obscured type 2 AGNs or transient events (see a dramatic example in, e.g., Yang et al. 2019).</footnote>
cadences, analyses of the light curves can produce artificial lags ("aliases"), often reaching the edges of the light curves, where there is limited overlap in data points of the pair of light curves.

There are various techniques to measure the time lag between two light curves. The most commonly adopted method is the interpolated cross-correlation function (ICCF; Gaskell & Peterson 1987; Peterson et al. 1998; see Section 3.2). The major advantage of ICCF is that it is an empirical and fast method that is model independent, and the simple linear interpolation scheme across light-curve gaps can recover some variability information lost. There have been some recent modifications on the ICCF method (e.g., Grier et al. 2017, 2019; I-Hsiu et al. 2019) that incorporate a weighting scheme in ICCF to down-weight time delays with less overlapping data points in the light curves in the lag search. This modified method (WCCF; see Section 3.3) proves to be effective in eliminating aliases near the edge of the light curves.

A more recently developed technique, implemented in the public code JAVELIN (Zu et al. 2011), improves the interpolation scheme within light-curve gaps by assuming a damped random walk (DRW) model for stochastic variability of AGNs. The DRW model has proven to be a reasonably good prescription to describe the optical continuum variability of quasars (e.g., Kelly et al. 2009; Kozłowski et al. 2010; MacLeod et al. 2010; Zu et al. 2013) on timescales of interest to most RM studies of quasars, offering a superior interpolation scheme for the light curves to that of ICCF. JAVELIN also implements more statistically rigorous procedures to estimate the uncertainties in the interpolations and measurements of the lags with Bayesian inference. There are other public codes such as CREAM (Starkey et al. 2016) that model the driving light curve with alternative statistical methods and perform equally well as JAVELIN.

I-Hsiu et al. (2019) performed a detailed comparison of lag measurement methodologies for ICCF, JAVELIN, and the z-transformed discrete correlation function (zDCF; Alexander 2013), in particular for light curves with moderate to low qualities. They found that JAVELIN provides overall the best performance in terms of recovering the correct lag and estimating the lag uncertainties, while ICCF is less effective for low to moderate light-curve quality, with zDCF being the least effective method. When the light-curve quality is sufficiently high, all methods converge to consistent lag measurements.

Since the light-curve quality for our quasars is moderate to low, we will rely on the more robust lag measurement method JAVELIN to provide our fiducial measurements. However, we will also use the WCCF method to guide the JAVELIN measurements and to provide additional criteria in eliminating false positives. We note that currently there is no perfect method to measure lags for light curves of different qualities. Thus, we will also impose a set of cuts assisted by visual inspection to select a final "cleaned" lag sample (Section 3.5), as often done in recent studies (e.g., Grier et al. 2017, 2019; Lyu et al. 2019).

### 3.1. Intrinsic Variability Cut

To robustly detect the lag, the variability characteristics of the light curves are of critical importance. We define a variability metric to ensure that variability is well detected in the MIR light curve. The observed rms variability includes both intrinsic variability and photometric uncertainties. To estimate the intrinsic rms magnitude (or "excess rms") of a light curve

and the uncertainty of the intrinsic rms, we utilize a maximum likelihood estimator detailed in Shen et al. (2019b) (their Equations (5)–(9)).$^{43}$ The estimate of the intrinsic rms in the MIR light curve, $\sigma_{MIR}$, and its uncertainty, $\Delta\sigma_{MIR}$, are defined by Equation (8) of Shen et al. (2019b). We then define the signal-to-noise ratio of the estimated intrinsic rms as $S/N_{MIR} = \sigma_{MIR}/\Delta\sigma_{MIR}$. We require significant variability detection in the MIR light curve as $S/N_{MIR} > 4$. Among the 8990 quasars with WISE light curves, 1588 quasars have $S/N_{MIR} > 4$ in the W1 band. Among the 7582 quasars covered by DES, 7384 are matched in WISE, and 1315 of them have $S/N_{MIR} > 4$. Figure 2 shows the distribution of $S/N_{MIR}$ (left) and the intrinsic rms magnitude in the W1 band versus the mean W1 magnitude. Requiring significant variability in the MIR light curve implies that the optical light curve generally also has large variations. Many of these highly variable quasars fall into the category of “extreme variability quasars” (EVQs; Rumbaugh et al. 2018) that have more than 1 mag maximum variations in g band over multiyear timescales. These gradual, large-amplitude, multiyear variability features greatly facilitate the measurement of the MIR lag, as already demonstrated with earlier WISE light curves (e.g., Sheng et al. 2017).

### 3.2. Interpolated Cross-correlation Function

The cross-correlation function (CCF) is commonly used to measure the time delay between two light curves. For two sets of variables (e.g., light curves) X and Y, the Pearson correlation coefficient, r, can be calculated by

$$r = \frac{\sum_{i=1}^{N}(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_{i=1}^{N}(x_i - \bar{x})^2 \sum_{i=1}^{N}(y_i - \bar{y})^2}}, \quad (1)$$

where $\bar{x}$ and $\bar{y}$ are the mean of X and Y, respectively, $x_i$ and $y_i$ are the ith members, and N is the number of $(x_i, y_i)$ pairs. For two time series, CCF is the Pearson correlation coefficient as a function of the time displacement ($\tau$)$^{44}$ of one signal relative to the other. The time delay between the two time series can be determined by the time displacement with the maximum CCF. Furthermore, the maximum CCF coefficient, $r_{max}$, is a useful value to evaluate the time lag significance; a larger $r_{max}$ indicates that the two time series are better correlated (e.g., Grier et al. 2017). For uneven-sampled light curves in essentially all RM studies, ICCF is deployed to linearly interpolate the shifted light curves. We used the PyCCF code (Sun et al. 2018) to perform the ICCF calculations.

### 3.3. Weighted Cross-correlation Function

Due to the relative sparse time sampling and the limited baseline of the light curves, there are often multiple peaks in the ICCF. The peaks near the boundaries of the lag search limited by the duration of the light curve, and sometimes within large gaps in the light curves, are usually false positives owing to the small number of overlapping data points in the light-curve pair (but see rare counterexamples in, e.g., Shen et al. 2019a). There are formal estimations of the statistical

<footnote>$^{43}$ We correct a typo in their Equation (9): $\text{Var}[\mu] = \sigma_0^2/\Sigma g_i$.</footnote>

<footnote>$^{44}$ We adopt the convention that positive values of $\tau$ correspond to the MIR light curve lagging behind the optical light curve.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left: histogram of the S/N<sub>MIR</sub> for the intrinsic rms magnitude in W1 calculated for the light curve. Right: distribution in the W1 intrinsic rms vs. average W1 magnitude space for the parent sample (gray points) and the S/N<sub>MIR</sub> > 4 subsample (red points) to illustrate the amplitudes and uncertainties of the intrinsic rms variability measurements. The intrinsic rms magnitude and its uncertainties are described in Section 3.1.</figcaption></figure>

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

significance of the correlation that depend on both the correlation coefficient *r* and the number of overlapping data points (e.g., Bevington 1969; Shen et al. 2015). However, in practice this statistical significance is almost never used because of the noisy ICCF for light curves with low to moderate quality. To remedy for these aliasing peaks due to less overlapped data, Grier et al. (2017, 2019) developed a quantitative scheme taking into account the number of overlapping data points in the ICCF lag search. They used a weight function defined as $P(\tau) = [N(\tau)/N(0)]^2$, where $N(\tau)$ is the number of overlapping epochs as a function of time delay $\tau$. The exponent is somewhat arbitrary as long as it is positive to down-weight time delays with fewer overlapping data points.

I-Hsiu et al. (2019) further tested the efficiency and robustness of this weighting scheme with simulated light curves that mimic the low to moderate quality of data from recent RM survey programs, but with known input lags. They found that for survey-quality light curves, it is essentially necessary to impose this weighting scheme because the raw ICCF will produce an overwhelmingly large number of aliasing peaks in the less overlapped regime, as well as lags with large uncertainties. The weighted CCF (WCCF) is a much more robust method to recover the true lag.

To account for flux uncertainties in each epoch, we adopt a modified weight function considering both the number of overlapping epochs and the flux uncertainty of each epoch:

$$P(\tau) = \frac{\sum_{j=1}^{n} (1/\sigma_j)^2}{\sum_{i=1}^{N} (1/\sigma_i)^2}, \qquad (2)$$

where $\sigma$ is the single-epoch flux measurement error, *n* is the number of overlapping epochs at time delay $\tau$, and *N* is the maximum number of overlapping epochs. In the case of equal flux uncertainties, our definition of the weight function is equivalent to $N(\tau)/N(0)$. Therefore, the suppression of correlations at large absolute values of $\tau$ is not as aggressive as the weight function used in Grier et al. (2017).

In principle, the weight function should be defined for the optical light curve ($P_{OPT}$) and the MIR light curve ($P_{MIR}$) separately, since both light curves have different lengths, and the surveys contributing to the optical light curves have different sampling densities. A combined weight function $P_{OM} \equiv P_{OPT}P_{MIR}$ can then be defined to impose stronger suppression of correlations in the less overlapped regime. However, since the MIR light curve is less well sampled than

the optical light curve and dominates the cross-correlation signal, we only use $P_{MIR}$ to down-weight the original ICCF. In addition, the weights computed for the optical light curve ($P_{OPT}$) are complicated by the very different sampling densities and flux uncertainties in various surveys contributing to the light curve. Using $P_{OM} \equiv P_{OPT}P_{MIR}$ thus could drastically change the shape of the ICCF, which is more than necessary to mitigate lag aliases in the less overlapped regime. However, in defining the maximum search range in JAVELIN as discussed in Section 3.4, using the more stringent weights $P_{OM}$ is more efficient for the fit to converge.

Multiplying the ICCF by $P_{MIR}(\tau)$, we obtain the weighted CCF (WCCF) as

$$r_{\text{WCCF}}(\tau) = r_{\text{ICCF}}(\tau) \times P_{\text{MIR}}(\tau). \qquad (3)$$

The maximum of $r_{WCCF}(\tau)$, $r_{wmax}$, is a parameter similar to $r_{max}$ that can be used to evaluate the time lag significance that incorporates information from the weight function.

Figure 3 shows an example of ICCF and WCCF of a quasar in our S82 sample, as well as the weight functions $P_{MIR}$ and $P_{OPT}$. The WCCF efficiently eliminated the fake peaks near the boundaries of the lag search window and enhanced the primary CCF peak. Given the benefits of the WCCF, we will use the WCCF and the weight functions defined here to guide our JAVELIN measurements and to refine our criteria of significant lag detection in Section 3.5. But we note that our fiducial lags are computed using the more robust JAVELIN method, and the WCCF is used only for cross-checks and to prevent JAVELIN lag searches near the boundaries of the light curves.

### 3.4. Lag Measurements with JAVELIN

In JAVELIN, the MIR light curve is a scaled and smoothed version of the driving optical/UV light curve owing to the extended structure of the dust torus. The responding torus IR light curve is the convolution of the optical continuum light curve, $f_{cont}(t)$, with a transfer function, $\Psi(\tau)$, determined by the geometry of the dust torus. At rest frame ~2 µm, the emission is dominated by the torus in most of our quasars. But in extreme cases, the accretion disk continuum may contribute a significant amount of IR emission, following the theoretical prediction of the $F_{\nu} \propto \nu^{1/3}$ law (Kishimoto et al. 2008), which is delayed relative to the optical light curve on much shorter timescales (the light-crossing time of the accretion disk) than the torus lag.
<figure><img src="image_7.png" /><figcaption>Figure 3. Example quasar in our S82 sample, J0023+0035, at z = 0.4219. The top two panels show the light curves in (a) optical and (b) MIR W1 band. We use combined optical data from multiple ground-based imaging surveys, including SDSS, CRTS, PS1, PTF, DES, and ZTF. In panel (c), the gray dashed line and the cyan dashed-dotted line are the normalized weight functions from overlapped MIR and optical data ($P_{MIR}(\tau)$ and $P_{OPT}(\tau)$), respectively. The purple dotted line denotes $P_{OM}$. The red solid line is the original ICCF, and the blue solid line is the weighted ICCF (WCCF = ICCF $\times P_{MIR}$). The WCCF efficiently eliminates aliasing peaks near the edges of the lag search range owing to limited overlap in the light curves. The purple histogram (with black solid outline) shows the JAVELIN posterior lag distribution in observed frame. The vertical green dashed line shows the peak of the JAVELIN posterior distribution. We measure an MIR (W1) lag as $306.5_{-25.8}^{+48.8}$ days using JAVELIN. Panel (d) shows the JAVELIN distribution of logarithmic probability vs. lag. The red dashed line marks the lag with the maximum probability. The two orange vertical lines in panel (d) show the lag search window in JAVELIN defined by $P_{OM} > 0.1$. The consistency between the maximum probability lag (panel (d)) and the peak of the lag posterior (panel (c)) is one of the criteria for our visual rejection of less secure lags (see Section 3.5).</figcaption></figure>

Neglecting the IR emission from the accretion disk, the responding IR light curve, $f_{IR}(t)$, can be written as

$$f_{IR}(t) = \alpha \int d\tau \Psi(\tau) f_{cont}(t - \tau), \quad (4)$$

where $\tau$ is the average light-travel time from the accretion disk to the dust torus and $\alpha$ is the ratio between the responding IR and optical variability amplitudes.

For our purposes, we use the two-band photometric RM model in JAVELIN (Zu et al. 2011), which models the quasar continuum variability with the DRW model and measures the lag between two photometric bands as described in Equation (4). For simplicity, we consider a top-hat transfer function for the dust torus:

$$\Psi(t) = \begin{cases} 1, & t_1 \le t \le t_2 \\ 0, & \text{otherwise} \end{cases}, \quad (5)$$

which has a mean lag of $\bar{\tau} = (t_1 + t_2)/2$ and a width of $\Delta\tau = (t_2 - t_1)$. We allow JAVELIN to fit $\Delta\tau$ as a free

<figure><img src="image_8.png" /><figcaption>Figure 4. Example of lag measurements using optical and WISE light curve for the nearby AGN Mrk 110. Notations are the same as in Figure 3. The optical fluxes are converted to ASAS-SN V band. Koshida et al. (2014) measured the K-band lag of Mrk 110 as $113.1_{-8.6}^{+8.8}$ and $124.1_{-7.1}^{+7.1}$ days in two campaigns in the observed frame. Using optical and WISE light curves, we successfully measure an MIR (W1) lag of Mrk 110 as $158.8_{-6.4}^{+8.0}$ days, where the longer lag is likely due to the longer-wavelength IR band. The WISE light curve is apparently smoother than the optical light curve because of the averaging effect of responses from different parts of the extended torus.</figcaption></figure>

parameter, and we find that there is a broad trend that the average $\Delta\tau$ increases with the best-fit lag.

For each of the quasars that pass the $S/N_{MIR} > 4$ cut, we run JAVELIN using 50,000 Markov Chain Monte Carlo (MCMC) chains, sufficient for convergence. To reduce the most time-consuming chains near the boundaries of the time delays that can be reasonably ruled out by WCCF, we allow JAVELIN to explore a range of lags where $P_{OM} > 0.1$ (typically $t_{min} \sim -2000$ and $t_{max} \sim 5500$ days). This maximum lag search range corresponds to ~70% of the maximum baseline defined by the optical or MIR light curve (or ~30% overlap in the optical+IR light curves), which is typically required as the bare minimum for robustly measuring a lag (e.g., Grier et al. 2017, 2019). In other words, the use of the WCCF weights does not impose stringent limits in the JAVELIN lag search other than eliminating lags (mostly false positives) near the boundaries of the light curves.

In the example shown in Figure 3, the black solid histogram in panel (c) indicates the JAVELIN posterior distribution of lags. We obtain the median lag from the posterior as $\tau_{JAVELIN}$ and the $1\sigma$ lag uncertainty as the semiamplitude of the range enclosing the 16th and 84th percentiles of the JAVELIN posterior distribution of the mean lag.

The moderate to low quality of our light curves generally does not allow us to explore more sophisticated torus transfer function forms, except for the best individual cases. However, we have tested varying the width of the top-hat transfer function and confirm that this somewhat arbitrary choice of the transfer function form in JAVELIN does not affect much the
measured mean lag and its uncertainties. This measured "average" lag simply reflects the inner boundary of the torus, as rest-frame K band traces the hottest dust near the sublimation temperature (e.g., Hönig & Kishimoto 2011).

To demonstrate the feasibility of our lag measurement approach, we show in Figure 4 a local AGN (Mrk 110) with well-measured K-band lags using earlier ground-based IR light curves. We successfully measured the lag with optical and WISE light curves, albeit with a slightly larger value owing to the different (3.5 µm) IR band used. We have also applied our methodology to the entire PG quasar sample studied in Lyu et al. (2019) and found consistent results in general (see discussion in Section 3.6), although the latter used a different method (a χ² minimization scheme) to measure the lag.

### 3.5. Criteria of Significant Lag Detection

To define a final “cleaned” sample of lags, we perform the following automatic cuts and manual rejections with visual inspection. Most importantly, these rejection criteria do not use any prior information of an anticipated R–L relation for the lags; therefore, they do not introduce any bias to the observed R–L relation.

As described in Section 3.1, we first cut the sample (matched with WISE and DES) by S/N<sub>MIR</sub> > 4, resulting in 1315 quasars. We then require r<sub>wmax</sub> > 0.5 in the WCCF (Section 3.3) to remove objects where the optical and MIR light curves are not well correlated at any time lag within the entire lag search window. This step results in 1283 quasars. To identify the primary peak of the lag distribution in the JAVELIN analysis, we bin the JAVELIN MCMC lag posterior with a bin size of 50 days. We identify the tallest peak of the binned distribution as the primary lag peak. From the primary peak, we search on each side to define the peak boundary where the binned number of lags falls below 5% of the number of lags in the tallest bin. We require the fraction of lags within the primary peak to all lags f<sub>peak</sub> > 0.5. This step using JAVELIN lag posterior is to ensure that there is clear evidence for a primary peak in the lag distribution. There are 1064 quasars remaining after these automatic cuts (autocut).

We then visually inspect the light curves, WCCF, and the JAVELIN lag posterior distribution in the remaining objects to remove unreliable lags (examples of each case are shown in the Appendix):

1. We reject objects with large gaps (≥2 yr) in the optical light curve, especially when the data do not overlap sufficiently to meaningfully measure the lag (see an example in Figure A1)—193 objects were removed.

2. We further require that the maximum probability lag from the MCMC chains in JAVELIN is consistent within the 1σ range of the primary peak of the lag posterior to ensure that the MCMC chains are well converged. Otherwise, the MIR light curve might align with the optical light curve at more than one potential lag (see Figure A1 for an example). In such cases, the MCMC chains are not well converged and the results are unreliable—171 objects with inconsistent primary peak and maximum probability lag were removed.

3. Very rarely the MIR variability is dominated by one season (e.g., a flare or a systematic outlier). Figure A1 shows an example of these objects that have extreme

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Bolometric luminosity vs. redshift for our parent quasar sample (blue points) and the final lag sample (red points). The four dashed lines correspond to observed-frame lags of 1, 5, 10, and 15 yr, respectively, using the best-fit R–L relation in Kishimoto et al. (2007). The blue circles are color-coded by the W1 intrinsic variability S/N<sub>MIR</sub>. Our final lag sample only includes those with S/N<sub>MIR</sub> > 4 and have lags that are measurable with the baselines of our light curves. The few objects scattered beyond the τ<sub>obs</sub> = 15 yr line have shorter actual lags (<15 yr) than predicted from the canonical R–L relation in Kishimoto et al. (2007).</figcaption></figure>

MIR variability in a single epoch and very noisy ICCF. We excluded six such objects.

4. Finally, we remove objects with multiple peaks in the JAVELIN lag posterior, or one extreme broad peak spanning over thousands of days, or with an obvious secondary peak (see an example in Figure A1). This is to ensure that we only select objects with one well-defined primary lag peak, in assist to the autocut criterion of f<sub>peak</sub> > 0.5 mentioned above. We removed 330 objects with multiple peaks.

Our finalcut sample includes 587 quasars. Although we do not require the lag to be positive, there are no negative lags remaining after our selection process.

Figure 5 displays the distribution of the S82 parent sample and the finalcut 587 quasars with high-fidelity MIR lag measurements in the quasar luminosity–redshift plane. Objects are removed from the parent sample mainly as a result of the S/N<sub>MIR</sub> > 4 cut, as well as the difficulty in measuring a long lag given the maximum baseline of the light curves. Nevertheless, the finalcut sample includes quasars over a broad range of redshift and luminosity that far extends the regime probed by previous samples with IR lag measurements at z < 0.3 with almost an order-of-magnitude increase in statistics.

Figure 6 displays the distribution of the parent quasar sample and the finalcut sample in the S/N<sub>MIR</sub> and WCCF peak space. Such a plot is often used in recent large-scale broad-line RM programs to demonstrate the statistical detection of lags using light curves of low to moderate quality (e.g., Shen et al. 2016; Grier et al. 2017, 2019). The lengths and overlap of our optical and MIR light curves allow negative (i.e., MIR leading optical) and positive lags in the approximate range of [-3000, 7000] days. If there is statistical evidence of the MIR light curve lagging behind the optical light curve, we would observe an asymmetry in the lag distribution toward more positive lags. Examining the range of ±3000 days in Figure 6, there is indeed a preference of positive lags, indicating that statistically the MIR light curve lags the optical, as expected from torus
<figure><img src="image_10.png" /><figcaption>**Figure 6.** S/N<sub>MIR</sub> of W1 intrinsic rms vs. the peak of WCCF for the parent sample (gray points) and the final lag sample (red points). The top panel shows the histograms of both samples. The asymmetry below the S/N = 4 cut suggests that most of these lags are genuine lags in a statistical sense, albeit with larger uncertainties compared with the high-fidelity lags shown in red.</figcaption></figure>

<table><caption>**Table 4**
Table Format for the Finalcut Lag Sample</caption><thead><tr><th>Column</th><th>Units</th><th>Description</th></tr></thead><tbody><tr><td>DBID</td><td></td><td>Object ID of SDSS S82 quasars</td></tr><tr><td>RA</td><td>deg</td><td>J2000 R.A.</td></tr><tr><td>DEC</td><td>deg</td><td>J2000 Decl.</td></tr><tr><td>REDSHIFT</td><td></td><td>Spectroscopic redshift</td></tr><tr><td>LOGLBOL</td><td>erg s<sup>-1</sup></td><td>Bolometric luminosity from Shen et al. (2011)</td></tr><tr><td>LOGLBOL_ERR</td><td>erg s<sup>-1</sup></td><td>Uncertainty in LOGLBOL</td></tr><tr><td>LOGBH</td><td>M<sub>☉</sub></td><td>Fiducial single-epoch BH mass from Shen et al. (2011)</td></tr><tr><td>LOGBH_ERR</td><td>M<sub>☉</sub></td><td>Uncertainty in LOGBH</td></tr><tr><td>LOGEDD Ratio</td><td></td><td>Eddington ratio based on fiducial single-epoch BH mass</td></tr><tr><td>N Wise</td><td></td><td>Number of WISE epochs</td></tr><tr><td>W1_AVG</td><td></td><td>W1-band weighted average magnitude</td></tr><tr><td>W1_RMS</td><td></td><td>W1-band intrinsic rms</td></tr><tr><td>W1_SIGRMS</td><td></td><td>Uncertainty in W1_rms</td></tr><tr><td>W1_S/N</td><td></td><td>S/N of W1-band intrinsic rms</td></tr><tr><td>RMAX</td><td></td><td>Peak correlation r<sub>max</sub> from ICCF</td></tr><tr><td>RWMAX</td><td></td><td>Peak correlation r<sub>wmax</sub> from WCCF</td></tr><tr><td>PEAK_WCCF</td><td>days</td><td>Peak location of WCCF</td></tr><tr><td>F PEAK</td><td></td><td>f<sub>peak</sub> (see Section 3.5)</td></tr><tr><td>TAU_JAVELIN</td><td>days</td><td>JAVELIN time delay in observed frame</td></tr><tr><td>TAU_JAVELIN_LOW</td><td>days</td><td>1σ lower limit of TAU_JAVELIN</td></tr><tr><td>TAU_JAVELIN Upper</td><td>days</td><td>1σ upper limit of TAU_JAVELIN</td></tr></tbody></table>

(This table is available in its entirety in machine-readable form.)

reprocessing. Figure 6 suggests that there are potentially thousands of MIR lags that are real below the S/N<sub>MIR</sub> = 4 cut. However, only the red points above the S/N<sub>MIR</sub> = 4 cut

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Comparison of IR lag measurements for the PG quasar sample studied in Lyu et al. (2019) using optical and WISE W1 light curves. Despite very different lag measurement methodologies, there is good agreement between the lags measured in both studies. The few outliers have low-quality WCCF (r<sub>max</sub> < 0.5) and so would not pass our selection criteria as reliable lags (Section 3.5).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Redshift distributions of different AGN/quasar samples with IR lag measurements.</figcaption></figure>

provide the most secure individual lag measurements to study the relation between dust lags and quasar luminosity.

We summarize our lag measurements for the 587 quasars in a machine-readable table, and we describe the columns in Table 4. In the Appendix (Figure A2) we show several additional examples of lag measurements in our finalcut sample.

### 3.6. Comparison with Earlier Work

Using our lag measurement methodology, we reanalyze the PG quasar sample studied in Lyu et al. (2019) using our compiled optical and WISE light curves. Figure 7 compares the lag measurements in both studies. We find good agreement between our measurements and those in Lyu et al. (2019). This result confirms that our lag measurement methodology is robust, even though it differs significantly from the χ<sup>2</sup>-fit method adopted by Lyu et al. (2019).

Figure 8 displays the redshift distribution for various samples with IR lag measurements. Our sample covers a much
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Correlation between the torus size (inferred from the rest-frame lag $\tau_{IR}$) and the rest-frame optical quasar luminosity. Our finalcut lag sample is shown with red filled circles with gray error bars, and the PG quasars from Lyu et al. (2019) are shown with open circles. The black diamonds and gray crosses are earlier lag measurements between K band and optical band for a small sample of local AGNs (Suganuma et al. 2006; Koshida et al. 2014). The magenta line is the $R-L$ relation in Kishimoto et al. (2007). The two blue dashed-dotted lines are the rest-frame timescales corresponding to the marked observed-frame time duration (7 and 14 yr) calculated at the 90th percentile redshift of our quasars in each luminosity bin. Regions above the $t_{obs} = 7$ yr line will suffer significant incompleteness in lag detection given the lengths of our optical and MIR light curves.</figcaption></figure>

broader redshift range and is the first statistical sample of IR lags at $z \gtrsim 0.3$.

## 4. Results

We show our main results using the finalcut sample of 587 quasars with high-fidelity MIR lags in Figure 9. We have converted all luminosities to V-band luminosity (assuming a bolometric correction of 10 in V band) to directly compare with earlier studies on low-redshift AGNs and PG quasars. Our sample well samples the high-luminosity end of the distribution and tightly follows the best-fit relation in Koshida et al. (2014) based on a local AGN sample and K-band lags. Our MIR lags also agree with the lags measured for tens of PG quasars in Lyu et al. (2019) over the same luminosity range, indicating negligible evolution from $z < 0.5$ to $z \sim 1$.

Because our sample dominates in number in the $R-L$ plot, a joint fit of all lags across the entire luminosity range will be heavily weighted by our sample. Furthermore, several selection effects will tend to bias the measurement of the $R-L$ relation, in particular the slope (see discussions in Section 5). For these reasons, we do not perform a formal fit to the IR $R-L$ relation in this work and defer a proper analysis of the $R-L$ relation to future work. Nevertheless, we can still estimate the intrinsic scatter of this relation using our large sample, which is much less biased than the slope owing to selection effects. Using the Bayesian regression algorithm developed by Kelly (2007), we estimate an intrinsic scatter of only $\sim 0.17$ dex for our IR lag sample.

Figure 9 demonstrates the power of combining decade-long optical and IR surveys in dust RM in distant quasars. The inferred $R-L$ relation over more than four orders of magnitude in AGN luminosity has profound implications for

torus structure and physics, as well as the utility of using this relation as a luminosity indicator for cosmology. We discuss potential selection biases in our sample and origins of the scatter in the relation in Section 5.

## 5. Discussion

Figure 9 reveals a slight trend of deviation from the local $R-L$ relation toward the high-luminosity end. This deviation is most likely caused by selection effects. First, higher-luminosity quasars on average are at higher redshifts (Figure 5), where the longer observed-frame lags due to higher luminosity and the $1+z$ cosmic time dilation are more difficult to measure given the fixed light-curve baseline. Modeling the detailed selection function given the light-curve duration, cadence, and variability S/N and for an underlying sample of quasars over broad redshift and luminosity ranges can be achieved with simulations (e.g., Shen et al. 2015; I-Hsiu et al. 2019) but is beyond the scope of this paper. Here we provide a qualitative understanding of the loss of long lags due to the baselines of our optical+MIR light curves.

In a given luminosity bin for our quasar sample in Figure 9, we determine the 90th percentile redshift, $z_{90}$, in that bin and calculate the rest-frame time at $z_{90}$ corresponding to a given observed-frame time duration. We use this rest-frame time as a rough estimate of the upper limit of measurable lags at that luminosity, given the survey length. Since higher-luminosity quasars on average have higher redshifts in our sample, this rest-frame timescale decreases with luminosity. Given the actual lengths of our optical and MIR light curves, we can reasonably assume that all lags shorter than 7 yr in the observed frame will be detected since the shifted MIR light curve is in complete overlap with the earlier optical light curve. However, the probability of detection will decrease toward longer lags in the observed frame. For an observed-frame lag of 14 yr, we only have $\sim 60\%$ of the full MIR light curve overlapping with the earlier optical light curve, and therefore the detection probability will be significantly reduced compared to the case of shorter lags. Therefore, we expect that our survey will start significantly losing lags beyond the line corresponding to 7 yr, shown in Figure 9 as the dashed-triple-dotted line. As a consequence, this selection effect will bias the average lags low toward the high-luminosity end.

A second selection effect is that as objects move toward the high-luminosity end (on average high-redshift end), the W1 3.5 µm band samples shorter rest-frame IR wavelength. We have found that the W1-band lags are systematically shorter than the W2-band lags in the same quasars, suggesting wavelength-dependent IR lags (Q. Yang et al. 2020, in preparation). This selection effect would bias the average lags lower toward higher redshift and therefore higher average luminosity in Figure 9, leading to slight flattening of the slope. Nonnegligible contamination from the accretion disk in the W1 band in extreme cases may also contribute to this bias.

The combination of these two selection effects can qualitatively explain the slight deviation of our sample from the local $R-L$ relation at the high-luminosity end. In future work, we will carefully model these selection effects with detailed simulations and to derive unbiased constraints on the IR $R-L$ relation.

To examine the scatter around the average $R-L$ relation, we show our finalcut lag sample in Figure 10, where we color-code the objects by various properties. We found that the scatter in
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Distributions of our final cut lag sample in the *R*−*L* plane, color-coded by different properties (top left: Eddington ratio; top right: rest-frame wavelength sampled by *W*1; bottom: intrinsic *W*1 variability in magnitude). At fixed luminosity, the scatter in lag does not depend on any of these additional parameters. The apparent trend of rest-frame wavelength with luminosity is driven by the apparent redshift−luminosity relation. The apparent trend of reducing average variability amplitude with luminosity is consistent with the variability−luminosity relation observed in the optical.</figcaption></figure>

MIR lag at fixed luminosity does not depend on the Eddington ratio or the variability amplitude of the quasar. This greatly simplifies the interpretation of the MIR *R*−*L* relation as due to a primary luminosity-driven effect. The apparent trend of rest-frame IR wavelength sampled by *W*1 as a function of luminosity (top right panel of Figure 10) is caused by the apparent luminosity−redshift relation. In addition, higher-luminosity quasars tend to have on average lower variability amplitude (bottom panel), consistent with previous studies of quasar variability in MIR (e.g., Kozłowski et al. 2016) and optical (e.g., MacLeod et al. 2012).

**6. Summary**

We have presented first results from our dust RM study of distant ($z \gtrsim 0.3$) quasars in the SDSS Stripe 82 region, using ground-based optical imaging and the MIR light curves from the WISE satellite. Our optical light curves span a baseline of ~20 yr, and the MIR light curves span almost a decade. We presented high-fidelity optical-WISE *W*1 (3.5 µm) lag measurements for 587 quasars at $0.3 \lesssim z \lesssim 2$, with a median redshift of 0.8. Our statistical analysis suggests that there are thousands more lags that are measurable (Figure 6), but their uncertainties are much larger than our high-fidelity sample, and therefore these less reliable lags are excluded from our analysis in this work.

The 587 quasars with high-fidelity MIR lags span more than two orders of magnitude in quasar luminosity and tightly follow the torus *R*−*L* relation observed for nearby AGNs observed in the optical and *K* band. The intrinsic scatter of the *R*−*L* relation defined by our sample alone is only ~0.17 dex. Furthermore, there is no apparent dependence of the scatter in IR lag at fixed luminosity on additional quasar properties such as the Eddington ratio and variability amplitude of the quasar, suggesting that luminosity is the primary driver to determine the IR lag. However, to robustly measure the *R*−*L* relation with the WISE sample, we must account for selection biases at the high-luminosity end owing to the duration of the light curves (Section 5). In addition, the WISE *W*1 band probes shorter rest-frame wavelengths toward higher redshift, and the wavelength dependence of torus lags must also be taken into consideration.

The observed global torus *R*−*L* relation over more than four orders of magnitude in AGN luminosity is remarkable, considering the small intrinsic scatter and the lack of dependence on additional quasar properties. The physics of setting the inner dust torus radius (e.g., dust sublimation) is also much simpler than photoionization in broad-line RM (e.g., Czerny & Hryniewicz 2011), making this dust *R*−*L* relation an attractive luminosity indicator to probe cosmology at high redshift (e.g., Yoshii et al. 2014; Koshida et al. 2017). However, the systematics in the dust *R*−*L* relation and selection biases affecting the measurement of this relation
must be thoroughly investigated before this technique can be applied to cosmology.

On the other hand, measuring reliable dust torus lags in AGNs with different physical properties across a broad range of redshift also has tremendous value to understanding the physics of AGNs. Our results based on current optical and IR imaging survey data provide strong endorsement for the joint analysis of the 10 yr optical light curves from the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST) and the planned 5 yr IR survey with the Nancy Grace Roman Space Telescope (NGRST, previously known as WFIRST). The LSST+NGRST data set resembles the light curves studied here with a similar shift in the onset of the optical and IR temporal coverage. However, LSST+NGRST will be able to measure light curves for AGNs that are up to 5 mag fainter than our S82 quasar sample. The expected dust lags from these lower-luminosity AGNs will be a factor of ~10 shorter and therefore easily measurable with the maximum ~10 yr LSST +NGRST baselines. Since NGRST only goes to ~2 µm, the dust torus lag measurements are necessarily limited to the $z \lesssim 0.5$ regime, where the variability at 2 µm (observed frame) from the torus emission can still be reasonably well measured on top of the constant host galaxy light. Importantly, by optimizing the cadence and overlap between NGRST and LSST, it is possible to compile a large sample of low-redshift and low-luminosity AGNs with well-measured torus lags to densely populate the low-luminosity regime in the torus $R-L$ relation. For reference, the cumulative sky density of $i < 24$ quasars at $z < 0.5$ is ~30 deg$^{-2}$ estimated using the extrapolation of the Hopkins et al. (2007) quasar luminosity function. The V-band luminosity for $z = 0.5$ AGNs at this flux limit is about $\log(L_v/\text{erg s}^{-1}) \sim 42.4$ using magnitude and luminosity conversions in Richards et al. (2006) and Shen et al. (2009), corresponding to an observed-frame torus lag of ~1 month. Thus, even the wide survey cadence of LSST (e.g., 15 days in $r$ and 3 days in merged multiband light curves; LSST Science Collaborations et al. 2017) would be sufficient for lag measurements. The LSST Deep Drill Fields will have much higher cadences over tens of square degrees, which could be the high-priority fields for NGRST to (partially) overlap in a medium-deep survey and to provide dense IR light curves.

We outline some ongoing work to further explore the utility of current optical+WISE light curves on AGN dust RM:

1. We are developing a more robust forward modeling scheme to measure the intrinsic torus $R-L$ relation from surveys that takes into account selection effects introduced by the duration (and, to a lesser extent, cadence) of light curves, using simulated light curves following the approaches described in, e.g., Shen et al. (2015) and I-Hsiu et al. (2019).

2. We are studying individual cases with sufficient light-curve quality to constrain the transfer function of dust RM and developing physical models for the structure of AGN dust torus that can be constrained with dust RM (e.g., Nenkova et al. 2008; Shen 2012; Almeyda et al. 2017; Yang et al. 2019).

3. We are investigating multiband IR lags ($W1$ and $W2$) for the S82 sample, and we will use these multiband lags to further constrain the temperature profile of the torus (e.g., Hönig & Kishimoto 2011) and the induced scatter in the $R-L$ relation. In the meantime, we are compiling NIR

light curves for a subset of our quasars to measure NIR lags to extend the wavelength coverage.

4. We are expanding our dust RM analysis to quasars outside the Stripe 82 region. The longest optical baseline then comes from the shallower CRTS survey, which necessarily limits our study to the rarer and brighter quasars. However, the much larger sky coverage compensates for the lower sky density, and therefore we expect a significant increase in the number of the brightest quasars with measurable lags from optical and WISE light curves, most of which will be at lower redshifts and lower luminosities than the S82 sample studied here.

We thank the anonymous referee for comments that improved the manuscript. Q.Y. and Y.S. acknowledge support from NSF grant AST-1715579 (Q.Y., Y.S.) and an Alfred P. Sloan Research Fellowship (Y.S.).

This publication makes use of data products from WISE, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This publication also makes use of data products from NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the Planetary Science Division of the National Aeronautics and Space Administration.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at The Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.
Based in part on observations at the Cerro Tololo Inter-
American Observatory, National Optical Astronomy Observa-
tory, which is operated by the Association of Universities for
Research in Astronomy (AURA) under a cooperative agree-
ment with the National Science Foundation.

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS website is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

The PS1 and the PS1 public science archive have been made
possible through contributions by the Institute for Astronomy,
the University of Hawaii, the Pan-STARRS Project Office, the
Max-Planck Society and its participating institutes, the Max
Planck Institute for Astronomy, Heidelberg and the Max
Planck Institute for Extraterrestrial Physics, Garching, Johns
Hopkins University, Durham University, the University of
Edinburgh, the Queen's University Belfast, the Harvard-
Smithsonian Center for Astrophysics, the Las Cumbres
Observatory Global Telescope Network Incorporated, the
National Central University of Taiwan, the Space Telescope
Science Institute, the National Aeronautics and Space Admin-
istration under grant No. NNX08AR22G issued through the
Planetary Science Division of the NASA Science Mission

Directorate, the National Science Foundation grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

The CSS survey is funded by the National Aeronautics and
Space Administration under grant No. NNG05GF22G issued
through the Science Mission Directorate Near-Earth Objects
Observations Program. The CRTS survey is supported by the
U.S. National Science Foundation under grants AST-0909182.

ZTF: Based on observations obtained with the Samuel Oschin 48-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under grant No. AST-1440341 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron and Humboldt University, Los Alamos National Laboratories, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW.

ASAS-SN is supported by the Gordon and Betty Moore Foundation through grant GBMF5490 to the Ohio State University and NSF grant AST-1515927. Development of ASAS-SN has been supported by NSF grant AST-0908816, the Mt. Cuba Astronomical Foundation, the Center for Cosmology and Astro-Particle Physics at The Ohio State University, the Chinese Academy of Sciences South America Center for Astronomy (CASSACA), the Villum Foundation, and George Skestos.

Appendix

In Figure A1 we show one example for each category of
objects excluded by visual inspection in our finalcut
(Section 3.5).

In Figure A2 we show four more examples of lag detections
in our finalcut sample that cover a range of observed-frame
lags. The full figure set for all 587 objects in the finalcut sample
is available online.
<figure><img src="image_15.png" /></figure>

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

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

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

**Figure A1.** Top left: this object was rejected by visual inspection owing to large gaps (>4 yr) in the optical light curve. Top right: this object was rejected by visual inspection because the JAVELIN posterior MCMC lag with the maximum probability (red dashed line in panel (d)) is inconsistent with the primary peak in the lag posterior (green dashed line in panel (c)). Bottom left: this object was rejected by visual inspection because its MIR light curve shows extreme variability in a single epoch and the ICCF is very noisy. Bottom right: this object was rejected by visual inspection because there is an obvious secondary peak. This is to ensure that we only select objects with one significant primary peak.
<figure><img src="image_19.png" /></figure>

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

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

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

**Figure A2.** The upper two panels show the light curves in optical (a) and MIR W1 band (b), respectively. In panel (c), the gray dashed line and the cyan dashed-dotted line are the normalized weight functions from overlapped MIR and optical data ($P_{\text{MIR}}(\tau)$ and $P_{\text{OPT}}(\tau)$), respectively. The purple dotted line denotes $P_{\text{OM}}$. The red solid line is the original ICCF, and the blue solid line is the weighted ICCF (WCCF = $\text{ICCF} \times P_{\text{MIR}}$). The histogram shows the JAVELIN posterior lag distribution in observed frame. The vertical green dashed line shows the peak of the JAVELIN posterior distribution. Panel (d) shows the JAVELIN distribution of logarithmic probability vs. lag. The red dashed line marks the lag with the maximum probability. The two orange vertical lines in panel (d) show the lag search window in JAVELIN defined by $P_{\text{OM}} > 0.1$. The upper two panels show the light curves in optical (a) and MIR W1 band (b), respectively. In panel (c), the gray dashed line and the cyan dashed-dotted line are the normalized weight functions from overlapped MIR and optical data ($P_{\text{MIR}}(\tau)$ and $P_{\text{OPT}}(\tau)$), respectively. The purple dotted line denotes $P_{\text{OM}}$. The red solid line is the original ICCF, and the blue solid line is the weighted ICCF (WCCF = $\text{ICCF} \times P_{\text{MIR}}$). The WCCF efficiently eliminates aliasing peaks near the edges of the lag search range owing to limited overlap in the light curves. The purple histogram (with black solid outline) shows the JAVELIN posterior lag distribution in observed frame. The vertical green dashed line shows the peak of the JAVELIN posterior distribution. Panel (d) shows the JAVELIN distribution of logarithmic probability vs. lag. The red dashed line marks the lag with the maximum probability. The two orange vertical lines in panel (d) show the lag search window in JAVELIN defined by $P_{\text{OM}} > 0.1$. Examples of lag detections in our finalcut sample that cover a range of observed-frame lags. Notations are the same as Figure 3. The complete figure set (587 images) is available in the online journal.

(The complete figure set (587 images) is available.)
ORCID iDs

Qian Yang https://orcid.org/0000-0002-6893-3742

Yue Shen https://orcid.org/0000-0003-1659-7035

Xin Liu https://orcid.org/0000-0003-0049-5210

James Annis https://orcid.org/0000-0002-0609-3987

Santiago Avila https://orcid.org/0000-0001-5047-3662

David Brooks https://orcid.org/0000-0002-8458-5047

David Burke https://orcid.org/0000-0003-1866-1950

Aurelio Carnero Rosell https://orcid.org/0000-0003-3044-5150

Matias Carrasco Kind https://orcid.org/0000-0002-4802-3194

Juan De Vicente https://orcid.org/0000-0001-8318-6813

Shantanu Desai https://orcid.org/0000-0002-0466-3288

H. Thomas Diehl https://orcid.org/0000-0002-8357-7467

Brenna Flaugher https://orcid.org/0000-0002-2367-5049

Josh Friedman https://orcid.org/0000-0003-4079-3263

Juan Garcia-Bellido https://orcid.org/0000-0002-9370-8360

David Gerdes https://orcid.org/0000-0001-6942-2736

Daniel Gruen https://orcid.org/0000-0003-3270-7644

Robert Gruendl https://orcid.org/0000-0002-4588-6517

Julia Gschwend https://orcid.org/0000-0003-3023-8362

Gaston Gutierrez https://orcid.org/0000-0003-0825-0517

Samuel Hinton https://orcid.org/0000-0003-2071-9349

Devon L. Hollowood https://orcid.org/0000-0002-9369-4157

Klaus Honscheid https://orcid.org/0000-0002-6550-2023

Nikolay Kuropatkin https://orcid.org/0000-0003-2511-0946

Marcio Maia https://orcid.org/0000-0001-9856-9307

Jennifer Marshall https://orcid.org/0000-0003-0710-9474

Paul Martini https://orcid.org/0000-0002-4279-4182

Peter Melchior https://orcid.org/0000-0002-8873-5065

Felipe Menanteau https://orcid.org/0000-0002-1372-2534

Ramon Miquel https://orcid.org/0000-0002-6610-4836

Andrés Plazas Malagón https://orcid.org/0000-0002-2598-0514

Kathy Romer https://orcid.org/0000-0002-9328-879X

Eusebio Sanchez https://orcid.org/0000-0002-9646-8198

Michael Schubnell https://orcid.org/0000-0001-9504-2059

Ignacio Sevilla https://orcid.org/0000-0002-1831-1953

Mathew Smith https://orcid.org/0000-0002-3321-1432

Eric Suchyta https://orcid.org/0000-0002-7047-9358

Gregory Tarle https://orcid.org/0000-0003-1704-0781

References

Abbott, T. M. C., Abdalla, F. B., Allam, S., et al. 2018, ApJS, 239, 18

Alexander, T. 2013, arXiv:1302.1508

Almeyda, T., Robinson, A., Richmond, M., Vazquez, B., & Nikutta, R. 2017, ApJ, 843, 3

Antonucci, R. 1993, ARA&A, 31, 473

Barvainis, R. 1987, ApJ, 320, 537

Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002

Bevington, P. R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York: McGraw-Hill)

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Birtscher, L., Meisenheimer, K., Tristram, K. R. W., et al. 2013, A&A, 558, A149

Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560

Clavel, J., Wamsteker, W., & Glass, I. S. 1989, ApJ, 337, 236

Czerny, B., & Hryniewicz, K. 2011, A&A, 525, L8

Done, C., Davis, S. W., Jin, C., Blaes, O., & Ward, M. 2012, MNRAS, 420, 1848

Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2009, ApJ, 696, 870

Elvis, M., Wilkes, B. J., McDowell, J. C., et al. 1994, ApJS, 95, 1

Gaskell, C. M., & Peterson, B. M. 1987, ApJS, 65, 1

Glass, I. S. 1992, MNRAS, 256, 23P

Glass, I. S. 2004, MNRAS, 350, 1049

GRAVITY Collaboration, Dexter, J., Shangguan, J., et al. 2019, arXiv:1910.00593

Grier, C. J., Shen, Y., Horne, K., et al. 2019, arXiv:1904.03199

Grier, C. J., Trump, J. R., Shen, Y., et al. 2017, ApJ, 851, 21

Haardt, F., & Maraschi, L. 1991, ApJL, 380, L51

Homayouni, Y., Trump, J. R., Grier, C. J., et al. 2019, ApJ, 880, 126

Hönig, S. F., Beckert, T., Ohnaka, K., & Weigelt, G. 2006, A&A, 452, 459

Hönig, S. F., & Kishimoto, M. 2011, A&A, 534, A121

Hopkins, P. F., Richards, G. T., & Hernquist, L. 2007, ApJ, 654, 731

I-Hsiu, Li, J., Shen, Y., Brandt, W. N., et al. 2019, arXiv:1909.03092

Jaffe, W., Meisenheimer, K., Röttgering, H. J. A., et al. 2004, Natur, 429, 47

Jarrett, T. H., Cohen, M., Masci, F., et al. 2011, ApJ, 735, 112

Jiang, Y.-F., Green, P. J., Greene, J. E., et al. 2017, ApJ, 836, 186

Kelly, B. C. 2007, ApJ, 665, 1489

Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895

Kishimoto, M., Antonucci, R., Blaes, O., et al. 2008, Natur, 454, 492

Kishimoto, M., Hönig, S. F., Antonucci, R., et al. 2009, A&A, 507, L57

Kishimoto, M., Hönig, S. F., Antonucci, R., et al. 2011, A&A, 527, A121

Kishimoto, M., Hönig, S. F., Beckert, T., & Weigelt, G. 2007, A&A, 476, 713

Koshida, S., Minezaki, T., Yoshii, Y., et al. 2014, ApJ, 788, 159

Koshida, S., Yoshii, Y., Kobayashi, Y., et al. 2009, ApJL, 700, L109

Koshida, S., Yoshii, Y., Kobayashi, Y., et al. 2017, ApJL, 842, L13

Kozłowski, S., Kochanek, C. S., Ashby, M. L. N., et al. 2016, ApJ, 817, 119

Kozłowski, S., Kochanek, C. S., Udalski, A., et al. 2010, ApJ, 708, 927

Law, N. M., Kulkarni, S. R., Dekany, R. G., et al. 2009, PASP, 121, 1395

LSST Science Collaborations, Marshall, P., Anguita, T., et al. 2017, arXiv:1708.04058

Lyu, J., Rieke, G. H., & Smith, P. S. 2019, arXiv:1909.11101

MacLeod, C. L., Ivezić, Ž., Kochanek, C. S., et al. 2010, ApJ, 721, 1014

MacLeod, C. L., Ivezić, Ž., Sesar, B., et al. 2012, ApJ, 753, 106

Mainzer, A., Bauer, J., Cutri, R. M., et al. 2014, ApJ, 792, 30

Mainzer, A., Grav, T., Bauer, J., et al. 2011, ApJ, 743, 156

Minezaki, T., Yoshii, Y., Kobayashi, Y., et al. 2004, ApJL, 600, L35

Minezaki, T., Yoshii, Y., Kobayashi, Y., et al. 2019, arXiv:1910.08722

Mudd, D., Martini, P., Zu, Y., et al. 2018, ApJ, 862, 123

Nelson, B. O. 1996, ApJL, 465, L87

Nenkova, M., Sirocky, M. M., Nikutta, R., Ivezić, Ž., & Elitzur, M. 2008, ApJ, 685, 160

Netzer, H. 2015, ARA&A, 53, 365

Oknyanskii, V. L. 1993, AstL, 19, 416

Oknyanskij, V. L., Lyuty, V. M., Taranova, O. G., & Shenavrin, V. I. 1999, AstL, 25, 483

Penston, M. V., Penston, M. J., Neugebauer, G., et al. 1971, MNRAS, 153, 29

Peterson, B. M., Wanders, I., Horne, K., et al. 1998, PASP, 110, 660

Richards, G. T., Strauss, M. A., Fan, X., et al. 2006, AJ, 131, 2766

Rieke, G. H. 1978, ApJ, 226, 550

Rumbaugh, N., Shen, Y., Morganson, E., et al. 2018, ApJ, 854, 160

Sanders, D. B., Phinney, E. S., Neugebauer, G., Soifer, B. T., & Matthews, K. 1989, ApJ, 347, 29

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48

Shen, Y. 2012, ApJ, 757, 152

Shen, Y., Brandt, W. N., Dawson, K. S., et al. 2015, ApJS, 216, 4

Shen, Y., Grier, C. J., Horne, K., et al. 2019a, arXiv:1908.00027

Shen, Y., Hall, P. B., Horne, K., et al. 2019b, ApJS, 241, 34

Shen, Y., Horne, K., Grier, C. J., et al. 2016, ApJ, 818, 30

Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45

Shen, Y., Strauss, M. A., Ross, N. P., et al. 2009, ApJ, 697, 1656

Sheng, Z., Wang, T., Jiang, N., et al. 2017, ApJL, 846, L7

Starkey, D. A., Horne, K., & Villforth, C. 2016, MNRAS, 456, 1960

Suganuma, M., Yoshii, Y., Kobayashi, Y., et al. 2006, ApJ, 639, 46

Sun, M., Grier, C. J., & Peterson, B. M. 2018, PyCCF: Python Cross-Correlation Function for Reverberation Mapping Studies, Astrophysics Source Code Library, ascl:1805.032

Swain, M., Vasisht, G., Akeson, R., et al. 2003, ApJL, 596, L163

Tristram, K. R. W., Meisenheimer, K., Jaffe, W., et al. 2007, A&A, 474, 837

Urry, C. M., & Padovani, P. 1995, PASP, 107, 803

Vazquez, B., Galianni, P., Richmond, M., et al. 2015, ApJ, 801, 127

Weigelt, G., Hofmann, K. H., Kishimoto, M., et al. 2012, A&A, 541, L9

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
Yang, Q., Shen, Y., Liu, X., et al. 2019, ApJ, 885, 110

York, D. G., Adelman, J., Anderson, J. E. J., et al. 2000, AJ, 120, 1579

Yoshii, Y., Kobayashi, Y., Minezaki, T., Koshida, S., & Peterson, B. A. 2014, ApJL, 784, L11

Yu, Z., Martini, P., Davis, T. M., et al. 2020, ApJS, 246, 16

Zu, Y., Kochanek, C. S., Kozłowski, S., & Udalski, A. 2013, ApJ, 765, 106

Zu, Y., Kochanek, C. S., & Peterson, B. M. 2011, ApJ, 735, 80