# Four hints and test candidates of the local cosmic expansion

Kazuhiro Agatsuma<EMAIL_ADDRESS>School of Physics and Astronomy
and Institute for Gravitational Wave Astronomy, University of Birmingham,
Edgbaston, Birmingham B15 2TT, United Kingdom

###### Abstract

The expansion of the universe on short distance scales is a new frontier to
investigate the dark energy. The excess orbital decay in binary pulsars may be
related to acceleration by the local cosmic expansion, called the cosmic drag.
Modern observations of two independent binaries (PSR B1534+12 and PSR
B1913+16) support this interpretation and result in a scale-independent
expansion with viscous uniformity, in which binary systems have a smaller
expansion rate than the Hubble constant. This paper shows additional
evidential binaries (PSR J1012–5307 and PSR J1906+0746), supporting the cosmic
drag picture. The total anomaly of the conventional model is about
$3.6\,\sigma$ including two evidential binaries reported before. In addition,
an observable range of the cosmic drag has been calculated for typical models
of both NS-NS binary and NS-WD binary. In this region, six test candidates are
listed with predictions of the excess orbital decay.

###### pacs:

##  I Introduction

A discovery of the expansion of the universe made a paradigm shift. The
recession velocities of distant galaxies follow the Hubble–Lemaître law Hubble
; Lemaitre . In these two decades, progress in observations revealed that the
cosmic expansion is accelerating Riess_1998 ; Perlmutter_1999 . The
acceleration mechanism, which is called the dark energy, is one of the most
intriguing mysteries in our universe. The cosmic expansion is visible on large
scale observations but not valid on small scales ($\lesssim 70$ Mpc) because
local peculiar velocities exceed the recession velocity Freedman_2001 . Since
the gravitationally bound systems maintain own size, the cosmic expansion is
separated from local dynamics as the current standard interpretation.

Another interpretation survives with following general relativity (GR)
Anderson1995 ; Cooperstock1998 ; Bonnor1999 ; Adkins2007 ; Carrera2010 ;
Price2012 ; Giulini2014 . They argue that the cosmic expansion can be
acceleration on local dynamics. The formulation was built using the geodesic
deviation equation Cooperstock1998 , in which Friedmann-Lemaître-Robertson-
Walker (FLRW) coordinates is transformed to Fermi normal coordinates as an
approximation of the locally inertial frame. The equations of motion to the
lowest order are expressed by

$\frac{d^{2}x_{\mathrm{F}}^{k}}{dt^{2}}-\left(\frac{\ddot{a}}{a}\right)x_{\mathrm{F}}^{k}=0$
(1)

where $x_{\mathrm{F}}^{k}$ is the geodesic distance in the Fermi normal
coordinates with $k$ from 1 to 3, and $a$ is the scale factor. The expansion
term is called the cosmic drag Anderson1995 , which represents a stretch of
space. The object’s position is not fixed on space but dragged by a flow of
the expanding space. In this framework, the Hubble–Lemaître law is related to
the cosmic drag by

$\ddot{R}_{\mathrm{H}}=\dot{H}R+H\dot{R}=(\ddot{a}/a)R=-qH^{2}R.$ (2)

The last term employs the deceleration parameter $q$, which is defined as

$q:=-\ddot{a}a/\dot{a}^{2}=-(\ddot{a}/\dot{a})H^{-1}=-(\ddot{a}/a)H^{-2}$ (3)

using $H=\dot{a}/a$ and $\dot{H}=(\ddot{a}/a)-H^{2}$. This expression has been
applied to orbital systems Carrera2010 ; Giulini2014 . The two body problem is
expressed by the pseudo-Newtonian approach:

$\ddot{R}=L^{2}/R^{3}-GM_{0}/R^{2}-qH^{2}R$ (4)

with $L=R^{2}\dot{\phi}$. Here, $R$ is the radial geodesic distance, $L$ the
conserved angular momentum per unit mass in planar polar coordinates
($R,\phi$), $G$ the gravitational constant, $M_{0}$ the mass of the central
body, and the present value of $H$ is the Hubble constant $H_{0}\approx 70$ km
$\mathrm{s}^{-1}\mathrm{Mpc}^{-1}$ (this value is used in this paper). The
same expression can be applied for the scales of galaxy clusters Nandra_1 ;
Nandra_2 . A critical radius can exist at

$R_{\mathrm{c}}=\left(\frac{GM_{0}}{-q_{0}H_{0}^{2}}\right)^{1/3}$ (5)

for $q_{0}<0$ as the present cosmological value of $q$, in which the
gravitational attraction force equals the local expansion ($L=0$ and
$\ddot{R}=0$). This model is consistent with observations because the
recession velocity asymptotes the Hubble–Lemaître law ($\dot{R}=H_{0}R$) on
large scales ($R\gg R_{\mathrm{c}}$). However, we cannot test this theory by
searching for deviations from standard Keplerian orbits because they are too
small to measure Carrera2010 ; Giulini2014 .

A unique approach has been proposed to observe this local expansion effect
using binary pulsar systems Agatsuma_2020 . By comparing with accurate
predictions of the orbital decay by gravitational wave (GW) emission, it is
shown that the braking effect of the local cosmic expansion can be observed by
precise measurements. The orbital decay of a binary system is given by

$\displaystyle\frac{dr}{dt}$ $\displaystyle=$ $\displaystyle
H_{\mathrm{s}}\cdot
r-\frac{64}{5}\frac{G^{3}}{c^{5}}\frac{(m_{\mathrm{1}}m_{\mathrm{2}})(m_{\mathrm{1}}+m_{\mathrm{2}})}{r^{3}}$
(6) $\displaystyle\equiv$ $\displaystyle H_{\mathrm{s}}\cdot
r-\frac{K_{\mathrm{gw}}}{r^{3}},$

including both the local expansion and GW emission Agatsuma_2020 ;
PetersMathews ; Peters . Here, $c$ is the speed of light, $r$ is the orbital
radius (separation) of binary stars, $m_{1}$ and $m_{2}$ are each mass of the
binary, and $H_{\mathrm{s}}$ is the Hubble parameter in this spiral decay
system. By differentiating Eq. (6) and using Eq. (2), effective force
expression is obtained:

$\ddot{r}=-q_{\mathrm{s}}H_{\mathrm{s}}^{2}r-3K_{\mathrm{gw}}^{2}r^{-7}.$ (7)

A purely Keplerian orbit is achieved at the critical radius $r_{\mathrm{c}}$
where the local expansion equals GW reaction. At this point, $dr/dt=0$ and
$\ddot{r}=0$ are imposed in Eq. (6) and Eq. (7), respectively. These provide

$r_{\mathrm{c}}=\left(\frac{K_{\mathrm{gw}}}{H_{\mathrm{s}}}\right)^{1/4}$ (8)

and

$r_{\mathrm{c}}=\left(\sqrt{-\frac{3}{q_{\mathrm{s}}}}\frac{K_{\mathrm{gw}}}{H_{\mathrm{s}}}\right)^{1/4}.$
(9)

Since both radii have to be equivalent, the deceleration parameter of this
system is determined to be $q_{\mathrm{s}}=-3$. This value brought a
characteristic of viscous uniformity, corresponding to $H_{\mathrm{s}}\approx
H_{0}/6$ Agatsuma_2020 .

According to GR, binary star systems dominated by GW reaction have to be the
slowest orbital decay. However, both PSR B1534+12 Stairs2002 ; B1534+12  and
PSR B1913+16 Weisberg2010 ; Weisberg2016  show slower decay than the theory,
called the excess $\dot{P_{b}}$, with total anomaly of about $2.5\,\sigma$
Agatsuma_2020 . A faster decay is possible by other energy losses, but slower
decay is forbidden by GR. The above model can explain this deviation from GR.
The cosmic drag picture is a promising model, featured by scale-independent
and uniform expansion.

This is still a hypothesis due to the consistency with few observations of the
$2.5\,\sigma$ anomaly. In order to validate it furthermore, it is significant
to scrutinize other observations. This paper shows additional two evidential
binary pulsars with future test candidates.

##  II Excess orbital decay

The estimation of the orbital decay in a binary star system is written by

$\dot{P_{b}}^{\mathrm{obs}}=\dot{P_{b}}^{\mathrm{Kin}}+\dot{P_{b}}^{\mathrm{GW}}+\dot{P_{b}}^{\mathrm{X}}.$
(10)

Here, $\dot{P_{b}}^{\mathrm{obs}}$ is the observed orbital decay (time
derivative of the orbital period), $\dot{P_{b}}^{\mathrm{Kin}}$ is the all
kinematic contributions in a galactic potential, $\dot{P_{b}}^{\mathrm{GW}}$
is the theoretical estimate from GW reaction (intrinsic decay),
$\dot{P_{b}}^{\mathrm{X}}$ is the excess $\dot{P_{b}}$. Other contributions
like the mass loss are omitted. The conventional model (no local expansion)
expects $\dot{P_{b}}^{\mathrm{X}}=0$ when GW emission is a dominant factor of
the energy loss.

In the cosmic drag picture, the local cosmic expansion brake the orbital
decay, which makes $\dot{P_{b}}^{\mathrm{X}}$ as predicted below. According to
Agatsuma_2020 , the orbital decay by GW emission is described by

$\frac{\dot{P_{b}}}{P_{b}}=\frac{3}{2}\left[H_{\mathrm{s}}-\frac{K_{\mathrm{gw}}}{r^{4}}\frac{1+(73/24)e^{2}+(37/96)e^{4}}{(1-e^{2})^{7/2}}\right],$
(11)

including the contribution from the cosmic drag. Here, $r$ is the orbital
radius (as the eccentricity $e=0$) or the semimajor axis (as $e\neq 0$). Thus,
the excess $\dot{P_{b}}$ by the cosmic drag ($\dot{P_{b}}^{\mathrm{CD}}$) is
the first term of it:

$\dot{P_{b}}^{\mathrm{CD}}=\frac{3}{2}H_{\mathrm{s}}P_{b}.$ (12)

This is applicable for $r\leq r_{\mathrm{c}}$. Since the system can be
regarded as the Keplerian orbit around $r_{\mathrm{c}}$, the orbital period is
written by

$\displaystyle P_{b}$
$\displaystyle=2\pi\left[\frac{r^{3}}{G(m_{1}+m_{2})}\right]^{1/2}$ (13)
$\displaystyle=2\pi\left[\frac{k^{3}}{G(m_{1}+m_{2})}\right]^{1/2}\left(\frac{K_{\mathrm{gw}}}{H_{\mathrm{s}}}\right)^{3/8}$
(14)

where the radius is normalized by $r_{\mathrm{c}}$ and $k=r/r_{\mathrm{c}}$ is
used. From the above equations, the expected excess is

$\dot{P_{b}}^{\mathrm{CD}}=3\pi
H_{\mathrm{s}}\left[\frac{k^{3}}{G(m_{1}+m_{2})}\right]^{1/2}\left(\frac{K_{\mathrm{gw}}}{H_{\mathrm{s}}}\right)^{3/8}$
(15)

