# Spectroscopic reverberation mapping of Quasar PKS 0736+017: Broad-Line
Region and Black-hole Mass

Shivangi Pandey1, Suvendu Rakshit1, Jong-Hak Woo2, and C. S. Stalin3
1Aryabhatta Research Institute of Observational Sciences, Nainital–263001,
Uttarakhand, India
2 Astronomy Program, Department of Physics and Astronomy, Seoul National
University, Seoul, 08826, Republic of Korea
3 Indian Institute of Astrophysics, Block II, Koramangala, Bangalore-560034,
India E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

To understand the mass distribution and co-evolution of supermassive black
holes with their host galaxy, it is crucial to measure the black hole mass of
AGN. Reverberation mapping is a unique tool to estimate the black hole masses
in AGN. We performed spectroscopic reverberation study using long-term
monitoring data with more than 100 spectra of a radio-loud quasar PKS 0736+017
to estimate the size of the broad line region (BLR) and black hole mass. The
optical spectrum shows strong H$\mathrm{\beta}$ and H$\mathrm{\gamma}$
emission lines. We generated the light curves of 5100Å continuum flux
($f_{5100}$), H$\beta$, and H$\gamma$. All the light curves are found to be
strongly variable with fractional variability of 69$\%$, 21$\%$, 30$\%$ for
V-band, H$\beta$, and H$\gamma$ light curves, respectively. Along with the
thermal contribution, non-thermal emission contributes to the estimated
continuum luminosity at 5100Å. Using different methods, e.g., CCF, JAVELIN,
von-neumann, we estimated the size of the BLR, which is found to be
66.4${}^{+6.0}_{-4.2}$ light days in the rest frame. The BLR size combined
with the line width of H$\beta$ provides a black hole mass of
7.32${}^{+0.89}_{-0.91}\times 10^{7}M_{\odot}$. The source closely follows the
BLR size-luminosity relation of AGN.

###### keywords:

galaxies:active – galaxies:individual: PKS0736+017 – quasars: supermassive
black holes – techniques: spectroscopic

††pubyear: 2022††pagerange: Spectroscopic reverberation mapping of Quasar PKS
0736+017: Broad-Line Region and Black-hole Mass–13

##  1 Introduction

Active galactic nuclei (AGNs) are powered by the accretion of matter onto a
central supermassive black hole (SMBH) of mass M > 10${}^{6}M_{\odot}$ (Woo &
Urry, 2002). These are one of the most bright and persistent objects in the
Universe ($\gtrapprox 10^{44}$ ergs s-1), that outshine their entire host
galaxy by emitting an abundant amount of radiation in a broad range of the
electromagnetic region. It shows flux variation over a time scale of months to
years (e.g., Wagner & Witzel, 1995; Ulrich et al., 1997). The emission from
AGN peaks at the UV/optical region of the spectral energy distribution (SED).
Due to the strong gravity of a black hole, matter attracted toward it, spirals
in and forms a disk-like structure known as an accretion disk. Surrounding
this is the broad-line region (BLR), which is made of gas clouds orbiting
around the SMBH with velocities of a few thousand km/s. The gas is ionized by
the continuum radiation and emits broad emission lines due to Doppler
broadening (Urry & Padovani, 1995).

The black hole mass is strongly correlated with the host galaxy properties
suggesting a co-evolution of black hole and host galaxy (Kormendy & Ho, 2013)
and subsequently studying these AGNs provides leverage to probe the growth and
evolution of black holes and their host galaxy across the Universe. The mass
of the black hole in AGN is challenging to measure because the bright central
core overpowers the host galaxy, and the spatial resolution needed to resolve
the central structure is beyond the capability of existing telescopes
(however, see Gravity Collaboration et al., 2018; GRAVITY Collaboration et
al., 2020). The black hole mass measurement of a radio-loud object is more
difficult due to the non-thermal contribution from jet, which affects the
optical continuum variation of AGN.

Reverberation Mapping (RM; Blandford et al., 1982; Peterson, 1993) uses the
flux variability to estimate black hole masses and constrain the geometry and
kinematics of the central engine. RM is based on the variation of line fluxes
due to the variation of optical/UV photons from the accretion disk. RM is a
well-known technique that has been applied so far in more than 100 objects to
provide the BLR size and black hole mass (e.g., Bahcall et al., 1972;
Blandford et al., 1982; Peterson, 1993; Peterson et al., 1998; Wandel et al.,
1999; Kaspi et al., 2000, 2007; Peterson et al., 2002, 2004, 2014; Greene &
Ho, 2005; Dietrich et al., 2012; Zu et al., 2013a; Bentz et al., 2009, 2010;
Barth et al., 2015; Du et al., 2014, 2015; Woo et al., 2015; Du et al., 2016a;
Pei et al., 2017; Fausnaugh et al., 2017; Du et al., 2016b, 2018; Shen et al.,
2016; Grier et al., 2017; Park et al., 2017; Rakshit et al., 2019; Zhang et
al., 2019a; Cho et al., 2020; Bontà et al., 2020; Rakshit, 2020; Bontà et al.,
2020; Williams et al., 2020; Amorim et al., 2021; Cackett et al., 2021;
Dehghanian et al., 2021; Bentz et al., 2022; Villafaña et al., 2022; U et al.,
2022). Remarkably it has provided a relation between H$\beta$ BLR size and
luminosity at 5100 Å ($L_{5100}$) that can be used to estimate BLR size and
black hole mass of any AGN having a single-epoch spectrum (Bentz et al.,
2013). Reverberation mapping study of high-accreting sources (Du et al.,
2016a, b, 2018) show significant deviation from the BLR size-luminosity
relation, that is found to be due to Eddington ratio/Fe II strength (Du &
Wang, 2019). Therefore, proper calibration of the size-luminosity relation is
crucial by increasing the sample size and objects with diverse properties.
Moreover, emission-line and continuum light curves can be used to estimate the
BLR size and constrain the geometry and kinematics of the BLR via geometrical
and dynamical models (e.g., Pancoast et al., 2014; Mandal et al., 2020).

RM of radio-loud AGN is challenging due to the contribution of non-thermal
emission, as a consequence RM have been successfully performed only for a few
radio-loud AGNs. For example, H$\beta$ lag of 3C120 and 3C273 has been well
measured by Peterson et al. (1998) and Kaspi et al. (2000), respectively (see
also Zhang et al., 2019b, for 3c273). Rakshit (2020) successfully measured
H$\beta$ lag of PKS1510-089. RM study based on Mg II line of radio-loud AGN
has also been carried out successfully for HE 0413-4031 by Zajaček et al.
(2020); however, in some cases strong non-thermal contribution prevented
reliable lag measurement (e.g., Nalewajko et al., 2019; Chavushyan et al.,
2020; Amaya-Almazán et al., 2022).

PKS 0736+017 (hereafter PKS0736) is a flat-spectrum radio quasar (FSRQ)
located at the redshift of 0.189 (Lynds, 1967; Tadhunter et al., 1993). The
SED of PKS0736 spans from the radio through the $\gamma$-ray wavelengths,
consisting of the typical double hump structure of Blazars where the low-
energy hump is due to synchrotron emission and high energy one is due to the
external Compton process (Impey & Neugebauer, 1988; Fossati et al., 1998;
Clements et al., 2003; Abdalla et al., 2020). The source shows a one-sided
parsec-scale jet (Kellermann et al., 1998; Bloom et al., 2009) and a compact
core (Gower et al., 1984; Romney et al., 1984) and is speculated to have been
hosted by an elliptical galaxy (Wright et al., 1998; McLure et al., 1999;
Falomo & Ulrich, 2000). The optical-UV spectrum is characterized by broad
emission lines and a big blue bump that is associated with thermal emission
from the accretion disk (Baldwin, 1975; Malkan & Moore, 1986). The previous
estimates of the black hole mass of PKS0736 based on the BLR size-luminosity
relation elucidates 108-108.73 M⊙ (Wandel et al., 1999; McLure & Dunlop, 2001;
Marchesini et al., 2004; Dai et al., 2007; Abdalla et al., 2020). Therefore,
the black hole mass of PKS0736 remains highly uncertain, and accurate
measurement of the black hole mass will help SED modeling. Moreover, the BLR
of PKS0736 is poorly studied due to the unavailability of long-duration
spectroscopic monitoring.

