# The SDSS-V Black Hole Mapper Reverberation Mapping Project: Unusual Broad-
Line Variability in a Luminous Quasar

Logan B. Fries Department of Physics, 196A Auditorium Road, Unit 3046, University of Connecticut, Storrs, CT 06269, USA  Jonathan R. Trump Department of Physics, 196A Auditorium Road, Unit 3046, University of Connecticut, Storrs, CT 06269, USA  Megan C. Davis Department of Physics, 196A Auditorium Road, Unit 3046, University of Connecticut, Storrs, CT 06269, USA  C. J. Grier Steward Observatory, The University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA  Department of Astronomy, University of Wisconsin-Madison, Madison, WI 53706, USA  Yue Shen Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  Scott F. Anderson Astronomy Department, University of Washington, Box 351580, Seattle, WA 98195, USA  Tom Dwelly Max-Planck-Institut für extraterrestrische Physik, Giessenbachstraße, 85748 Garching, Germany  Michael Eracleous Department of Astronomy & Astrophysics and Institute for Gravitation and the Cosmos, The Pennsylvania State University, 525 Davey Lab, University Park, PA 16802, USA  Y. Homayouni Department of Astronomy & Astrophysics and Institute for Gravitation and the Cosmos, The Pennsylvania State University, 525 Davey Lab, University Park, PA 16802, USA  Keith Horne SUPA Physics and Astronomy, University of St. Andrews, Fife, KY16 9SS, Scotland, UK  Mirko Krumpe  Leibniz-Institut für Astrophysik (AIP), An der Sternwarte 16, 14482 Potsdam, Germany  Sean Morrison Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  Jessie C. Runnoe Department of Physics and Astronomy, Vanderbilt University, Nashville, TN 37235, USA  Benny Trakhtenbrot School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel  Roberto J. Assef Instituto de Estudios Astrofísicos, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile  W. N. Brandt Department of Astronomy & Astrophysics, 525 Davey Lab, The Pennsylvania State University, University Park, PA 16802, USA  Institute for Gravitation and the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA  Department of Physics, 104 Davey Lab, The Pennsylvania State University, University Park, PA 16802, USA  Joel Brownstein Department of Physics and Astronomy, University of Utah, 115 S. 1400 E., Salt Lake City, UT 84112, USA  Collin Dabbieri Department of Physics and Astronomy, Vanderbilt University, Nashville, TN 37235, USA  Alexander Fix Department of Astrophysical & Planetary Sciences, University of Colorado, 2000 Colorado Ave, Boulder, CO 80309, USA  Gloria Fonseca Alvarez Department of Physics, 196A Auditorium Road, Unit 3046, University of Connecticut, Storrs, CT 06269, USA  NSF’s National Optical-Infrared Astronomy Research Laboratory, 950 N. Cherry Ave., Tucson, AZ 85719, USA  Sara Frederick Department of Physics and Astronomy, Vanderbilt University, Nashville, TN 37235, USA  P. B. Hall Department of Physics and Astronomy, York University, Toronto, ON M3J 1P3, Canada  Anton M. Koekemoer Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA  Jennifer I-Hsiu Li Department of Astronomy, University of Michigan, Ann Arbor, MI, 48109, USA  Xin Liu Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  Mary Loli Martínez-Aldama Instituto de Física y Astronomía, Facultad de Ciencias, Universidad de Valparaíso, Gran Bretaña 1111, Valparaíso, Chile  Departamento de Astronomía, Universidad de Chile, Casilla 36D, Santiago, Chile  Claudio Ricci Instituto de Estudios Astrofísicos, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile  Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People’s Republic of China  Donald P. Schneider Department of Astronomy & Astrophysics and Institute for Gravitation and the Cosmos, The Pennsylvania State University, 525 Davey Lab, University Park, PA 16802, USA  Hugh W. Sharp Department of Physics, 196A Auditorium Road, Unit 3046, University of Connecticut, Storrs, CT 06269, USA  Matthew J. Temple Instituto de Estudios Astrofísicos, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército Libertador 441, Santiago, Chile  Qian Yang Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Grisha Zeltyn School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel  Dmitry Bizyaev Apache Point Observatory and New Mexico State University, P.O. Box 59, Sunspot, NM, 88349-0059, USA  Sternberg Astronomical Institute, Moscow State University, Moscow

###### Abstract

We present a high-cadence multi-epoch analysis of dramatic variability of
three broad emission lines (MgII, H$\beta$, and H$\alpha$) in the spectra of
the luminous quasar ($\lambda L_{\lambda}\rm(5100\mbox{\AA})$ = $4.7\times
10^{44}$ erg s-1) SDSS J141041.25+531849.0 at $z=0.359$ with 127 spectroscopic
epochs over 9 years of monitoring (2013-2022). We observe anti-correlations
between the broad emission-line widths and flux in all three emission lines,
indicating that all three broad emission lines “breathe” in response to
stochastic continuum variations. We also observe dramatic radial velocity
shifts in all three broad emission lines, ranging from $\Delta{v}$ $\sim$400
km s-1 to $\sim$800 km s-1, that vary over the course of the monitoring
period. Our preferred explanation for the broad-line variability is complex
kinematics in the broad-line region gas. We suggest a model for the broad-line
variability that includes a combination of gas inflow with a radial gradient,
an azimuthal asymmetry (e.g., a hot spot), superimposed on the stochastic
flux-driven changes to the optimal emission region (“line breathing”). Similar
instances of line-profile variability due to complex gas kinematics around
quasars are likely to represent an important source of false positives in
radial velocity searches for binary black holes, which typically lack the kind
of high-cadence data we analyze here. The long-duration, wide-field, and many-
epoch spectroscopic monitoring of SDSS-V BHM-RM provides an excellent
opportunity for identifying and characterizing broad emission-line
variability, and the inferred nature of the inner gas environment, of luminous
quasars.

††software: AstroPy (Astropy Collaboration et al. 2013, 2018), Matplotlib
(Hunter 2007), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020)

##  1 Introduction

Active Galactic Nuclei (AGN) are powered by supermassive black holes that are
actively accreting matter at the centers of massive galaxies (Lynden-Bell
1969). A hallmark signature of many AGN is the existence of broad emission
lines, as first described by Seyfert (1943). Such lines arise from
photoionization of the gas of the broad-line region (BLR), which is a
distribution of gas that is thought to be moving in virialized orbits (Bentz
et al. 2009; Shapovalova, A. I. et al. 2010; Barth et al. 2011a, b; Grier et
al. 2013) close to the central black hole. As the radius of the BLR is on the
order of light-days, it is difficult to spatially resolve with current
technology. To date, there have been only a handful of studies that have been
able to spatially resolve the BLR using near-infrared interferometry (Gravity
Collaboration et al. 2020). Investigating the structure of the BLR generally
requires indirect techniques like reverberation mapping (Blandford & McKee
1982; Peterson 1993; Cackett et al. 2021).

Reverberation mapping (RM) is a technique that utilizes the fact that
variations in the continuum flux of the AGN are followed by variations in the
broad-emission lines from the BLR, with a time delay, $\tau$, that corresponds
to the light-travel time between the continuum-emitting region and the BLR.
The time delay is thus related to the typical radius of the BLR by the
relation $R_{\rm BLR}=c\tau$. If we assume that the BLR is virialized, then
the mass of the central black hole is determined by the virial product,
namely,

$M_{\rm{RM}}=f\frac{v^{2}R_{\rm BLR}}{G}$ (1)

where $v$ is the velocity of the BLR gas, $R_{\rm BLR}$ is the BLR radius, and
$G$ is the gravitational constant. A dimensionless factor $f$ is introduced
into the equation to parameterize the inclination angle, unknown geometry, and
orientation of the BLR. The dimensionless factor $f$ is calibrated by
comparison with spatially resolved BLR measurements (Sturm et al. 2018), the
kinematics of stars and gas near the AGN (e.g., Grier et al. 2013; Woo et al.
2015), and/or velocity-resolved reverberation mapping (e.g., Pancoast et al.
2014; Grier et al. 2017a).

Both observations (Wilhite et al. 2006; Park et al. 2012; Barth et al. 2015;
Dexter et al. 2019; Homan et al. 2020; Wang et al. 2020) and photoionization
modeling (Korista & Goad 2004; Cackett & Horne 2006; Goad & Korista 2014) have
demonstrated that an increase in central flux from the AGN accretion disk
results in an increase in the radius of the BLR “optimal emitting region” and
a corresponding decrease in the orbital velocity of the gas since the line
emission originates further out. This phenomenon is known as “line breathing”
and manifests itself as an anti-correlation between the broad-line flux and
width of the BLR over time. Quantifying the breathing behaviors of broad
emission lines provides valuable constraints on the geometry, kinematics, and
photoionzation of the BLR (Wang et al. 2020).

The continuum flux of quasars exhibits variability of order $\sim$10% on
timescales of weeks to years that is thought to be driven by thermal
fluctuations in the accretion disk (Ulrich et al. 1997; Vanden Berk et al.
2004; Kelly et al. 2009; MacLeod et al. 2010). Periodicity in quasar
photometric lightcurves has been a popular method for searching for SMBH
binary candidates (Valtonen et al. 2008; Graham et al. 2015; Ackermann et al.
2015; Li et al. 2016; Charisi et al. 2016; Sandrinelli et al. 2018; Severgnini
et al. 2018; Li et al. 2019; Chen et al. 2020; Liao et al. 2021; Zhang 2022).
However, Vaughan et al. (2016) demonstrated that the stochastic variability of
single (non-binary) quasars can resemble a periodic signal that can span a few
periods, leading to false-positive identification of SMBH binaries. An
alternative method for identifying binary SMBHs is searching for periodic
broad emission-line radial-velocity shifts (Gaskell 1983; Loeb 2010; Eracleous
et al. 2012; Shen et al. 2013; Liu et al. 2014; Runnoe et al. 2017; Guo et al.
2019) akin to spectroscopic binary stars. This method requires high-cadence,
long-term spectroscopic monitoring as the binary orbits are on the order of
decades (Yu 2002).

Figure 1: Spectroscopic light curves for the rest-frame 5100Å continuum and
the three emission lines (H$\beta$, H$\alpha$, MgII). All 3 of the lines are
brightest in the 2021-2022 monitoring period, although H$\beta$ and H$\alpha$
have more extreme brightening than MgII. The emission-line variability
generally appears to follow the continuum variability with lags consistent
with the previously measured $\tau\sim 30$ days for the Balmer lines (Grier et
al. 2017b) and $\tau\sim 145$ days for MgII (Homayouni et al. 2020).

Radial velocity shifts of broad emission lines can also result from recoiling
BHs (Herrmann et al. 2007; Eracleous et al. 2012), tidal disruption events
(Gezari 2021), and gas outflows/inflows (Brotherton et al. 1994; Storchi-
Bergmann 2010; Rakshit & Woo 2018; Kovačević-Dojčinović et al. 2022). In
addition, Barth et al. (2015) found that H$\beta$ velocity centroids can
undergo dramatic changes, on timescales of a month, in response to continuum
flux variations. The radial velocity shifts in this case are a product of
asymmetric reverberation by the BLR, and can appear as false-positive
detections in binary black hole searches.