with viscous uniformity ($H_{\mathrm{s}}=H_{0}/6$) Agatsuma_2020 .

For a region of $r_{\mathrm{c}}\leq r\ll R_{\mathrm{c}}$, as discussed in the
prior paper Agatsuma_2020 , the cosmic drag completely cancel the GW reaction
(orbital decay), which indicates

$\dot{P_{b}}^{\mathrm{CD}}=-\dot{P_{b}}^{\mathrm{GW}}.$ (16)

Here, $\dot{P_{b}}^{\mathrm{GW}}$ is obtained from Eq. (11) as
$H_{\mathrm{s}}=0$. A fixed orbit is expected in spite of GW emission.

If the observed $\dot{P_{b}}^{\mathrm{X}}$ is consistent with the prediction
$\dot{P_{b}}^{\mathrm{CD}}$, it supports the cosmic drag picture.

##  III The third and fourth evidential objects

PSR J1012+5307 has become observable of the cosmic drag by recent refinements
of parameters Ding_2020 . J1012+5307 is a binary system, which consist of a
pulsar of $1.82\,M_{\odot}$ and a white dwarf of $0.174\,M_{\odot}$ (NS-WD
system) with a small eccentricity ($e=1.2\times 10^{-6}$). The orbital radius
slightly exceeds the critical radius ($r=1.24r_{\mathrm{c}}$) according to the
orbital period $P_{b}$ of 0.6 days. It means that the local expansion cancels
the orbital decay of GW emission under the cosmic drag picture, and hence, the
orbit should asymptote Eq. (4) (almost Keplerian) without orbital decay
Agatsuma_2020 . When the system follows Keplerian orbit without orbital decay,
$\dot{P_{b}}^{\mathrm{obs}}=\dot{P_{b}}^{\mathrm{Kin}}$ should be observed in
Eq. (10). In other words,
$\dot{P_{b}}^{\mathrm{X}}=-\dot{P_{b}}^{\mathrm{GW}}$ (Eq. (16)) is expected
when $\dot{P_{b}}^{\mathrm{GW}}$ is included in the estimate. According to the
latest observation of PSR J1012+5307 (see Table 5 in Ding_2020 ), the relevant
parameters are

$\displaystyle\dot{P_{b}}^{\mathrm{GW}}$ $\displaystyle=$
$\displaystyle-13\times 10^{-15},$ (17)
$\displaystyle\dot{P_{b}}^{\mathrm{X}}$ $\displaystyle=$
$\displaystyle(10.6\pm 6.1)\times 10^{-15},$ (18)

which is consistent with the above interpretation. This system looks a fixed
orbit as $\dot{P_{b}}^{\mathrm{X}}$ cancels $\dot{P_{b}}^{\mathrm{GW}}$. It
may be explained by the local cosmic expansion. This is the first hint to
support the cosmic drag picture for both a region of $r>r_{\mathrm{c}}$ and an
NS-WD binary system.

Besides, a binary neutron star PSR J1906+0746 is also marginally evidential.
The predicted $\dot{P_{b}}^{\mathrm{CD}}$ is $0.8\times 10^{-14}$, which is
consistent with the observed $\dot{P_{b}}^{\mathrm{X}}$ of $(3\pm 3)\times
10^{-14}$ (see Table 3 in VanLeeuwen_2015 ). The deviation from the
conventional model is only $1\,\sigma$.

By adding anomalies in the above two evidential objects to the previous value
($2.5\,\sigma$), the total significance beyond the conventional model is about
$3.6\,\sigma$ (99.97%). Their relevant parameters are listed in Table 1. The
observed $\dot{P_{b}}^{\mathrm{X}}$ in J1012+5307, B1534+12, and J1906+0746
have a good agreement with the cosmic drag prediction
($\dot{P_{b}}^{\mathrm{X}}\approx\dot{P_{b}}^{\mathrm{CD}}$), and B1913+16 has
a close value. All four independent binaries in Table 1 have a positive
$\dot{P_{b}}^{\mathrm{X}}$ (slower decay) with $1.0-1.8\,\sigma$ significance,
which is a forbidden region of GR without the cosmic drag.

##  IV Test region and candidates

There is an observable range of the cosmic drag in the binary pulsars. For a
region of $r<r_{\mathrm{c}}$, the shorter the radius, the smaller the
expansion effect. For a region of $r>r_{\mathrm{c}}$, the expected difference
($\dot{P_{b}}^{\mathrm{CD}}=-\dot{P_{b}}^{\mathrm{GW}}$) is getting smaller
and smaller with increasing the radius. For the comparison,
$\dot{P_{b}}^{\mathrm{CD}}$ has to be larger than the observational
uncertainty. The binary pulsars with too small separation or too large
separation are out of scope.

Figure 1 (top) shows $\dot{P_{b}}^{\mathrm{GW}}$ with/without
$\dot{P_{b}}^{\mathrm{CD}}$ for NS-NS binary calculated by the above equations
(11), (15), and (16). Fig. 1 (bottom) shows the expected
$\dot{P_{b}}^{\mathrm{CD}}$ for both NS-NS binary and NS-WD binary. Modern
observations have an uncertainty of $\dot{P_{b}}^{\mathrm{X}}$ of about
$0.4-1.1\times 10^{-14}$. Therefore, the observable range requires
$\dot{P_{b}}^{\mathrm{CD}}\gtrsim 1.0\times 10^{-14}$, which is roughly
$0.4\,r_{\mathrm{c}}\lesssim r\lesssim 1.8\,r_{\mathrm{c}}$ for NS-NS binary,
corresponding to

$0.2\lesssim P_{b}\lesssim 2.0\,\,\mathrm{[days]}$ (19)

as $P_{b}=0.8$ days at $r=r_{\mathrm{c}}$. In the case of NS-WD binary, the
observable range is shallower ($0.6\,r_{\mathrm{c}}\lesssim r\lesssim
1.4\,r_{\mathrm{c}}$), corresponding to

$0.2\lesssim P_{b}\lesssim 0.6\,\,\mathrm{[days]}$ (20)

as $P_{b}=0.4$ days at $r=r_{\mathrm{c}}$. Note that $H_{\mathrm{s}}$ is
observable only for $r<r_{\mathrm{c}}$ by Eq. (12). Increase of mass and/or
eccentricity expand the observable range because the increased
$r_{\mathrm{c}}$ makes the peak excess higher. The He star companion,
redbacks, and black-widow pulsars would be difficult for this test because the
binary evolution is not GW dominant, meaning increased uncertainties.

Figure 1: Rough estimates of $\dot{P_{b}}$ and the excess $\dot{P_{b}}$ due to
the cosmic drag ($\dot{P_{b}}^{\mathrm{CD}}$). For simplicity, the same mass
condition ($1.4\,M_{\odot}$) for neutron star (NS) and the 1/10 ratio
($0.14\,M_{\odot}$) for white dwarf (WD) are assumed with $e=0$. Viscous
uniformity ($H_{\mathrm{s}}=H_{0}/6$) is adopted for $r<r_{\mathrm{c}}$.

In addition to the evidential four binary pulsars, there are six test
candidates, which can estimate $\dot{P_{b}}^{\mathrm{CD}}$ in this observable
range. Their parameters are listed in Table 2. In this list, the first two
(J0509+3801 and J0751–1807) have no information about
$\dot{P_{b}}^{\mathrm{X}}$. The last two (J1829+2456 and B2127–11C) are
uncertain values for $\dot{P_{b}}^{\mathrm{X}}$. Except for B2127–11C, only
J1738–0333 meets the conventional model ($\dot{P_{b}}^{\mathrm{X}}=0$) in this
region. J1756–2251 has a negative $\dot{P_{b}}^{\mathrm{X}}$ (faster decay),
which fit neither the conventional model nor cosmic drag picture.

[htb]  Binary pulsar Pulsar mass Companion mass Eccentricity Orbital period
Separation $\dot{P_{b}}^{\mathrm{X}}$ $\dot{P_{b}}^{\mathrm{CD}}$ Reference
$m_{1}(M_{\odot})$ $m_{2}(M_{\odot})$ (e)a $P_{b}$ (days)a $(r_{\mathrm{c}})$b
($10^{-14}$) ($10^{-14}$)c J1012+5307 1.817d 0.174(11)d $1.2\times 10^{-6}$
0.605 1.24(3) 1.1(6) 1.3(1) Ding_2020  B1534+12 1.3330(2) 1.3455(2) 0.274
0.421 0.57 2.0(1.1) 2.1 B1534+12  J1906+0746 1.291(11) 1.322(11) 0.085 0.166
0.35 3(3) 0.8 VanLeeuwen_2015  B1913+16 1.438(1) 1.390(1) 0.617 0.323 0.28
0.5(4), 0.8(5) 1.6 Weisberg2016 , Weisberg2010

Table 1:  Evidential binary pulsars.

[htb]  Binary pulsar Pulsar mass Companion mass Eccentricity Orbital period
Separation $\dot{P_{b}}^{\mathrm{X}}$ $\dot{P_{b}}^{\mathrm{CD}}$ Reference
$m_{1}(M_{\odot})$ $m_{2}(M_{\odot})$ (e)a $P_{b}$ (days)a $(r_{\mathrm{c}})$b
($10^{-14}$) ($10^{-14}$)c J0509+3801 1.34(8) 1.46(8) 0.586 0.380 0.34(1) —
1.9 Lynch_2018  J0751–1807 1.64(15)e 0.16(1)e $3.3\times 10^{-6}$ 0.263
0.74(2) — 1.3 Desvignes_2016  J1738–0333 1.46(6) 0.181(8) $3.4\times 10^{-7}$
0.355 0.89(1) 0.2(4) 1.7 Freire_2012  J1756–2251 1.341(7) 1.230(7) 0.181 0.320
0.52 $-1.7^{+0.9}_{-0.6}$ 1.6 Ferdman_2014  J1829+2456 1.306(4) 1.299(4) 0.139
1.176 1.25 $-2.1(1.1)$f 2.3 Haniewicz_2021  B2127–11C 1.358(10) 1.354(10)
0.681 0.335 0.25 $-1(13)$ 1.6 Jacoby_2006

Table 2:  Test candidates.

  * a

The number of digits are rounded for simplicity.

  * b

Semimajor axis when $e\neq 0$. The uncertainty is due to mass estimation but
it ($\sim 2\%$) from the Hubble tension is omitted.

  * c

The main uncertainty is due to mass estimation for $r>r_{\mathrm{c}}$, or the
Hubble tension of $\sim 7\%$ for $r<r_{\mathrm{c}}$. The latter is omitted.

  * d

The pulsar mass is derived from the companion mass Antoniadis_2016  and mass
ratio of 10.44 Sanchez_2020 .

  * e

For the mass estimation, $\dot{P_{b}}^{\mathrm{GW}}$ is used as the true value
($\dot{P_{b}}^{\mathrm{X}}=0$).

  * f

Not claimed value. This is derived from DDH solution but there is a
discrepancy with DDGR solution.

##  V Discussion