In this paper, we performed a spectroscopic variability study of PKS0736 to
estimate the size of the BLR and black hole mass through reverberation
mapping. The optical data is taken from the Steward Observatory (SO) as a part
of the spectropolarimetric monitoring project
111http://james.as.arizona.edu/~psmith/Fermi/, which is a support program for
the Fermi Gamma-Ray Space Telescope. Apart from the 3C273 reported by Zhang et
al. (2019a) and PKS 1510-089 reported by Rakshit (2020), PKS0736 have good
quality data with more than 100 spectra suitable for reverberation mapping
based on H$\beta$ emission line. The optical spectrum shows a blue continuum
and the presence of strong Balmer lines (H$\beta$ and H$\gamma$) and Fe II
emission. We performed several methods to calculate the size of the BLR and
the black hole mass of PKS0736 from time-series analysis. In section 2 we
describe the data analysis, and in section 3 we present the result of the
spectral analysis and time-lag measurements using various methods.
Subsequently, in section 4 we discuss the result and conclude in section 5.

##  2 Data

###  2.1 Optical data

Spectroscopic monitoring of PKS0736 was performed in SO as a part of
spectropolarimetric monitoring project (Smith et al., 2009), which is a
support program for the Fermi Gamma-Ray Space Telescope. Observations were
carried out using the 2.3m Bok Telescope on Kitt Peak and the 1.54m Kuiper
Telescope on Mount Bigelow in Arizona using the spectrophotometric instrument
SPOL (Schmidt et al., 1992) with a 600mm-1 grating, which provides the
spectral range of 4000–7550Å with a dispersion of 4 Å/pixel. Depending on the
width of the slit used for the observation, the resolution is typically
between 16-24 Å. Details regarding observations and data reduction are given
in Smith et al. (2009). In short, differential photometry using a standard
field star was performed to calibrate photometric magnitudes. The instrumental
magnitudes of the AGN and the comparison star are determined by using a
synthetic Johnson V bandpass for spectroscopic data. A total of 133 V-band
photometric observations carried out between November 2014 and May 2018 were
used. Spectra were flux-calibrated using the average sensitivity function
derived from multiple observations of several spectrophotometric standard
stars throughout an observing campaign. Final flux calibrations were
performed, re-scaling the spectrum from a given night to match the synthetic
V-band photometry of that night (see Smith et al., 2009). Therefore, a total
of 127 photometrically calibrated spectra obtained between November 2014 and
May 2018 were downloaded from the SO database and used in this work. Among the
127 spectra, 107 spectra were observed with a slit width of 4.1′′.

###  2.2 $\gamma$-ray and radio data

The $\gamma$-ray data within the energy range of 100MeV to 300GeV were
collected from the publicly available database of the Large Area Telescope
(LAT) on board the Fermi Gamma-Ray Space Telescope (Abdo et al., 2009) between
December 2014 and November 2017. The weekly-binned reduced light curve was
downloaded and used without any further processing from the light curve
Repository of
LAT222https://fermi.gsfc.nasa.gov/ssc/data/access/lat/LightCurveRepository/source.html?source_name=4FGL_J0739.2+0137
(Fermi Large Area Telescope Collaboratin, 2021).