In this paper, we present observations of an AGN with unusual broad emission-
line variability, SDSS J141041.25+531849.0 (hereafter RM160), found within a
broader search for variability in broad emission-line profiles in the recently
started Sloan Digital Sky Survey V (SDSS-V, Kollmeier et al. 2017, Almeida et
al. 2023). Section 2 describes the sample selection, our criteria to identify
quasars with unusual variability in their broad emission-line shapes, and
provide the general characteristics of the object of interest. Section 3
describes the methods we use to quantify the broad emission-line profiles.
Section 4 describes the broad emission-line profile changes and presents a
physical model to explain the observations. Section 5 summarizes our results.

Throughout this work, we assume a $\Lambda$CDM cosmology with
$\Omega_{\Lambda}$ = 0.7, $\Omega_{M}$ = 0.3, and $H_{0}$ = 70 km s-1 Mpc-1.

##  2 Observations and Parent Sample

###  2.1 Data

The data are from the third (Eisenstein et al. 2011), fourth (Blanton et al.
2017), and the ongoing fifth generation (Kollmeier et al. 2017, Almeida et al.
2023) of the Sloan Digital Sky Survey (SDSS, York et al. 2000). The data were
obtained using the plate-based, fiber-fed BOSS spectrograph (Smee et al. 2013)
mounted on the 2.5m SDSS telescope (Gunn et al. 2006) at the Apache Point
Observatory. The spectrograph has a dual-channel design, a blue channel (3000Å
< $\lambda$ < 6350Å) and a red channel (5650Å < $\lambda$ < 10,400Å), both
with a spectral resolution of $R\sim 2000$. The SDSS-III and SDSS IV data
(2013-2020) were reduced with the v5_13_0 version of idlspec2d and the SDSS-V
data (2021-present) were reduced with the v6_0_9 version of idlspec2d, the
SDSS BOSS spectroscopic reduction pipeline (Bolton et al. 2012).

The spectroscopic monitoring spans a range of 9 years (2013 - 2022) with 127
epochs. An ‘epoch’ generally represents observations taken in a single night,
but in some cases epochs will include observations coadded from up to 3
consecutive nights in order to pass the ‘epoch-completion threshold’ defined
as a blue-channel based signal-to-noise ratio (SNR) threshold of ${\rm
SNR}^{2}(g)>20$ for a target of fiducial point spread function (PSF) magnitude
$g=22$. This observing design aims to maintain a constant and useful SNR for
all epochs, although some epochs have lower SNR because they could not be
completed (i.e., pass the epoch completion threshold) within 3 nights.

Figure 1 shows the spectroscopic light curves 111The observables shown here
are explained further in Section 2.4 for our quasar of interest, RM160, as a
demonstration of the SDSS Reverberation Mapping (SDSS-RM) and SDSS-V Black
Hole Mapper Reverberation Mapping (BHM-RM) data sets. The highest density of
monitoring from SDSS-III and SDSS-IV occurred in 2014 (30 epochs) and at the
start of SDSS-V in 2021 (39 epochs). These light curves include a second-order
calibration in flux and wavelength using the [Oiii]$\lambda$5007 narrow
emission line (see Section 3.3 for details).

###  2.2 Parent Sample

The parent sample for our broader unusual variability search consists of all
320 quasars that have been monitored by both the SDSS-V BHM-RM program (for
details: see Trump et al. in prep.) and SDSS-RM (for details: see Shen et al.
2015, 2019). These targets lie within the SDSS-RM field, which is a 7 deg2
field that has been observed as a part of SDSS-RM in SDSS-III and SDSS-IV from
2013-2020, and then by BHM-RM in SDSS-V since 2021 (with monitoring scheduled
to continue through at least 2026). The parent sample spans a redshift range
of $0.1~{}<~{}z~{}<~{}4.34$ and is magnitude-limited by $i_{\rm PSF}<21.7$.
The median redshift of the parent sample is $z_{\rm med}=1.52$ and the median
i-band magnitude of the parent sample is $\textit{i}_{\rm med}=21.09$.

###  2.3 Identifying Unusual Line Profile Changes in SDSS-RM / BHM-RM Quasars

Our object of interest, RM160, was found during a broad search for quasars
with variability in their broad emission-line profiles. To measure
variability, we quantified changes to the broad emission-line profiles in each
epoch using the reduced chi-squared
$\chi_{\nu}^{2}=\frac{1}{\nu}\sum{\frac{\left(O_{i}-M_{i}\right)^{2}}{U_{i}^{2}}}$,
with $\nu$ the degrees of freedom, $O_{i}$ the observed flux density, $M_{i}$
the model flux density (described below), and $U_{i}$ the observational
uncertainty summed across a wavelength range indexed by $i$.

We use the median spectrum across all 127 epochs as our model flux density in
our chi-squared calculation. We seek to identify changes in the line-profile
shape, rather than the (commonly observed) brightening or dimming of the
overall broad line. To accomplish this, we allow the median spectrum in the
line-profile region to scale up and down by a multiplicative factor computed
from the ratio of the median flux density of the line-profile region at each
epoch divided by the average flux density of the line-profile region across
all epochs. Figure 2 shows an example spectrum for RM160, where the blue line
shows the median spectrum centered on the emission-line region and spanning a
total width of 1.5 $\times$ FWHM reported in Shen et al. (2019), which is the
extent to which we measure variability of the broad emission-line profile.

We identified interesting candidates using the 90th percentile of the
distribution of the $\chi_{\nu}^{2}$ values from each epoch, choosing
$\chi_{\nu}^{2}(90\%)$ > 6 as a threshold associated with the tail of line-
profile variability. Figure 3 shows the distribution of H$\beta$
$\chi_{\nu}^{2}(90\%)$ of our sample.

There are 15 objects that fit our $\chi_{\nu}^{2}(90\%)$ criterion. We
visually inspected the variability of the 15 objects by creating animations of
the time-variable spectra for all 127 epochs. Out of those 15 objects, we
identified the subject of this study, RM160. The $\chi_{\nu}^{2}(90\%)$ for
RM160 is 8.44. We noticed from visual inspection that the H$\beta$, H$\alpha$,
and MgII broad emission lines of RM160 all appear to undergo significant
velocity shifts over time. The apparent variability of the remaining objects
were largely the result of noise.

Figure 2: The H$\beta$ spectral region for RM160 in Epoch 111 (May 18, 2021
corresponding to an MJD of 59352). The dark blue line is the median spectrum
computed for the H$\beta$ line centered on the line center and has a total
width of 1.5 $\times$ FWHM reported by Shen et al. (2019). In this particular
epoch the H$\beta$ emission-line profile is redder than the median profile.
The online version of this figure is an animation. The animation is 1:00
minutes long and shows the time evolution of RM160 for all 127 epochs of
spectra beginning in April 11, 2013 and ending on June 03, 2022.  Figure 3:
Distribution of the $\chi_{\nu}^{2}(90\%)$ measured from the multi-epoch
H$\beta$ region spectra for the 320 objects in our parent SDSS-RM sample. The
vertical, dashed purple line indicates our criterion ($\chi_{\nu}^{2}(90\%)$
$\geq$ 6) for identifying unusual line-profile variability. The
$\chi_{\nu}^{2}(90\%)$ for RM160 is 8.44 and was chosen as it had the highest
SNR and the most apparent broad emission-line profile variability, while the
other ones that satisfied our $\chi_{\nu}^{2}(90\%)$ criterion were either
noisy and/or had less apparent variability in their broad emission-line
profiles.

###  2.4 The Source of Interest: RM160

The subject of this study is RM160, a luminous quasar ($\lambda
L_{\lambda}\rm(5100\mbox{\AA})$= $4.7\times 10^{44}$ erg s-1) in the SDSS-RM
field. It has a redshift of $z=0.359$ and an i-band magnitude of $i=19.68$.

There are published H$\beta$, H$\alpha$, and MgII RM time lags for this
object. The H$\beta$ and H$\alpha$ lags were measured by Grier et al. (2017b)
using only the 2014 data. The observed-frame H$\beta$ lag is
$\tau_{\rm{H\beta,}\rm{obs}}=31.3^{+8.1}_{-4.1}$ days and the observed-frame
H$\alpha$ lag is $\tau_{\rm{H\alpha,}\rm{obs}}=27.7^{+5.3}_{-4.7}$ days. The
observed-frame MgII lags were measured in Homayouni et al. (2020) using 4
years of data (2014-2017). The observed-frame MgII lag is
$\tau_{\rm{MgII,}\rm{obs}}=144.7^{+24.7}_{-22.6}$ days. We note that the MgII
lag for this object has a false-positive rate of 16%, which is not in the
"gold-sample" (false-positive rate of $\leq$10%) of Homayouni et al. (2020).
As such, the MgII lag for RM160 may be unreliable. In general, we assume that
the MgII lag is longer than the lags of the Balmer lines (i.e., we assume the
BLR is stratified; see Clavel et al. 1991; Reichert et al. 1994). The black
hole mass of RM160 was computed in Grier et al. (2017b) to be $(M_{\rm
BH}/10^{7}M_{\odot})=7.0^{+1.7}_{-1.3}$.

From HST imaging (taken on 2020 September 28), the host-galaxy contribution
for RM160 is 14% in F606W (Li et al. 2023). This measurement does not take
into account the 2″ SDSS fiber which would make the host-fraction even
smaller, thus we do not account for the marginal host contribution in our
analysis and we assume that the spectrum is dominated by the quasar.

Figure 1 shows the continuum and emission-line light curves for RM160. The
light curve behavior appears to be qualitatively consistent with the
previously measured lags, with the Balmer-line variability appearing to follow
the same pattern as the continuum after a lag of $\sim$30 days and the MgII
variability following the same pattern after an additional $\sim$145 days.

Figure 4: Median spectra for MgII, H$\beta$, and H$\alpha$ of RM160. The top
row shows the full spectrum in black along with the narrow emission-line fits
in red and the continuum fits in orange. The bottom row shows the continuum
and narrow-line subtracted median spectra for MgII, H$\beta$, and H$\alpha$,
respectively. Our continuum and narrow-line subtracted spectra show clear
broad-line profiles (the residuals from non-Gaussian narrow lines do not
affect the broad-line measurements).

##  3 Quantifying the Emission-Line Profiles

###  3.1 Fitting the Continuum

To isolate and model the broad emission-line regions, we first need to
subtract the continuum from the spectra. To subtract the continuum, we fit a
first-order polynomial to the spectrum based on the median of the continuum
over 50 pixels from the line-free regions blueward and redward of the broad
lines (and nearby narrow lines for H$\beta$ and H$\alpha$). The continuum fits
for the median spectrum are shown in yellow in the top panels in Figure 4.
Appendix A presents alternative modeling with PyQSOFit (Guo et al. 2018; Shen
et al. 2019) that separately fits the continuum and iron psuedo-continuum and
finds consistent variability patterns in the fitted line properties.

###  3.2 Fitting the Narrow Emission Lines