In the observable range of the cosmic drag, there are ten binary pulsars at
least. At the present observations, four binaries of them have not reached the
precision of interest. Four of the others (six) are almost consistent with the
cosmic drag prediction. From their results, anomalies of the conventional (no
expansion) model appear in the forbidden region with a total significance of
about $3.6\,\sigma$. It is worth considering although the statistical
significance has not reached a criteria of the “discovery ($5\,\sigma$)”.

As a possible reason of $\dot{P_{b}}^{\mathrm{X}}$, a distance estimate by the
dispersion measure is doubted in PSR B1534+12 B1534+12 ; Stairs2002 . While,
the distance uncertainty cannot be a reason of $\dot{P_{b}}^{\mathrm{X}}$ in
PSR J1012+5307 Ding_2020  because this observation employs a parallax-based
distance with less uncertainty.

PSR J1738–0333 and PSR J1756–2251 have a discrepancy with the cosmic drag
picture. Their results of $\dot{P_{b}}^{\mathrm{X}}<\dot{P_{b}}^{\mathrm{CD}}$
(faster decay) are possible when uncounted energy losses contribute to
$\dot{P_{b}}^{\mathrm{obs}}$ (i.e. not GW dominant). If the cosmic drag
picture is correct, J1738–0333 and J1756–2251 may be such situation or any
systematic issue. It would be likely to reach $\dot{P_{b}}^{\mathrm{CD}}$
observable by refining parameters as well as the case of J1012+5307 Ding_2020
.

There is another interpretation of $\dot{P_{b}}^{\mathrm{X}}$ as alternative
theory of gravity to change $G$ with time and/or gravitational dipole
radiation (e.g. tested in Lazaridis_2009 ; Ding_2020 ). The main difference of
the cosmic drag from them is the limited observable range around
$r_{\mathrm{c}}$ as shown in Fig. 1.

##  VI Acknowledgements

I am grateful to Dr. Alberto Vecchio and Dr. Andreas Freise for their great
support to continue this study. I thank Dr. Kazuhiro Yamamoto for a useful
discussion. I appreciate Dr. Paulo Freire and his team for confirming
measurement values in J1829+2456.

## References

  * (1) E. Hubble, A relation between distance and radial velocity among extra-galactic nebulae, Proc. Natl. Acad. Sci. 15 (1929) 168–173.
  * (2) G. Lemaître, Un Univers homogène de masse constante et de rayon croissant rendant compte de la vitesse radiale des nébuleuses extra-galactiques, Ann. La Société Sci. Bruxelles. 47 (1927) 49–59.
  * (3) A.G. Riess, et. al., Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant, Astron. J. 116 (1998) 1009–1038.
  * (4) S. Perlmutter, et. al. [T.S.C. Project], Measurements of $\Omega$ and $\Lambda$ from 42 High‐Redshift Supernovae, Astrophys. J. 517 (1999) 565–586.
  * (5) W.L. Freedman, et. al., Final Results from the Hubble Space Telescope Key Project to Measure the Hubble Constant, Astrophys. J. 553 (2001) 47–72.
  * (6) J.L. Anderson, Multiparticle dynamics in an expanding universe, Phys. Rev. Lett. 75 (1995) 3602–3604.
  * (7) F.I. Cooperstock, V. Faraoni, D.N. Vollick, The Influence of the Cosmological Expansion on Local Systems, Astrophys. J. 503 (1998) 61–66.
  * (8) W.B. Bonnor, Size of a hydrogen atom in the expanding universe, Class. Quantum Gravity. 16 (1999) 1313–1321.
  * (9) G.S. Adkins, J. McDonnell, R.N. Fell, Cosmological perturbations on local systems, Phys. Rev. D - Part. Fields, Gravit. Cosmol. 75 (2007) 64011.
  * (10) M. Carrera, D. Giulini, Influence of global cosmological expansion on local dynamics and kinematics, Rev. Mod. Phys. 82 (2010) 169–208.
  * (11) R.H. Price, J.D. Romano, In an expanding universe, what doesn’t expand?, Am. J. Phys. 80 (2012) 376–381.
  * (12) D. Giulini, Does cosmological expansion affect local physics?, Stud. Hist. Philos. Sci. Part B - Stud. Hist. Philos. Mod. Phys. 46 (2014) 24–37.
  * (13) R. Nandra, A.N. Lasenby, M.P. Hobson, The effect of a massive object on an expanding universe, Mon. Not. R. Astron. Soc. 422 (2012) 2931–2944.
  * (14) R. Nandra, A.N. Lasenby, M.P. Hobson, The effect of an expanding universe on massive objects, Mon. Not. R. Astron. Soc. 422 (2012) 2945–2959.
  * (15) K. Agatsuma, The expansion of the universe in binary star systems, Physics of the Dark Universe. 30 (2020) 100732.
  * (16) P.C. Peters, J. Mathews, Gravitational radiation from point masses in a keplerian orbit, Phys. Rev. 131 (1963) 435–440.
  * (17) P.C. Peters, Gravitational radiation and the motion of two point masses, Phys. Rev. 136 (1964) B1224–B1232.
  * (18) I.H. Stairs, S.E. Thorsett, J.H. Taylor, A. Wolszczan, Studies of the Relativistic Binary Pulsar PSR B1534+12. I. Timing Analysis, Astrophys. J. 581 (2002) 501–508.
  * (19) E. Fonseca, I.H. Stairs, S.E. Thorsett, A comprehensive study of relativistic gravity using PSR B1534+12, Astrophys. J. 787 (2014) 82.
  * (20) J.M. Weisberg, D.J. Nice, J.H. Taylor, Timing measurements of the relativistic binary pulsar PSR B1913+16, Astrophys. J. 722 (2010) 1030-1034.
  * (21) J.M. Weisberg, Y. Huang, Relativistic Measurements From Timing the Binary Pulsar Psr B1913+16, Astrophys. J. 829 (2016) 55.
  * (22) H. Ding, A.T. Deller, P. Freire, D.L. Kaplan, T.J. W. Lazio, R. Shannon, B. Stappers, Very Long Baseline Astrometry of PSR J1012+5307 and its Implications on Alternative Theories of Gravity, Astrophys. J. 896 (2020) 85.
  * (23) J. Van Leeuwen, L. Kasian, I.H. Stairs, D.R. Lorimer, F. Camilo, S. Chatterjee, I. Cognard, THE BINARY COMPANION OF YOUNG , RELATIVISTIC PULSAR J1906 + 0746, 118 (2015).
  * (24) J. Antoniadis, T.M. Tauris, F. Ozel, E. Barr, D.J. Champion, P.C.C. Freire, The millisecond pulsar mass distribution: Evidence for bimodality and constraints on the maximum neutron star mass, (2016). http://arxiv.org/abs/1605.01665.
  * (25) D.M. Sánchez, A.G. Istrate, M.H. Van Kerkwijk, R.P. Breton, D.L. Kaplan, PSR J1012+5307: A millisecond pulsar with an extremely low-mass white dwarf companion, Mon. Not. R. Astron. Soc. 494 (2020) 4031–4042.
  * (26) R.S. Lynch, J.K. Swiggum, V.I. Kondratiev, D.L. Kaplan, K. Stovall, E. Fonseca, M.S.E. Roberts, L. Levin, M.E. DeCesar, B. Cui, S.B. Cenko, P. Gatkine, A.M. Archibald, S. Banaszak, C.M. Biwer, J. Boyles, P. Chawla, L.P. Dartez, D. Day, A.J. Ford, J. Flanigan, J.W.T. Hessels, J. Hinojosa, F.A. Jenet, C. Karako-Argaman, V.M. Kaspi, S. Leake, G. Lunsford, J.G. Martinez, A. Mata, M.A. McLaughlin, H. Al Noori, S.M. Ransom, M.D. Rohr, X. Siemens, R. Spiewak, I.H. Stairs, J. van Leeuwen, A.N. Walker, B.L. Wells, The Green Bank North Celestial Cap Pulsar Survey. III. 45 New Pulsar Timing Solutions, Astrophys. J. 859 (2018) 93.
  * (27) G. Desvignes, R.N. Caballero, L. Lentati, J.P.W. Verbiest, D.J. Champion, B.W. Stappers, G.H. Janssen, P. Lazarus, S. Osłowski, S. Babak, C.G. Bassa, P. Brem, M. Burgay, I. Cognard, J.R. Gair, E. Graikou, L. Guillemot, J.W.T. Hessels, A. Jessner, C. Jordan, R. Karuppusamy, M. Kramer, A. Lassus, K. Lazaridis, K.J. Lee, K. Liu, A.G. Lyne, J. McKee, C.M.F. Mingarelli, D. Perrodin, A. Petiteau, A. Possenti, M.B. Purver, P.A. Rosado, S. Sanidas, A. Sesana, G. Shaifullah, R. Smits, S.R. Taylor, G. Theureau, C. Tiburzi, R. Van Haasteren, A. Vecchio, High-precision timing of 42 millisecond pulsars with the European Pulsar Timing Array, Mon. Not. R. Astron. Soc. 458 (2016) 3341–3380.
  * (28) P.C.C. Freire, N. Wex, G. Esposito-Farèse, J.P.W. Verbiest, M. Bailes, B.A. Jacoby, M. Kramer, I.H. Stairs, J. Antoniadis, G.H. Janssen, The relativistic pulsar-white dwarf binary PSR J1738+0333 - II. The most stringent test of scalar-tensor gravity, Mon. Not. R. Astron. Soc. 423 (2012) 3328–3343.
  * (29) R.D. Ferdman, I.H. Stairs, M. Kramer, G.H. Janssen, C.G. Bassa, B.W. Stappers, P.B. Demorest, I. Cognard, G. Desvignes, G. Theureau, M. Burgay, A.G. Lyne, R.N. Manchester, A. Possenti, PSR J1756-2251: A pulsar with a low-mass neutron star companion, Mon. Not. R. Astron. Soc. 443 (2014) 2183–2196.
  * (30) H.T. Haniewicz, R.D. Ferdman, P.C.C. Freire, D.J. Champion, K.A. Bunting, D.R. Lorimer, M.A. McLaughlin, Precise mass measurements for the double neutron star system J1829+2456, Mon. Not. R. Astron. Soc. 500 (2021) 4620–4627.
  * (31) B.A. Jacoby, P.B. Cameron, F.A. Jenet, S.B. Anderson, R.N. Murty, S.R. Kulkarni, Measurement of Orbital Decay in the Double Neutron Star Binary PSR B2127+11C, Astrophys. J. 644 (2006) L113–L116.
  * (32) K. Lazaridis, N. Wex, A. Jessner, M. Kramer, B.W. Stappers, G.H. Janssen, G. Desvignes, M.B. Purver, I. Cognard, G. Theureau, A.G. Lyne, C.A. Jordan, J.A. Zensus, Generic tests of the existence of the gravitational dipole radiation and the variation of the gravitational constant, Mon. Not. R. Astron. Soc. 400 (2009) 805–814.


# Multi-wavelength Variation Phenomena of PKS 0735+178 on Diverse Timescale

Yue Fang  Qihang Chen  Yan Zhang  Jianghua Wu† Department of Astronomy,
Beijing Normal University, 100875, Beijing, China Jianghua Wu<EMAIL_ADDRESS>
###### Abstract