The 15GHz radio data observed using the 40m Telescope at the Owens Valley
Radio Observatory (OVRO333https://sites.astro.caltech.edu/ovroblazars/) were
also collected (Richards et al., 2011).

##  3 Results and Analysis

###  3.1 Optical spectral analysis and lightcurves

First, the spectra are brought to the rest frame, then emission line fluxes
are measured. Two methods are usually adopted in previous studies for flux
measurement (1) direct integration of emission lines within a given window
after subtracting the power-law continuum (e.g., Kaspi et al., 2000; Grier et
al., 2012; Fausnaugh et al., 2017; Du et al., 2018) and (2) detailed multi-
component spectral decomposition (e.g., Barth et al., 2015; Rakshit et al.,
2019) which also isolates Fe II emission. However, the latter method demands
high S/N spectra, while the former works well even with moderate-low S/N
spectra. Both methods are found to provide consistent lag estimates (Du et
al., 2016a). Here we adopted the direct integration method due to its better
robustness and applicability to the moderate S/N spectra like the one used in
this work.

To calculate the emission line flux, we have fitted a power-law i.e.,
P($\lambda$) = $\alpha\lambda^{\beta}$ where $\alpha$ is a normalization
constant, and $\beta$ is the spectral index, locally for H$\beta$ plus [OIII]
region in the wavelength window 4760-4790Å and 5080-5120Å, and separately for
H$\gamma$ region in the wavelength window of 4270-4290Å and 4400-4405Å as
shown in Fig. LABEL:fig:spectrum_plot. Then, we subtracted the best-fit power-
law model from each spectrum for both regions separately. Then using the
power-law subtracted spectrum, we have integrated the flux values of each
emission line region within given wavelength windows as mentioned in Table 1.
It provides us with the flux for each line in each respective spectrum. Note
that we consider the region 4820-4900 Å to calculate H$\beta$ line flux
because beyond 4900 Å there is an excess contribution from the Fe II 4924 Å
multiplets and [OIII]4959 . This eventually creates a red asymmetry in the
line profile of H$\beta$ (see e.g., Joly, 1988). Since the Fe II emission is
also variable in AGNs it will also be present in the RMS spectrum and hinder
the calculation of width of the line emission of H$\beta$.

We have generated the light curves using the integrated flux values and the
photometric V-band light curve obtained from SO. Using the best-fit power-law
used for H$\beta$ and the [OIII] region, we have calculated the continuum flux
at 5100Å. As mentioned in the previous section, the flux calibration of SO
spectra is performed using standard stars observed throughout the campaign;
therefore, seeing effect, slit loss, and systematic uncertainty contributes to
the high uncertainty in the flux calibration. Hence, we re-scaled the flux
values based on the [O III] $\lambda$5007 emission lines considering it non-
variable during the time-scale of the monitoring program. The final light
curves are shown in Fig. 2. Uncertainties in the calculated emission line flux
and continuum flux include systematic uncertainty and Poisson noise added in
quadrature. The systematic uncertainties are estimated by the median filter
method as described in Zhang et al. (2019a), first, smoothing the light curve
by a median filter with five points and then subtracting the smoothed light
curve from the original one. The standard deviation of the residual is used as
the estimate of the systematic uncertainty. The flux values are provided in
Table 2.

Table 1: Emission line region enlisted to calculate their respective flux. Line | Wavelength range (Å)
---|---
H$\gamma$ | 4300-4400
Fe II | 4435-4685
H$\beta$ | 4820-4900
OIII | 4980-5020

Figure 2: light curves of PKS0736. Top to bottom: $\gamma$ ray flux in photons
s-1 cm-2, photometric V-band magnitude, spectroscopic flux at 5100Å in ergs
cm-2s-1 Å-1, H$\beta$ and H$\gamma$ flux in ergs cm-2s-1, radio flux in
Jansky(Jy), Non thermal dominance parameter (NTD). The vertical line at 57285
and 58250 is the MJD range which is used for time analysis and horizontal
lines represent NTD=1 and 2.

###  3.2 Variability

To quantify the variability amplitude of the light curves, we estimated the
fractional variability amplitude (Fvar) for V-band, $f_{5100}$, H$\beta$, and
H$\gamma$ light curves using the following equation (Rodriguez-Pascual et al.,
1997).

$\displaystyle
F_{\mathrm{var}}=\dfrac{\sqrt{(\sigma^{2}-<\sigma^{2}_{\mathrm{err}}>)}}{<f>}$
(1)

where $\sigma^{2}$ is the variance, $\sigma^{2}_{\mathrm{err}}$ is the mean
square error, and <$f$> is the arithmetic mean of the light curves. The ratio
of maximum to minimum flux variation ($R_{\mathrm{max}}$) was also calculated
for photometric and spectroscopic light curves. The values are given in Table
3. All the continuum and line light curves show significant variations with
Fvar=21% for H$\beta$ and 30% for H$\gamma$. The $R_{\mathrm{max}}$ is found
to 3 and 5 for H$\beta$ and H$\gamma$, respectively. The higher variability
found in H$\gamma$ compared to H$\beta$ agrees well with the previous finding
in the literature (see Bentz et al., 2010; Zhang et al., 2019a) and as
predicated from photo-ionization calculation (Korista & Goad, 2004).

Between MJD=57000 - 57285, the $\gamma$\- ray flux is in an active or flaring
state, and the same is reflected in the optical flux showing correlated
variation. It suggests a strong non-thermal contribution from the jet
dominates this part of the light curve. Although, the same effect is not
visible in the emission line fluxes. Considering the source is FSRQ, a non-
thermal contribution from the jet to the optical flux is expected. To
determine the dominance of non-thermal contribution over the thermal emission
from the accretion disk, we calculated the non-thermal dominance (NTD; Shaw et
al., 2012) using

$\displaystyle\mathrm{NTD}=\dfrac{L_{o}}{L_{p}}\quad\mathrm{and}\quad
L_{o}=L_{\mathrm{disk}}+L_{\mathrm{jet}}$ (2)

where $L_{o}$ is the observed continuum luminosity having a combination of
luminosity emitted from the accretion disk ($L_{\mathrm{disk}}$) and the jet
($L_{\mathrm{jet}}$). The $L_{p}$ is the predicted disk continuum luminosity
that can be assumed to be responsible for line luminosity coming out from the
BLR. Therefore, if the thermal emission from the disk only ionizes the broad
line clouds and no effect of non-thermal emission is present for the same,
then, $L_{p}$ = $L_{\mathrm{disk}}$ and NTD will equalize to 1 +
$L_{\mathrm{jet}}$ /$L_{\mathrm{disk}}$. If only thermal contribution from the
disk is responsible for continuum luminosity, NTD becomes unity. NTD can be
larger than 2 if the jet contribution is greater than the disk contribution.

To measure $L_{p}$, we used the correlation of
$L_{\mathrm{H\beta}}$-$L_{5100{\AA}}$ obtained by Rakshit et al. (2020) for
SDSS DR14 quasars and the $L_{\mathrm{H\beta}}$ estimated in this work via
$\log L_{\mathrm{H\beta}}=(1.057\pm 0.002)\log L_{5100{\AA}}+(-4.41\pm 0.10)$.
The variation of NTD with time is shown in the last panel of Fig. 2. We note
the following points: (1) The NTD varies between $1.1-11.2$ with an average of
2.7, suggesting that the non-thermal emission from the jet contributes to the
continuum variation, and (2) From MJD = 57000 - 57285; the NTD shows strong
spikes, which are correlated with the flaring event in the $\gamma$-ray light
curve, increasing up to NTD=11.

The correlation between the continuum and emission-line luminosity of PKS0736
is studied. In Fig. LABEL:fig:NTD_plot, H$\beta$ luminosity (upper panel),
$\gamma$-ray (middle), NTD (bottom panel) are plotted against 5100Å continuum
luminosity. The NTD is larger than 2 when $\log L_{5100}$ > 44.6. Therefore, a
significant non-thermal contribution in the optical emission is found above
this luminosity and disk contribution found to dominate below this range. This
is also clear from the correlation of line and continuum luminosity plot,
which shows a large scatter above $\log L_{5100}$ > 44.6 (see upper panel of
Fig. LABEL:fig:NTD_plot). The Spearman correlation coefficient calculated for
$L_{\mathrm{H\beta}}$ vs. $L_{5100}$ is noted as $r_{s}=0.61$ with a p-value
of no correlation to be $p=10^{-15}$ whereas a much stronger and positive
correlation is seen for data points where NTD<2 with $r_{s}=0.71$ and
$p=10^{-7}$. As shown in the Fig.2, the $\gamma$-ray flux changes by a factor
of 20 between MJD=57000-57285 with the peak flux of 1.1$\times 10^{-7}$
photons s-1 cm-2 on MJD=57124.0. Apart from the flaring event during
MJD=57000-57285, the light curve is mostly quiescent. The $F_{\mathrm{var}}$
is found to be 99.5% for full light curve and 51% for MJD>57285\. Therefore,
the variations in $\gamma$-ray is much larger than the optical flux variation.
Optical continuum flux is found to be correlated with the $\gamma$-ray flux
with $r_{s}=0.79$ and $p=10^{-15}$ (see Fig.LABEL:fig:NTD_plot).

MJD | V | f(5100) | f(H$\beta$) | f(H$\gamma$)
---|---|---|---|---
(1) | (2) | (3) | (4) | (5)
56989 | 15.09 | 5.28$\pm$0.59 | 29.15$\pm$3.20 | 13.63 $\pm$1.62
56990 | 15.21 | 5.71$\pm$0.49 | 28.45$\pm$4.14 | 11.22 $\pm$1.63
56991 | 15.48 | 6.72$\pm$1.18 | 25.32$\pm$2.64 | 15.98 $\pm$2.95
57015 | 16.42 | 1.06$\pm$0.12 | 28.09$\pm$4.19 | 13.05 $\pm$1.55
Table 2: Columns as: (1) Modified Julian date (MJD), (2) photometric V-band
magnitude, (3) spectroscopic flux at 5100$\AA$ in units of 10-15 ergs cm-2s-1
Å-1, (4) and (5) are H$\beta$ and H$\gamma$ line flux in units of 10-15 ergs
cm-2s-1, respectively. This table is available in its entirety in machine-
readable form. A portion is shown here for guidance.  Table 3: Variability
statistics.

light curve | Median flux |  $F_{\mathrm{var}}$ ($\%$) | $R_{\mathrm{max}}$
---|---|---|---
(1) | (2) | (3) | (4)
V-band | 0.76$\pm$0.07 | 69.86$\pm$4.30 | 12.94$\pm$0.67
f5100 | 0.82$\pm$0.11 | 85.43$\pm$5.57 | 18.87$\pm$4.11
H$\beta$ | 21.95$\pm$0.66 | 21.10$\pm$1.79 | 2.92$\pm$0.63
H$\gamma$ | 10.11$\pm$0.41 | 29.73$\pm$2.25 | 5.41$\pm$0.88

  * •

Notes. Columns are: (1) light curve and (2) median flux of the light curve in
units of 10-15 ergs s-1cm-2 Å-1 for $f_{5100}$, and 10-15 ergs s-1cm-2 for
emission lines; (3) fractional rms variability in percentage; and (4) the
ratio of maximum to minimum flux variation.

###  3.3 Time delay measurement

As mentioned in the previous section, a strong flare in the $\gamma$-ray light
curves is seen between 57000-57285 MJD due to the dominance of non-thermal
contribution over the disk contribution. Therefore, we excluded this region
and used the time window between 57285 and 58250 for further analysis.

####  3.3.1 Interpolated cross Correlation Analysis

The cross-correlation technique quantifies the degree of similarity between
two sets of numbers. The process involves two light curves to assess
information between peak/centroid values. One can shift either continuum
concerning emission-line light curve or vice versa and calculate the
correlation coefficient ($r$) at each shifted time. Then the peak or the
centroid ($\tau_{\mathrm{cent}}$) covering 0.8$\times r_{\mathrm{peak}}$ is
considered as the lag (Peterson et al., 2004). Following Peterson et al.
(1998), we performed interpolated cross-correlation (ICCF)
444https://bitbucket.org/cgrier/python_ccf_code/src/master/. First, the cross-
correlation function (CCF) is measured by interpolating the continuum series
keeping the emission-line light curve unchanged, and second interpolating the
emission-line light curve keeping the continuum light curve unchanged. The
final ICCF is determined by averaging these two results. The CCF is
characterized by (1) its peak value $r_{\mathrm{max}}$, (2) the time delay
corresponding to this value $\tau_{\mathrm{peak}}$, and (3) the centroid
$\tau_{\mathrm{cent}}$ of the CCF. The Monte Carlo method of flux
randomization (FR) and random subset selection (RSS) is used to estimate the
uncertainty in the lag (Peterson et al., 1998; Peterson et al., 2004).

We have performed the auto-correlation (ACF) between the two continuum light
curves, the 5100Å and photometric V-band in Fig. 4. The ACF shows a lag of
zero-day as expected. We also estimated the time lag between photometric
V-band flux vs. H$\beta$ and H$\gamma$. ICCF shows a strong peak around 80
days for H$\beta$ and H$\gamma$, and a minor peak at $\sim$180 days that could
be due to the seasonal gaps (see discussion in section 4.1). Therefore, we
restricted our analysis to 150 days covering the major peak and estimated lag
within that. The ICCF method provides lags between V-band vs. H$\beta$ of
79.0${}^{+7.1}_{-5.0}$ days and between V-band and H$\gamma$ of
71.8${}^{+10.9}_{-18.1}$ days in the observer-frame. We obtained a shorter lag
for H$\gamma$ compared to H$\beta$, which agrees with the previous studies on
different objects (e.g., Bentz et al., 2010).

Figure 4: Cross-correlation function and lag distribution for V-band vs
$f_{5100}$ (left), H$\beta$ (middle) and H$\gamma$ (right). The ICCF (line)
and centroid probability distribution from ICCF (filled histogram) are shown.

####  3.3.2 JAVELIN

We have also used JAVELIN555https://github.com/nye17/JAVELIN developed by Zu
et al. (2011, 2013b) which estimates lag by modeling the continuum and line
light curves. According to Kelly et al. (2009), quasar variability can be well
described by a damped random walk (DRW) process. JAVELIN models the continuum
light curve using a DRW having two parameters: amplitude and time scale of
variability. Then the emission line light curve is generated, which is a
shifted, scaled, and smoothed version of the continuum light curve. It
interpolates between data points and self-consistently estimates and includes
the uncertainties in the interpolation. Error estimation uses the Markov Chain
Monte Carlo (MCMC) approach to calculate the statistical confidence limits on
each best-fit parameter and accurately model the continuum and emission-line
light curve simultaneously.

In Fig. 5, the probability distribution of the observed frame lag is plotted,
as computed by JAVELIN. The left panel of the Figure shows the probability
distribution of JAVELIN lag between the V-band light curve vs. H$\beta$; here,
too, we noticed a significant peak at $\sim$80 days and multiple minor peaks.
Therefore, we calculated a lag between 0 to 150 days, similar to ICCF. We also
plotted the probability distribution of lag between the V-band light curve vs.
H$\gamma$ in the right panel for the lag range of 0 to 150 days. The plots
show prominent peaks at $\sim$ 78 and 75 days for H$\beta$ and H$\gamma$,
respectively. These results are consistent with the ICCF. Time delays obtained
from JAVELIN are given in Table 4.

Figure 5: Probability distribution of observed frame lag computed by JAVELIN
between V-band vs H$\beta$ (left), and with H$\gamma$ (right).

####  3.3.3 The von Neumann and Bartels estimators

This method is based on a mean-square successive-difference estimator (see von
Neumann, 1941) which handles sparsely and irregularly sampled data and is
relevant for the process underlying the variability that cannot be adequately
modeled. It does not rely on data interpolation, binning in correlation space,
ergodicity arguments, and stochastic models for quasar variability. Hence
time-delay measurements are restricted to the information embedded in the
light curves being processed (Chelouche et al., 2017).

Figure 6: Probability distribution of the observed frame time lag based on the
von Neumann estimator (top panels) and the Bartels estimator (bottom panels)
for H$\beta$ (left) and H$\gamma$ (right).

The time delay distribution obtained from the von Neumann method after a Monte
Carlo simulation of FR-RSS, as done for the CCF analysis, is shown in the
upper panels of Fig. 6. The H$\beta$ and H$\gamma$ show firm peaks at $\sim$68
and 60 days, respectively. The Bartels estimator is the modified method of the
von Neumann estimator (Bartels, 1982) which uses a ranked light curve instead
of a sorted one, which can also be used to measure time delay based on the
regularity or randomness of data. The lag results are given in Table 4 between
time interval taken from 0-300 days. The results obtained using von Neumann
and Bartels methods are consistent with the lags obtained using ICCF and
JAVELIN.

Method | Lag (days) |
---|---|---
| V vs. H$\beta$ | V vs. H$\gamma$
ICCF | 79.0${}^{+7.1}_{-5.0}$ | 71.8${}^{+10.9}_{-18.1}$
JAVELIN | 78.2${}^{+6.8}_{-3.5}$ | 74.7${}^{+22.7}_{-2.8}$
von Neumann | 67.6${}^{+37.6}_{-31.4}$ | 60.0${}^{+33.5}_{-31.6}$
Bartels | 73.7${}^{+37.5}_{-28.3}$ | 62.2${}^{+31.4}_{-31.1}$
Table 4: Columns as follows: (1) method used to calculate the lag(days), (2)
depicts lag in observed frame between photometric V-Band vs. H$\beta$ and (3)
lag in observed frame between photometric V-Band vs. H$\gamma$

The above methods suggest a time lag of $\sim$80 days between the photometric
V-band continuum and the H$\beta$ line light curve. We visually analyze the
consistency of the measured lag; the V-band flux light curve along with the
H$\beta$ light curve back-shifted by 80 days are plotted in Fig. 7. The
continuum and back-shifted line light curves agreeably match. Therefore, we
considered the lag of 79.0${}^{+7.1}_{-5.0}$ days in the observed frame
obtained by the ICCF method as the best lag measurement for PKS 0736+017.

Figure 7: Normalized V-band light curve plotted along with the H$\beta$ light
curve back-shifted by 80 days.

####  3.3.4 Mean and RMS spectra

We constructed mean spectrum and root-mean-square (rms) spectrum as follows.

$\displaystyle\bar{F(\lambda)}=\dfrac{1}{N}\sum_{i=0}^{N-1}F_{i}(\lambda),$
(3)

and

$\displaystyle
S(\lambda)=\sqrt{\Big{[}\dfrac{1}{N-1}\sum_{i=0}^{N-1}(F_{i}(\lambda)-\bar{F(\lambda)}^{2})\Big{]}}$
(4)

where F${}_{i}(\lambda)$ is the ith spectrum of the N(127) spectra that
comprise the database. The mean spectrum is the average spectrum of all the
spectra obtained during the campaign, whereas the rms spectrum is based on the
variation around this mean. The rms spectrum isolates the constant features
such as narrow emission line and the host galaxy contribution of the spectrum
and hinders accurate line width measurements. However, the rms spectrum is
often too noisy therefore measuring the line width from the rms spectrum is
challenging.

In Fig. 8, we showed the mean and rms spectra constructed from all the flux-
calibrated spectrum without subtracting from the fitted continuum. The lower
panel depicting the rms spectra with and without subtraction from fitted
power-law continuum. Both the mean and rms spectra show the presence of strong
Balmer lines. Most importantly, the narrow [O III]5007, present in the mean
spectrum, disappeared in the rms spectrum suggesting the flux calibration is
proper. As noted earlier, although the SO monitoring campaign has used
different slits throughout the campaign, the majority of them (107 out of 127)
are obtained using 4.1′′ slit width. To find out any effect of slit width in
the construction of the mean and rms spectra and subsequent analysis, we re-
constructed them using only spectra taken with 4.1′′ slit. The resultant mean
and rms spectra (see Fig. A) is remarkably similar to that obtained using the
entire spectral data.

Figure 8: Mean and RMS spectra of PKS0736 constructed from the nightly
spectrum after rescaling based on [O III]$\lambda$5007\. Top: mean spectrum.
Bottom: RMS spectra before (upper plot) and after (lower plot) the continuum
subtraction in the region of wavelength window 4760-4790Å and 5080-5120Å.

####  3.3.5 Detection significance of cross-correlation lag

Recent reverberation mapping studies on large samples perform simulation to
access the lag significance (Penton et al., 2022; U et al., 2022; Homayouni et
al., 2022), however, there are no universal criteria to define the lag
reliability. We used the publicly available PyIICCF
code666https://github.com/legolason/PyIICCF/ developed by Guo et al. (2022)
which uses the method described in U et al. (2022) to assess the lag
significance of our measurement. In this method, two uncorrelated red-noise
light curves are simulated with the same S/N and cadence as the observed data
to determine the probability of finding the peak of correlation
($r_{\mathrm{max}}$) larger than the observed ICCF peak via null hypothesis
test.

For this purpose, we start by simulating the continuum and line light curves
using the DRW model. First, we fitted the continuum V-band and H$\beta$ line
light curves using the DRW model and used the best-fit parameters to simulate
the corresponding light curves. Second, we simulated a 100 times longer light
curve than the observed data and randomly selected a portion of the light
curve. Third, we added Gaussian noise to the mock light curve based on the
uncertainty in the observed light curve. Forth, we down-sampled the mock light
curve to have the same cadence and length as the observed light curve. We
simulated 1000 such sets of continuum and line light curves. Finally, we
cross-correlated mock continuum light curve with the observed line light curve
and vice versa within a lag search range of 0-300 days. From these CCF, we
obtained a distribution of $\tau_{\mathrm{cent}}$ and $r_{\mathrm{max}}$ (see
Fig. 9). Then we derived the p($r_{\mathrm{max}}$) which is the ratio of the
number of positive lag with $r_{\mathrm{max}}$ larger than the observed
$r_{\mathrm{max}}$ and $\tau_{\mathrm{cent}}>0$ to the total number of
positive lag. This p($r_{\mathrm{max}}$) provides a way to access the
robustness of our lag estimation, a smaller p($r_{\mathrm{max}}$) gives more
robust and reliable lag detection. We found p($r_{\mathrm{max}}$) is 0.05
which is $<0.2$ adopted by U et al. (2022) as the limiting condition for
reliable lag determination. This suggests that our lag measurement is robust
and reliable.

Figure 9:  Left: The plot shows the CCF coefficient(r) Vs. lag between V-band
continuum and H$\beta$ line light curve taken in 0-300 days, for all 1000
simulations and the red curve shows the r Vs. lag for the observed light
curves with 3-$\sigma$ confidence interval. Right: CCF peak and centroid
distribution of the coefficient and the no of simulations showing r>0.5.

####  3.3.6 Line width and black hole mass measurement

To estimate the black hole masses, we measured the FWHM and the line
dispersion ($\sigma_{\mathrm{line}}$) of the H$\beta$ emission line from both
the mean and rms spectra after fitting the continuum and subtracting with the
same. The full width at half maximum (FWHM) is calculated by finding
0.5$\times$F$(\lambda)_{\mathrm{max}}$ both from the left ($\lambda_{l}$) and
right side ($\lambda_{r}$) of the curve and subtracting the two wavelengths
i.e, $\lambda_{r}-\lambda_{l}$ (Peterson et al., 2004). Whereas, for
calculation of $\sigma_{\mathrm{line}}$, the flux weighted line center was
first determined as follows:

$\displaystyle\lambda_{0}=\dfrac{\int\lambda f_{\lambda}\,d\lambda}{\int
f_{\lambda}\,d\lambda}$ (5)

and the line dispersion is

$\displaystyle\sigma_{\mathrm{line}}^{2}=\dfrac{\int\lambda^{2}f_{\lambda}\,d\lambda}{\int
f_{\lambda}\,d\lambda}-\lambda_{0}^{2}$ (6)

The values are mentioned in Table 5. The same for spectra taken with 4.1′′
slit is given in Table 6, which shows consistent results justifying the line
width measurements from the entire spectral data.

Assuming that the black hole’s gravitational potential governs the gas motion
in the BLR, black hole mass can be estimated by combining the size of the BLR
($R_{\mathrm{BLR}}$) and the velocity width of broad emission lines
($\Delta$V) based on the virial relation:

$\displaystyle M_{\mathrm{BH}}=\dfrac{f\times R_{\mathrm{BLR}}(\Delta
V)^{2}}{G}$ (7)

where $f$ is a dimensionless scale factor that depends on the kinematics and
geometry of BLR gas clouds. The BLR size (calculated from time lag
$R_{\mathrm{BLR}}$=c$\tau$) and line width (from mean and rms spectra) with
$f$ of 4.47($\sigma$) or 1.12(FWHM) taken from (Woo et al., 2015) is used to
calculate black hole mass. The Monte Carlo bootstrap method (Peterson et al.,
2004) was used to estimate uncertainty in the line width measurements and
hence the black hole mass. N spectra were chosen randomly from a set of N
spectra without replacement for each realization, and line width was measured
from the mean and rms spectra. We varied the endpoints randomly within
$\pm$10Å of the initially selected H$\beta$ region (i.e., 4800-4920Å).
Finally, a total of 5000 realizations were performed, providing a distribution
of FWHM and $\sigma_{\mathrm{line}}$. The mean of the distribution is the
final line width, and the standard deviation of the distribution ( 34
percentile of both sides of the mean i.e. the 16th and the 84th percentiles of
the distribution) was considered the measurement uncertainty. The resultant
line width measurement along with black hole mass is given in Table 5.

Using the cross-correlation method, we obtained a lag ($\tau_{\mathrm{cent}}$)
between optical V-band and H$\beta$ (H$\gamma$) of 66.4${}^{+6.0}_{-4.2}$
(60.4${}^{+9.2}_{-15.3}$) days in rest frame. We used the virial relation with
$f$ of 4.47(1.12) for $\sigma_{\mathrm{line}}$(FWHM) providing four different
black hole masses based on the four different choices of line widths. We note
that $\sigma_{\mathrm{line}}$ is less sensitive to the line peak and that FWHM
is less sensitive to the line wing; therefore, black hole masses based on the
$\sigma_{\mathrm{line}}$ are widely adopted as the best mass measurement
(e.g., Peterson et al., 2004; Peterson, 2014). Therefore, a black hole mass of
7.32${}^{+0.89}_{-0.91}\times 10^{7}M_{\odot}$ from the
$\sigma_{\mathrm{line}}$ of the rms spectrum is found for PKS0736.

Table 5: Rest-frame resolution uncorrected line width and black hole mass
measurements from mean and rms spectra.

| Spectrum | Type | $\Delta V(\mathrm{km\,s^{-1}})$ |  $M_{\mathrm{BH}}$($\times 10^{7}M_{\odot}$)
---|---|---|---|---
| Mean | FWHM | 2160.14${}^{+198.65}_{-210.48}$ | 6.78${}^{+1.30}_{-1.26}$
|  | $\sigma_{\mathrm{line}}$ | 1230.92${}^{+79.46}_{-74.03}$ | 8.78${}^{+1.17}_{-1.02}$
| RMS | FWHM | 1946.48${}^{+284.58}_{-257.41}$ | 5.50${}^{+1.72}_{-1.36}$
|  | $\sigma_{\mathrm{line}}$ | 1123.04${}^{+66.58}_{-72.48}$ | 7.32${}^{+0.89}_{-0.91}$

  * •

Notes. Columns: Spectrum, line width type, line width, black hole mass.

##  4 Discussion

###  4.1 Effect of seasonal gaps on lag estimate

Due to the low declination, the monitoring observations of PKS0736 are
affected by the seasonal gaps, as can be seen from the light curve plot. We
have simulated light curves to investigate the impact of seasonal gaps on the
lag measurement. We have constructed a mock continuum light curve using the
DRW model of JAVELIN having the same properties (i.e., same amplitude and
time-scale of variation) as the observed V-band continuum light curve. Then we
simulated the emission line light curve from the mock continuum light curve
with a shift of 80 days based on our calculated lag from the observed data,
smoothing and scaling the mock light curve to mimic the observed light curve.
We then down-sampled the light curves to have the same observational gaps.

We then performed all the time-series analysis methods to recover the time lag
of 80 days from the simulated data. The results are shown in Fig. 12 and
tabulated in Table 7 for a single set of mock lightcurves. We noticed a strong
peak at $\sim$80 days in all the methods, i.e., ICCF, JAVELIN, von Neumann,
and Bartels, that well-recovered the input lag suggesting the observed cadence
and seasonal gaps do not affect our lag estimation. Although a secondary peak
of around 180 days has been seen from CCF and JAVELIN, it is well-separated
from the primary peak and does not affect our lag measurements. Moreover, the
peak of the ICCF ($r_{\mathrm{max}}\sim 0.5$) matches well with that obtained
from the observed data (see middle panel of Fig.4). We further simulated such
500 sets of continuum and line lightcurves and repeated the above process of
measuring lag. In Fig. LABEL:fig:simulated_tlaghist we showed the distribution
of lag ratio recovered from the ICCF to the input lag of 80 days, which
provides a ratio of unity. Based on these simulations, we conclude a H$\beta$
lag of $\sim$80 days from the observed data of PKS0736.

###  4.2 Size-luminosity relation

The PKS0736 is a radio-loud AGN having a combination of thermal emission from
the disk and non-thermal emission from the jet elucidated by Fig.
LABEL:fig:NTD_plot. Consequently, the measurement of $L_{5100}$ is affected by
non-thermal emission. However, we considered the quiescent state by removing
the flaring region as visible in the $\gamma$-ray light curve (Fig. 2).

We plotted PKS0736 in the size-Luminosity diagram as shown in Fig. 10 along
with various objects from the literature. The PKS0736 follows the best-fit
size luminosity relation given by Bentz et al. (2013). Hence, its position is
found to be consistent with the size$-$luminosity relation of other AGN. Du &
Wang (2019) has provided a new scaling relation taking into account the Fe II
emission contamination in the spectrum of quasars and for especially high
accreting AGN. This new scaling relation can explain the deviation seen in
high-accreting sources from the size-luminosity relation. We estimate Fe II
strength from the continuum subtracted mean spectrum. The strength of Fe II
(RFeII), which is the flux ratio of Fe II (4435-4685 Å) to H$\beta$ is found
to be 0.58$\pm$0.08. Therefore, PKS0736 is moderately strong Fe II emitter.
The time-lag estimated using scaling relation of Du & Wang (2019) is around 63
light days, which is slightly lower than the 87 light days estimated using
Bentz et al. (2013).

The disk luminosity of PKS0736 was measured by several authors and found to be
in the range of 1044.6 to 1045.7 ergs s-1 (see Abdalla et al., 2020; Dai et
al., 2007; McLure & Dunlop, 2001). The scaling relation RBLR = 1017
$\sqrt{L_{\mathrm{disk}}/10^{45}}$cm (Ghisellini & Tavecchio, 2009), provides
RBLR = 24-86 light-days, our measured BLR size of 66 light-days is within this
range. Alternatively, the NTD in the quiescent state can be used as an
estimator for removing non-thermal contribution from the observed L5100. The
median of NTD which is measured as 2.51${}^{+1.52}_{-0.65}$ indicates
$\sim$33$\%$ of disk contribution in L5100. If we use this revised L5100 value
of 1.44 $\times$ 1044, the expected BLR size based on Bentz et al. (2013) is
found to be 41 light-days, which is slightly smaller than our estimated rest
frame BLR size of 66 light-days.

We have also performed the ICCF analysis between radio (see Fig. 2) and V-band
light curves. The correlation between radio and optical is found to be very
weak with $r_{\mathrm{max}}$ $\sim$ 0.3. The measured lag is $\sim$ 255 light
days which is much longer than the lag of H$\beta$. The weak correlation
between radio and optical and a much longer radio lag than H$\beta$ lag
suggest the radio emitting region is different than the H$\beta$ emitting
region and radio is not the primary contributor in the V-band fluxes within
the time-range studied.

Figure 10: The plot is between the BLR size vs. $L_{5100}$ relation of AGNs.
The source PKS0736 is following the relation very well. The best-fit relation
of Bentz et al. (2013) is shown along with various RM results from the
literature.

###  4.3 The effect of diffuse continuum

Continuum reverberation mapping study of AGNs showed the UV-optical lags
generally follow the $\tau\propto\lambda^{4/3}$ in agreement with the simple
irradiated disk models; however, larger by a factor of 2-4 than predicated by
standard accretion disk (see, e.g., Fausnaugh et al., 2016; Edelson et al.,
2019; Guo et al., 2022, and references therein). Diffuse continuum (DC)
emission from the BLR has been suggested as the possible origin of this
larger-than-expected lag (Korista & Goad, 2001; Lawther et al., 2018;
Chelouche et al., 2019; Korista & Goad, 2019). Recently, Netzer (2022) showed
that time-dependent emission of the diffuse BLR gas could explain the observed
large lag in the continuum RM. The DC could affect our NTD calculation and the
BLR size-luminosity relation. Assuming the radiation pressure confined (RPC)
cloud models and a covering factor of 0.2, the Netzer (2022) showed that the
DC fraction could be 23% of the total continuum at 5100A for a source with
$L_{\mathrm{total,5100}}=L_{\mathrm{disk,5100}}+L_{\mathrm{DC,5100}}=10^{44}$
erg s-1. The DC has two effects on the BLR size-luminosity relation. First,
the true $L_{\mathrm{disk,5100}}$ will be lower than the observed
$L_{\mathrm{total,5100}}$ and this will reduce the expected BLR size by 12%
($\sim$4 days) for a source with $L_{\mathrm{total,5100}}=10^{44}$ erg s-1
(see their equation 3). The NTD values are calculated based on the observed to
predicted luminosity ratio. The $L_{5100,disk}$ should be lower if DC is taken
into account; however, the predicted luminosity will also be affected as the
correlation between L(H$\beta$) vs. $L_{\mathrm{total,5100}}$ for the SDSS
DR14 quasars is not corrected for DC. Second, since the lag is measured with
respect to V-band, the reference point will be shifted toward a higher lag due
to DC. Therefore, the true lag should be smaller by a few days, $\sim$4 days
for $L_{\mathrm{total,5100}}=10^{44}$ erg s-1 (see equation 6 of Netzer,
2022). Both the size of the BLR and the luminosity of the continuum will be
smaller by a similar extent for all the objects depending on the luminosity.
Hence, the size-luminosity relation may not be affected much. A detailed
quantitative analysis of this issue is beyond the scope of this paper.

###  4.4 Black hole mass measurement

Previous studies have estimated a range of black hole masses of PKS0736 from
108 to 10${}^{8.73}M_{\odot}$ (Schmidt et al., 1992; Abdalla et al., 2020; Dai
et al., 2007; Marchesini et al., 2004; McLure & Dunlop, 2001). These estimates
are based on the scaling relationship (RBLR $\propto$ $L_{5100}^{\alpha}$)
(e.g., Kaspi et al., 2000) measured from the single-epoch spectrum. Since
PKS0736 is highly variable, the measurement of black hole mass from the
single-epoch spectrum is highly uncertain. We note that single-epoch mass
measurements are affected by choice of line width used to estimate the black
hole. Generally, most of the single-epoch masses are calculated using FWHM;
however, it has been found from the reverberation mapping study of AGNs with
multiple emission lines that $\sigma_{\mathrm{line}}$ provides a better
measure of the black hole mass (Peterson et al., 2004). The rms spectrum
isolates the non-varying components (e.g., narrow emission lines, host galaxy)
and provides a robust estimation of black hole masses. Our reverberation-based
black hole mass is found to be between 5.50-8.78 $\times 10^{7}M_{\odot}$
depending on the line width estimator and type of spectrum. However, we
consider the black hole masses of 7.32${}^{+0.89}_{-0.91}\times
10^{7}M_{\odot}$ calculated from the $\sigma_{\mathrm{line}}$ of the rms
spectrum as the best black hole mass measurement.

The bolometric luminosity of PKS0736 is calculated as 5.52 $\times$ 1045 erg
s-1 using mean $L_{5100}$ and the relation $L_{\mathrm{BOL}}$ = 9 $\times$
$L_{5100}$ (Kaspi et al., 2000). The Eddington luminosity ($L_{\mathrm{EDD}}$)
estimated as 9.22 $\times$ 1045 erg s-1 whereas the Eddington ratio
($\lambda_{\mathrm{EDD}}$) is 0.60 calculated by $L_{\mathrm{EDD}}$ = 1.26
$\times$ 10${}^{38}\,M_{\mathrm{BH}}$ and the black hole mass based on the
$\sigma_{\mathrm{line}}$ of rms spectrum. This suggest that PKS0736 is
accreting at a sub-Eddington rate.

##  5 Conclusions

We analyzed the SO spectro-photometric data for PKS0736 obtained between
November 2014 and May 2018, rendering more than 100 spectra and photometric
data points. We found strong variability in the photometric light curve with
$F_{\mathrm{var}}$=69.86$\pm$4.30%, which is also reflected in the
spectroscopic continuum with zero lag. Both the H$\beta$ and H$\gamma$ show
high variability with $F_{\mathrm{var}}$ of 21.10$\pm$1.79% and 29.73$\pm$
2.25%, respectively. The estimated NTD is found to be $2.51$ suggesting about
67% contribution from the non-thermal emission towards the observed mean
$L_{5100}$ of 6.13$\times 10^{44}$ ergs s-1. Using the cross-correlation
method, we obtained a lag between optical V-band and H$\beta$ (H$\gamma$) of
66.4${}^{+6.0}_{-4.2}$ (60.4${}^{+9.2}_{-15.3}$) in the rest frame. Using
virial relation and a scale factor of $f=1.12$ ($f=4.47$) for FWHM
($\sigma_{\mathrm{line}}$) a black hole mass of 7.32${}^{+0.89}_{-0.91}\times
10^{7}M_{\odot}$ is obtained from the rms spectrum and
$\sigma_{\mathrm{line}}$ of H$\beta$ line profile. The position of PKS0736 is
consistent with the size$-$luminosity relation of AGNs.

## Acknowledgements

We thank the referee for comments and suggestions that helped to improve the
quality of the manuscript. SR thanks Hagai Netzer for useful discussion about
diffuse continuum. This publication makes use of data products from the Fermi
Gamma-ray Space Telescope and accessed from the Fermi Science Support Center.
Data from the Steward Observatory spectro-polarimetric monitoring project were
used. This research has made use of data from the OVRO 40-m monitoring program
(Richards et al., 2011) which is supported in part by NASA grants NNX08AW31G,
NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911. SR
acknowledges the partial support of SRG-SERB, DST, New Delhi through grant no.
SRG/2021/001334. J.H.W. acknowledges funding from the Basic Science Research
Program through the National Research Foundation of the Korean government
(NRF- 2021R1A2C3008486).

## Data Availability

The data underlying this article were accessed from the Steward Observatory
database with the following link: http://james.as.arizona.edu/~psmith/Fermi/.
The derived data generated in this research are available in the article and
in its online supplementary material.

## References

  * Abdalla et al. (2020) Abdalla H., et al., 2020, Astronomy and Astrophysics, 633, A162
  * Abdo et al. (2009) Abdo A. A., et al., 2009, ApJ, 707, 55
  * Amaya-Almazán et al. (2022) Amaya-Almazán R. A., Chavushyan V., Patiño-Álvarez V. M., 2022, The Astrophysical Journal, 929, 14
  * Amorim et al. (2021) Amorim A., et al., 2021, Astronomy and Astrophysics, 654, A85
  * Bahcall et al. (1972) Bahcall J. N., Kozlovsky B.-Z., Salpeter E. E., Bahcall J. N., Kozlovsky B.-Z., Salpeter E. E., 1972, ApJ, 171, 467
  * Baldwin (1975) Baldwin J. A., 1975, ApJ, 201, 26
  * Bartels (1982) Bartels R., 1982, Journal of the American Statistical Association, 77, 40
  * Barth et al. (2015) Barth A. J., et al., 2015, The Astrophysical Journal Supplement Series, 217, 26
  * Bentz et al. (2009) Bentz M. C., et al., 2009, Astrophysical Journal, 705, 199
  * Bentz et al. (2010) Bentz M. C., et al., 2010, Astrophysical Journal Letters, 720, L46
  * Bentz et al. (2013) Bentz M. C., et al., 2013, Astrophysical Journal, 767
  * Bentz et al. (2022) Bentz M. C., Williams P. R., Treu T., 2022, ] 10.48550/arxiv.2206.03513
  * Blandford et al. (1982) Blandford R. D., McKee C. F., Blandford R. D., McKee C. F., 1982, ApJ, 255, 419
  * Bloom et al. (2009) Bloom S., Marscher A., Moore E., Gear W., Teräsranta H., Valtaoja E., Aller H., Aller M., 2009, The Astrophysical Journal Supplement Series, 122, 1
  * Bontà et al. (2020) Bontà E. D., et al., 2020, The Astrophysical Journal, 903, 112
  * Cackett et al. (2021) Cackett E. M., Bentz M. C., Kara E., 2021, iScience, 24, 102557
  * Chavushyan et al. (2020) Chavushyan V., Patiño-Álvarez V. M., Amaya-Almazán R. A., Carrasco L., 2020, ApJ, 891, 68
  * Chelouche et al. (2017) Chelouche D., Pozo-Nuñez F., Zucker S., 2017, The Astrophysical Journal, 844, 146
  * Chelouche et al. (2019) Chelouche D., Pozo Nuñez F., Kaspi S., 2019, Nature Astronomy, 3, 251
  * Cho et al. (2020) Cho H., et al., 2020, ApJ, 892, 93
  * Clements et al. (2003) Clements S. D., Jenks A., Torres Y., 2003, The Astronomical Journal, Volume 126, Issue 1, pp. 37-46., 126, 37
  * Dai et al. (2007) Dai H., Xie G. Z., Zhou S. B., Li H. Z., Chen L. E., Ma L., 2007, The Astronomical Journal, Volume 133, Issue 5, pp. 2187-2191., 133, 2187
  * Dehghanian et al. (2021) Dehghanian M., et al., 2021, The Astrophysical Journal, 906, 14
  * Dietrich et al. (2012) Dietrich M., et al., 2012, Astrophysical Journal, 757, 53
  * Du & Wang (2019) Du P., Wang J.-M., 2019, The Astrophysical Journal, 886, 42
  * Du et al. (2014) Du P., et al., 2014, The Astrophysical Journal, 782, 45
  * Du et al. (2015) Du P., et al., 2015, The Astrophysical Journal, 806, 22
  * Du et al. (2016a) Du P., et al., 2016a, ApJ, 820, 27
  * Du et al. (2016b) Du P., et al., 2016b, ApJ, 825, 126
  * Du et al. (2018) Du P., et al., 2018, ApJ, 856, 6
  * Edelson et al. (2019) Edelson R., et al., 2019, The Astrophysical Journal, 870, 123
  * Falomo & Ulrich (2000) Falomo R., Ulrich M. H., 2000, A&A, 357, 91
  * Fausnaugh et al. (2016) Fausnaugh M. M., et al., 2016, The Astrophysical Journal, 821, 56
  * Fausnaugh et al. (2017) Fausnaugh M. M., et al., 2017, The Astrophysical Journal, 840, 97
  * Fermi Large Area Telescope Collaboratin (2021) Fermi Large Area Telescope Collaboratin 2021, The Astronomer’s Telegram, 15110, 1
  * Fossati et al. (1998) Fossati G., Maraschi L., Celotti A., Comastri A., Ghisellini G., 1998, MNRAS, 299, 433
  * GRAVITY Collaboration et al. (2020) GRAVITY Collaboration et al., 2020, Astronomy and Astrophysics, 643
  * Ghisellini & Tavecchio (2009) Ghisellini G., Tavecchio F., 2009, MNRAS, 397, 985
  * Gower et al. (1984) Gower A. C., Hutchings J. B., Gower A. C., Hutchings J. B., 1984, AJ, 89, 1658
  * Gravity Collaboration et al. (2018) Gravity Collaboration et al., 2018, Nature, 563, 657
  * Greene & Ho (2005) Greene J. E., Ho L. C., 2005, The Astrophysical Journal, 630, 122
  * Grier et al. (2012) Grier C. J., et al., 2012, The Astrophysical Journal, 755, 60
  * Grier et al. (2017) Grier C. J., et al., 2017, ApJ, 851, 21
  * Guo et al. (2022) Guo H., Barth A. J., Wang S., 2022, arXiv e-prints, p. arXiv:2207.06432
  * Homayouni et al. (2022) Homayouni Y., et al., 2022, ApJ, 926, 225
  * Impey & Neugebauer (1988) Impey C. D., Neugebauer G., 1988, AJ, 95, 307
  * Joly (1988) Joly M., 1988, A&A, 192, 87
  * Kaspi et al. (2000) Kaspi S., Smith P. S., Netzer H., Maoz D., Jannuzi B. T., Giveon U., 2000, The Astrophysical Journal, 533, 631
  * Kaspi et al. (2007) Kaspi S., Brandt W. N., Maoz D., Netzer H., Schneider D. P., Shemmer O., 2007, The Astrophysical Journal, 659, 997
  * Kellermann et al. (1998) Kellermann K. I., Vermeulen R. C., Zensus J. A., Cohen M. H., 1998, AJ, 115, 1295
  * Kelly et al. (2009) Kelly B. C., Bechtold J., Siemiginowska A., 2009, ApJ, 698, 895
  * Korista & Goad (2001) Korista K. T., Goad M. R., 2001, ApJ, 553, 695
  * Korista & Goad (2004) Korista K. T., Goad M. R., 2004, ApJ, 606, 749
  * Korista & Goad (2019) Korista K. T., Goad M. R., 2019, MNRAS, 489, 5284
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Lawther et al. (2018) Lawther D., Goad M. R., Korista K. T., Ulrich O., Vestergaard M., 2018, MNRAS, 481, 533
  * Lynds (1967) Lynds C. R., 1967, ApJ, 147, 837
  * Malkan & Moore (1986) Malkan M. A., Moore R. L., 1986, ApJ, 300, 216
  * Mandal et al. (2020) Mandal A. K., Rakshit S., Stalin C. S., et al., 2020, Monthly Notices of the Royal Astronomical Society, 501, 3905
  * Marchesini et al. (2004) Marchesini D., Celotti A., Ferrarese L., 2004, Monthly Notices of the Royal Astronomical Society, 351, 733
  * McLure & Dunlop (2001) McLure R. J., Dunlop J. S., 2001, Monthly Notices of the Royal Astronomical Society, Volume 327, Issue 1, pp. 199-207., 327, 199
  * McLure et al. (1999) McLure R. J., Kukula M. J., Dunlop J. S., Baum S. A., O’Dea C. P., Hughes D. H., 1999, MNRAS, 308, 377
  * Nalewajko et al. (2019) Nalewajko K., Gupta A. C., Liao M., Hryniewicz K., Gupta M., Gu M., 2019, A&A, 631, A4
  * Netzer (2022) Netzer H., 2022, MNRAS, 509, 2637
  * Pancoast et al. (2014) Pancoast A., Brewer B. J., Treu T., et al., 2014, Monthly Notices of the Royal Astronomical Society, 445, 3073
  * Park et al. (2017) Park S., et al., 2017, ApJ, 847, 125
  * Pei et al. (2017) Pei L., et al., 2017, The Astrophysical Journal, 837, 131
  * Penton et al. (2022) Penton A., et al., 2022, MNRAS, 509, 4008
  * Peterson (1993) Peterson B. M., 1993, Publications of the Astronomical Society of the Pacific, 105, 247
  * Peterson (2014) Peterson B. M., 2014, Space Sci. Rev., 183, 253
  * Peterson et al. (1998) Peterson B. M., Wanders I., Bertram R., Hunley J. F., Pogge R. W., Wagner R. M., 1998, The Astrophysical Journal, 501, 82
  * Peterson et al. (2002) Peterson B. M., et al., 2002, The Astrophysical Journal, 581, 197
  * Peterson et al. (2004) Peterson B. M., et al., 2004, The Astrophysical Journal, 613, 682
  * Peterson et al. (2014) Peterson B. M., et al., 2014, Astrophysical Journal, 795, 149
  * Rakshit (2020) Rakshit S., 2020, A&A, 642, A59
  * Rakshit et al. (2019) Rakshit S., Woo J.-H., Gallo E., et al., 2019, ApJ, 886, 93
  * Rakshit et al. (2020) Rakshit S., Stalin C. S., Kotilainen J., 2020, ApJS, 249, 17
  * Richards et al. (2011) Richards J. L., et al., 2011, ApJS, 194, 29
  * Rodriguez-Pascual et al. (1997) Rodriguez-Pascual P. M., et al., 1997, The Astrophysical Journal Supplement Series, Volume 110, Issue 1, pp. 9-20., 110, 9
  * Romney et al. (1984) Romney J., et al., 1984, A&A, 135, 289
  * Schmidt et al. (1992) Schmidt G. D., Stockman H. S., Smith P. S., 1992, ApJ, 398, L57
  * Shaw et al. (2012) Shaw M. S., Romani R. W., Cotter G., et al., 2012, The Astrophysical Journal, 748, 49
  * Shen et al. (2016) Shen Y., et al., 2016, ApJ, 818, 30
  * Smith et al. (2009) Smith P. S., Montiel E., Rightley S., Turner J., Schmidt G. D., Jannuzi B. T., 2009, arXiv e-prints, p. arXiv:0912.3621
  * Tadhunter et al. (1993) Tadhunter C. N., Morganti R., di Serego Alighieri S., Fosbury R. A. E., Danziger I. J., 1993, MNRAS, 263, 999
  * U et al. (2022) U V., et al., 2022, The Astrophysical Journal, 925, 52
  * Ulrich et al. (1997) Ulrich M.-H., Maraschi L., Urry C. M., 1997, Annual Review of Astronomy and Astrophysics, 35, 445
  * Urry & Padovani (1995) Urry C. M., Padovani P., 1995, Publications of the Astronomical Society of the Pacific, 107, 803
  * Villafaña et al. (2022) Villafaña L., et al., 2022, The Astrophysical Journal, 930, 52
  * Wagner & Witzel (1995) Wagner S. J., Witzel A., 1995, Annual Review of Astronomy and Astrophysics, 33, 163
  * Wandel et al. (1999) Wandel A., Peterson B. M., Malkan M. A., 1999, ApJ, 526, 579
  * Williams et al. (2020) Williams P. R., et al., 2020, The Astrophysical Journal, 902, 74
  * Woo & Urry (2002) Woo J.-H., Urry C. M., 2002, ApJ, 579, 530
  * Woo et al. (2015) Woo J. H., Yoon Y., Park S., Park D., Kim S. C., 2015, Astrophysical Journal, 801, 1
  * Wright et al. (1998) Wright S. C., McHardy I. M., Abraham R. G., 1998, MNRAS, 295, 799
  * Zajaček et al. (2020) Zajaček M., et al., 2020, ApJ, 896, 146
  * Zhang et al. (2019a) Zhang Z.-X., et al., 2019a, The Astrophysical Journal, 876, 49
  * Zhang et al. (2019b) Zhang Z.-X., Du P., Smith P. S., et al., 2019b, ApJ, 876, 49
  * Zu et al. (2011) Zu Y., Kochanek C. S., Peterson B. M., 2011, ApJ, 735, 80
  * Zu et al. (2013a) Zu Y., Kochanek C. S., Kozłowski S., Udalski A., 2013a, Astrophysical Journal, 765
  * Zu et al. (2013b) Zu Y., Kochanek C. S., Kozłowski S., et al., 2013b, ApJ, 765, 106
  * von Neumann (1941) von Neumann J., 1941, The Annals of Mathematical Statistics, 12, 367