We fit each narrow emission line in the H$\alpha$ and H$\beta$ regions with a
single Gaussian. We assume that the narrow lines are constant over the course
of our monitoring period of $\sim$9 years. In Section 3.3, we confirm this
assumption by demonstrating that the $\hbox{{\rm[O}\kern 1.00006pt{\sc
iii}{\rm]}}\lambda$5007 narrow-line flux is constant over the monitoring
period.

We used the $\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007
parameters to constrain the other narrow-line fits. Specifically, we tied the
line centers of the narrow-line Gaussian fits, for each epoch, to
$\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007 using the narrow-
line wavelength centers from Vanden Berk et al. (2001). We also tied the line
widths of the narrow-line Gaussian fits, for each epoch, to the line width
from the $\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007 fit. The
$\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$4959 line flux was also
constrained to be 1/3 that of the $\hbox{{\rm[O}\kern 1.00006pt{\sc
iii}{\rm]}}\lambda$5007 line flux (Storey & Zeippen 2000).

For each individual fit for the H$\alpha$ region, we fit the H$\alpha$ narrow
line without constraints on the amplitude and found that this resulted in a
clean-looking residual (i.e., a clear broad emission line profile with the
narrow emission lines cleanly subtracted out). However, this unconstrained
approach to the narrow H$\beta$ line resulted in poor fits with large apparent
residuals. Therefore, we constrained the amplitude of the narrow H$\beta$ line
to a value that produced a smooth broad-line residual that lacked a cuspy
narrow-line peak in the median spectrum fit. We then applied that narrow
H$\beta$ amplitude to the fits for all epochs.

Figure 5: The [Oiii]$\lambda$5007 flux (top), line center (middle), and line
width $\sigma$ (bottom) measured from single-Gaussian fits of the continuum
subtracted spectra at each epoch. In all panels, the black line indicates the
median and the gray shading represents the NMAD. We use the apparent changes
in $\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007 flux and line
center for a second-order flux and wavelength calibration for the spectra of
each epoch.

Figure 4 illustrates the process of fitting and subtracting both the continuum
and the narrow emission lines. There are limitations to modeling complex
spectra with a single Gaussian. For example, the $\hbox{{\rm[O}\kern
1.00006pt{\sc iii}{\rm]}}\lambda$5007 line appears to have a blue wing that
cannot be modeled using a single Gaussian and would be better fit with two or
more Gaussians. This asymmetric, non-Gaussian profile has been shown to be due
to outflowing, ionized gas in the narrow-line region (NLR) (Rojas et al. 2020;
Ayubinia et al. 2022; Molina et al. 2022). However, these small residuals
caused by non-Gaussian shapes in the narrow emission lines do not affect the
measured broad-line profiles. We visually inspected the fits for each broad
emission line in each epoch and confirmed that our method produces clean
broad-line profiles that have effective subtractions of the continuum and
narrow emission lines.

###  3.3 Second-Order Calibrations Based on Narrow Emission Lines

We investigated the stability of the narrow emission lines to examine and
improve the flux and wavelength calibration. The $\hbox{{\rm[O}\kern
1.00006pt{\sc iii}{\rm]}}\lambda$5007 line is observed to be stable over
timescales of a few years (Foltz et al. 1981; Peterson et al. 1982) and so it
is often used as a flux-calibration standard in AGN spectra. We examine the
possibility of applying a spectrophotometric calibration by fitting the
$\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007 line with a single
Gaussian in each epoch, which is separate from the aforementioned fitting
procedure in Section 3.2. In Figure 5, we show the Gaussian fit parameters for
$\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007 line flux, center,
and width, as well as the respective median values (black) and normalized
median absolute deviations (NMAD; gray regions). Figure 5 confirms our
assumption of a non-variable $\hbox{{\rm[O}\kern 1.00006pt{\sc
iii}{\rm]}}\lambda$5007 line throughout the course of our monitoring period.
We note that in the case of NGC 5548, the $\hbox{{\rm[O}\kern 1.00006pt{\sc
iii}{\rm]}}\lambda$5007 line has been shown to vary over a timescale of
$\sim$30 years by $\sim$10% (Peterson et al. 2013).

Since the $\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007 narrow
line does not vary over these timescales, the changes in $\hbox{{\rm[O}\kern
1.00006pt{\sc iii}{\rm]}}\lambda$5007 flux represent epoch-dependent changes
in the spectrophotometric calibration, with 8 epochs in 2021 that fall well
below the median. The $\hbox{{\rm[O}\kern 1.00006pt{\sc
iii}{\rm]}}\lambda$5007 line center shifts by exactly 1 or 2 pixels in 3
epochs, indicating a shift in the wavelength calibration. The
$\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007 line width is
constant within its uncertainties, indicating that the spectral resolution is
stable throughout the observations.

We observe similar changes in the fitted line fluxes and centers for the [OII]
and [Sii] emission lines (with larger uncertainties for these weaker lines).
This suggests that the flux and wavelength changes are gray (not wavelength-
dependent) and systematic (not intrinsic to RM160). In other words, the flux
and wavelength changes observed for $\hbox{{\rm[O}\kern 1.00006pt{\sc
iii}{\rm]}}\lambda$5007 represent gray calibration issues for the entire
spectrum.

We perform a second-order flux and wavelength calibration that forces the
$\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007 flux and wavelength
to be constant across all epochs and apply it to the spectra at each epoch. We
scale the spectrum at each epoch by a factor of $\tilde{f}\hbox{{\rm[O}\kern
1.00006pt{\sc iii}{\rm]}}/f\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}$ and
we correct the wavelength of each spectrum by a factor of
$\tilde{\mu}\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}/\mu\hbox{{\rm[O}\kern
1.00006pt{\sc iii}{\rm]}}$ , where $\tilde{f}\hbox{{\rm[O}\kern 1.00006pt{\sc
iii}{\rm]}}$ and $\tilde{\mu}\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}$ are
the median $\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007 flux and
line center, respectively across all epochs and $f\hbox{{\rm[O}\kern
1.00006pt{\sc iii}{\rm]}}$ and $\mu\hbox{{\rm[O}\kern 1.00006pt{\sc
iii}{\rm]}}$ are the $\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007
flux and line center at each epoch.

###  3.4 Quantifying the Broad Emission-Line Profile Variability

We measure the broad emission-line properties of the continuum and narrow-line
subtracted spectra using the moments of a distribution:

  * •

Line Flux:
$f=\sum_{\lambda_{1}}^{\lambda_{2}}{f_{\lambda}(\lambda)\Delta\lambda}$

  * •

Line Center: $\lambda_{\rm{C}}=\frac{\sum_{\lambda_{1}}^{\lambda_{2}}{\lambda
f_{\lambda}(\lambda)\Delta\lambda}}{\sum_{\lambda_{1}}^{\lambda_{2}}{f_{\lambda}(\lambda)\Delta\lambda}}$

  * •

Line Width:
$\rm{\sigma}=\sqrt{\frac{\sum_{\lambda_{1}}^{\lambda_{2}}{(\lambda_{\rm{C}}-\mu)^{2}{\textit{f}_{\lambda}(\lambda)\Delta\lambda}}}{\sum_{\lambda_{1}}^{\lambda_{2}}{\textit{f}_{\lambda}(\lambda)\Delta\lambda}}}$

  * •

FWHM: We measured FWHM by first applying a median boxcar smoothing of 5 pixels
to our spectra at each epoch. We then used the FWHM routine from the specutils
package (Earl et al. 2022) to compute the FWHM.

Here $f_{\lambda}$ is the flux density, $\lambda$ is the wavelength,
$\Delta\lambda$ is the wavelength per pixel, and $\lambda_{1}$ and
$\lambda_{2}$ for each line are from Vanden Berk et al. (2001). These non-
parametric measurements were chosen as they make no assumptions about the
underlying shape of the broad line profiles. We used the non-parametric model
on the continuum and narrow-line subtracted spectra (bottom row panel in
Figure 4).

We employ a Monte Carlo resampling approach to estimate the uncertainties in
line flux, line-center, and line width (both $\sigma$ and FWHM). For each
epoch, we measured these quantities for 200 resampled spectra using the
corresponding error spectrum. We adopted the standard deviations in the
measured parameters from the 200 re-sampled spectra as our uncertainties for
each epoch.

Figure 6: Variability of the MgII broad emission-line profile, as quantified
by our non-parametric measurements of the line flux (top row), FWHM (2nd row),
line-width $\sigma$ (3rd row), and the line-center velocity shift $\Delta{v}$
(bottom row). The units for FWHM, line-width $\sigma$, and line-center
velocity shift $\Delta{v}$ are all in km s-1. The MgII line gets brighter and
narrower, fainter and broader and then brighter and narrower again. Meanwhile,
the line-center starts very red (compared to the systemic redshift from the
narrow emission lines), gets bluer and then gets as red as it was at the
beginning. The velocity shifts in line-center do not appear to be
simultaneously correlated with the changes in flux and line-width.  Figure 7:
Variability of the H$\beta$ broad emission-line profile, as quantified by our
non-parametric measurements of the line flux (top row), FWHM (2nd row), line-
width $\sigma$ (3rd row), and the line-center velocity shifts $\Delta{v}$
(bottom row). The units for FWHM, line-width $\sigma$, and line-center
velocity shift $\Delta{v}$ are all in km s-1. The H$\beta$ line gets brighter
and narrower, fainter and broader and then brighter and narrower again.
Meanwhile, the line-center starts very red (compared to the systemic redshift
from the narrow emission lines), gets bluer and then gets more redder and
seems to plateau after the 2019 monitoring period. The velocity shifts in
line-center do not appear to be simultaneously correlated with the changes in
flux and line-width.  Figure 8: Variability of the H$\alpha$ broad emission-
line profile, as quantified by our non-parametric measurements of the line
flux (top row), FWHM (2nd row), line-width $\sigma$ (3rd row), and the line-
center velocity shifts $\Delta{v}$ (bottom row). The units for FWHM, line-
width $\sigma$, and line-center velocity shift $\Delta{v}$ are all in km s-1.
The H$\alpha$ line gets brighter and narrower, fainter and broader and then
brighter and narrower again. Meanwhile, the line-center starts very red
(compared to the systemic redshift from the narrow emission lines), gets bluer
and then gets more redder, but not as red as it was in the 2014 monitoring
period. The velocity shifts in line-center do not appear to be simultaneously
correlated with the changes in flux and line-width.

To test the robustness of our non-parametric measurements, we employed two
alternative parametric fitting procedures: (1) using single Gaussians to model
the broad emission lines and (2) using PyQSOFit which models the Fe ii psuedo-
continuum and fits multiple Gaussians to the broad and narrow emission lines.
The different methods all show the same relative variability of the broad-line
profiles, with the exception of the $\sigma$ measurements from PyQSOFit: for
details, see Appendix A. For clarity, we use the nonparametric measurements
for our analysis throughout the text.

##  4 Results and Discussion

###  4.1 Broad Emission-Line Profile Variability