BL Lac object PKS 0735+178 showed some complex multi-wavelength variation
phenomena in the previous studies, especially for its color behavior. Bluer-
when-brighter, redder-when-brighter and achromatic behavior were all found to
be possible long-term trends of PKS 0735+178. In this work, we collected the
long-term multi-wavelength data of PKS 0735+178, and also performed a multi-
color optical monitoring on intraday timescale. The intraday variability was
detected on one night. On the long timescale, a possible 22-day time lag was
found between the $R$ and $\gamma$-ray bands. The results of the cross-
correlation analysis exhibited strong correlations between various optical
bands on both intraday and long timescale. However, only a mild correlation
was found between the long-term $\gamma$-ray and $R$-band light curves, which
could be interpreted by the different emission mechanisms of the $\gamma$-ray
and optical emissions. PKS 0735+178 showed a significant harder-when-brighter
in the $\gamma$-ray band, which is consistent with the observed optical bluer-
when-brighter trend on both long-term and intraday timescales. We found that
the HWB and BWB trends will be enhanced during the active states, especially
for the historical low state. Such phenomenon indicates a special activity-
dependent color behavior of PKS 0735+178, and it could be well interpreted by
the jet emission model.

BL Lacertae objects (158), Active galaxies (17), Galaxy photometry (611)

††software: IRAF (Tody, 1986, 1993), BCES (Akritas & Bershady, 1996; Nemmen et
al., 2012), PyCCF (Peterson et al., 1998b; Sun et al., 2018), JAVELIN (Zu et
al., 2011).

##  1 INTRODUCTION

Blazars are the most violently variable class of active galactic nuclei (AGN)
(Urry & Padovani, 1995), which have their relativistic jets directed towards
the earth. They display violent variability of flux and polarization of non-
thermal radiation across the entire electromagnetic spectrum. According to
whether there are some broad emission lines, blazars can be further divided
into flat-spectrum radio quasars (FSRQs) and BL Lac objects. The broadband
spectral energy distributions (SEDs) of blazars exhibit two characteristic
bumps. One locates at low frequencies (from radio to the UV or X-ray bands)
and is dominated by synchrotron radiation. The other locates at high
frequencies (the X-ray and $\gamma$-ray bands) and is interpreted as the
inverse Compton scattering processes (e.g., Ulrich et al., 1997; Böttcher,
2007; Sikora et al., 2009). The color behavior can serve as a tool to
investigate the nature of emission of blazars and to probe the physical
processes in them. The analysis of color-magnitude trends in the monitoring
data of blazars has revealed three common patterns, i.e., bluer-when-brighter
(BWB), redder-when-brighter (RWB), and achromatic. Through some comprehensive
statistical analysis, BL Lac objects frequently show BWB trend, while FSRQs
tend to show RWB trend (e.g., Gu et al., 2006; Wu et al., 2011; Gaur et al.,
2012; Isler et al., 2017).

PKS 0735+178 (also known as OI 158, S3 0735+17, 4FGL J0738.1+1742) was
identified as a BL Lac object (Carswell et al., 1974). Its J2000 coordinates
are RA= 07h38m07.3937s, DEC = +17d42m18.998s, and the redshift is $z=0.45$
(Rector & Stocke, 2001). In 1999, Hartman et al. (1999) reported $\gamma$-ray
detection of this blazar using Energetic Gamma Ray Experiment Telescope, and
it was found to be quite steady in X-ray and $\gamma$-ray bands with the
continuous monitoring (Bregman et al., 1984; Nolan et al., 2003). On the
contrary, in the optical and radio bands, PKS 0735+178 showed violent long-
term variability in the past decades (Webb et al., 1988; Ciprini et al.,
2007). Also, some intraday variabilities were reported on its optical history
(Sagar et al., 2004; Goyal et al., 2009). Through a century long optical light
curve, Fan et al. (1997) reported a possible period of $\sim 14$ years.
Ciprini et al. (2007) found 3 characteristic timescales of about 4.5, 8.5, and
11–13 years.