##  Appendix A Analysis of total 107 spectra taken with 4.1′′ slit

Although monitoring observations were carried out using multiple slits, the
majority of them (84% spectra) are taken with 4.1′′ slit width. To investigate
whether the different slit widths used during spectroscopic observations of
PKS0736 have any effect on the line width measurement and subsequently the
black hole mass estimation, we constructed the mean and rms spectra only with
the 107 spectra taken with 4.1′′ slit width.

Figure 11: Mean and RMS spectra of PKS0736 from the 107 spectra obtained using
4.1′′ slit after rescaling based on [O III]$\lambda$5007 without continuum
subtracted from nightly spectra. Table 6: Rest-frame resolution uncorrected
line width and black hole mass measurements from mean and rms spectra.

| Spectrum | Type | $\Delta V(\mathrm{km\,s^{-1}})$ |  $M_{\mathrm{BH}}$($\times 10^{7}M_{\odot}$)
---|---|---|---|---
| Mean | FWHM | 2147.08${}^{+196.95}_{-222.35}$ | 6.60${}^{+1.26}_{-1.29}$
|  | $\sigma_{\mathrm{line}}$ | 1223.05${}^{+82.62}_{-73.72}$ | 8.55${}^{+1.19}_{-1.00}$
| RMS | FWHM | 1934.60${}^{+257.63}_{-258.21}$ | 5.36${}^{+1.52}_{-1.33}$
|  | $\sigma_{\mathrm{line}}$ | 1064.81${}^{+80.42}_{-79.15}$ | 6.48${}^{+1.01}_{-0.92}$

  * •