Figures 6, 7, and 8 show the variability of the MgII, H$\beta$, and H$\alpha$
broad emission-line profiles as quantified by our non-parametric measures of
flux, line width (FWHM and $\sigma$), and line center. In general, in each of
the figures, the broad emission lines get brighter and narrower, fainter and
broader, and then brighter and narrower again, with an anti-correlation
between the flux and line width. This phenomenon is known as ‘line breathing’
and we discuss it in more detail in Section 4.2. We note that the flux and
line-width (FWHM and $\sigma$) of MgII appear to vary as a lagged version of
the Balmer lines, consistent with the observed lags in each emission line of
this quasar presented in Grier et al. (2017b) and Homayouni et al. (2020).
This suggests that the MgII-emitting region is further away from the central
engine, compared with the Balmer-emitting regions.

All three broad emission lines have line centers that are much redder than the
systemic redshift, as determined from the narrow emission lines. The broad-
line centers all have a similar general variability pattern of starting red,
shifting bluer over a few years, and then getting redder near the end of our
monitoring period. The broad H$\beta$ line is redder than the other broad
lines and has the largest radial velocity shifts (maximum change in
$\Delta{v}$ of $\sim$800 km s-1, compared to $\sim$400 km s-1 for H$\alpha$
and MgII). All three lines become bluest (but still redder than the systemic
narrow-line redshift) in $\sim$2017 and become red again by $\sim$2020, with
MgII returning to its initial red center while H$\alpha$ and H$\beta$ do not
become as red as when they started.

Unlike the light curves of the emission-line flux, the comparison between the
radial velocity shifts in the Balmer lines (H$\beta$ and H$\alpha$) and MgII
is inconsistent with the measured lags reported by Grier et al. (2017b) and
Homayouni et al. (2020). Specifically, the MgII radial velocity shifts do not
appear to mirror the Balmer-line shifts after a lag of $\sim$100 days.
Instead, the MgII radial velocity shifts appear to be a smoother version of
the Balmer-line shifts with no apparent lag between them.

The radial velocity shifts of RM160 are qualitatively similar to, but more
extreme than, what is observed for the lower luminosity Seyfert 1 AGN
discussed by Barth et al. (2015). For example, the largest shifts reported by
Barth et al. (2015) are for NGC 4593 with radial velocity shifts of 266 $\pm$
11 km s-1 for the broad H$\beta$ emission line. RM160 has a much more dramatic
H$\beta$ velocity shift of $\sim$800 km s-1 that occurs over $\sim$4 years.
Sergeev et al. (2007) studied NGC 5548 and found radial velocity shifts of
$\sim$1,000 km s-1 over a 30 year period.

###  4.2 Line Breathing

Line breathing is an anti-correlation between the flux and the width of a
broad emission line over time that has long been predicted by photoionization
modeling (Korista & Goad 2004). Gas in a Keplerian orbit, as is the likely
case for the BLR, has higher velocity at small radii and lower velocity at
large radii. For BLR gas ionized by a central continuum, an increase in the
continuum flux will over-ionize the gas nearest to the black hole and will
thus increase the emissivity-weighted BLR radius. The increase in the optimal
emitting radius will result in a decrease in the line width, since gas
orbiting at a larger radius has a lower orbital velocity. Observations have
shown that BLR line breathing can occur on timescales of days to weeks (Barth
et al. 2015). Line breathing is described by the following relation:

$\Delta\,\log\,W=\beta\,\Delta\,\log\,L$ (2)

where $W$ is the line-width, $L$ is the luminosity and $\beta$ is some
constant of proportionality. The broad-line luminosity is generally used for
$L$ in Equation 2 as a (lagged) representation of the continuum luminosity
that is responsible for driving the change in optimal emitting radius.
Assuming a typical radius-luminosity relationship of $R_{\rm BLR}\propto
L^{0.5}$ (Bentz et al. 2013) and virial orbits ($v\propto R_{\rm
BLR}^{-0.5}$), the expected relation between the changes in line-width and
luminosity is
$\Delta\,\textrm{log}\,W=\textrm{-0.25}\,\Delta\,\textrm{log}\,L$.

Reverberation mapping studies over time have shown that line breathing is
observed for H$\beta$ (Cackett & Horne 2006; Park et al. 2012; Barth et al.
2015; Wang et al. 2020), typically consistent with the expected constant of
proportionality of $\beta=-0.25$. On the other hand, observations of MgII have
found weak or no anti-correlations between line width and flux (Dexter et al.
2019; Homan et al. 2020; Yang et al. 2020; Wang et al. 2020). In CIV, some
observations have found a “reverse breathing” effect whereby the line width of
CIV increases with increasing flux (Wilhite et al. 2006; Wang et al. 2020).
Furthermore, Wang et al. (2020) found that H$\alpha$ shows much less breathing
than H$\beta$ on average. The differences in line breathing patterns between
broad emission lines indicate that there could be a difference in the
distribution of the gas around a particular “optimal emitting region” and also
can indicate that the structure of the BLR is not uniform at all radii.
Studying line breathing behavior for different emission lines can give clues
about the multi-scale structure of the BLR.

Figure 9: Relationships between line width and line flux for RM160. The green
points are the data over our 9 year monitoring period. The maroon line is the
best-fit linear relation and the gray shading indicates the distribution of
best-fit lines from the Markov Chains used for the fitting. The top panels
show the line FWHM vs. flux for MgII, H$\beta$, and H$\alpha$ (left to right).
The bottom panels show the line width computed from the second moment vs. flux
for MgII, H$\beta$, and H$\alpha$ (left to right). Line widths (FWHM and
$\sigma$) are in units of km s-1 and line fluxes are in units of 10-17 erg s-1
cm-2.

We analyze line breathing for the three major broad lines accessible in RM160:
MgII, H$\beta$, and H$\alpha$ over the entire 9 year monitoring period studied
here (2013 - 2022). Figure 9 shows the relationship between broad-line width
(FWHM and $\sigma$) and broad-line flux for MgII, H$\beta$, and H$\alpha$. We
fit the lines with linear models, motivated by Equation 2, using the Bayesian
linear regression package linmix (Kelly 2007). The $\beta$ slopes determined
from linmix are shown in Fig. 9 on the bottom left. All three emission lines
exhibit an anti-correlation between broad-line width and broad-line flux, with
a slope that is steeper for FWHM and shallower for $\sigma$. The H$\beta$ line
has the steepest anti-correlation and MgII has the shallowest anti-
correlation, in general agreement with previous work (Wang et al. 2020).

The line breathing properties of RM160 were previously measured by Wang et al.
(2020) using the continuum flux and broad-line widths from 29 epochs of
spectra observed over 2014-2017 (choosing only epochs that are 2$\sigma$ above
the mean S/N of each season from the 64 total epochs in 2014-2017). Our study
measures the line breathing properties over a longer period with 127 epochs of
spectra over 2013-2022 using the respective broad-line flux, which is a lagged
representation of the continuum flux, and broad-line widths. Table 1 compares
the slopes of the relationship between broad-line width (both FWHM and
$\sigma$) and broad-line flux from the present study and from Wang et al.
(2020). We find almost identical results for $\beta_{\sigma}$ (slope for
$\sigma$ line-width) for MgII and H$\alpha$ and $\beta_{\sigma}$ is fairly
close for H$\beta$. However, we see large disparities between $\beta_{\rm
FWHM}$ for all three broad emission lines, likely due to the differences in
flux (continuum vs. broad-line) used in the analyses.

Figure 10: Ratio of FWHM and $\sigma$ vs. MJD as quantified by our non-
parametric measurements of FWHM and $\sigma$ for MgII (top row), H$\beta$
(middle row), and H$\alpha$ (bottom row). The dashed red line indicates the
median $\frac{\rm FWHM}{\sigma}$ value for the 2021 monitoring period and the
colorbar represents the integrated flux for each emission line in units of
$10^{-17}$ erg s-1 cm-2. The ratio of FWHM and $\sigma$ indicates the
peakiness of the distribution. A $\frac{\rm FWHM}{\sigma}$ = 2.355 would
indicate a Gaussian distribution, while a value higher than 2.355 would
indicate a boxier line and a value lower than 2.355 would indicate a cuspier
line. We see that the peakiness of MgII is relatively stable over our 9 year
monitoring period, while H$\beta$ and H$\alpha$ become more peakier with time
and as they get brighter.

In Figure 9, there are distinct differences in the line-breathing slopes
measured from line-width $\sigma$ and FWHM. We investigate the differences in
FWHM and $\sigma$ further in Figure 10, which shows the ratio of FWHM to
$\sigma$ as a function of time with a colorbar indicating the relative
emission-line flux at each epoch. The FWHM/$\sigma$ ratio measures the
boxiness/peakiness of an emission line, with
$\frac{\mathrm{FWHM}}{\sigma}=2.355$ for a Gaussian, a higher ratio
($\frac{\mathrm{FWHM}}{\sigma}>2.355$) for a boxy line (i.e., low kurtosis),
and a lower ratio ($\frac{\mathrm{FWHM}}{\sigma}<2.355$) for a peaky line
(i.e., high kurtosis). The MgII line in RM160 becomes slightly peakier with
time (and with increasing flux). The Balmer lines change much more
dramatically, going from Gaussian or boxy profiles at a low-flux at the
beginning of the monitoring to much peakier profiles at late times and high
fluxes.

The change from boxier to peakier broad-line profiles explains the difference
in line-breathing slopes measured for FWHM and $\sigma$. The changes in
received continuum (as measured in the varying emission-line flux) do not
drive monolithic changes in the observed BLR orbits of RM160, and instead
result in changes to the shape of the velocity distribution. This likely
indicates asymmetry in the radial distribution of the broad-line gas where
more contribution from distant, low velocity gas increases the line core and
makes it appear more peakier in the luminous phase. We return to a discussion
of the complex kinematics of the BLR of RM160 in Section 4.4.

| MgII | H$\beta$ | H$\alpha$
---|---|---|---
Study | $\beta_{\rm FWHM}$ | $\beta_{\sigma}$ | $\beta_{\rm FWHM}$ | $\beta_{\sigma}$ | $\beta_{\rm FWHM}$ | $\beta_{\sigma}$
Wang et al. 2020 | -0.06 $\pm$ 0.01 | -0.08 $\pm$ 0.01 | -0.08 $\pm$ 0.01 | -0.08 $\pm$ 0.01 | -0.09 $\pm$ 0.01 | -0.16 $\pm$ 0.01
This work | -0.19 $\pm$ 0.01 | -0.11 $\pm$ 0.01 | -0.47 $\pm$ 0.02 | -0.18 $\pm$ 0.01 | -0.32 $\pm$ 0.01 | -0.17 $\pm$ 0.01
Table 1: Table comparing the best-fit slopes of the relationship between
broad-line width (both FWHM and $\sigma$) and broad-line flux. We note that
the results match very well for the MgII and H$\alpha$ $\sigma$ slopes, but
the FWHM slopes for MgII, H$\beta$, and H$\alpha$ have large disparities. The
differences between slopes measured by Wang et al. (2020) and this study are
likely due to line-profile changes between 2014-2017 and 2018-2022,
differences in fitting methods, and the difference in flux (continuum vs.
broad-line) used.

###  4.3 The Binary Supermassive Black Hole Hypothesis