A bronze neutrino event, IC211208A, with a $>30$% probability of being of
astrophysical origin, was localised to the vicinity of PKS 0735+178 (2.1 deg
separation, GCN 31191). Follow-up observations showed that it flared in radio
(ATel $\\#$15105), optical (ATel $\\#$15098), X-rays (ATel $\\#$15102) and
$\gamma$-rays (ATel $\\#$15099) nearly simultaneously, with a probability of
$>30\%$ being the astrophysical origin of this neutrino event. During the
IC211208A event, PKS 0735+178 reached its brightest state in the $R$ band
($\sim 14.2$
mag111http://herculesii.astro.berkeley.edu/kait/agn/lightcurve_all/lightcurve_J0738+1742_all_psf_natural_group_kait.png,
ATel $\\#$15021), and showed a peak $\gamma$-ray daily flux value of
$5^{+2}_{-2}\times 10^{-7}\ ph\ cm^{-2}\ s^{-1}$, which was 10 times greater
than its average 4FGL-DR2 flux (GCN 31194). Sahakyan et al. (2022) presented a
multi-wavelength analysis of PKS 0735+178, for which a radio context was also
provided by high angular resolution imaging (Nanci et al., 2022; Weaver et
al., 2022).

Specially, it seems that PKS 0735+178 didn’t follow a simple BWB or RWB trend,
and researchers have been arguing about its color behavior for more than two
decades. As early as in 2000, Fan & Lin (2000) collected 10-year historical
optical data of PKS 0735+178 and found that it showed a BWB trend, and the
studies of Rani et al. (2010) and Meng et al. (2018) also showed this trend.
However, Sandrinelli et al. (2014) found a positive correlation between the
$R-H$ color vs. the $H$-band magnitude, with a general trend indicating bluer
color for decreasing flux, i.e., a RWB trend. In the $\gamma$-ray band, a
comprehensive study of blazars (including PKS 0735+178) illustrated a
significant harder-when-brighter (HWB) trend for FSRQs, while the BL Lacs
showed no propensity toward a HWB or softer-when-brighter trend (SWB)
(Williamson et al., 2014). Achromatic behavior or weak correlation in optical
bands also appeared in the previous studies on PKS 0735+178 (Gu et al., 2006;
Ciprini et al., 2007). Recently, Yuan & Fan (2021) found that PKS 0735+178
showed two opposite spectral behaviors among different bands, i.e., BWB for
spectral indices (SIs) vs. $g$-band flux but RWB for SIs vs. $r$\- and
$i$-band flux. All these results suggested a complex pattern of its color
behavior.

In order to explore the specific color/spectral behaviors of PKS 0735+178, we
collected its long-term optical and $\gamma$-ray data, and we also carried out
a multi-color optical monitoring program with the 85cm telescope at Xinglong
observatory. The descriptions of observations and data reductions could be
found in Section 2. In Section 3, we studied the intraday variability (IDV) of
this object and searched for the possible inter-band time lags through the
interpolated cross correlation function (ICCF). In Section 4, we analyzed its
color/spectral behavior in detail, including the behaviors under different
timescales and different activity states, to figure out the primary variation
mechanism. Finally, a summary is given in Section 5.

##  2 OBSERVATION AND DATA REDUCTION

###  2.1 Optical: Intraday

The intraday monitoring was performed with the 85cm telescope at Xinglong
Station of the National Astronomical Observatories Chinese Academy of Sciences
(NAOC). This telescope uses the prime focus optical design with a focal ratio
$f/3.3$. The CCD was a 2048 × 2048 chip with a field of view of $\sim
32.8\times 32.8$ $arcmin^{2}$. We observed PKS 0735+178 for 3 nights in Dec.
2020, and $\sim$400 data points were collected. We have listed the relevant
observation information in Table 1.

Table 1: Intraday Results of PKS 0735+178

Julian Date | Date | Filter | No. Of Expo. | Enhanced $F$-test | ANOVA test | Var? | Amp
---|---|---|---|---|---|---|---
(MJD) | (ISO) |  |  | $\nu_{1}$ | $\nu_{2}$ | $F$ | $F_{c}$ | $\nu_{1}$ | $\nu_{2}$ | $F$ | $F_{star}$ | $F_{c}$ |  | (%)
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11) | (12) | (13) | (14) | (15)
59188 | 2020-12-04 | B | 69 | 68 | 340 | 1.57 | 1.51 | 12 | 56 | 7.87 | 0.33 | 2.52 | V | $10.6\pm 4.7$
|  | V | 69 | 68 | 340 | 1.80 | 1.51 | 12 | 56 | 6.09 | 0.38 | 2.52 | V | $10.1\pm 4.8$
|  | R | 69 | 68 | 340 | 1.99 | 1.51 | 12 | 56 | 5.66 | 0.74 | 2.52 | V | $10.1\pm 3.8$
59190 | 2020-12-06 | B | 29 | 28 | 140 | 1.42 | 1.87 | 4 | 24 | 4.07 | 1.21 | 4.22 | N |
|  | V | 29 | 28 | 140 | 0.43 | 1.87 | 4 | 24 | 0.79 | 0.64 | 4.22 | N |
|  | R | 29 | 28 | 140 | 0.35 | 1.87 | 4 | 24 | 1.41 | 1.33 | 4.22 | N |
59191 | 2020-12-07 | B | 32 | 31 | 155 | 0.11 | 1.81 | 5 | 26 | 0.70 | 0.58 | 3.82 | N |
|  | V | 32 | 31 | 155 | 0.16 | 1.81 | 5 | 26 | 1.32 | 0.60 | 3.82 | N |
|  | R | 32 | 31 | 155 | 0.23 | 1.81 | 5 | 26 | 0.38 | 0.36 | 3.82 | N |

Note. — The columns are (1) observational date (MJD), (2) observational date
(International Organization for Standardization, ISO), (3) filter, (4) number
of exposures, (5) the degree of freedom within group in the enhanced $F$-test
(the number of values within each group that are free to vary, and it is
calculated by $\nu_{1}=N_{b}-1$, where $N_{b}$ is the number of data points in
the blazar’s light curve), (6) the degree of freedom between groups in the
enhanced $F$-test (the number of values between all the groups that are free
to vary, and it is calculated by $\nu_{2}=(\sum_{j=1}^{k}N_{j})-k$, where
$N_{j}$ is the number of data points in the $j$th comparison star’s light
curve and $k$ is the number of comparison stars), (7-8) $F$ and the critical
value $F_{c}$ in the enhanced $F$-test, (9-10) two degrees of freedom in the
ANOVA test (their meanings are the same as those in the enhanced $F$-test, but
calculated by $\nu_{1}=g-1$ and $\nu_{2}=N_{b}-g$, respectively, where $g$ is
the number of groups), (11-13) $F$ of PKS 0735+178, $F$ of star C and $F_{c}$
in the ANOVA test, (14) variable or not, (15) variability amplitude and
associated uncertainty, respectively.

As shown in Fig. 1, PKS 0735+178 and six field stars are labeled. The data
reduction followed the standard process including bias-subtraction, flat-
fielding and instrument magnitudes extraction with the Image Reduction and
Analysis Facility (IRAF)222IRAF is distributed by the National Optical
Astronomy Observatories, which are operated by the Association of Universities
for Research in Astronomy, Inc., under cooperative agreement with the National
Science Foundation. software. Through aperture photometry, we extracted the
instrument magnitudes of PKS 0735+178 and six field stars. In order to obtain
the best aperture, the photometry was carried out with ten different aperture
radii range from 1 to 3 times of the full width at half-maximum (FWHM), and
the inner and outer radii of the sky annulus were accordingly set as 6 and 8
times of the FWHMs, respectively. We finally selected 1.6$\times$FWHM, which
gave minimum standard deviations of the differential magnitudes between the
check star and the comparison star. In order to minimize the error of
calibration, star 1 was chosen as the comparison star, which is the brightest
one of all six. Moreover, the closest star 2 was chosen as the check star. The
fluxes of PKS 0735+178 and check star were calibrated with respect to a
comparison star. The calibrated magnitudes of our target and star 2 were
presented in Fig. 2. For clarity, the $B$\- and $V$-band light curves are
shifted by $-0.3$ mag and $-0.4$ mag, respectively. Similar shifts of $-1.1$
mag and $-0.5$ mag were also applied respectively to the $B$\- and $R$-band
light curves in the small panels.

Figure 1: Finding chart of PKS 0735+178 in the $R$ band. The labels “T” and
“1-6” represent PKS 0735+178, the comparison star, the check star, and four
field stars, respectively.

###  2.2 Optical: Long-term

We collected the historical optical data from various archives, including
those of the Steward Observatory (Smith et al.,
2009)333http://james.as.arizona.edu/~psmith/Fermi, the Small and Moderate
Aperture Research Telescope System (SMARTS)444www.astro.yale.edu/smarts, the
Catalina Real-Time Transient Survey (CRTS) (Drake et al.,
2009)555http://crts.caltech.edu/, the Katzman Automatic Imaging Telescope
(KAIT)666http://herculesii.astro.berkeley.edu/kait/agn/ and the All-Sky
Automated Survey for Supernovae (Jayasinghe et al.,
2019)777http://www.astronomy.ohio-state.edu/asassn. Notice that the data from
SMARTS are actually raw images after bias subtracted and flat field
calibrated, which can be queried from NOIRLab Astro Data
Archive888https://astroarchive.noao.edu, and the subsequent processing is the
same as in Sec. 2.1. The collected data are in the $V$-, $R$-, and $I$-bands,
but we will not consider the $I$-band due to the low sampling. If more than
one data point is present in a night, the average magnitudes are taken.

###  2.3 $\gamma$-ray

The Fermi large area telescope (LAT) data of PKS 0735+178 from 2009 January 1
to 2019 December 31 were downloaded, and a circular region of interest (ROI)
of $20^{\circ}$ was chosen. The energy range was restricted to 0.1-200 GeV. We
performed the unbinned likelihood analysis using version v2.0.8 of the Fermi
Science Tools, the background models from the iso_P8R3_SOURCE_V3_v1.txt
isotropic template, and the gll_iem_v07.fits Galactic diffuse emission model.
We utilized analysis cuts of evclass = 128, evtype = 3, and zmax=90 of the
photon data. We generated the XML model files using make4FGLxml.py, and all
the spectral indices and fluxes of the sources were fixed, other than the
target, at their 4FGL values. We considered one measurement to be a successful
detection when the test statistics (TS) value exceeded 10, which corresponds
to a $3\sigma$ criteria. The light curve was generated with a 15-days bin, and
we used a power law (PL) function, i.e.,
$dN(E)/dE=N_{0}(E/E_{0})^{-\Gamma_{\gamma}}$, to fit the $\gamma$-ray fluxes
in eight logarithmically equal energy bins on each time bin. The variation
behavior of spectra will be studied by pairing the simultaneous flux and SI
$\Gamma_{\gamma}$. The long-term light curves in all three bands were
presented in Fig. 3. Moreover, the measurements with TS $<10$ were labeled by
“upper limits” (red triangles).

##  3 Variations

###  3.1 Light Curves

From visual inspection of Figure 2 and 3, we can find that the optical light
curves are well correlated with each other on both intraday and long
timescales, and a significant IDV is exhibited on MJD 59188. However, on the
long timescale, the correlation between the optical band and the $\gamma$-ray
band is clearly weaker than that between two optical bands. A prominent flare
was observed during MJD 56500-56800 in the $V$ and $R$ bands. A historical low
state (17.23 mags) was observed in the $V$ band during MJD 57350-57450, which
was $\sim$0.2 mag fainter than the previous lowest state ($\sim$17 mags during
MJD 50400-50600, see Ciprini et al. (2007)). During MJD 59188-59191, the
source was the faintest on MJD 59190 with 15.95 mags in the $R$ band and
reached the brightest state on MJD 59188 with 15.78 mags. During our long-term
optical observations, PKS 0735+178 reached the brightest state on MJD 56649 in
the $R$ band, which was $\sim 0.3$ magnitude fainter than its historical
brightest state on MJD 59551, i.e., six days before the neutrino event
IC211208A. We noticed that all the optical flares over the past decades have
far lower intensity than the most recent flare possibly associated with the
neutrino event IC211208A. It brightened by $\sim$1 mag in 50 days since MJD
59500, and then dimmed at nearly the same rate (GCN 31529), which was about
three times the rate of that flare during MJD 56500-56800. In the $\gamma$-ray
band, the object is 44 times brighter in the highest state than in the lowest
state.

Figure 2: The intraday light curves of PKS 0735+178 and the check star are
given in the large and small panels, respectively. For clarity, the $B$\- and
$V$-band light curves of PKS 0735+178 are shifted by $-0.4$ and $-0.3$ mag,
respectively, while those of the check star are shifted by $-1.1$ and $-0.5$
mag, respectively.  Figure 3: Top panel: the light curve of $\gamma$-ray in
the range of 0.1-200 GeV, and the red triangles represent the upper limits of
the measurements with TS $<10$. The dashed vertical lines divide the light
curve into five phases according to the activities (see Section 3.1 for
detail). The second and third panels: the $V$\- and $R$-band light curves,
which were segmented according to the phases in the $\gamma$-ray band. The
fourth panel and bottom panel are the SI of the $\gamma$-ray and $V-R$ CI,
respectively. Red dots represent those CIs obtained with the original
magnitude pairs without damped random walk interpolation.

In order to determine the variability on intraday timescale, we used two
robust statistical tests, namely, the power-enhanced version of $F$-test (de
Diego, 2014) and analysis of variance (ANOVA) (de Diego et al., 1998), which
were widely used in the past IDV studies (e.g. Gaur et al., 2015; Polednikova
et al., 2016; Meng et al., 2017; Zhang et al., 2018; Pandey et al., 2020). The
detailed algorithms can be found in the corresponding references. If the $F$
exceeds the critical value $F_{c}$ at the $99\%$ confidence level, the null
hypothesis will be rejected and the existence of variability (V) will be
confirmed, otherwise it is non-variable (N). To avoid detecting spurious
variability, an additional ANOVA test was performed to the check star. Only
both tests of PKS 0735+178 show “V” and the test of check star is “N”, PKS
0735+178 is genuinely variable (“V”), otherwise it is non-variable (“N”). The
final results were listed in Table. 1, including two corresponding degrees of
freedom, the $F$ value and the critical value $F_{c}$. PKS 0735+178 only
showed IDV on MJD 59188. We found that the observing duration on MJD 59188 was
the longest among three nights, and most of the variabilities on this night
occurred in the second half night. Longer sampling would naturally have higher
probability to detect IDV.

We also calculated the variation amplitudes on this night, by

$Amp=100\%\times\sqrt{(A_{max}-A_{min})^{2}-2\sigma^{2}},$ (1)

where $A_{max}$ and $A_{min}$ are the maximum and minimum magnitudes,
respectively, and $\sigma$ is the measurement error (Heidt & Wagner, 1996).
The propagating uncertainty of $Amp$ is given by
$\sigma_{Amp}=(\frac{A_{max}-A_{min}}{Amp})\cdot\sqrt{(\sigma_{A_{max}}^{2}+\sigma_{A_{min}}^{2})}$,
where the $\sigma_{A_{max}}$ and $\sigma_{A_{min}}$ are the measurement errors
corresponding to $A_{max}$ and $A_{min}$, respectively. The amplitudes and
associated uncertainties were also given in Table 1. We can find that the
variation amplitude tends to decrease with decreasing frequency, which is
consistent with the common tendency, i.e., the higher the energy band, the
greater the IDV amplitude (Webb et al., 1998; Meng et al., 2017).

On the long-term trend, three distinct $\gamma$-ray flares were observed
during MJD 55950-56420, 56420-56780 and 57840-58040, respectively, and the
second flare was observed quasi-simultaneously in the $V$ and $R$ bands.
According to the activities, we divided the $\gamma$-ray light curve into five
phases, namely, phase $A$ (the 1st quiescent, and from MJD 54839 to 55720),
phase $B$ (the 1st active state, and from MJD 55720 to 57280), phase $C$ (the
2nd quiescent, and from MJD 57280 to 57850), phase $D$ (the 2nd active state,
and from MJD 57850 to 58230) and phase $E$ (the 3rd quiescent, and from MJD
58230 to the end), which were distinguished by vertical dashed lines in Fig.
3. The optical light curves were segmented according to the $\gamma$-ray
phase.

Figure 4: The rows from top to bottom give the lag measurements on MJD
54839-58821, MJD 59188 and MJD 59190, respectively. The measurements on MJD
59191 were excluded due to the too small $r_{max}$ (see Table 2). The solid
and dotted lines indicate the measured lags and uncertainties, respectively.
The corresponding lag calculation results are listed in each panel. The dashed
lines indicate the smoothing of the weighted CCCDs and the histograms are the
unweighted CCCDs. The gray regions highlight the ranges of lags included in
the final calculation, and the removal efficiency $f_{out}$ were also given.

###  3.2 Cross-Correlation Analysis

We performed the correlation analysis to search for the possible inter-band
time lags by using ICCF999The public available algorithm PyCCF can be found in
http://ascl.net/code/v/1868 (Peterson et al., 1998a, b). Through shift and
linear interpolation, the Pearson coefficient $r$ between the two light curves
will be calculated at each given time shift $\tau$. The optimal time lag and
its uncertainties were obtained via the flux randomization/random subset
sampling (FR/RSS) method, using Monte Carlo (MC) realizations, as discussed in
Peterson et al. (2004). Those realisation with a peak correlation coefficient
$r_{peak}<0.5$ will be labeled as a failed measurement, and the failure rate
$f_{fail}$ will be recorded. The centroid lag will be estimated using points
surrounding the $r_{peak}$, and we can obtain the cross correlation centroid
distribution (CCCD) after 10000 MC realizations. The search range of time lag
was empirically fixed as $[-200,200]$ days and $[-60,60]$ minutes for long-
term and intraday timescales, respectively.

However, the alias of multiple significant peaks has been reported in the
previous studies (Grier et al., 2017; Zhang et al., 2018; Homayouni et al.,
2019; Li et al., 2019). Thus we followed Grier et al. (2017) and introduced
the weighted lag identification. The weight of each data point in CCCD were
defined as $P=[N(\tau)/N(0)]^{2}$, where $N(\tau)$ corresponds to the number
of overlapping epochs at a time shift $\tau$. After a Gaussian kernel
estimation (KDE), we can obtain the peak with the largest area in CCCD, and
the data points beyond the range of this primary peak will be removed. The
final lag was estimated as the median of the unweighted & unsmoothed CCCD
within the range of the primary peak, and the 16th and 84th percentiles were
chosen as the lower and upper uncertainties, respectively. Empirically, a
reliable result should satisfy $f_{fail}\leq 30\%$, $f_{out}\leq 50\%$ (the
ratio of the removed data points to the total data points) and $r_{max}\geq
0.5$. All the measurement results that satisfied these criteria were presented
in Fig. 4, and corresponding $f_{fail}$, $f_{out}$, the maximum correlation
coefficient $r_{max}$ and time lags were listed in Table 2. The negative lags
indicate that the former band leads the latter one.

On the intraday timescale, all the time lags were near-zero within the error
range. A similar result was obtained between the $V$ and $R$ bands on the long
timescale. The non-detection of the lag between variations in different
optical bands may be because the emitting regions of these optical bands are
too close to each other to result in the time lag (Carini et al., 2011; Wu et
al., 2012; Meng et al., 2017). Moreover, according to the data compiled from
Sandrinelli et al. (2014), a near-zero lag ($-3.6^{+4.2}_{-4.2}$ days, see
Appendix B) was also found between the $R$ and $J$ bands. The only one non-
zero result was the 22-day lag between the $R$ and $\gamma$-ray bands. On the
long timescale, the $r_{max}$ between the $\gamma$-ray and $R$ bands is only
0.5, which is much less than the $r_{max}=0.82$ between two optical bands.
Several studies also found that the correlation between the optical and
$\gamma$-ray bands was much weaker than that between various optical bands
(e.g., Chatterjee et al., 2013; Cohen et al., 2014; Liodakis et al., 2019;
Rajput et al., 2020). This is mainly due to the different origins of the
optical and $\gamma$-ray emissions, the former is dominated by the synchrotron
radiation while the latter is by the inverse Compton process. Alternatively,
Rajput et al. (2020) demonstrated that the changes of the magnetic field
strength could lead to this phenomenon.

Table 2: Results of cross-correlation analysis.

Date | Bands | $f_{fail}$ | $f_{out}$ | $r_{max}$ | Lag
---|---|---|---|---|---
|  | ($\%$) | ($\%$) |  |
54839-58821 |  $\gamma$-ray vs. $R$ | 0.5 | 18.7 | 0.52 |  $-22.1^{+18.3}_{-17.1}$ d
|  $V$ vs. $R$ | 0.2 | 0 | 0.82 |  $-1.6^{+7.4}_{-8.9}$ d
59188 |  $V$ vs. $R$ | 5.6 | 0 | 0.86 |  $-10.2^{+17.5}_{-17.5}$ min
|  $B$ vs. $R$ | 2.8 | 0 | 0.90 |  $2.6^{+15.1}_{-14.4}$ min
59190 |  $V$ vs. $R$ | 5.5 | 7.7 | 0.62 |  $1.3^{+13.0}_{-12.7}$ min
|  $B$ vs. $R$ | 1.8 | 1.9 | 0.84 |  $-2.1^{+13.0}_{-14.5}$ min
59191 |  $V$ vs. $R$ | 27.5 | 44.9 | 0.35 | /
|  $B$ vs. $R$ | 17.5 | 50.2 | 0.18 | /

Note. — The columns stand for observational date in the unit of MJD, bands,
failure rate, removal efficiency, the maximum correlation coefficient and the
measured lag through PyCCF. The negative lags indicate that the former band
leads the latter one.

##  4 Color/Spectral Behavior

###  4.1 Observational Findings

In principle, the color index (CI) should be obtained by pairing the
simultaneous data points. Hence, our quasi-simultaneous intraday light curves
were linearly interpolated to get the intraday CIs. The intraday color
behaviors on MJD 59188 were shown in Fig. 5. In order to comprehensively
consider the errors in both magnitudes and CIs, we introduced the BCES
estimator (Akritas & Bershady, 1996). Only when the absolute value of the
correlation coefficient exceeded 0.2 and the statistical significance of the
correlation exceeded 99$\%$ (i.e., $p<0.01$), the target was considered to
show a confident color-magnitude correlation. It is clear that all the
diagrams in Fig. 5 show significant BWB trends.

Figure 5: Intraday color-magnitude diagrams of PKS 0735+178 on MJD 59188.

However, since our long-term optical data were collected from multiple
datasets, we can only pair two data points on the same night. Even so, all the
obtained CIs clustered in the period of MJD 55000-55400. Thus, we interpolated
the long-term light curves with the damped random walk (DRW) method, which was
proven to model quasar light curve behavior quite well on time scales of
months to years (Gaskell & Peterson, 1987; MacLeod et al., 2010; Zu et al.,
2013). More details can be found in Appendix A. We interpolated the $R$-band
light curve to obtain the simultaneous data points with those of the $V$ band
and calculated the CIs. For comparison, the original CIs without interpolation
were also calculated. The variation of CIs was presented in Fig. 3. The long-
term color trends of PKS 0735+178 were presented in Fig. 6, including the
overall trend and four individual trends101010Due to the small sample of data,
we did not give the diagram of phase $E$..

The BWB behaviors were exhibited on the overall phase, phases $A$, $B$ and
$C$. Such behaviors were also proven by the color-magnitude correlation
without interpolation (red dots and dashed line in the top-left panel). Our
results are consistent with the common color trend of BL Lac objects (Fan &
Lin, 2000). We noticed that the BWB trend exhibited in phase $B$, which
contained a significant flare component, was stronger than that of phase $A$.
Several studies have demonstrated that the BWB trend of some blazars will be
enhanced during active states, e.g., S5 0716+714 (Wu et al., 2007), PG
1553+113 (Agarwal et al., 2021) and AO 0235+164 (Wang & Jiang, 2020). As a
comparison, the phase $A$ that contained only one incomplete flare presented a
mild BWB trend, and an achromatic trend was exhibited in phase $D$ that lacked
a flare component.

The most complex color behavior was presented in phase $C$. From visual
inspection, two distinct components were found in the color-magnitude
correlation, i.e., a strong BWB trend and an achromatic one. To further
investigate the color behavior, we divided the variations in phase $C$ into
low and high states by the mean magnitude ($\sim 16.23$ mags) of this period,
as shown in Fig. 7. Particularly, the low state corresponds to the historical
low state mentioned above. PKS 0735+178 showed no trend in the high state and
a strong BWB trend in the low state. One can find that the correlation
coefficient of the low state reached 0.664, which means the most significant
BWB trend was observed in the historical low state.

Figure 6: Long-term optical color-magnitude correlations of PKS 0735+178. The
panels from top to bottom and left to right represent the overall phase and
four individual phases, respectively, and the corresponding formula,
correlation coefficient $r$ and $p$ value were given in the top-left corner of
each panel. We also fitted those data points without interpolation, which are
labeled by the red dots and dashed line in the top-left panel, and the
specific parameters were given by the red legends.  Figure 7: The fittings of
the low and high states in phase $C$ divided by the vertical dash-dotted line.
The solid and dashed lines represent the fittings of the low and high states,
respectively. The corresponding Spearman correlation coefficients are given in
the upper left corner.

The $\gamma$-ray spectral behaviors were presented in Fig. 8. In the past
decade, PKS 0735+178 showed a significant HWB trend, which is consistent to
the BWB trend in optical bands. Moreover, the overall HWB trend tends to
saturate toward higher flux (see the red dashed lines in the first panel).
This is a common phenomenon in the X-ray and $\gamma$-ray bands (e.g., Xue &
Cui, 2005; Giebels et al., 2007; Weng et al., 2020; Acciari et al., 2021),
which suggests a more efficient particle acceleration mechanism in the high
state (Giebels et al., 2007). As for the individual trends, two significant
HWB behaviors appeared in phases $A$ and $B$, and two low significant HWB
behavior appeared in phases $C$ and $D$. We noticed that the HWB trend of
phase $B$ with flare components was obviously stronger than those of phases
$A$, $C$ and $D$. Such phenomenon is consistent with the situations in the
optical bands. The flare in phase $D$ was the shortest one and most of the
data points were concentrated in the highest state, so the exhibited HWB trend
was not so significant. The spectral hardening during the flare state
indicates the increasing of detected high energy photons, which is expected if
the inverse Compton peak shifts to higher energy (Shah et al., 2019, 2021).
Moreover, a SWB trend was found with a low significance level in phase $E$.

Figure 8: Long-term $\gamma$-ray SI-flux correlations of PKS 0735+178. The
panels from top to bottom and from left to right represent the overall phase
and five individual phases, respectively. The change of the spectral behavior
from a HWB trend to a saturation in the overall phase is labeled by the red
dashed lines.

###  4.2 Interpretation and Discussions on Color Behavior

Since the lack of a radiative efficiently accretion disk, the color behavior
of BL Lac objects were always interpreted by the jet emission model (Marscher
& Gear, 1985). In the frame of the leptonic model of blazars, the electrons
are accelerated near the root of the relativistic jet and the magnetic field
is compressed therein. This process leads to the observed flux and spectral
variabilities (Marscher & Gear, 1985; Qian et al., 1991). Additionally, the
weak intraday BWB behavior also could be explained by the superposition of
many distinct new variable components (Gaur et al., 2015).

Furthermore, the significant BWB trend in the historical low state suggests
that the acceleration process in flare and low activity states may have
different origins. In principle, we should consider carefully the
contamination from the host galaxy in such a low state. However, the
brightness of the host galaxy ($>20.8$ mag, Scarpa et al. (2000)) was too
faint with respect to the nucleus, and therefore it was unlikely to affect the
color behavior. Yan et al. (2013) found that the shock acceleration is
dominant in the low state, while stochastic turbulence acceleration is
dominant in the flare state.

According to the scenario proposed by Virtanen & Vainio (2005), the electrons
are accelerated at shock front and move along the jet. The shock acceleration
will continuously affect electrons, both in the low and high states. However,
the turbulence appears rather strong in the downstream jet, and the stochastic
turbulence can have considerable contributions to the acceleration process of
electrons (Virtanen & Vainio, 2005; Marscher, 2013). The shock-accelerated
electrons are efficiently re-accelerated in the downstream jet through a
stochastic turbulence, and the increasing radiation contributes to the total
emission and the blazar will turn to a flare state (Yan et al., 2013; Feng et
al., 2020). Moreover, the high-energy $\gamma$-ray photons were also thought
to be originated in the downstream region (Jorstad et al., 2001). Both shock
acceleration and stochastic turbulence acceleration lead to variations in the
magnetic field inside the jet, which further change the brightness of the
blazar (Kirk et al., 1998). If the inter-band amplitude difference in the low
state is larger than that of the flare state, the corresponding color behavior
will be more significant (Dai et al., 2011).

Based on a review of previous studies, Ciprini et al. (2007) proposed an
achromatic trend of PKS 0735+178. Unfortunately, they lacked the variation
data during the 2001 outburst. So they weren’t able to check the spectral
behavior. Thus, their data suggested a rather constant SI and presented an
achromatic trend. This is not in conflict with our conclusions. A similar
achromatic trend caused by the incomplete data was also reported by Gu et al.
(2006). The RWB trend proposed by Sandrinelli et al. (2014) needs an in-depth
discussion. We found that their data also did not contain any complete flare.
In our analysis, such properties should be the features of an achromatic or
weak BWB trend rather than a clear RWB trend. We noticed that their variation
analysis involved the infrared $J$-band. Isler et al. (2017) mentioned that
the color behavior depended on the ratio of the contributions from the
accretion disk and relativistic jet, i.e., the ratio of thermal emission to
non-thermal emission. In the infrared band, the thermal emission is thought to
partly come from the dusty molecular torus that reprocesses the photons from
the broad line region and the accretion disk into the infrared region
(Antonucci, 1993; Perlman et al., 2008). This torus component increases the
proportion of the thermal emission and causes the source to exhibit a RWB
trend. Moreover, Safna et al. (2020) also showed that the RWB seems a common
trend among blazars when an infrared band was involved.

The two opposite color trends reported by Yuan & Fan (2021) seem to interfere
by the clustering of the data points (see Fig. 9 therein). We noticed that
they didn’t take a nightly average for those data points, which would have
affected the final result by superimposing the intraday trend on the long-term
trend. To verify our idea, we reconstructed the color-magnitude correlations
of their data with our method, as shown in Fig. 9. Unfortunately, all the
three diagrams didn’t show any significant correlations, which were mainly
caused by the small sample of data points and associated large uncertainties.
So we couldn’t conclude which color-magnitude correlation exhibited during
their observations.

Figure 9: Reconstructed color-magnitude correlations. The nightly averaged
magnitudes were used.

##  5 Summary

In this work, we monitored PKS 0735+178 with the 85cm telescope at Xinglong
observatory for three nights, and we also collected 11-years multi-wavelength
light curves of this source. The inter-band time lags were calculated through
ICCF and simply discussed the correlations between various bands. We analyzed
the color and spectral behavior to investigate the variation mechanism. Our
conclusions are summarized as follows:

  * •

IDV was observed on MJD 59188. Based on the results of ICCF, the only one non-
zero lag of $22.1$d was detected between the $R$ and $\gamma$-ray bands. Too
close emitting regions prevent a non-zero lag been detected between two
optical bands.

  * •

The correlation between the $\gamma$-ray and $R$ bands is much weaker than
that between the $V$ and $R$ bands. This is mainly due to the different
origins of the optical and $\gamma$-ray emission.

  * •

We found that PKS 0735+178 showed a BWB trend in the optical bands on both
long-term and intraday timescales, and a HWB trend was found in the
$\gamma$-ray band. These color or spectral behaviors will be enhanced during
the active states. Such color/spectral behaviors could be naturally explained
by the jet emission model.

  * •

The most significant BWB trend was found during the historical low state. It
could be explained by the larger difference of the inter-band variation
amplitude of the trough state compared to that of the flare state, which may
be caused by the acceleration process with different origins. Moreover, the
saturation of HWB trend in the high state suggests an efficient particle
acceleration mechanism.

Either as a blazar with special color behaviors or as a potential or even
promising neutrino emitter, PKS 0735+178 is an interesting target and can be
monitored intensively at multiple electromagnetic wavelengths and, especially,
with IceCube or future neutrino telescope (e.g., Cubic Kilometre Neutrino
Telescope), in order to disclose the mechanisms and regions responsible for
the neutrino and electromagnetic emissions. Improved localization and
sensitivity on neutrino detection is required, and a joint electromagnetic-
neutrino observations will be a key in identifying flaring neutrino sources.

The authors thank the anonymous referee for the valuable comments which
improved the quality and clarity of the manuscript. This work has been
supported by the Chinese National Natural Science Foundation nts 11973017.
Data from the Steward Observatory spectropolarimetric monitoring project were
used. This program is supported by Fermi Guest Investigator grants NNX08AW56G,
NNX09AU10G, NNX12AO93G, and NNX15AU81G. The CRTS survey is supported by the
U.S. National Science Foundation under grants AST-0909182.

##  Appendix A DRW Interpolation

The covariance function of a DRW has an exponential form

$S(\Delta t)=\sigma_{DRW}^{2}\exp(-|\Delta t/\tau_{DRW}|)$

where $\Delta t$ is the time interval between two epochs, $\sigma_{DRW}$ is
the DRW amplitude, and $\tau_{DRW}$ is the DRW timescale. The DRW
interpolation was implemented through
JAVELIN111111https://github.com/nye17/javelin algorithm (Zu et al., 2011).
First we fitted the $R$-band light curve to constrain the $\sigma_{DRW}$ and
$\tau_{DRW}$. The $V$-band light curve was considered as a shifted, smoothed
and scaled version of the $R$-band light curve and shared the same
$\tau_{DRW}$ and $\sigma_{DRW}$. In this process, four additional parameters
are added, i.e., the time lag, the kernel width, the transfer function
amplitude, and the ratio between the two bands. Through the Markov chain Monte
Carlo (MCMC) method, we can obtain the best fitting parameters and use them to
predict the two light curves. Furthermore, we restricted the time lags between
$V$ and $R$ bands in $[-10,10]$ days, and the DRW timescale have been
restricted in $[100,300]$ days, which is found for a larger sample of quasars
from the Sloan Digital Sky Survey (MacLeod et al., 2010). The final
predictions of the two light curves were presented in Fig. A1.

Figure A1: The predicted light curves of $V$ and $R$ bands. The shade region
indicates the error bound and the black points indicate the observations.

To check whether the interpolated light curve preserves the properties of the
original light curve, its time-domain power density spectrum (PSD), i.e.,
structure function (SF), was calculated. The SF describes the growth of
variability with time and is widely used in the analysis of AGN variabilities
(e.g., Ulrich et al., 1997; Abdo et al., 2010; Ackermann et al., 2011).
Compared to the PSD in the frequency domain, the SF is less subject to
sampling problems in the presence of very irregular time series (Simonetti et
al., 1985; Paltani et al., 1997). One of the most common form of SF is defined
as

$SF(\tau)=\frac{1}{N(\tau)}\sum_{i=1}\\{[m(t_{i})-m(t_{i}+\tau)]^{2}-[\sigma_{err}(t_{i})^{2}+\sigma_{err}(t_{i}+\tau)^{2}]\\},$

where $\tau$ is the time difference between the data pair, $m(t_{i})$ and
$\sigma_{err}(t_{i})$ are the magnitude and associated uncertainty at time
$t_{i}$, respectively (Bauer et al., 2009). The uncertainty of SF is defined
as

$\sigma_{SF}(\tau)=\frac{1}{N_{{i}}\cdot
SF}\sqrt{\sum_{i=1}\left\\{[m(t_{i}+\tau)-m(t_{i})]^{2}\cdot\left[\sigma_{{err}}^{2}(t_{i}+\tau)+\sigma_{{err}}^{2}(t_{i})\right]\right\\}},$

where $N_{i}$ is the number of data point pairs in the bin. In our settings, a
linear bin $\tau=4$ days was adopted, and the time scales $\tau>200$ days was
discarded.

We calculated the SFs of the $R$-band interpolated light curve ($SF_{i}$) and
original light curve ($SF_{o}$), and the comparison between the two SFs were
presented in Fig. A2. Notice that our CIs were obtained through the
interpolated $R$-band light curve and the original $V$-band light curve, and
thus the comparison between the two SFs in the $V$ band was not given. Except
for the deviation in the low $\tau$ case ($\tau<32$ days), the two SFs are
well consistent with each other, which proves that the interpolation preserves
the properties of the original light curve. One possible reason for the
low-$\tau$ deviation is the poor fit of the DRW model to the short timescale
light curves.

Figure A2: The SFs of the $R$-band interpolated light curve and original light
curve, which are denoted by the black crosses and red dots, respectively.

As stated in Sec. 4.1, we obtained the CIs without interpolation during MJD
55000-55400, and a comparison between the two CIs would also be a good
visualization to the impacts of interpolation. As shown in Fig. A3, the CIs
obtained with and without interpolation are presented. The time baseline of
the two CIs was restricted to the same period, i.e., MJD 55000-55400. It’s
clear that the color-magnitude correlation with interpolation is only sightly
weaker than that without interpolation. It suggests that DRW interpolation
does not produce a spurious correlation. Unfortunately, we don’t have more
data to extend the time baseline to assess its performance on longer
timescale.

Figure A3: The comparison between the CIs obtained with and without
interpolation, which were denoted by the black and red dots, respectively. The
corresponding formula, correlation coefficient $r$ and $p$ value are given.

##  Appendix B ICCF of $R$-$J$ bands

We used the raw data from Sandrinelli et al. (2014) to implement the ICCF. We
restricted the search range of time lags to $[-200,200]$ days and 10000 MC
realizations were applied. The final result showed a $0\%$ failure rate, $0\%$
removal efficiency and a maximum correlation coefficient $\sim 0.8$, which
indicate a reliable measurement, and a $-3.6^{+4.2}_{-4.2}$ days lag was
detected (see Fig. B4).

Figure B4: The same as Fig. 4, but for $R-J$ bands.

## References

  * Abdo et al. (2010) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 722, 520.  https://ui.adsabs.harvard.edu/abs/2010ApJ...722..520A
  * Acciari et al. (2021) Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2021, MNRAS, 504, 1427
  * Ackermann et al. (2011) Ackermann, M., Ajello, M., Allafort, A., et al. 2011, ApJ, 743, 171.  https://ui.adsabs.harvard.edu/abs/2011ApJ...743..171A
  * Agarwal et al. (2021) Agarwal, A., Mihov, B., Andruchow, I., et al. 2021, A&A, 645, A137
  * Akritas & Bershady (1996) Akritas, M. G., & Bershady, M. A. 1996, ApJ, 470, 706.  https://ui.adsabs.harvard.edu/abs/1996ApJ...470..706A
  * Antonucci (1993) Antonucci, R. 1993, ARA&A, 31, 473
  * Bauer et al. (2009) Bauer, A., Baltay, C., Coppi, P., et al. 2009, ApJ, 696, 1241.  https://ui.adsabs.harvard.edu/abs/2009ApJ...696.1241B
  * Böttcher (2007) Böttcher, M. 2007, Ap&SS, 309, 95
  * Bregman et al. (1984) Bregman, J. N., Glassgold, A. E., Huggins, P. J., et al. 1984, ApJ, 276, 454
  * Carini et al. (2011) Carini, M. T., Walters, R., & Hopper, L. 2011, AJ, 141, 49
  * Carswell et al. (1974) Carswell, R. F., Strittmatter, P. A., Williams, R. E., Kinman, T. D., & Serkowski, K. 1974, ApJ, 190, L101
  * Chatterjee et al. (2013) Chatterjee, R., Fossati, G., Urry, C. M., et al. 2013, ApJ, 763, L11
  * Ciprini et al. (2007) Ciprini, S., Takalo, L. O., Tosti, G., et al. 2007, A&A, 467, 465
  * Cohen et al. (2014) Cohen, D. P., Romani, R. W., Filippenko, A. V., et al. 2014, ApJ, 797, 137
  * Dai et al. (2011) Dai, Y., Wu, J., Zhu, Z.-H., Zhou, X., & Ma, J. 2011, AJ, 141, 65
  * de Diego (2014) de Diego, J. A. 2014, AJ, 148, 93
  * de Diego et al. (1998) de Diego, J. A., Dultzin-Hacyan, D., Ramírez, A., & Benítez, E. 1998, ApJ, 501, 69
  * Drake et al. (2009) Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2009, ApJ, 696, 870
  * Fan & Lin (2000) Fan, J. H., & Lin, R. G. 2000, ApJ, 537, 101
  * Fan et al. (1997) Fan, J. H., Xie, G. Z., Lin, R. G., et al. 1997, A&AS, 125, 525
  * Feng et al. (2020) Feng, H.-C., Yang, S., Yang, Z.-X., et al. 2020, ApJ, 902, 42
  * Gaskell & Peterson (1987) Gaskell, C. M., & Peterson, B. M. 1987, ApJS, 65, 1
  * Gaur et al. (2012) Gaur, H., Gupta, A. C., Strigachev, A., et al. 2012, MNRAS, 425, 3002
  * Gaur et al. (2015) Gaur, H., Gupta, A. C., Bachev, R., et al. 2015, MNRAS, 452, 4263
  * Giebels et al. (2007) Giebels, B., Dubus, G., & Khélifi, B. 2007, A&A, 462, 29
  * Goyal et al. (2009) Goyal, A., Gopal-Krishna, Anupama, G. C., et al. 2009, MNRAS, 399, 1622
  * Grier et al. (2017) Grier, C. J., Trump, J. R., Shen, Y., et al. 2017, ApJ, 851, 21
  * Gu et al. (2006) Gu, M. F., Lee, C. U., Pak, S., Yim, H. S., & Fletcher, A. B. 2006, A&A, 450, 39
  * Hartman et al. (1999) Hartman, R. C., Bertsch, D. L., Bloom, S. D., et al. 1999, ApJS, 123, 79
  * Heidt & Wagner (1996) Heidt, J., & Wagner, S. J. 1996, A&A, 305, 42
  * Homayouni et al. (2019) Homayouni, Y., Trump, J. R., Grier, C. J., et al. 2019, ApJ, 880, 126
  * Isler et al. (2017) Isler, J. C., Urry, C. M., Coppi, P., et al. 2017, ApJ, 844, 107
  * Jayasinghe et al. (2019) Jayasinghe, T., Stanek, K. Z., Kochanek, C. S., et al. 2019, MNRAS, 485, 961
  * Jorstad et al. (2001) Jorstad, S. G., Marscher, A. P., Mattox, J. R., et al. 2001, ApJS, 134, 181.  https://ui.adsabs.harvard.edu/abs/2001ApJS..134..181J
  * Kirk et al. (1998) Kirk, J. G., Rieger, F. M., & Mastichiadis, A. 1998, A&A, 333, 452
  * Li et al. (2019) Li, I-Hsiu, J., Shen, Y., Brandt, W. N., et al. 2019, ApJ, 884, 119
  * Liodakis et al. (2019) Liodakis, I., Romani, R. W., Filippenko, A. V., Kocevski, D., & Zheng, W. 2019, ApJ, 880, 32
  * MacLeod et al. (2010) MacLeod, C. L., Ivezić, Ž., Kochanek, C. S., et al. 2010, ApJ, 721, 1014
  * Marscher (2013) Marscher, A. P. 2013, in European Physical Journal Web of Conferences, Vol. 61, European Physical Journal Web of Conferences, 04001.  https://ui.adsabs.harvard.edu/abs/2013EPJWC..6104001M
  * Marscher & Gear (1985) Marscher, A. P., & Gear, W. K. 1985, ApJ, 298, 114
  * Meng et al. (2017) Meng, N., Wu, J., Webb, J. R., Zhang, X., & Dai, Y. 2017, MNRAS, 469, 3588
  * Meng et al. (2018) Meng, N., Zhang, X., Wu, J., Ma, J., & Zhou, X. 2018, ApJS, 237, 30
  * Nanci et al. (2022) Nanci, C., Giroletti, M., Orienti, M., et al. 2022, arXiv e-prints, arXiv:2203.13268.  https://ui.adsabs.harvard.edu/abs/2022arXiv220313268N
  * Nemmen et al. (2012) Nemmen, R. S., Georganopoulos, M., Guiriec, S., et al. 2012, Science, 338, 1445\.  https://ui.adsabs.harvard.edu/abs/2012Sci...338.1445N
  * Nolan et al. (2003) Nolan, P. L., Tompkins, W. F., Grenier, I. A., & Michelson, P. F. 2003, ApJ, 597, 615
  * Paltani et al. (1997) Paltani, S., Courvoisier, T. J. L., Blecha, A., & Bratschi, P. 1997, A&A, 327, 539.  https://ui.adsabs.harvard.edu/abs/1997A&A...327..539P
  * Pandey et al. (2020) Pandey, A., Gupta, A. C., Kurtanidze, S. O., et al. 2020, ApJ, 890, 72
  * Perlman et al. (2008) Perlman, E., Addison, B., Georganopoulos, M., Wingert, B., & Graff, P. 2008, in Blazar Variability across the Electromagnetic Spectrum, 9
  * Peterson et al. (1998a) Peterson, B. M., Wanders, I., Bertram, R., et al. 1998a, ApJ, 501, 82
  * Peterson et al. (1998b) Peterson, B. M., Wanders, I., Horne, K., et al. 1998b, PASP, 110, 660
  * Peterson et al. (2004) Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682
  * Polednikova et al. (2016) Polednikova, J., Ederoclite, A., de Diego, J. A., et al. 2016, MNRAS, 460, 3950
  * Qian et al. (1991) Qian, S. J., Quirrenbach, A., Witzel, A., et al. 1991, A&A, 241, 15
  * Rajput et al. (2020) Rajput, B., Stalin, C. S., & Sahayanathan, S. 2020, MNRAS, 498, 5128
  * Rani et al. (2010) Rani, B., Gupta, A. C., Strigachev, A., et al. 2010, MNRAS, 404, 1992
  * Rector & Stocke (2001) Rector, T. A., & Stocke, J. T. 2001, AJ, 122, 565
  * Safna et al. (2020) Safna, P. Z., Stalin, C. S., Rakshit, S., & Mathew, B. 2020, MNRAS, 498, 3578
  * Sagar et al. (2004) Sagar, R., Stalin, C. S., Gopal-Krishna, & Wiita, P. J. 2004, MNRAS, 348, 176
  * Sahakyan et al. (2022) Sahakyan, N., Giommi, P., Padovani, P., et al. 2022, arXiv e-prints, arXiv:2204.05060.  https://ui.adsabs.harvard.edu/abs/2022arXiv220405060S
  * Sandrinelli et al. (2014) Sandrinelli, A., Covino, S., & Treves, A. 2014, A&A, 562, A79
  * Scarpa et al. (2000) Scarpa, R., Urry, C. M., Padovani, P., Calzetti, D., & O’Dowd, M. 2000, ApJ, 544, 258
  * Shah et al. (2021) Shah, Z., Jithesh, V., Sahayanathan, S., & Iqbal, N. 2021, MNRAS, 2103.13657
  * Shah et al. (2019) Shah, Z., Jithesh, V., Sahayanathan, S., Misra, R., & Iqbal, N. 2019, MNRAS, 484, 3168
  * Sikora et al. (2009) Sikora, M., Stawarz, Ł., Moderski, R., Nalewajko, K., & Madejski, G. M. 2009, ApJ, 704, 38
  * Simonetti et al. (1985) Simonetti, J. H., Cordes, J. M., & Heeschen, D. S. 1985, ApJ, 296, 46.  https://ui.adsabs.harvard.edu/abs/1985ApJ...296...46S
  * Smith et al. (2009) Smith, P. S., Montiel, E., Rightley, S., et al. 2009, arXiv e-prints, arXiv:0912.3621
  * Sun et al. (2018) Sun, M., Grier, C. J., & Peterson, B. M. 2018, PyCCF: Python Cross Correlation Function for reverberation mapping studies, , , ascl:1805.032
  * Tody (1986) Tody, D. 1986, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 627, Instrumentation in astronomy VI, ed. D. L. Crawford, 733
  * Tody (1993) Tody, D. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 173
  * Ulrich et al. (1997) Ulrich, M.-H., Maraschi, L., & Urry, C. M. 1997, ARA&A, 35, 445
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803
  * Virtanen & Vainio (2005) Virtanen, J. J. P., & Vainio, R. 2005, ApJ, 621, 313
  * Wang & Jiang (2020) Wang, Y.-F., & Jiang, Y.-G. 2020, ApJ, 902, 41
  * Weaver et al. (2022) Weaver, Z. R., Jorstad, S. G., Marscher, A. P., et al. 2022, ApJS, 260, 12.  https://ui.adsabs.harvard.edu/abs/2022ApJS..260...12W
  * Webb et al. (1988) Webb, J. R., Smith, A. G., Leacock, R. J., et al. 1988, AJ, 95, 374
  * Webb et al. (1998) Webb, J. R., Freedman, I., Howard, E., et al. 1998, AJ, 115, 2244
  * Weng et al. (2020) Weng, S.-S., Chen, Y., Wang, T.-T., et al. 2020, MNRAS, 491, 2576
  * Williamson et al. (2014) Williamson, K. E., Jorstad, S. G., Marscher, A. P., et al. 2014, ApJ, 789, 135
  * Wu et al. (2012) Wu, J., Böttcher, M., Zhou, X., et al. 2012, AJ, 143, 108
  * Wu et al. (2011) Wu, J., Zhou, X., Ma, J., & Jiang, Z. 2011, MNRAS, 418, 1640
  * Wu et al. (2007) Wu, J., Zhou, X., Ma, J., et al. 2007, AJ, 133, 1599
  * Xue & Cui (2005) Xue, Y., & Cui, W. 2005, ApJ, 622, 160
  * Yan et al. (2013) Yan, D., Zhang, L., Yuan, Q., Fan, Z., & Zeng, H. 2013, ApJ, 765, 122
  * Yuan & Fan (2021) Yuan, Y. H., & Fan, J. H. 2021, PASP, 133, 074101
  * Zhang et al. (2018) Zhang, X., Wu, J., & Meng, N. 2018, MNRAS, 478, 3513
  * Zu et al. (2013) Zu, Y., Kochanek, C. S., Kozłowski, S., & Udalski, A. 2013, ApJ, 765, 106
  * Zu et al. (2011) Zu, Y., Kochanek, C. S., & Peterson, B. M. 2011, ApJ, 735, 80