Notes. Columns: Spectrum, line width type, line width, black hole mass.

Fig. 11 is remarkably similar to the mean and rms spectra obtained from entire
spectroscopic data (see Fig. 8) provided from SO irrespective of the slit
width used during observation. Moreover, we estimated FWHM and $\sigma$ from
these mean and rms spectra. The results are given in Table 6, which shows
consistent results with that obtained with entire spectroscopic data (see
Table 5). It suggests that the use of entire spectral data do not affect black
hole mass measurement.

##  Appendix B Time series analysis of simulated light curves

We have performed time-series analysis on the simulated data as mentioned in
section 4.1. The results of the time-series analysis are shown in Fig. 12 for
ICCF, JAVELIN, Von-Neumann, and Bartels. Table 7 summarizes the results.

Figure 12: Lag estimation from the simulated light curves. Top: Cross-correlation function and lag distribution for simulated data. The ICCF (line) and centroid probability distribution from ICCF (filled histogram) are shown. Middle: Lag probability distribution obtained from JAVELIN for simulated data. Bottom: Lag probability distribution obtained from Von-Neumann (left) and Bartels (right) for simulated data. Method | Lag (days)
---|---
ICCF | 84.0${}^{+7.1}_{-11.0}$
JAVELIN | 97.7${}^{+2.8}_{-17.7}$
von Neumann | 90.4${}^{+13.4}_{-13.1}$
Bartels | 90.5${}^{+12.6}_{-11.5}$
Table 7: Columns as follows: (1) the method used to calculate the lag (days),
(2) lag recovered from the simulated light curves with an input lag of 80
days.