A periodic Doppler shift of the broad emission lines could be indicative of a
SMBH binary (Eracleous et al. 2012). However, in the binary scenario,
$\Delta{v}$ should be consistent through the different broad emission lines as
the systemic velocities should all trace the center of gravity of the active
SMBH. This contradicts the observations of RM160 (see: Figure 11), where
$\Delta{v}$ is different for the different broad lines. The differences in
$\Delta{v}$ for different lines provide evidence to disfavor a SMBH binary as
the cause of the radial velocity variations of RM160. The simplest binary
scenario (i.e., a scenario in which only one SMBH is active), should also have
only bulk radial velocity shifts of the broad emission lines with no
variability in their shapes (Fig. 1 in Guo et al. 2019), contradicting the
${\rm FWHM}/\sigma$ variations observed in our observations (Figure 10). A
binary scenario with velocity shifts driven entirely by periodic Doppler
motion would also have a $\Delta{v}$ curve that is symmetric about zero,
opposite in shape (concave down) to our observed (concave up) $\Delta{v}$
curves.

We provide further evidence against the binary hypothesis by considering the
relationship between binary separation and orbital period (Kepler’s third
law):

$a^{3}=\frac{GM}{4\pi^{2}}P^{2}$ (3)

where $a$ is the binary separation, $G$ is the gravitational constant, $M$ is
the total binary mass, and $P$ is the period of the binary.

There are two possibilities for a binary SMBH system that has a BLR undergoing
radial velocity shifts: a close-pair binary with a circumbinary BLR ($a\ll
R_{\rm BLR}$ and $P\ll P_{\rm BLR}$) or a wide-separation binary in which each
SMBH has its own BLR ($a\gg R_{\rm BLR}$ and $P\gg P_{\rm BLR}$). A binary
SMBH with a separation similar to the BLR ($a\approx R_{\rm BLR}$) would
disrupt the orbiting gas and such a system would not have observable broad
emission-lines (unlike RM160).

In Figure 11, we fit the radial velocity shifts of all three broad emission
lines using a sine function with the LMFIT package (Newville et al. 2014). The
best-fit sine functions have characteristic periods of $P(\hbox{{\rm Mg}\kern
1.00006pt{\sc II}})=16.5\pm 7.1$ yr, $P(\hbox{{\rm H}$\beta$})=5.8\pm 4.2$ yr
(or $P(\hbox{{\rm H}$\beta$})=14.4\pm 3.4$ yr if fit to only the 2014-2020
data), and $P(\hbox{{\rm H}$\alpha$})=22.1\pm 1.0$ yr.

Using Equation 3 and assuming that the RM mass ($\log{M_{\rm
BH}/M_{\odot}}=7.8$) for RM160 is the total binary mass (which is the case for
the binary scenario with a circumbinary BLR) and that the radial velocity
shifts correspond to the period of the binary, the best-fit periods imply
binary separations of $a(\hbox{{\rm Mg}\kern 1.00006pt{\sc II}})=14.9\pm 8.5$
light-days, $a(\hbox{{\rm H}$\beta$})=7.4\pm 6.0$ light-days (or $a(\hbox{{\rm
H}$\beta$})=13.6\pm 5.2$ light-days if fit to only the 2014-2020 data), and
$a(\hbox{{\rm H}$\alpha$})=18.1\pm 2.3$ light-days. We also use Equation 3 to
calculate BLR periods using BLR radii equal to the measured reverberation lags
from Grier et al. (2017b) for the Balmer lines and from Homayouni et al.
(2020) for the MgII line. These periods are shown on the top right of each
panel in Figure 11.

The best-fit sine periods range from 6 to 22 yr, with semi-major axes of 7 to
18 light-days for the implied binary orbit. Comparing these values to the
observed lags from Grier et al. (2017b) and Homayouni et al. (2020), we find
that the range of binary semi-major axis is similar to the observed lag for
H$\alpha$ (i.e., $a_{\rm BLR,H\alpha}=27.7^{+5.3}_{-4.7}$ light-days). In the
binary scenario this would place the black hole orbit within the BLR, which
would cause the BLR to become unstable due to the gravitational interactions
with the orbiting black hole. We therefore disfavor a binary explanation for
the radial velocity shifts in RM160.

If we instead consider the possibility of a wide separation binary system with
a BLR around one (active) black hole, then the RM mass represents a minimum
for the total binary mass and the semi-major axes from Equation 3 are also
minimum values. Even for a maximum total binary mass of $5\times
10^{10}M_{\odot}$ (the upper limit for an AGN, e.g., King 2016), and an
extreme binary mass ratio of $\sim$1000:1, the binary separation would be
$a(\rm\hbox{{\rm Mg}\kern 1.00006pt{\sc II}})=166.8\pm 5.1$ light-days, again
overlapping with the observed BLR size of $a_{\rm MgII}=106$ light-days
(Homayouni et al. 2020). In other words, the binary hypothesis is ruled out
for both a circumbinary BLR and a single-AGN BLR due to overlap between the
putative binary orbit and the observed BLR size.

Figure 11: Sine fits to the radial velocity shifts of MgII (top panel),
H$\beta$ (middle panel), and H$\alpha$ (bottom panel). The tangerine points
represent the observations and the black line is the best-fit sine function
using the fitting routine LMFIT (Newville et al. 2014). The purple line in the
middle panel is the best-fit sine function for the H$\beta$ radial velocity
shifts excluding the 2021 and 2022 data (MJD - 56393 < 2555). For each
emission line, we present the period of the best-fit sine function and the
period of the BLR at the top right of each panel. We disfavor a binary
explanation for the radial velocity variations because the binary period is
similar to the broad-line orbital period, and because the three lines have
best fit sine functions with different periods.

###  4.4 BLR Kinematics and Geometry

We begin by calculating the dynamical timescale for the BLR, $\tau_{\rm
dyn}\simeq(R_{\rm BLR}^{3}/GM_{\rm BH})^{1/2}$ = 7.7 yr, using the H$\beta$
based parameters for $R_{\rm BLR}$ and $M_{\rm BH}$. The best-fit variability
periods in Figure 11 are longer than this dynamical timescale, indicating that
the observed variability is consistent with dynamical changes in the BLR. We
present an example of a phenomenological model for these dynamical changes
below.

RM160 has large ($\sim$400-800 km s-1) radial velocity shifts in its broad
emission lines. Although the source exhibits line breathing in a fashion that
is similar to other quasars, its line shifts do not follow the same pattern as
the changes in line brightness and width. One plausible explanation of the
radial velocity shifts is a BLR with azimuthal asymmetry and a gradient of
inflow velocity in the radial direction. An inflow model is motivated by
previous velocity-resolved reverberation mapping observations (e.g., Bentz et
al. 2010a; Grier et al. 2017a; Bentz et al. 2021; U et al. 2022; Villafaña et
al. 2022) that frequently find evidence for inflowing BLRs in nearby Seyfert
AGN. Alternatively, we could be seeing an outflow due to BLR emission being
preferentially emitted back toward the ionizing source (Ferland et al. 1992),
with a gradient of outflow velocity decreasing with radial distance from the
quasar. We refer to an inflow hereafter, but acknowledge that a decelerating
outflow with the right azimuthal asymmetry might also explain the
observations.

Our model is illustrated by Figure 12, which includes an animation of the BLR
and its kinematics connected to the observed broad-line light curves and
widths.

We begin by noting that all three emission lines have large redshifts
($\sim$500-1600 km s-1) with respect to the narrow lines (i.e., the systemic
redshift) at all epochs. This likely indicates bulk inflow of the BLR gas in
the line of sight. At all epochs the H$\beta$ line is most redshifted
($\sim$900-1600 km s-1) and the MgII line is the least redshifted
($\sim$500-800 km s-1). The reverberation lags indicate that the MgII-emitting
gas (rest-frame $R_{\rm BLR}=106$ light-days; Homayouni et al. 2020) is much
further from the continuum emission than the H$\beta$-emitting gas (rest-frame
$R_{\rm BLR}=23$ days; Grier et al. 2017b). This further implies that the
inflow of the BLR gas has a radial gradient, with higher inflow velocity for
gas closer to the SMBH (like the H$\beta$ emission region) and lower inflow
for more distant gas (like the MgII emission region). The H$\alpha$ line is a
bit puzzling in this picture because it has a similar reverberation lag (rest-
frame $R_{\rm BLR}=20$ days; Grier et al. 2017b) to H$\beta$ but has a smaller
inflow velocity. In general the H$\alpha$ line is expected to be emitted from
slightly larger radii than H$\beta$ due to radial stratification and optical
depth effects (Netzer 1975; Rees et al. 1989; Korista & Goad 2004; Bentz et
al. 2010b) and we assume that this is also the case here, despite the
similarity in measured reverberation lags for the two lines.

Alternatively, the redshifted broad emission lines might be explained by
gravitational redshift (Tremaine et al. 2014). This scenario similarly
predicts that the H$\beta$ line would be more redshifted than the H$\alpha$
and MgII lines, due to H$\beta$ being emitted from gas closer to the black
hole that has larger orbital velocities. However the broad emission line
widths ($\sigma<3000$ km s-1 for all three lines) do not imply relativistic
orbits unless the BLR is observed at a nearly face-on inclination. We thus
prefer an inflow as the explanation for the redshifted broad emission lines
rather than gravitational redshift.

An inflowing BLR explains the redshifted lines, and a radial gradient inflow
explains the difference in redshift from H$\beta$ to H$\alpha$ to MgII. But it
does not explain the variability of the broad-line centers. Figure 11 shows
that the observed line-center variations, especially for H$\alpha$ and MgII,
are best-fit by a sine function with a period that is similar to the BLR
orbital period implied from the H$\beta$ lag (i.e., 22.1 years for H$\alpha$
and 16.5 years for MgII). This suggests that the line-center shifts might be
related to an azimuthal asymmetry in the broad-line emission that orbits the
central SMBH. The asymmetric BLR emission might be associated with a higher
density in the gas, a hot spot (or hot “wedge” or spiral arm), higher
responsivity of the gas on one side, or could be associated with asymmetric
illumination from the accretion disk.

The detailed radial velocity shifts can be explained by a combination of
inflowing gas onto the BLR with a radial gradient, orbiting asymmetric gas
emission, and flux-driven changes to the optimal emission region (line
breathing). At the start of our monitoring in $\sim$2014, the BLR receives low
continuum flux and the emission region is close to the SMBH, with high bulk
inflow velocity and an asymmetric gas region that is additionally on the
receding (redshifted) side of its orbit. The asymmetric gas region reaches the
approaching part of the orbit in $\sim$2017, although the line is still
redshifted due to the bulk inflow. The modest ($\sim$2$\times$) brightening of
the line emission and coordinated decrease in line-width over 2014-2017 causes
the line emission region to move slightly outward, also resulting in slightly
lower redshift of the line. After 2017 the asymmetric region of the BLR begins
receding again, and together with the line flux decrease (and line width
increase, and emission radius decrease) the line once again becomes redder. We
hypothesize that, in $\sim$2020, the orbit of the asymmetric region would have
caused the lines to reach the high redshifts they began with in 2014, if not
for the dramatic brightening of the Balmer lines observed in 2020. The
significant increase in flux results in the line emission region becoming
larger (with observed lower line widths), corresponding to lower bulk inflow
velocities. Combined with the redder emission from the orbiting azimuthal
asymmetry, the decreased bulk inflow velocity results in a lower integrated
redshift velocity for the line in 2020 as compared to 2014.

Figure 12:  Our geometric and kinematic model of the BLR in RM160 (top),
connected to the observed broad-line flux light curve and radial velocity
shifts (bottom). The viewing angle of the observer is indicated in the top
left. The radii of H$\beta$, H$\alpha$, and MgII are indicated by green blue,
and orange in both the top and bottom panels, scaled to match the measured
reverberation lags (Grier et al. 2017b; Homayouni et al. 2020) and with gray
shading to indicate the radial extent of each emission line region. The broad-
line regions “breathe,” moving inward when the quasar is fainter (e.g., in
2014 and 2019) and outward when the quasar is brighter (e.g., in 2021). The
red dot represents an azimuthal asymmetry in the BLR that orbits the black
hole during our spectroscopic monitoring. In the animation, darker circles
indicate the flux and radial velocity in the bottom panels at the same time as
the model in the top panels. In the still images, the labels a, b, c, and d
indicate representative time windows (also shaded purple) for the BLR model in
the top panels and corresponding flux and line-center in the bottom panels. A
combination of radial inflow (vectors), azimuthal asymmetry (red spot), and
line breathing (green, blue, and red lines and gray shading) can explain the
unusual broad-line variability of RM160. The online version of this figure is
an animation. The animation is 0:24 minutes long and shows the time-evolution
of our BLR model from 2013 to 2022 (top) along with the corresponding time-
evolution of the flux and line-center, respectively, in the bottom panels.

The combined effects of bulk inflow, azimuthal asymmetry, and line breathing
are strongest for H$\beta$, as the line is at a smaller $R_{\rm BLR}$, and
weakest for MgII, as the line is at a much larger $R_{\rm BLR}$. The pattern
over 2014-2020 of the MgII line-center shifts are much more symmetric than the
H$\beta$ line-center shifts. The H$\beta$ line also has much larger redshift
velocities than the other lines due to its emission region occupying smaller
BLR radii.

Velocity-resolved reverberation mapping would provide a test of our geometric
and kinematic model for the BLR in RM160 shown in Figure 12. Velocity-resolved
RM involves measuring how different segments of an emission line reverberate
in response to continuum variations (Denney et al. 2009; Bentz et al. 2009,
2010b; U et al. 2022; Li et al. 2022; Villafaña et al. 2022). If the BLR is
virialized, the lags at the center of the emission line are the longest since
they preferentially correspond to gas further from the black hole, with
shorter lags measured for the line wings. For an inflowing BLR (like our
phenomenological model), the lags of the blue wing would be longest and the
lags of the red wing would be shortest. Velocity-resolved RM could also
isolate the putative azimuthal asymmetry that orbits around the black hole in
our model. We anticipate performing velocity-resolved RM and further testing
our model for RM160 in future work.

##  5 Summary

We presented multi-epoch optical spectroscopy of RM160, a luminous quasar
which exhibits unusual broad emission-line variability in the SDSS-RM field.
This object was identified from a broad search for extreme variability in
quasar broad emission-line profiles and has been observed 127 times over the 9
year monitoring period with plans to continue observations through 2026,
within the SDSS-V project.

We find that RM160 exhibits normal line-breathing behavior consistent with
many previously studied AGN (e.g., Barth et al. 2015 and Wang et al. 2020) in
the variations of flux and line-width $\sigma$, but there is a sizable
discrepancy between line-width $\sigma$ and FWHM for the Balmer series
(H$\beta$ and H$\alpha$). We find that the shape of the line-profile (as
indicated by FWHM/$\sigma$) changes over time for the Balmer series from a
boxy line profile (FWHM/$\sigma$ > 2.355) to a cuspy profile (FWHM/$\sigma$ <
2.355) with increasing flux throughout the 9 year monitoring period. This
likely indicates asymmetry in the radial distribution of the broad-line gas,
where more contribution from distant, low-velocity gas increases the line core
and makes it appear more peakier in the luminous phase.

Dramatic radial velocity variations occur in each of the three broad emission
lines (MgII, H$\beta$, and H$\alpha$) that all follow the same qualitative
trend of starting red, shifting bluer over a few years, and then getting
redder near the end of the monitoring. The radial velocity shifts are not well
explained by a black hole binary because the best-fit period corresponds to
the inner BLR orbits, such that a putative binary would have disrupted the BLR
gas.

Our explanation for the large radial velocity shifts in the broad emission
lines of RM160 is a BLR with azimuthal asymmetry and a gradient of inflow
velocity in the radial direction coupled with flux-driven changes to the
optimal emission region (line breathing), as illustrated in Figure 12. Similar
instances of line-profile variability due to complex gas kinematics in the BLR
are likely to represent an important source of false positives in radial
velocity searches for binary black holes. The long-duration, wide-field, and
many-epoch spectroscopic monitoring of SDSS-V/BHM-RM will be excellent for
studying such systems and helping understand the various mechanisms driving
BLR dynamics.

##  6 Acknowledgements

LBF, JRT, and MCD acknowledge support from NSF grant CAREER-1945546, and with
CJG acknowledge support from NSF grant AST-2108668. JRT, CJG, and YS also
acknowledge support from NSF grant AST-2009539. MK acknowledges support by DFG
grant KR 3338/4-1. B.T. acknowledges support from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
program (grant agreement 950533) and from the Israel Science Foundation (grant
1849/19). X.L. acknowledges support from NSF grant AST-2206499. CR
acknowledges support from the Fondecyt Iniciación grant 11190831 and ANID
BASAL project FB210003. RJA acknowledges support by ANID BASAL FB 210003 and
by FONDECYT grant number 1191124. M. L. M.-A. acknowledges financial support
from Millenium Nucleus NCN$19\\_058$ (TITANs).

Funding for the Sloan Digital Sky Survey V has been provided by the Alfred P.
Sloan Foundation, the Heising-Simons Foundation, the National Science
Foundation, and the Participating Institutions. SDSS acknowledges support and
resources from the Center for High-Performance Computing at the University of
Utah. The SDSS web site is www.sdss5.org.

SDSS is managed by the Astrophysical Research Consortium for the Participating
Institutions of the SDSS Collaboration, including the Carnegie Institution for
Science, Chilean National Time Allocation Committee (CNTAC) ratified
researchers, the Gotham Participation Group, Harvard University, Heidelberg
University, The Johns Hopkins University, L’Ecole polytechnique fédérale de
Lausanne (EPFL), Leibniz-Institut für Astrophysik Potsdam (AIP), Max-Planck-
Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für
Extraterrestrische Physik (MPE), Nanjing University, National Astronomical
Observatories of China (NAOC), New Mexico State University, The Ohio State
University, Pennsylvania State University, Smithsonian Astrophysical
Observatory, Space Telescope Science Institute (STScI), the Stellar
Astrophysics Participation Group, Universidad Nacional Autónoma de México,
University of Arizona, University of Colorado Boulder, University of Illinois
at Urbana-Champaign, University of Toronto, University of Utah, University of
Virginia, Yale University, and Yunnan University.

, linmix (Kelly 2007), PyQSOFit (Guo et al. 2018; Shen et al. 2019), specutils
(Earl et al. 2022)

##  Appendix A Alternative Fitting Procedures

Our paper uses a non-parametric approach to measuring the broad emission lines
of RM160. In this Appendix, we fit the emission-line properties with
alternative methods in order to ensure that our conclusions about the quasar’s
dramatic changes in line flux, width and radial velocity are robust to the
choice of fitting method. Section A.1 describes the results of single-Gaussian
fits for the broad emission lines. Section A.2 instead uses PyQSOFit (Guo et
al. 2018; Shen et al. 2019), a spectral fitting code that fits the quasar
continuum, iron psuedo-continuum, and broad and narrow emission lines with
multiple Gaussians. Section A.3 compares the three fitting methods and
concludes that our conclusions about the quasar’s dramatic broad-line
breathing and radial velocity shifts are robust to different fitting methods.

###  A.1 Single Gaussian Fits

The first alternative fitting method was to fit the broad emission lines with
a single Gaussian. Similarly to our non-parametric measurements, we begin with
the same continuum-subtracted spectra describe in Section 3.1. We then used
the Gaussian1D method from the AstroPy package (Astropy Collaboration et al.
2013, 2018) in order to fit a single Gaussian to each of the narrow and broad
emission lines.

Figure 13: Variability of the MgII broad emission-line profile, as quantified
by our three fitting methods for the line flux (top row), FWHM (2nd row),
line-width $\sigma$ (3rd row), and the line-center velocity shifts $\Delta{v}$
(bottom row). The units for FWHM, line-width $\sigma$, and line-center
velocity shift $\Delta{v}$ are all in km s-1. The measured flux and FWHM are
consistent across all methods. There are differences in the measured $\sigma$
and $\Delta{v}$ due to differences in fitting the wings of the line profile,
but all three methods find similar qualitative trends of line breathing and
shifts in the line center.  Figure 14: Variability of the H$\beta$ broad
emission-line profile, as quantified by the three fitting methods for the line
flux (top row), FWHM (2nd row), line-width $\sigma$ (3rd row), and the line-
center velocity shifts $\Delta{v}$ (bottom row). The units for FWHM, line-
width $\sigma$, and line-center velocity shift $\Delta{v}$ are all in km s-1.
The H$\beta$ line has a complex profile (see Figure 2) and the three methods
result in different line fluxes, $\sigma$ widths, and line centers due to
including different amounts of the line wings in the fit. Despite these
differences, all three methods find similar qualitative trends for the
variability of the line.  Figure 15: Variability of the H$\alpha$ broad
emission-line profile, as quantified by our three fitting methods for the line
flux (top row), FWHM (2nd row), line-width $\sigma$ (3rd row), and the line-
center velocity shifts $\Delta{v}$ (bottom row). The units for FWHM, line-
width $\sigma$, and line-center velocity shift $\Delta{v}$ are all in km s-1.
As for the fit to MgII, the measured flux and FWHM is consistent across all
methods. The measured flux and FWHM are consistent across all methods.There
are differences in the measured $\sigma$ and $\Delta{v}$ due to differences in
fitting the wings of the line profile, but all three methods find similar
qualitative trends of line breathing and shifts in the line center.

For each individual emission-line region, we followed the initial procedure in
Section 3.2, in that we tied the line-centers and line-widths of the narrow
emission lines to $\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007
and constrained the amplitude of $\hbox{{\rm[O}\kern 1.00006pt{\sc
iii}{\rm]}}\lambda$4959 to be 1/3 that of $\hbox{{\rm[O}\kern 1.00006pt{\sc
iii}{\rm]}}\lambda$5007\. In what follows, we briefly describe the steps to
obtain the broad emission-line component in each emission-line region.

For MgII, we fit the individual spectra with a single Gaussian as there are no
narrow emission lines within the MgII emission-line region.

Figure 16: Line-width vs. flux plots for RM160 for all fitting methods. The
purple points are from the non-parametric fitting method, the pink points are
from QSOfit, and the green points are from the single Gaussian fitting method.
The colors of the lines of best-fit correspond to the colors of the data
points. The top panels show the FWHM vs. flux for MgII, H$\beta$, and
H$\alpha$ (left to right). The bottom panels show the line-width $\sigma$ vs.
flux for MgII, H$\beta$, and H$\alpha$ (left to right). The individual slopes
are shown at the bottom left of each panel where $\beta_{QF}$ is the slope for
the QSOFit fitting method, $\beta_{SG}$ is the slope for the single Gaussian
fiting method, and $\beta_{NP}$ is the slope for our non-parametric
measurements. The best-fit slopes for FWHM are broadly consistent, while we
see large discrepancies between the slopes for line-width $\sigma$. This seems
consistent with what we see for the FWHM values and line-width $\sigma$ values
in Figures 13, 14, and 15 where we see a general consensus for FWHM values
across all schemes, and differences between values in line-width $\sigma$.

For H$\alpha$, we simultaneously fit the narrow-lines ($\hbox{{\rm[N}\kern
1.00006pt{\sc ii}{\rm]}}\lambda$6548, narrow H$\alpha$, $\hbox{{\rm[N}\kern
1.00006pt{\sc ii}{\rm]}}\lambda$6584, $\hbox{{\rm[S}\kern 1.00006pt{\sc
ii}{\rm]}}\lambda$6718, and $\hbox{{\rm[S}\kern 1.00006pt{\sc
ii}{\rm]}}\lambda$6732) and the broad H$\alpha$ emission line with no
constraints on the amplitude of any emission line. This simultaneous fitting
of all of the narrow lines and the broad H$\alpha$ component produced reliable
fits.

For H$\beta$, degeneracies between the lines made it difficult to
simultaneously fit the narrow lines (narrow H$\beta$, $\hbox{{\rm[O}\kern
1.00006pt{\sc iii}{\rm]}}\lambda$4959, and $\hbox{{\rm[O}\kern 1.00006pt{\sc
iii}{\rm]}}\lambda$5007) while also fitting the broad H$\beta$ component
reliably well. Instead, we constrained the amplitude of the narrow H$\beta$
line to that which produced the smoothest residual in the median spectrum fit.
We then applied that narrow H$\beta$ amplitude to the fits for all epochs
along with the $\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$4959 and
$\hbox{{\rm[O}\kern 1.00006pt{\sc iii}{\rm]}}\lambda$5007 narrow lines and
produced a narrow line subtracted spectrum. From that, we fit a single
Gaussian to the continuum and narrow-line subtracted spectra in all epochs to
fit the broad H$\beta$ component.

For all emission lines in this fitting scheme, we calculated the FWHM using
the relationship between FWHM and line-width ($\sigma$) for a Gaussian,
namely, FWHM = 2.355$\sigma$.

We employed bootstrap resampling to estimate the uncertainties in line flux,
line center, line width ($\sigma$), and FWHM similar to our uncertainty
estimation approach in Section 3.4. We resampled each individual spectra 200
times by sampling the fluxes within their Gaussian uncertainties in the
respective epoch. We adopted the standard deviations of the 200 re-sampled
spectra as our uncertainties for each epoch.

###  A.2 PyQSOFit

The second alternative fitting method was to fit the entire spectrum with the
robust quasar-fitting algorithm PyQSOFit. The considerations for the continuum
fit include a power law, an optical and UV FeII template, and a fifth-order
polynomial to account for possible dust reddening. PyQSOFit allows the user to
select the emission-line and continuum components that are included in the
overall fit as well as selecting the range of fitting for the emission lines.
We show our fitting parameters in Table 2.

We left most of the default parameters unchanged, except for turning off the
dereddening and host decomposition since our fits within the relatively narrow
wavelength regions around the broad lines will be unaffected by reddening and
contributions from the host galaxy. We decided to fit the narrow emission
lines with a single Gaussian, consistent with the fitting methods in Section
3.4 and Appendix A.1. We allowed 2 Gaussians to fit the broad lines in order
to better capture the complex shapes of the line profiles.

line | type | n_gauss | $\lambda_{\rm rest}$
---|---|---|---
H$\alpha$ | broad | 2.0 | 6564.61
H$\alpha$ | narrow | 1.0 | 6564.61
[Nii]$\lambda$6549 | narrow | 1.0 | 6549.85
[Nii]$\lambda$6585 | narrow | 1.0 | 6585.28
[Sii]$\lambda$6718 | narrow | 1.0 | 6718.29
[Sii]$\lambda$6732 | narrow | 1.0 | 6732.67
H$\beta$ | broad | 2.0 | 4862.68
H$\beta$ | narrow | 1.0 | 4862.68
[Oiii]$\lambda$4959 | narrow | 1.0 | 4960.30
[Oiii]$\lambda$5007 | narrow | 1.0 | 5008.24
MgII | broad | 2.0 | 2798.75
MgII | narrow | 1.0 | 2798.75
Table 2: Table of our parameters used for PyQSOFit. The first column is the
emission-line name, the second column is the type of emission line (broad or
narrow), the third column is the number of Gaussians used in the fit, and the
fourth column is the rest-frame central wavelength.

PyQSOFit was used to obtain the flux, line width ($\sigma$), FWHM, and line
center of the MgII, H$\beta$, and H$\alpha$ emission lines for each epoch. WE
use the uncertainties reported by PyQSOFit for each of the fitted quantities.

Note that our PyQSOFit analysis is not exactly analogous to Wang et al. (2020)
because they used 3 Gaussians to fit each broad emission line. Furthermore,
they rejected epochs that are 2$\sigma$ below the mean S/N for each season.
Wang et al. (2020) also use a window of [-2.5 $\times$ MAD, 2.5 $\times$ MAD]
to compute $\sigma_{line}$ so as to eliminate the effects of noise and
blending in the line wings.

###  A.3 Comparison of the Three Fitting Methods

Figures 13, 14, and 15 show the variability of the MgII, H$\beta$, and
H$\alpha$ broad emission-line profiles as quantified by the three different
fitting methods.

The three methods result in very similar flux light curves for the H$\alpha$
and MgIIlines. The different methods also find very similar FWHM measurements
for the three lines, implying that the FWHM measurement is robust to the
details of the emission-line models (Dalla Bontà et al. 2020). On the other
hand, the three methods have large differences in the $\sigma$ measurement.
The $\sigma$ line width is highly sensitive to the wings of the line and so is
susceptible to small differences in the fitted line profile. Among the three
lines, H$\beta$ has the largest differences in measured quantities between the
three lines, owing to its complex and asymmetric line profile (shown in Figure
2) that is highly sensitive to small differences in fitting.

Despite differences in the details of the line fits, the three methods result
in similar qualitative trends for all three emission lines. That is, the
differences between methods are generally systematic offsets: e.g., QSOFit
measures broader $\sigma$ and lower $\Delta v$ for all lines, as well as
brighter H$\beta$, than the other methods. The similar overall qualitative
trends mean that our general conclusions about the relative changes in the
line profiles are not dependent on the details of the fitting method.

Figure 16 shows the line breathing behavior measured by the three fitting
methods. We fit the relationships in each panel using the Bayesian linear
regression package linmix (Kelly 2007). The slopes are shown at the bottom
left of each panel where $\beta_{QF}$ corresponds to the slope for the QSOFit
fitting method, $\beta_{SG}$ corresponds to the slope for the single Gaussian
fitting method, and $\beta_{NP}$ corresponds to the slope for our non-
parametric measurements. The slopes for FWHM are broadly similar for all three
fits to each emission line, while the slopes for line-width $\sigma$ show
large differences. As noted above in the discussion of Figures 13, 14, and 15,
$\sigma$ line-width measurements are more sensitive to the details of the line
profile and so are more dependent on the differences in best-fit models from
each method. From this investigation we note that characterizing quasar line
breathing with $\sigma$ is likely to depend significantly on the details of
the fitting method especially for emission lines that change shape like RM160
(see Figure 10). FWHM line-width measurements, on the other hand, are less
sensitive to details of the model fitting.

## References

  * Ackermann et al. (2015) Ackermann, M., Ajello, M., Albert, A., et al. 2015, The Astrophysical Journal, 813, L41.  https://doi.org/10.1088/2041-8205/813/2/l41
  * Almeida et al. (2023) Almeida, A., Anderson, S. F., Argudo-Fernández, M., et al. 2023, arXiv e-prints, arXiv:2301.07688
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Ayubinia et al. (2022) Ayubinia, A., Xue, Y., Woo, J.-H., et al. 2022, Universe, 8, 559
  * Barth et al. (2011a) Barth, A. J., Nguyen, M. L., Malkan, M. A., et al. 2011a, ApJ, 732, 121
  * Barth et al. (2011b) Barth, A. J., Pancoast, A., Thorman, S. J., et al. 2011b, ApJ, 743, L4
  * Barth et al. (2015) Barth, A. J., Bennert, V. N., Canalizo, G., et al. 2015, ApJS, 217, 26
  * Bentz et al. (2021) Bentz, M. C., Williams, P. R., Street, R., et al. 2021, ApJ, 920, 112
  * Bentz et al. (2009) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2009, ApJ, 705, 199
  * Bentz et al. (2010a) Bentz, M. C., Horne, K., Barth, A. J., et al. 2010a, ApJ, 720, L46
  * Bentz et al. (2010b) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2010b, ApJ, 716, 993
  * Bentz et al. (2013) Bentz, M. C., Denney, K. D., Grier, C. J., et al. 2013, The Astrophysical Journal, 767, 149.  https://doi.org/10.1088/0004-637x/767/2/149
  * Blandford & McKee (1982) Blandford, R. D., & McKee, C. F. 1982, ApJ, 255, 419
  * Blanton et al. (2017) Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28
  * Bolton et al. (2012) Bolton, A. S., Schlegel, D. J., Aubourg, É., et al. 2012, AJ, 144, 144
  * Brotherton et al. (1994) Brotherton, M. S., Wills, B. J., Steidel, C. C., & Sargent, W. L. W. 1994, ApJ, 423, 131
  * Cackett et al. (2021) Cackett, E. M., Bentz, M. C., & Kara, E. 2021, iScience, 24, 102557
  * Cackett & Horne (2006) Cackett, E. M., & Horne, K. 2006, Monthly Notices of the Royal Astronomical Society, 365, 1180.  https://doi.org/10.1111/j.1365-2966.2005.09795.x
  * Charisi et al. (2016) Charisi, M., Bartos, I., Haiman, Z., et al. 2016, MNRAS, 463, 2145
  * Chen et al. (2020) Chen, Y.-C., Liu, X., Liao, W.-T., et al. 2020, MNRAS, 499, 2245
  * Clavel et al. (1991) Clavel, J., Reichert, G. A., Alloin, D., et al. 1991, ApJ, 366, 64
  * Dalla Bontà et al. (2020) Dalla Bontà, E., Peterson, B. M., Bentz, M. C., et al. 2020, ApJ, 903, 112
  * Denney et al. (2009) Denney, K. D., Watson, L. C., Peterson, B. M., et al. 2009, ApJ, 702, 1353
  * Dexter et al. (2019) Dexter, J., Xin, S., Shen, Y., et al. 2019, The Astrophysical Journal, 885, 44\.  https://doi.org/10.3847/1538-4357/ab4354
  * Earl et al. (2022) Earl, N., Tollerud, E., Jones, C., et al. 2022, astropy/specutils: V1.7.0, vv1.7.0, Zenodo, doi:10.5281/zenodo.6207491.  https://doi.org/10.5281/zenodo.6207491
  * Eisenstein et al. (2011) Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, AJ, 142, 72
  * Eracleous et al. (2012) Eracleous, M., Boroson, T. A., Halpern, J. P., & Liu, J. 2012, ApJS, 201, 23
  * Ferland et al. (1992) Ferland, G. J., Peterson, B. M., Horne, K., Welsh, W. F., & Nahar, S. N. 1992, ApJ, 387, 95
  * Foltz et al. (1981) Foltz, C. B., Peterson, B. M., Capriotti, E. R., et al. 1981, ApJ, 250, 508
  * Gaskell (1983) Gaskell, C. M. 1983, in Liege International Astrophysical Colloquia, Vol. 24, Liege International Astrophysical Colloquia, ed. J.-P. Swings, 473–477
  * Gezari (2021) Gezari, S. 2021, ARA&A, 59, arXiv:2104.14580
  * Goad & Korista (2014) Goad, M. R., & Korista, K. T. 2014, Monthly Notices of the Royal Astronomical Society, 444, 43.  https://doi.org/10.1093/mnras/stu1456
  * Graham et al. (2015) Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015, Nature, 518, 74
  * Gravity Collaboration et al. (2020) Gravity Collaboration, Amorim, A., Bauböck, M., et al. 2020, A&A, 643, A154
  * Grier et al. (2017a) Grier, C. J., Pancoast, A., Barth, A. J., et al. 2017a, ApJ, 849, 146
  * Grier et al. (2013) Grier, C. J., Peterson, B. M., Horne, K., et al. 2013, The Astrophysical Journal, 764, 47.  https://doi.org/10.1088/0004-637x/764/1/47
  * Grier et al. (2013) Grier, C. J., Martini, P., Watson, L. C., et al. 2013, ApJ, 773, 90
  * Grier et al. (2017b) Grier, C. J., Trump, J. R., Shen, Y., et al. 2017b, ApJ, 851, 21
  * Gunn et al. (2006) Gunn, J. E., Siegmund, W. A., Mannery, E. J., et al. 2006, AJ, 131, 2332
  * Guo et al. (2019) Guo, H., Liu, X., Shen, Y., et al. 2019, MNRAS, 482, 3288
  * Guo et al. (2018) Guo, H., Shen, Y., & Wang, S. 2018, PyQSOFit: Python code to fit the spectrum of quasars, Astrophysics Source Code Library, , , ascl:1809.008
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357\.  https://doi.org/10.1038/s41586-020-2649-2
  * Herrmann et al. (2007) Herrmann, F., Hinder, I., Shoemaker, D., Laguna, P., & Matzner, R. A. 2007, ApJ, 661, 430
  * Homan et al. (2020) Homan, D., MacLeod, C. L., Lawrence, A., Ross, N. P., & Bruce, A. 2020, MNRAS, 496, 309
  * Homayouni et al. (2020) Homayouni, Y., Trump, J. R., Grier, C. J., et al. 2020, ApJ, 901, 55
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Kelly (2007) Kelly, B. C. 2007, ApJ, 665, 1489
  * Kelly et al. (2009) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895
  * King (2016) King, A. 2016, MNRAS, 456, L109
  * Kollmeier et al. (2017) Kollmeier, J. A., Zasowski, G., Rix, H.-W., et al. 2017, arXiv e-prints, arXiv:1711.03234
  * Korista & Goad (2004) Korista, K. T., & Goad, M. R. 2004, The Astrophysical Journal, 606, 749.  https://doi.org/10.1086/383193
  * Kovačević-Dojčinović et al. (2022) Kovačević-Dojčinović, J., Dojčinović, I., Lakićević, M., & Popović, L. Č. 2022, A&A, 659, A130
  * Li et al. (2023) Li, J. I.-H., Shen, Y., Ho, L. C., et al. 2023, arXiv e-prints, arXiv:2301.04177
  * Li et al. (2022) Li, S.-S., Feng, H.-C., Liu, H. T., et al. 2022, ApJ, 936, 75
  * Li et al. (2016) Li, Y.-R., Wang, J.-M., Ho, L. C., et al. 2016, The Astrophysical Journal, 822, 4.  https://doi.org/10.3847/0004-637x/822/1/4
  * Li et al. (2019) Li, Y.-R., Wang, J.-M., Zhang, Z.-X., et al. 2019, ApJS, 241, 33
  * Liao et al. (2021) Liao, W.-T., Chen, Y.-C., Liu, X., et al. 2021, MNRAS, 500, 4025
  * Liu et al. (2014) Liu, X., Shen, Y., Bian, F., Loeb, A., & Tremaine, S. 2014, ApJ, 789, 140
  * Loeb (2010) Loeb, A. 2010, Phys. Rev. D, 81, 047503
  * Lynden-Bell (1969) Lynden-Bell, D. 1969, Nature, 223, 690
  * MacLeod et al. (2010) MacLeod, C. L., Ivezić, Ž., Kochanek, C. S., et al. 2010, The Astrophysical Journal, 721, 1014.  https://doi.org/10.1088/0004-637x/721/2/1014
  * Molina et al. (2022) Molina, J., Ho, L. C., Wang, R., et al. 2022, ApJ, 935, 72
  * Netzer (1975) Netzer, H. 1975, MNRAS, 171, 395
  * Newville et al. (2014) Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, v0.8.0, Zenodo, doi:10.5281/zenodo.11813.  https://doi.org/10.5281/zenodo.11813
  * Pancoast et al. (2014) Pancoast, A., Brewer, B. J., & Treu, T. 2014, MNRAS, 445, 3055
  * Park et al. (2012) Park, D., Woo, J.-H., Treu, T., et al. 2012, The Astrophysical Journal, 747, 30\.  https://doi.org/10.1088/0004-637x/747/1/30
  * Peterson (1993) Peterson, B. M. 1993, Publications of the Astronomical Society of the Pacific, 105, 247.  https://doi.org/10.1086/133140
  * Peterson et al. (1982) Peterson, B. M., Foltz, C. B., Byard, P. L., & Wagner, R. M. 1982, ApJS, 49, 469
  * Peterson et al. (2013) Peterson, B. M., Denney, K. D., De Rosa, G., et al. 2013, ApJ, 779, 109
  * Rakshit & Woo (2018) Rakshit, S., & Woo, J.-H. 2018, ApJ, 865, 5
  * Rees et al. (1989) Rees, M. J., Netzer, H., & Ferland, G. J. 1989, ApJ, 347, 640
  * Reichert et al. (1994) Reichert, G. A., Rodriguez-Pascual, P. M., Alloin, D., et al. 1994, ApJ, 425, 582
  * Rojas et al. (2020) Rojas, A. F., Sani, E., Gavignaud, I., et al. 2020, MNRAS, 491, 5867
  * Runnoe et al. (2017) Runnoe, J. C., Eracleous, M., Pennell, A., et al. 2017, MNRAS, 468, 1683
  * Sandrinelli et al. (2018) Sandrinelli, A., Covino, S., Treves, A., et al. 2018, A&A, 615, A118
  * Sergeev et al. (2007) Sergeev, S. G., Doroshenko, V. T., Dzyuba, S. A., et al. 2007, ApJ, 668, 708
  * Severgnini et al. (2018) Severgnini, P., Cicone, C., Della Ceca, R., et al. 2018, MNRAS, 479, 3804
  * Seyfert (1943) Seyfert, C. K. 1943, ApJ, 97, 28
  * Shapovalova, A. I. et al. (2010) Shapovalova, A. I., Popovi´c, L. C., Burenkov, A. N., et al. 2010, A&A, 509, A106.  https://doi.org/10.1051/0004-6361/200912311
  * Shen et al. (2013) Shen, Y., Liu, X., Loeb, A., & Tremaine, S. 2013, ApJ, 775, 49
  * Shen et al. (2015) Shen, Y., Brandt, W. N., Dawson, K. S., et al. 2015, ApJS, 216, 4
  * Shen et al. (2019) Shen, Y., Hall, P. B., Horne, K., et al. 2019, ApJS, 241, 34
  * Smee et al. (2013) Smee, S. A., Gunn, J. E., Uomoto, A., et al. 2013, AJ, 146, 32
  * Storchi-Bergmann (2010) Storchi-Bergmann, T. 2010, in Co-Evolution of Central Black Holes and Galaxies, ed. B. M. Peterson, R. S. Somerville, & T. Storchi-Bergmann, Vol. 267, 290–298
  * Storey & Zeippen (2000) Storey, P. J., & Zeippen, C. J. 2000, MNRAS, 312, 813
  * Sturm et al. (2018) Sturm, E., Dexter, J., Pfuhl, O., et al. 2018, Nature, 563, 657.  https://doi.org/10.1038/s41586-018-0731-9
  * Tremaine et al. (2014) Tremaine, S., Shen, Y., Liu, X., & Loeb, A. 2014, The Astrophysical Journal, 794, 49.  https://doi.org/10.1088/0004-637x/794/1/49
  * U et al. (2022) U, V., Barth, A. J., Vogler, H. A., et al. 2022, ApJ, 925, 52
  * Ulrich et al. (1997) Ulrich, M.-H., Maraschi, L., & Urry, C. M. 1997, Annual Review of Astronomy and Astrophysics, 35, 445.  https://doi.org/10.1146/annurev.astro.35.1.445
  * Valtonen et al. (2008) Valtonen, M. J., Lehto, H. J., Nilsson, K., et al. 2008, Nature, 452, 851.  https://doi.org/10.1038/nature06896
  * Vanden Berk et al. (2001) Vanden Berk, D. E., Richards, G. T., Bauer, A., et al. 2001, AJ, 122, 549
  * Vanden Berk et al. (2004) Vanden Berk, D. E., Wilhite, B. C., Kron, R. G., et al. 2004, ApJ, 601, 692
  * Vaughan et al. (2016) Vaughan, S., Uttley, P., Markowitz, A. G., et al. 2016, MNRAS, 461, 3145
  * Villafaña et al. (2022) Villafaña, L., Williams, P. R., Treu, T., et al. 2022, ApJ, 930, 52
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Wang et al. (2020) Wang, S., Shen, Y., Jiang, L., et al. 2020, ApJ, 903, 51
  * Wilhite et al. (2006) Wilhite, B. C., Berk, D. E. V., Brunner, R. J., & Brinkmann, J. V. 2006, The Astrophysical Journal, 641, 78.  https://doi.org/10.1086/500421
  * Woo et al. (2015) Woo, J.-H., Yoon, Y., Park, S., Park, D., & Kim, S. C. 2015, ApJ, 801, 38
  * Yang et al. (2020) Yang, Q., Shen, Y., Chen, Y.-C., et al. 2020, Monthly Notices of the Royal Astronomical Society, 493, 5773.  https://doi.org/10.1093/mnras/staa645
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579
  * Yu (2002) Yu, Q. 2002, MNRAS, 331, 935
  * Zhang (2022) Zhang, X. 2022, MNRAS, 512, 1003