# Hyperbolic Audio Source Separation

###### Abstract

We introduce a framework for audio source separation using embeddings on a
hyperbolic manifold that compactly represent the hierarchical relationship
between sound sources and time-frequency features. Inspired by recent
successes modeling hierarchical relationships in text and images with
hyperbolic embeddings, our algorithm obtains a hyperbolic embedding for each
time-frequency bin of a mixture signal and estimates masks using hyperbolic
softmax layers. On a synthetic dataset containing mixtures of multiple people
talking and musical instruments playing, our hyperbolic model performed
comparably to a Euclidean baseline in terms of source to distortion ratio,
with stronger performance at low embedding dimensions. Furthermore, we find
that time-frequency regions containing multiple overlapping sources are
embedded towards the center (i.e., the most uncertain region) of the
hyperbolic space, and we can use this certainty estimate to efficiently trade-
off between artifact introduction and interference reduction when isolating
individual sounds.

Index Terms—  audio source separation, hyperbolic space, speech, music, sound
hierarchy

##  1 Introduction

A fundamental paradigm in deep learning-based audio source separation
algorithms is the idea of applying a mask to a feature representation (e.g., a
magnitude spectrogram or learned basis) of an audio mixture signal [1, 2, 3,
4]. By inverting or decoding the masked feature representation, we obtain the
isolated time-domain source signals. While techniques that learn feature
encoders and decoders directly based on waveform signals have achieved
impressive performance [3, 5], they lack interpretability compared to
techniques based on time-frequency (T-F) representations such as the short-
time Fourier transform (STFT) spectrogram [6, 7]. Among these, algorithms such
as deep clustering [2] and deep attractor networks [8] learn an embedding
vector for each T-F bin, and create masks using classifiers and/or clustering
algorithms to assign embeddings to sources. A fundamental problem for these
approaches then becomes how to best learn a discriminative embedding for each
T-F bin.

In this work, we take inspiration from recent advances in modeling language
[9, 10, 11], graphs [12, 13, 11], and images [14, 15, 16, 17, 18, 19] in
hyperbolic space, and explore their relevance for audio source separation.
Unlike Euclidean spaces, hyperbolic spaces have an inherent ability to infer
hierarchical representations from data with very little distortion [20, 21].
Hierarchical and tree-like structures are ubiquitous in many types of audio
processing problems such as musical instrument recognition [22, 23] and
separation [24], speaker identification [25], and sound event detection [26,
27]. However, all of these approaches model the hierarchical information
globally by computing a single embedding vector for an entire audio clip.
Recent work in image segmentation [17] learns a hyperbolic embedding at the
pixel level, and we take a similar approach by computing a hyperbolic
embedding for each T-F bin of an audio mixture spectrogram, as illustrated in
Fig. 1.

Fig. 1: Illustration of hyperbolic source separation. Left plot demonstrates
the process of taking a T-F bin from a mixture spectrogram $X_{t,f}$ and
projecting its Euclidean embedding to $z^{h}_{t,f}$ on the Poincaré ball (red
squares); the dashed purple arrow represents the distance of the resulting
embedding to the hyperplane $\text{H}^{c}_{p_{k},a_{k}}$. Right plot shows the
Poincaré ball with the resulting decision boundaries from a trained model
(top). The mixture embeddings (scatters) are colored-coded according to the
maximum softmax layer output. The predicted T-F masks for each of the leaf
sources (bottom) are plotted similarly, and silent bins (-40 dB and below) are
colored in white.

As has been shown by many recent computer vision studies [14, 19], hyperbolic
network layers [28] can be added to the end of existing network architectures
to learn and classify hyperbolic embeddings. In addition to benefits in terms
of hierarchical modeling, hyperbolic layers can learn strong embeddings at low
dimensions [25], represent uncertainty [14, 19], and estimate features such as
edges in images based on their distance to the origin in the hyperbolic space
[17]. Our main contribution is to explore whether similar observations are
applicable to hyperbolic embeddings of T-F bins, and we propose new avenues
for interpreting and interacting with hyperbolic embeddings in audio source
separation applications.

Specifically, we use a simulated dataset containing multiple overlapping
speakers and musical instruments to illustrate our approach as shown in Fig.
1. We classify each T-F bin embedding using hyperbolic softmax layers [28],
and each T-F bin is further assigned to both a parent and child class using a
hierarchical softmax [17]. We compare various widely used source separation
loss functions, but find that energy-weighted cross-entropy performs best.
Corroborating work in application domains such as computer vision and speaker
identification, we find that performance with low embedding dimensions is
relatively strong in hyperbolic space. In terms of uncertainty, we show how
the distance to the origin can be used to approximate computationally
expensive certainty estimates obtained with Monte-Carlo dropout [29] using
only a single forward pass. We also observe that T-F bins containing multiple
overlapping sources are consistently embedded near the origin (i.e., they are
considered the most uncertain embeddings). We then provide experiments
demonstrating how this uncertainty can be exploited to navigate the trade-off
between interference reduction and artifact introduction, which is fundamental
to audio source separation.

##  2 Hyperbolic Audio Source Separation

###  2.1 Embedding-based source separation

We consider the problem of separating a mixture audio signal into $K$ sources
of interest. While we later focus more precisely on class-based separation,
where the sources belong to $K$ different sound classes, we first give a
general presentation of separation approaches relying on time-frequency (T-F)
embeddings.

Source separation is often formulated as a mask-inference problem [1] in a
time-frequency domain, such as that obtained by a short-time Fourier transform
(STFT), where one seeks to obtain the complex spectrogram
$S^{k}\in\mathbb{C}^{T\times F}$ of each of the $K$ sources by multiplying
element-wise the spectrogram of the mixture $X\in\mathbb{C}^{T\times F}$ with
a real-valued mask $M^{k}\in\mathbb{R}^{T\times F}$: $\hat{S}^{k}=M^{k}\odot
X.$ A common way to obtain the mask $M^{k}_{t,f}$ at time-frequency bin
$(t,f)$ relies on first computing an $L$-dimensional Euclidean embedding
$z^{\text{e}}_{t,f}\in{\mathbb{R}}^{L}$ from $X$ using a deep neural network
$f_{\theta}(.)$ parameterized by $\theta$, such that
$Z^{\text{e}}=f_{\theta}(X)\in\mathbb{R}^{T\times F\times L}$. These
embeddings can be projected to dimension $K$ and passed through mask-output
non-linearities such as sigmoid or softmax, or clustered in ${\mathbb{R}}^{L}$
to obtain a clustering of the T-F bins, such as in deep clustering [2]. In
this work, we are interested in replacing these Euclidean embeddings and the
way they are used to derive masks by equivalents in a hyperbolic space, as we
hope that this leads to better embeddings that encompass hierarchical
relationships between sounds. As a starting point, we focus here on the case
of class-based separation with a softmax output layer, as this has been the
primary non-linearity used in hyperbolic neural networks in computer vision
[14, 19, 17], where the Euclidean embeddings and the multinomial logistic
regression (MLR) classification approach applied to them are replaced by
hyperbolic embeddings and a hyperbolic MLR. We leave the application to other
non-linearities and to deep clustering to future work.

###  2.2 Hyperbolic embedding-based source separation

We first briefly review several notions pertaining to Riemannian manifolds and
hyperbolic spaces that are useful in introducing hyperbolic embeddings and
MLR, following [28, 30, 10].

A Riemannian manifold is defined as a pair consisting of a manifold
$\mathcal{M}$ and a Riemannian metric $g$, where $g=(g_{x})_{x\in\mathcal{M}}$
defines the local geometry $g_{x}$ (i.e., an inner product) in the tangent
space $T_{x}\mathcal{M}$ at each point $x\in\mathcal{M}$. While $g$ defines
the geometry locally on $\mathcal{M}$, it also defines the global shortest
path, or _geodesic_ (analogous to a straight line in Euclidean space), between
two given points on $\mathcal{M}$. One can define an _exponential map_
$\text{exp}_{x}$ which projects any vector $v$ of the tangent space
$T_{x}\mathcal{M}$ onto $\mathcal{M}$, such that
$\text{exp}_{x}(v)\in\mathcal{M}$, and inversely a logarithmic map which
projects any point in $\mathcal{M}$ back onto the tangent space at $x$.

The $L$-dimensional hyperbolic space is an $L$-dimensional Riemannian manifold
of constant negative curvature $-c$. It can be described using several
isometric models, among which we focus here on the Poincaré unit ball model
$(\mathbb{D}_{c}^{L},g^{\mathbb{D}}_{c})$, defined in the space
$\mathbb{D}_{c}^{L}=\\{x\in\mathbb{R}^{L}\ |\ c||x||^{2}<1\\}\\}$. We assume
$c>0$, such that $\mathbb{D}_{c}^{L}$ corresponds to a ball of radius
$\frac{1}{\sqrt{c}}$ in Euclidean space. Its Riemannian metric is given by
$g^{\mathbb{D}}_{c}(x)=(\lambda^{c}_{x})^{2}g^{E}$, where
$\lambda^{c}_{x}=2/(1-c\|x\|^{2})$ is a so-called _conformal factor_ and
$g^{E}$ the Euclidean metric. Given two points $x,y\in\mathbb{D}_{c}^{L}$,
their induced hyperbolic distance $d_{c}$ is obtained as

$d_{c}(x,y)=\frac{2}{\sqrt{c}}\text{tanh}^{-1}(\sqrt{c}\|-x\oplus_{c}y\|),$
(1)

where $\oplus_{c}$ denotes the Möbius addition in $\mathbb{D}_{c}^{L}$,
defined as

$x\oplus_{c}y=\frac{(1+2c\langle
x,y\rangle+c\|y\|^{2})x+(1-c\|x\|^{2})y}{1+2c\langle
x,y\rangle+c^{2}\|x\|^{2}\|y\|^{2}}.$ (2)

One way to go back and forth between the Euclidean space ${\mathbb{R}}^{L}$
and the hyperbolic space $\mathbb{D}_{c}^{L}$ is to use the exponential and
logarithmic maps at the origin $0$, as
$T_{0}\mathbb{D}_{c}^{L}={\mathbb{R}}^{L}$, which can be obtained for
$v\in{\mathbb{R}}^{L}\setminus\\{0\\}$ and
$y\in\mathbb{D}_{c}^{L}\setminus\\{0\\}$ as:

$\exp^{c}_{0}(v)=\frac{\text{tanh}(\sqrt{c}\|v\|)}{\sqrt{c}\|v\|}v,\quad\log^{c}_{0}(y)=\frac{\text{tanh}^{-1}(\sqrt{c}\|y\|)}{\sqrt{c}\|y\|}y.$
(3)

We can thus obtain hyperbolic embeddings $z^{h}_{t,f}\in\mathbb{D}_{c}^{L}$ at
the output of a classical neural network such as $f_{\theta}$ by simply
projecting the usual Euclidean embeddings in that way:

$z^{h}_{t,f}=\exp_{0}(z^{e}_{t,f})=\exp_{0}(f_{\theta}(X)_{t,f}).$ (4)

To define a hyperbolic softmax based on these hyperbolic embeddings, the
Euclidean MLR can be generalized to the Poincaré ball as in [28]. In the
Euclidean space, MLR is performed by considering the logits obtained by
calculating the distance of an input’s embedding $z\in{\mathbb{R}}^{L}$ (such
as $z=z^{e}_{t,f}$) to each of $K$ class hyperplanes, where the $k$-th class
hyperplane is determined by a normal vector $a_{k}\in\mathbb{R}^{L}$ and a
point $p_{k}\in\mathbb{R}^{L}$ on that hyperplane. Analogously, one can define
a Poincaré hyperplane $\text{H}^{c}_{p_{k},a_{k}}$ by considering the union of
all geodesics passing by a point $p_{k}$ and orthogonal to a normal vector
$a_{k}$ in the tangent space $T_{p_{k}}\mathbb{D}_{c}^{L}$ at $p_{k}$.
Hyperbolic MLR can then be defined by considering the distance from a
hyperbolic embedding $z=z^{d}_{t,f}\in\mathbb{D}_{c}^{L}$ to each
$\text{H}^{c}_{p_{k},a_{k}}$, leading to the following formulation as shown in
[28]:

$p(\kappa\\!=\\!k|z)\\!\propto\exp\Bigl{(}\frac{\lambda^{c}_{p_{k}}\|a_{k}\|}{\sqrt{c}}\text{sinh}^{-1}\\!\Bigl{(}\frac{2\sqrt{c}|\langle-
p_{k}\\!\oplus_{c}\\!z,a_{k}\rangle|}{(1\\!-\\!c\|\\!-\\!p_{k}\\!\oplus_{c}\\!z\|^{2})\|a_{k}\|}\Bigr{)}\Bigr{)}.$
(5)

We can use the probability $p(\kappa=k|z^{d}_{t,f})$ that T-F bin $(t,f)$ is
dominated by the $k$-th source to obtain $K$ source-specific mask values for
each T-F bin in the input spectrogram $X$. This procedure is illustrated in
Fig. 1. Note that $p_{k}$ and $a_{k}$ both parameterize the $k$-th hyperbolic
hyperplane and are trainable. All parameters of the network $f_{\theta}$ and
of the hyperplanes can be optimized using classical source separation
objective functions, either on the masks or on the reconstructed signals.

###  2.3 Uncertainty in Audio Source Separation

It has been shown repeatedly [17, 30, 28, 14] that the distance
$d_{c}(0,z^{h}_{t,f})$ of the projected embeddings to the center of the
Poincaré ball can serve as a reliable measure of classification certainty.
This hyperbolic distance can be computed from the $L_{2}$ norm of
$z^{h}_{t,f}$ as considered in the Euclidean space using the monotonous
relationship $\log((1+c\|z^{h}_{t,f}\|^{2})/(1-c\|z^{h}_{t,f}\|^{2}))$, so
that the $L_{2}$ norm can equivalently be considered as a measure of
certainty. We here aim at validating this hypothesis with audio embeddings and
bringing further light onto the notion of certainty in the context of audio
source separation. While hyperbolic uncertainty has previously been exploited
towards the network optimization stage [25], the direct impact of this notion
on quality criteria such as amount of artifacts, distortion, or interferences
present in resulting signal has yet to be explored.

Furthermore, following the observation presented in [17], we argue that the
notion of certainty observed on audio embeddings provides meaningful insights
of signal content, such as the complexity of the scene or the number of active
sources it may contain.

##  3 Experimental Setup

Slakh+Speech Dataset: To test our proposed method we built a simple
hierarchical source separation dataset containing mixtures from two “parent”
classes - _music_ and _speech_ , and five “leaf” classes - _bass_ , _drums_ ,
_guitar_ , _speech-male_ , and _speech-female_. As building blocks, we use the
clean subset of LibriSpeech [31], and Slakh2100 [32], which is a dataset of
2,100 synthetic musical mixtures, each containing bass, drums, and guitar
stems in addition to various other instruments. We built a dataset consisting
of 1947 mixtures, each 60 s in length, for a total of about 32 hours. The data
splits are 70%, 20%, and 10% for training, validation, and testing sets,
respectively. The _speech-male_ source target consists of male speech
utterances randomly picked (without replacement) and concatenated
consecutively (without overlap) until the 60 s track length is reached. Any
signal from the last concatenated utterance exceeding that length was
discarded. We used the same procedure for _speech-female_. For Slakh2100, we
only selected the first 60 s of the _bass_ , _drums_ , and _guitar_ stems for
each track. Any tracks with a duration less than 60 s were discarded. All
sources were summed without applying additional gains to make the overall
mixture along with the _speech_ and _music_ submixes. This lead to challenging
input SDR values (“No Proc.” in Table 1), with standard deviation values
ranging from 2-13 dB depending on the class.

Network architecture and training setup: Our model consists of four BLSTM
layers with 600 units in each direction, followed by a dense layer to obtain
an $L$-dimensional Euclidean embedding for each T-F bin. A dropout of $0.3$ is
applied on the output of each BLSTM layer, except the last. For the hyperbolic
models ($c\\!>\\!0)$, an exponential projection layer is placed after the
dense layer, mapping the Euclidean embeddings onto the Poincaré ball with
curvature $-c$. As discussed in Section 2, MLR layers, either Euclidean or
hyperbolic with softmax activation functions are then used to obtain masks for
each of the source classes. In practice, we follow the hierarchical softmax
approach from [17], and have two MLR layers: one with $K\\!=\\!2$ for the
parent (speech/music) sources, and a second with $K\\!=\\!5$ for the leaf
classes. We use the mixture phase for resynthesis and compare multiple
training objectives in Section 4.

We use the ADAM optimizer for Euclidean parameters, and the Riemannian ADAM
[33] implementation from _geoopt_ [34] for hyperbolic parameters. All models
are trained using chunks of 3.2 s and a batch size of 10 for 300 epochs using
an initial learning rate of $10^{-3}$, which is halved if the validation loss
does not improve for 10 epochs. We use an STFT size of 32 ms with 50% overlap
and square-root Hann window.

##  4 Experimental Analysis

Table 1: SI-SDR in dB on the Slakh+Speech test set with different loss
functions and embedding dimensions $L$. The “Hyp.” column denotes whether the
model was trained using hyperbolic ($c=1$) or Euclidean embeddings.

|  |  | Parents | Music Leaves | Speech Leaves |
---|---|---|---|---|---|---
Loss | $L$ | Hyp. | Music | Speech | Bass | Drums | Guitar | Male | Female | Avg.
No Proc. |  |  | $-2.94$ | $2.99$ | $-8.84$ | $-11.8$ | $-8.73$ | $-3.17$ | $-4.16$ | $-5.24$
Oracle${}_{\text{PSF}}$ |  |  | $10.25$ | $13.80$ | $5.56$ | $9.79$ | $6.98$ | $10.10$ | $11.05$ | $9.65$
PSA | 2 | $\times$ | $7.73$ | $11.28$ | $2.70$ | $4.49$ | $2.98$ | $6.15$ | $6.48$ | $5.97$
WA | 2 | $\times$ | $7.08$ | $10.83$ | $2.80$ | $3.18$ | $2.33$ | $5.49$ | $5.64$ | $5.34$
CE${}_{\text{IBM}}$ | 2 | $\times$ | $5.06$ | $8.76$ | $-2.29$ | $1.84$ | $-1.03$ | $2.76$ | $2.93$ | $2.57$
CE${}_{\text{IBM, W.}}$ | 2 | $\times$ | $8.03$ | $11.47$ | $2.65$ | $4.86$ | $2.48$ | $6.44$ | $6.64$ | $6.08$
PSA | 2 | ✓ | $7.46$ | $10.87$ | $-4.34$ | $-6.88$ | $-4.61$ | $5.97$ | $6.28$ | $2.11$
WA | 2 | ✓ | $7.26$ | $10.96$ | $2.50$ | $3.72$ | $2.32$ | $5.34$ | $5.43$ | $5.36$
CE${}_{\text{IBM}}$ | 2 | ✓ | $5.42$ | $9.04$ | $-0.26$ | $2.56$ | $0.78$ | $3.40$ | $3.83$ | $3.54$
CE${}_{\text{IBM, W.}}$ | 2 | ✓ | $7.90$ | $11.36$ | $2.93$ | $5.10$ | $3.19$ | $6.52$ | $6.99$ | $6.29$
PSA | 128 | $\times$ | $7.73$ | $11.16$ | $3.25$ | $4.95$ | $3.28$ | $6.10$ | $6.37$ | $6.12$
WA | 128 | $\times$ | $7.40$ | $11.20$ | $3.30$ | $4.80$ | $3.01$ | $5.89$ | $6.11$ | $5.96$
CE${}_{\text{IBM}}$ | 128 | $\times$ | $6.54$ | $10.09$ | $1.73$ | $6.06$ | $2.48$ | $5.07$ | $5.72$ | $5.39$
CE${}_{\text{IBM, W.}}$ | 128 | $\times$ | $7.93$ | $11.47$ | $3.31$ | $6.84$ | $3.83$ | $6.78$ | $7.30$ | $6.78$
PSA | 128 | ✓ | $7.45$ | $10.89$ | $2.89$ | $4.88$ | $3.13$ | $5.78$ | $5.99$ | $5.86$
WA | 128 | ✓ | $6.97$ | $11.23$ | $3.26$ | $4.76$ | $2.87$ | $5.87$ | $5.99$ | $5.85$
CE${}_{\text{IBM}}$ | 128 | ✓ | $6.02$ | $9.72$ | $1.27$ | $5.78$ | $2.01$ | $4.50$ | $5.18$ | $4.93$
CE${}_{\text{IBM, W.}}$ | 128 | ✓ | $7.80$ | $11.34$ | $3.28$ | $6.13$ | $3.72$ | $6.66$ | $7.09$ | $6.57$

Model comparisons: Table 1 presents the scale-invariant signal-to-distortion
ratio (SI-SDR) [35] on the Slakh+Speech dataset described in Sec. 3. We
include the no processing condition (lower bound, using the mixture as
estimate) and oracle phase sensitive mask [1] (upper bound). For loss
functions, we use the phase-sensitive approximation (PSA) [1] with $L_{1}$
loss [36], and the waveform approximation (WA) loss training through the iSTFT
[37]. Additionally, to compare with related image segmentation approaches
[17], we explore using a cross entropy (CE) loss with the ideal binary mask
(IBM) as a training target [4], denoted by CE${}_{\text{IBM}}$ in Table 1, and
a magnitude-ratio-weighted CE loss inspired by the deep clustering weighting
function from [36], denoted by CE${}_{\text{IBM, W.}}$, where the weight for
each T-F bin is set to the ratio of the magnitude of the mixture at that bin
to the sum over all bins.

In Table 1, CE${}_{\text{IBM}}$ consistently performs poorly, likely because
low-energy T-F bins are weighted equally with high energy bins. However,
CE${}_{\text{IBM, W.}}$ consistently outperforms PSA and WA, which we
hypothesize is because the CE objective is particularly well-matched to the
softmax mask nonlinearity we use in this work. We also note that, consistent
with previous hyperbolic work [38, 28, 17, 25], performance of hyperbolic
models is relatively stronger than Euclidean ones at low ($L=2$) embedding
dimension. This result is further confirmed in Fig. 2, where we explore
multiple embedding dimensions and curvature parameters. Unless otherwise
stated, the hyperbolic configuration used in all subsequent experiments uses
the CE${}_{\text{IBM, W.}}$ loss, curvature parameter $c=0.1$, and embedding
dimension $L=2$.

Fig. 2: SI-SDR [dB] averaged across all (parent and leaf) sources for various
curvature values and embedding dimensions.  Fig. 3: Histograms denoting the
distribution of $L_{2}$ norms to the origin of the Poincaré Ball for all
embeddings in our test-set. We define the number of active sources for each
embedding based on the ground-truth mask values of their associated T-F bin.
Fig. 4: Separation metrics as a function of certainty thresholding (here
applied to the $L_{2}$ norm $\|z^{h}_{t,f}\|_{2}$). A threshold of $0$ means
all hyperbolic embeddings are taken into account during inference, and a
threshold of $0.95$ means only embeddings with high certainty (i.e., high
distance from the origin) are considered, while bins below the threshold are
silenced.

Hyperbolic Certainty and Audio Semantics: In Fig. 3, we show that embeddings
associated with T-F bins containing many active sources (i.e., 4+) tend to be
positioned close to the center of the Poincaré Ball, where the certainty is
low. In contrast, embeddings with a single active source tend to be located
closer to the edge of the Poincaré Ball, where the certainty is high.
Intuitively, the higher the number of sound classes (i.e., interferences)
active, the more uncertain as for what specific single source the bin
embedding may belong to. We consider a source active for a given T-F bin if it
satisfies the following conditions: (1) it is non-silent, i.e., it is not more
than $20.0$ dB below the maximum T-F bin energy in a given sound file, and (2)
the ratio of the ground-truth T-F bin magnitude for that source divided by the
sum of ground-truth magnitudes for all sources (i.e., the ideal ratio mask) is
greater than $0.1$.

Hyperbolic Distance - Interference and Artifacts: We also explore a simple
approach for using hyperbolic certainty estimates to modify and interact with
mask values. Specifically, we embed all T-F bins for a given mixture in the
Poincaré ball, and any bins that are close to the origin (i.e., low certainty)
have their mask values set to zero. Figure 4 shows the impact of varying this
certainty threshold on the SI-SDR, SIR, SAR [35, 39] metrics. Each data point
denotes an evaluation pass on the entire test set given a certainty threshold
(x-axis). We see that if we silence uncertain T-F bins in the mask
predictions, the resulting signals will contain less interference at the
expense of an increase in artifacts. Using hyperbolic certainty to control the
trade-off between artifact introduction and interference reduction could be an
exciting area to further explore in future work.

Comparison with Bayesian certainty: To further validate our hyperbolic
certainty measurements, we lead an experiment comparing hyperbolic and
Bayesian certainty. The latter approach is generally achieved by means of
Monte-Carlo dropout [40, 41] where multiple stochastic forward passes using a
trained network with different dropout realizations are performed to obtain
the posterior distribution for a single input example. The network prediction
and its associated certainty map can then be inferred from the mean and
variance computed over all stochastic passes. In order to deal with the fact
that our network makes predictions over multiple sound classes, we follow the
approach presented in [29] and formulate the Bayesian certainty $\zeta_{t,f}$
at T-F bin $(t,f)$ as the negative predictive entropy

$\zeta_{t,f}\\!=\\!\\!\sum_{k=1}^{K_{\text{leaf}}}\Bigl{(}\frac{1}{N}\\!\sum_{n}p(\kappa\\!=\\!k|\hat{z}^{(n)}_{t,f})\Bigr{)}\text{log}\Bigl{(}\frac{1}{N}\\!\sum_{n}p(\kappa\\!=\\!k|\hat{z}^{(n)}_{t,f})\Bigr{)},\\!\\!$
(6)

where $K_{\text{leaf}}$ denote the number of leaf classes (i.e., without the
parent classes), $N$ the total number of stochastic passes, and
$p(\kappa=k|\hat{z}^{(n)}_{t,f})$ the probability, expressed as the softmax
output, that T-F bin $(t,f)$ belongs to source class $k$, and
$\hat{z}^{(n)}_{t,f}$ is the embedding obtained with the dropped out network
parameters for stochastic pass $n$. Following [29], we set the dropout rate to
$0.5$ at the output of each BLSTM layer of our network. A low entropy (i.e.,
high certainty) means that the predictability is high, while a high entropy
(low certainty) means that mask values fluctuate significantly across passes.

Fig. 5: Comparison between Bayesian and hyperbolic certainty. The Bayesian
certainty is computed by running 1000 passes with Monte-Carlo dropout during
inference and computing the bin-wise negative entropy across all passes on the
five leaf nodes (i.e., excluding the parent nodes). The hyperbolic certainty
is obtained by computing the hyperbolic distance from the Poincaré ball origin
$d_{c}(0,z^{h}_{t,f})$ for each T-F bin embedding $z^{h}_{t,f}$.

Fig. 5 provides a visual contrast of the certainty map computed using 1000
Monte-Carlo dropout iterations and the one obtained using the hyperbolic
distance from the Poincaré ball origin for each T-F embedding in a single
forward pass. The colorscale limits on both maps are set using the 30th and
95th percentiles to ensure a fair qualitative comparison, even though the
scales are different. We observe a clear resemblance between both maps. The
correlation coefficient between the two certainty maps is $\rho=0.75$. This
observation tells us that the hyperbolic certainty map, which comes for free
(i.e., $N=1$), can be as interpretable as its Bayesian counterpart, which
requires multiple forward passes at inference time.

A demo of a hyperbolic separation interface is available
111https://darius522.github.io/hyperbolic-audio-sep/.

##  5 Conclusion

We have investigated the use of the Poincaré Ball model to perform audio
source separation in the hyperbolic space. Our hyperbolic model operates and
computes T-F embeddings in the Euclidean space and projects them onto the
hyperbolic space. Masks are obtained by hyperbolic multinomial logistic
regression considering the distance from hyperbolic embeddings to hyperbolic
hyperplanes. Through our experimental setup, we have demonstrated that
hyperbolic audio embeddings can convey useful information in regards to
uncertainty and underlying hierarchical sound structures. We associated these
notions to known audio concepts such as artifacts and interferences. In the
future, we aim at exploring how the hyperbolic space, and especially the
notion of hierarchy, could benefit additional audio-related tasks, such as
audio tagging and sound-event detection. In the context of audio source
separation, we believe that further connections to deep clustering, and more
complex and deeper sound taxonomies are worth exploring.

## References

  * [1] H. Erdogan, J. R. Hershey, S. Watanabe, and J. Le Roux, “Phase-sensitive and recognition-boosted speech separation using deep recurrent neural networks,” in _Proc. ICASSP_ , 2015, pp. 708–712.
  * [2] J. R. Hershey, Z. Chen, J. Le Roux, and S. Watanabe, “Deep clustering: Discriminative embeddings for segmentation and separation,” in _Proc. ICASSP_ , 2016, pp. 31–35.
  * [3] Y. Luo and N. Mesgarani, “Conv-TasNet: Surpassing ideal time–frequency magnitude masking for speech separation,” _IEEE/ACM Trans. Audio, Speech, Lang. Process._ , vol. 27, no. 8, pp. 1256–1266, 2019.
  * [4] D. Wang and J. Chen, “Supervised speech separation based on deep learning: An overview,” _IEEE/ACM Trans. Audio, Speech, Lang. Process._ , vol. 26, no. 10, pp. 1702–1726, 2018.
  * [5] A. Défossez, “Hybrid spectrogram and waveform source separation,” in _Proc. ISMIR Workshop on Music Source Separation_ , 2021.
  * [6] J. Heitkaemper, D. Jakobeit, C. Boeddeker, L. Drude _et al._ , “Demystifying TasNet: A dissecting approach,” in _Proc. ICASSP_ , 2020, pp. 6359–6363.
  * [7] R. Parikh, I. Kavalerov, C. Espy-Wilson, and S. Shamma, “Harmonicity plays a critical role in DNN based versus in biologically-inspired monaural speech segregation systems,” in _Proc. ICASSP_ , 2022, pp. 536–540.
  * [8] Y. Luo, Z. Chen, and N. Mesgarani, “Speaker-independent speech separation with deep attractor network,” _IEEE/ACM Trans. Audio, Speech, Lang. Process._ , vol. 26, no. 4, pp. 787–796, 2018.
  * [9] A. Suzuki, A. Nitanda, J. Wang, L. Xu _et al._ , “Generalization error bound for hyperbolic ordinal embedding,” in _Proc. ICML_ , 2021, pp. 10 011–10 021.
  * [10] R. Shimizu, Y. Mukuta, and T. Harada, “Hyperbolic neural networks++,” in _Proc. ICLR_ , 2021.
  * [11] J. Sun, Z. Cheng, S. Zuberi, F. Perez _et al._ , “HGCF: Hyperbolic graph convolution networks for collaborative filtering,” in _Proc. WWW_ , 2021\.
  * [12] Q. Liu, M. Nickel, and D. Kiela, “Hyperbolic graph neural networks,” in _Proc. NeurIPS_ , 2019.
  * [13] G. Bachmann, G. B’ecigneul, and O.-E. Ganea, “Constant curvature graph convolutional networks,” in _Proc. ICML_ , 2020.
  * [14] V. Khrulkov, L. Mirvakhabova, E. Ustinova, I. Oseledets _et al._ , “Hyperbolic image embeddings,” in _Proc. CVPR_ , Jun. 2020.
  * [15] I. Chami, A. Gu, V. Chatziafratis, and C. Ré, “From trees to continuous embeddings and back: Hyperbolic hierarchical clustering,” in _Proc. NeurIPS_ , Dec. 2020.
  * [16] Y. Guo, X. Wang, Y. Chen, and S. X. Yu, “Free hyperbolic neural networks with limited radii,” in _Proc. CVPR_ , 2022.
  * [17] M. G. Atigh, J. Schoep, E. Acar, N. van Noord _et al._ , “Hyperbolic image segmentation,” in _Proc. CVPR_ , Jun. 2022, pp. 4453–4462.
  * [18] Z. Weng, M. G. Ogut, S. Limonchik, and S. Yeung, “Unsupervised discovery of the long-tail in instance segmentation using hierarchical self-supervision,” in _Proc. CVPR_ , 2021, pp. 2603–2612.
  * [19] D. Surís, R. Liu, and C. Vondrick, “Learning the predictability of the future,” in _Proc. CVPR_ , 2021, pp. 12 607–12 617.
  * [20] F. Sala, C. De Sa, A. Gu, and C. Re, “Representation tradeoffs for hyperbolic embeddings,” in _Proc. ICML_ , Jul. 2018, pp. 4460–4469.
  * [21] R. Sarkar, “Low distortion delaunay embedding of trees in hyperbolic plane,” in _Proc. International Symposium On Graph Drawing_ , 2012, pp. 355–366.
  * [22] S. Essid, G. Richard, and B. David, “Hierarchical classification of musical instruments on solo recordings,” in _Proc. ICASSP_ , 2006, pp. 817–820.
  * [23] H. F. Garcia, A. Aguilar, E. Manilow, and B. Pardo, “Leveraging hierarchical structures for few-shot musical instrument recognition,” in _Proc. ISMIR_ , 2021, pp. 220–228.
  * [24] E. Manilow, G. Wichern, and J. Le Roux, “Hierarchical musical instrument separation.” in _Proc. ISMIR_ , 2020, pp. 376–383.
  * [25] J. Lee, K. Sung-Bin, S. Kang, and T.-H. Oh, “Lightweight speaker recognition in Poincaré spaces,” _IEEE Signal Process. Lett._ , vol. 29, pp. 224–228, 2022.
  * [26] A. Jati, N. Kumar, R. Chen, and P. Georgiou, “Hierarchy-aware loss function on a tree structured label space for audio event detection,” in _Proc. ICASSP_ , 2019, pp. 6–10.
  * [27] J. Cramer, V. Lostanlen, A. Farnsworth, J. Salamon _et al._ , “Chirping up the right tree: Incorporating biological taxonomies into deep bioacoustic classifiers,” in _Proc. ICASSP_ , May 2020, pp. 901–905.
  * [28] O. Ganea, G. Bécigneul, and T. Hofmann, “Hyperbolic neural networks,” in _Proc. NeurIPS_ , vol. 31, 2018.
  * [29] J. Mukhoti and Y. Gal, “Evaluating bayesian deep learning methods for semantic segmentation,” _arXiv preprint arXiv:1811.12709_ , 2018.
  * [30] W. Peng, T. Varanka, A. Mostafa, H. Shi _et al._ , “Hyperbolic deep neural networks: A survey,” _IEEE Trans. Pattern Anal. Mach. Intell._ , 12 2021\.
  * [31] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, “LibriSpeech: An ASR corpus based on public domain audio books,” in _Proc. ICASSP_ , Apr. 2015, pp. 5206–5210.
  * [32] E. Manilow, G. Wichern, P. Seetharaman, and J. Le Roux, “Cutting music source separation some Slakh: A dataset to study the impact of training data quality and quantity,” in _Proc. WASPAA_ , 2019.
  * [33] H. Kasai, P. Jawanpuria, and B. Mishra, “Riemannian adaptive stochastic gradient algorithms on matrix manifolds,” in _Proc. ICML_ , Jun. 2019, pp. 3262–3271.
  * [34] M. Kochurov, R. Karimov, and S. Kozlukov, “Geoopt: Riemannian optimization in PyTorch,” _arXiv preprint arXiv:2005.02819_ , 2020.
  * [35] J. Le Roux, S. Wisdom, H. Erdogan, and J. R. Hershey, “SDR–half-baked or well done?” in _Proc. ICASSP_ , May 2019.
  * [36] Z.-Q. Wang, J. Le Roux, and J. R. Hershey, “Alternative objective functions for deep clustering,” in _Proc. ICASSP_ , 2018, pp. 686–690.
  * [37] Z.-Q. Wang, J. Le Roux, D. Wang, and J. R. Hershey, “End-to-end speech separation with unfolded iterative phase reconstruction,” in _Proc. Interspeech_ , Sep. 2018.
  * [38] M. Nickel and D. Kiela, “Poincaré embeddings for learning hierarchical representations,” in _Proc. NeurIPS_ , vol. 30, 2017.
  * [39] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” _IEEE Trans. Audio, Speech, Lang. Process._ , vol. 14, pp. 1462 – 1469, Aug. 2006.
  * [40] Y. Gal and Z. Ghahramani, “Dropout as a Bayesian approximation: Representing model uncertainty in deep learning,” in _Proc. ICML_ , Jun. 2015.
  * [41] P. Nazreen and A. G. Ramakrishnan, “DNN based speech enhancement for unseen noises using Monte Carlo dropout,” in _Proc. ICSPCS_ , 2018, pp. 1–6.


# Three Cases of Optical Periodic Modulation in Active Galactic Nuclei

Jie Li,1 Zhongxiang Wang,1,2 and Dong Zheng1
1Department of Astronomy, School of Physics and Astronomy, Yunnan University,
Kunming 650091, China
2Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan
Road, Shanghai 200030, China E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We report on the case of optical periodic modulation discovered in two Active
Galactic Nuclei (AGN) and one candidate AGN. Analyzing the archival optical
data obtained from large transient surveys, namely the Catalina Real-Transient
Survey (CRTS) and the Zwicky Transient Facility (ZTF), we find periodicities
of 2169.7, 2103.1, and 1462.6 day in sources J0122+1032, J1007+1248 (or PG
1004+1248), and J2131$-$1127, respectively. The optical spectra of the first
two indicate that the first is likely a blazar and the second a type 1 Seyfert
galaxy, and while no spectroscopic information is available for the third one,
its overall properties suggest that it is likely an AGN. In addition, mid-
infrared (MIR) light curve data of the three sources, taken by the Wide-field
Infrared Survey Explorer (WISE), are also analyzed. The light curves show
significant variations, but not appearing related to the optical
periodicities. Based on the widely-discussed super-massive black hole binary
(SMBHB) scenario, we discuss the origin of the optical modulation. Two
possible interesting features, an additional 162-day short optical periodicity
in J2131$-$1127 and the consistency of the X-ray flux variations of J1007+1248
with its optical periodicity, are also discussed within the SMBHB scenario.

###### keywords:

Quasars:supermassive black holes

††pubyear: 2022††pagerange: Three Cases of Optical Periodic Modulation in
Active Galactic Nuclei–A

##  1 Introduction

Enabled by the availability of huge amounts of data collected from large
survey programs conducted in recent years, searches and studies of long-term
variabilities of different types of sources become possible. The Active
Galactic nuclei (AGN) are naturally the focus of these studies, as their
variability timescales can be long, driven by activities of the super-mass
black holes (SMBHs) at their center. Related phenomena include the so-called
quasi-periodic oscillations (QPOs), a fraction of which have year-long periods
(e.g., see Zhang & Wang 2021 and references therein) and have been discussed
to reflect the binary SMBH nature for the reported AGN systems (see, e.g.,
Valtonen et al. 2008; Ackermann et al. 2015; Sobacchi et al. 2017). Another
similar and probably more interesting phenomenon is the optical periodicities,
revealed in large part by recent large optical transient surveys. In the
Catalina Real-time Transient Survey (CRTS; Drake et al. 2009), 111 quasars
were reported to show candidate periodicity signals in 9-yr long light curves
(Graham et al., 2015a), in addition to the quasar PG 1302$-$102 that was found
to show a $\sim$1884 day periodicity (Graham et al., 2015b). By systematically
analyzing the Medium Deep Survey data obtained from the Panoramic Survey
Telescope and Rapid Response System (Pan-STARRS), Liu et al. (2015) and Liu et
al. (2019) respectively reported a possible periodic signal of $\sim$542 day
from a luminous radio-loud quasar (PSO J334.2028+01.4075) and 26 candidates of
similar periodic patterns. Also, by searching through data from the Palomar
Transient Factory (PTF), Charisi et al. (2016) reported 50 candidates with
significant periodicities, 33 of which remained significant when the CRTS data
were added in the analysis.

These long-term periodicities have been widely discussed as signals indicating
the existence of SMBH binaries (SMBHBs) situated at the center of AGN.
According to hierarchical models of galaxy formation and evolution, galaxy
mergers are common in their formation history (e.g., De Lucia et al. 2006),
and as a result, a significant fraction of galaxies would contain SMBHBs at
their center (e.g., Volonteri et al. 2003). Such a binary system would go
through complicated evolutionary processes, ending with coalescence driven by
gravitational radiation (see, e.g., Colpi 2014 for a review). Over the course
of the evolution, the dominant time period ($\sim$Gyr) the SMBHB spends at
would be when the two SMBHs have a sub-parsec separation distance (i.e.,
hardening phase; e.g., Yu 2002; Haiman et al. 2009; Colpi 2014). Depending on
the masses of the binary system, the corresponding orbital period could be
several years. Therefore, the periodicity signals detected from the surveys
match the theoretical expectations, although how the signals are produced is
not totally clear (Graham et al., 2015b). Different numerical simulations for
the SMBHB systems have shown that each black hole in such a binary, surrounded
by a circumbinary disc, would have a so-called mini-disc due to the accretion
supplied by the mass transfer from the circumbinary disc. The accretion rate
of the black holes would be modulated at the orbital period (e.g., Hayasaki et
al. 2007; MacFadyen & Milosavljević 2008; Cuadra et al. 2009; Farris et al.
2014). The presence of periodic flux variations would thus be related to the
orbital modulation of the accretion rate (Haiman et al., 2009).

The Zwicky Transient Facility (ZTF) survey, starting from 2018 March,
represents a more powerful optical transient survey due to its high cadence in
covering a given sky and fast data releases (Bellm et al., 2019). We explored
its released data and tested to search for periodic signals from AGN. In
essence, we queried the ZTF public database for AGN light curves using the
Automatic Learning for the Rapid Classification of Events (ALeRCE; Förster et
al. 2021; Sánchez-Sáez et al. 2021), for which we required $\geq$100 data
points in any light curve to ensure a clear view of any variation patterns. We
went over the output light curves by eye and selected those showing long-term
trends in the $\sim$4 yr ZTF-survey time period. We then checked the CRTS
light curve data for possible connectable patterns to the ZTF long-term
trends. In this way, we were able to find approximately 30 AGN sources that
showed potential long-term periodic variations in an approximately 6000-day
time period, set by the CRTS and ZTF data. Among them, there are three good
cases. Possibly related variations for the three cases were also seen in the
mid-infrared (MIR) data obtained with the Wide-field Infrared Survey Explorer
(WISE; Wright et al. 2010). Here we report on our findings of the three cases.
One is likely a blazar and the other two are a known AGN and a candidate AGN,
previously not known with the periodicity phenomenon. In Section 2, we
describe the known properties of the three sources and the available data for
them. In Section 3, the periodicity analysis method is described, and the
analyzing results are presented. We discuss and summarize the results in
Section 4.

Figure 1: Optical spectrum of J0122+1032 from the SDSS database. The absence
of any line features strongly suggests that the source is a blazar.

##  2 Sources and data

###  2.1 J0122+1032

Little information is available for J0122+1032 (hereafter J0122), but there is
an optical spectrum, taken on 2018 Nov. 29 (MJD 58451), in the Sloan Digital
Sky Survey (SDSS) database (Figure 1). The spectrum shows no obvious features,
which indicates the source is likely a BL Lac type blazar.

The source’s coordinates in ZTF, as well as the object ID (oid), are provided
in Table 1. We used ZTF’s light curve data (Masci et al., 2019). In order to
ensure the cleanness and goodness of the data, we set catflags=0 and chi$<$4
when querying (the same requirements were set for the ZTF data of the other
two sources), where the first parameter flags the photometric/image quality
and the second is the root-mean-square metric from point-spread-function fit
(see The ZTF Science Data System Explanatory
Supplement111https://www.ztf.caltech.edu/ztf-public-releases.html). The light
curve data at the ZTF bands $zg$ and $zr$ (Bellm et al., 2019) are in a
magnitude range of 19–20, and the time range of the data is given in Table 1.

We also obtained the source’s CRTS $V$ band data, for which the object IDs in
CRTS and time range are given in Table 1. The magnitude range is between
19–20.

The source was detected with WISE at the W1 (3.4 $\mu$m) and W2 (4.6 $\mu$m)
bands in the NEOWISE Post-Cryo survey phase. We downloaded the magnitude
measurements in the two bands from the NEOWISE-R Single-exposure Source
database.

We searched for X-ray observations that might cover the source, but did not
find any in the database of the major X-ray telescopes. However, we note a
$\gamma$-ray source 4FGL J0122.4+1034, detected with the Large Area Telescope
(LAT) onboard the Fermi Gamma-ray Space Telescope (Fermi), has a position
$\simeq$2.2$\arcmin$ away from that of J0122, while its 95%-confidence error
circle is $\simeq$2.5$\arcmin$. Thus this $\gamma$-ray source could be
associated with J0122. It is faint at $\gamma$-rays, having a test statistic
(TS) value of 32.2 (implies a $\sim 5\sigma$ detection significance). Its
emission is described with a power law (the photon index $\Gamma=1.91\pm
0.21$). This source is not identified with a type or associated with any known
source in the Fermi LAT catalog (Abdollahi et al., 2022). We tested to
construct a 90-day binned light curve, from which flare-like events might be
found and the blazar nature would be identified (e.g., Wang et al. 2020) since
the dominant sources detected with Fermi LAT are blazars (Abdollahi et al.,
2022). However no obvious flaring-type events were seen. Whether the
$\gamma$-ray source is associated with J0122 remains to be investigated.

Figure 2: Left panels: Optical light curve of J0122+1032. The WISE W1 (down-
shifted by 2 mag for clarity) and W2 band measurements are also shown in the
top panel. The normalized optical light curve is shown in the bottom panel,
which can be well described with a sinusoidal fit (gray curve). Right panel:
Sinusoidal fits to the WISE W1 and W2 band light curves (where the W1 band is
down-shifted by 2 mag for clarity), of which the green curves are the ones
with the period fixed at the optical value. The blue dashed curve is the model
fit to the optical light curve, shown for comparison.

Figure 3: Same as Figure 2 for J1007+1248. The W2 band light curve in the
right panel is down-shifted by 0.5 mag for clarity.

Figure 4: Same as Figure 2 for J2131$-$1127\. The W1 band light curve in the right panel is down-shifted by 0.6 mag for clarity.  Table 1: CRTS, ZTF, and WISE data for the three AGN Source | Coordinate | CRTS ID | ZTF oid | CRTS | ZTF | WISE
---|---|---|---|---|---|---
| (RA, Dec) |  |  | (MJD) | (MJD) | (MJD)
J0122+1032 | (01:22:23.63, +10:32:13.30) | MLS_J012223.5+103213 | 502108400001209 | 53680–56636 | 58288–59550 | 56668–59411
|  |  | 502208400002496 |  |  |
J1007+1248 | (10:07:26.11, +12:48:56.20) | CSS_J100726.1+124856 | 1566102100003309 | 53527–56596 | 58202–59550 | 56794–59545
|  | MLS_J100726.1+124856 | 520110100002864 |  |  |
|  |  | 520210100007454 |  |  |
|  |  | 1566202100004632 |  |  |
J2131$-$1127 | (21:31:06.96, $-$11:27:25.20) | CSS_J213107.1$-$112724 | 390107400004429 | 53528–56570 | 58257–59543 | 56790–59511
|  | MLS_J213107.1$-$112725 | 390207400009360 |  |  |
|  | SSS_J213107.0$-$112724 |  |  |  |

###  2.2 J1007+1248

The source J1007+1248 (hereafter J1007), known as PG 1004+130 (or 4C 13.41),
is a radio loud quasar (Wills et al., 1999) and has been relatively well
studied. There is an optical spectrum of it, taken on 2004 Feb. 20, in the
SDSS database (SDSS J100726.10+124856.2) and it has been identified as a type
1 Seyfert galaxy (e.g., Toba et al. 2014). It has redshift $z\simeq 0.24$, and
based on the full width at half maximum of the broad H$\beta$ line and the
optical continuum luminosity, was estimated to contain a SMBH with mass
$M_{t}=1.87\times 10^{9}\,M_{\odot}$ (Vestergaard & Peterson, 2006). Scott et
al. (2015) conducted detailed multi-epoch studies of the source at X-ray and
ultraviolet (UV) wavelengths, and found significant variations for the X-ray
flux and broad UV absorption lines.

This source is relatively bright at the optical and MIR wavelengths, with ZTF
and CRTS magnitudes around 15 and WISE W1 and W2 magnitudes between 10–12.
Same as the above for J0122, the light curve data were obtained from these
surveys’ database.

There were seven X-ray observations of the source field, two each with
Chandra, XMM-Newton, and Neil Gehrels Swift Observatory (Swift), and one with
Nuclear Spectroscopic Telescope Array (NuSTAR). Scott et al. (2015) analyzed
the Chandra and XMM-Newton data in detail, and Luo et al. (2013) reported the
analysis of the NuSTAR data. We used their results and derived the unabsorbed
fluxes in 0.5–10 keV for the source, in which the tool PIMMS was used. We
analyzed the Swift data and obtained fluxes in the same energy range; the
details of the analysis are provided in Section A. The X-ray flux results are
summarized in Table 5.

###  2.3 J2131$-$1127

This source was detected by CRTS, ZTF, and WISE. The CRTS light curve is
between 16–17 mag, and the ZTF $zr$ is within the same range but the $zg$ is
notably 1-mag larger. The WISE W1 and W2 magnitudes are around 13 and 12
respectively. The information for the data is given in Table 1.

There is no identification information for this source, but it was detected by
the Two Micron All Sky Survey in the infrared $JHK_{s}$ (2MASS; Skrutskie et
al. 2006) and by WISE in the W1–W4 bands. Both its MIR colors
(W1$-$W2$\simeq$1.1, W2$-$W3$\simeq$2.9; see, e.g., Wright et al. 2010) and
infrared colors (such as $J-K_{s}\simeq 0.77$, $K_{s}-$W3$\simeq 5.5$; see,
e.g., Tu & Wang 2013) put the source in the quasar category.

There were six usable Swift XRT observations conducted in 2019 covering the
field of J2131$-$1127 (hereafter J2131). The exposure times were short, from
178 sec to 1373 sec, and the source was not detected in any of them. We used
the online tools222https://www.swift.ac.uk/user_objects/ for the data
analysis. Using the data with the longest exposure (1373 sec), we derived an
unabsorbed flux upper limit of 6.2$\times 10^{-13}$ erg s-1 cm-2 in 0.3–10
keV, for which a power law was assumed with the hydrogen column density and
photon index fixed at the Galactic value 4.6$\times 10^{20}$ cm2 (HI4PI
Collaboration et al., 2016) and 1.0, respectively.

There is also a $\gamma$-ray source, 4FGL J2131.4$-$1124 in the Fermi LAT
catalog, possibly associated with J2131. The latter is $\simeq 5\arcmin$ away
from the former and within the 95%-confidence error circle of the former
($\sim 9\arcmin$). The $\gamma$-ray source (unidentified or unassociated in
the Fermi LAT catalog; Abdollahi et al. 2022) is faint, having a total TS
value of 46.4, and has emission described with a power law ($\Gamma=2.81\pm
0.15$). Similar to that for J0122, we constructed a 90-day binned light curve
of the source, but no obvious variability was seen in the light curve.

Based on the source’s infrared colors and its possible association with a
$\gamma$-ray source, we suggest it as a candidate AGN. Observations for
verification are needed. We note that given its optical brightness, if it is
an AGN, we would detect its X-ray emission around a flux of $10^{-13}$ erg s-1
cm-2 (e.g., He et al. 2019).

##  3 Periodicity analysis and results

The optical light curves of J0122, J1007, and J2131, as well as their WISE W1
and W2 band magnitudes, are shown in Figures 2, 3, and 4, respectively. As can
be seen, the sources all appear to have clearly visible modulation at the
optical wavelengths, and even their WISE light curves show possibly related
variations. We studied the modulation in each optical light curve and
investigated the variations seen in the MIR bands. Below, we first describe
our periodicity analysis method (Section 3.1), and then provide the results
from our studies of the optical light curves (Section 3.2) and MIR light
curves (Section 3.3).

###  3.1 Periodicity analysis method

As the observed light curves of the three sources all show sinusoidal-like
modulation (Figures 2, 3, 4) and the previous cases, i.e., those reported such
as in Graham et al. (2015b) and Liu et al. (2015), were well described by a
simple sinusoidal function, we adopted the same model form to study the
periodicities. A sinusoidal function, $A\sin 2\pi(t-B)/P+C$, was used, where
$t$ is time, $P$ is the period, and $A$, $B$, and $C$ are the amplitude, the
starting time (or the phase), and the constant magnitude, respectively.
However, because the CRTS data are in $V$ band and ZTF data are in $zg$ and
$zr$ bands, possibly having different $A$ and $C$ values, we first normalized
each set of the light-curve data in a band with the sinusoidal function (to
obtain $A$ and $C$). The $\chi^{2}$ fitting was used to find the best fits.
The obtained results are given in Table 2, in which the reduced $\chi^{2}$
values (=$\chi^{2}$/DoF, where DoF is the degrees of freedom) are included to
indicate the goodness of the fits. As a check by comparison, we also provide
the $\chi^{2}/{\rm DoF}$ values when fitting each set of data with a constant
(given in the lines of only parameter $C$ in Table 2). The comparison shows
significant improvements. For example, we may use the Akaike Information
Criterion (AIC) to check the improvements quantitatively, $AIC=N\ln({\rm
RSS}/N)+(2k+1)$, where RSS is the residual sum of squares, $N$ is the number
of the data points, and $k$ is the parameter number (Banks & Joyner, 2017).
The $\Delta AIC$ value is $\simeq-126$ when we compare the fitting results
from the sinusoid model to those of the constant model for the CRTS data of
J0122, and similar results can be obtained for the other data sets given the
large numbers of the data points. The large $\Delta AIC$ values indicate high
signficances for the presence of the sinusoidal variations in the light curves
(e.g., Rueda et al. 2022 and references therein).

For the sinusoidal fitting, the $\chi^{2}/{\rm DoF}$ values are mostly in a
range of $\simeq$0.6–2.9. Here it should be noted that AGN generally show
stochastic flux variations at a level of $<10\%$ (e.g., Vanden Berk et al.
2004), which implies that a systematic uncertainty may be added to account for
the variations in our fitting and thus the $\chi^{2}/{\rm DoF}$ values could
be further reduced. However for J2131, the $\chi^{2}/{\rm DoF}$ values from
fitting the ZTF $zg$ and $zr$ data sets are large, $\simeq 8.8$ and 25
respectively (Table 2), the reason for which is discussed below in Section
3.2.

With the obtained $A$ and $C$, we normalized the light-curve sets of a source
by subtracting $C$ from it and dividing it by $A$. In this way, a normalized
light curve of a source was constructed, with an average of $\sim$0 and a
range of 1 to $-$1 (see the left bottom panels of Figures 2, 3, and 4). These
light-curve data were then fitted with a sinusoidal function again, to mainly
determine $P$ and $B$. Here we set $B$ as one parameter, since the three
optical bands have close wavelength values and we would not expect phase
shifts in these optical bands based on the current studies of SMBHB systems.
In this last step, the Markov Chain Monte Carlo (MCMC) code emcee (Foreman-
Mackey et al., 2013) was used.

Table 2: Initial fitting results to individual CRTS and ZTF light curves

Source | Data | $\chi^{2}/{\rm DoF}$ | DoF | $A$ | $C^{\dagger}$
---|---|---|---|---|---
J0122+1032 | CRTS-$V$ | 2.39 | 129 | 0.321 | 19.314
|  | 6.46 | 132 |  | 19.262
| ZTF-$zg$ | 1.86 | 264 | 0.489 | 19.540
|  | 8.06 | 267 |  | 19.306
| ZTF-$zr$ | 2.22 | 408 | 0.394 | 19.051
|  | 13.22 | 411 |  | 18.876
J1007+1248 | CRTS-$V$ | 0.59 | 174 | 0.168 | 15.061
|  | 3.57 | 177 |  | 15.145
| ZTF-$zg$ | 2.75 | 178 | 0.210 | 15.159
|  | 147.21 | 181 |  | 15.101
| ZTF-$zr$ | 2.92 | 258 | 0.192 | 15.042
|  | 128.33 | 261 |  | 15.007
J2131-1127 | CRTS-$V$ | 2.87 | 383 | 0.151 | 16.898
|  | 4.68 | 386 |  | 16.880
| ZTF-$zg$ | 8.77 | 148 | 0.366 | 17.541
|  | 117.93 | 151 |  | 17.485
| ZTF-$zr$ | 25.02 | 168 | 0.381 | 16.622
|  | 500.97 | 171 |  | 16.616

  * •

† Lines of only $C$ parameter are the fits to a constant.

###  3.2 Optical light curve fitting results

From fitting the long-term light curves, which have a total time length of
approximately 6,000 days, we obtained the periods of 2169.7, 2103.1, and
1462.6 day for J0122, J1007, and J2131, respectively. The other results are
given in Table 3. The modulation of each source is reasonably well described
by each sinusoidal fit (shown in the left bottom panels of Figures 2, 3, and
4), as the CRTS and ZTF light-curve parts for each source are smoothly
connected by the fit (note the two sets of data have a gap of $\sim$1600 day).
The values of $A$ and $C$ are close to be 1 and 0 respectively, indicating
there were no significant mis-fit problems in the light-curve normalization
step; there could be such problems since the ZTF data do not cover a full
cycle of a periodicity.

We note that the $\chi^{2}/{\rm DoF}$ value for J2131 is again large (=10.5).
Carefully examining the ZTF $zg$ and $zr$ light curve parts, wiggling around
the model fit can be seen (Figure 5). Interestingly, the residuals of the
light curve parts subtracted from the best fit suggest the existence of short-
term periodic modulation. We tested to fit the residuals with a sinusoidal
function, and obtained $P=162.17\pm 0.18$ day and $A=0.188\pm 0.004$ mag with
$\chi^{2}/{\rm DoF}=13.2$. This large value reflects the significant
deviations of some data points from the short-term model fit (particularly
visible at $\sim$MJD 59000). In any case, the reduced $\chi^{2}$ from the
long-term ($P=1462.6$ day) fit for this part of the data is 21.2, indicating
significant improvement for describing the data variations. As the total $zg$
and $zr$ data points are 324, $\Delta AIC\simeq-126$ is obtained when
comparing the model fit having the short periodicity included with that
without it. The residuals and the fitting results strongly suggest that there
is another short-term periodicity, in addition to the 1462.6-day long
periodicity.

Table 3: Results from fitting the normalized light curves of the three AGN. Source | $\chi^{2}/{\rm DoF}$ | $A$ | $B$ | $C$ |  $P$ (day)
---|---|---|---|---|---
J0122+1032 | 2.8 | $1.12_{-0.02}^{+0.02}$ | $693.2_{-209.8}^{+211.3}$ | $0.19_{-0.01}^{+0.01}$ | $2169.7_{-7.9}^{+7.7}$
J1007+1248 | 2.9 | $1.036_{-0.005}^{+0.005}$ | $954.3_{-159.8}^{+157.0}$ | $0.060_{-0.005}^{+0.005}$ | $2103.1_{-5.7}^{+5.8}$
J2131$-$1127 | 10.5 | $1.015_{-0.004}^{+0.004}$ | $1098.9_{-95.5}^{+93.3}$ | $-0.019_{-0.003}^{+0.003}$ | $1462.6_{-2.4}^{+2.4}$
Figure 5: Sinusoidal fit (with $P=1462.6$ day) to the ZTF $zg$ (green) and
$zr$ (red) light curve parts of J2131 (top panel) and the residuals (bottom
panel). Modulation with a period of $\simeq$162.2 day and an amplitude of
$\simeq$0.188 mag possibly exists in the light curve part, which is indicated
by the black curve in the bottom panel.

###  3.3 Studies of MIR variations

The WISE light curves of the three sources show variations, particularly clear
and visible for J1007 and J2131 (Figures 2, 3, and 4). We calculated the
reduced $\chi^{2}$ values of the W1 and W2 light curves of the sources by
fitting to a constant, and provided them in Table 4. The $\chi^{2}/{\rm DoF}$
values for J0122 are $\geq 2.0$, and for J1007 and J2131, the values are $\geq
5$. In particular, the light curves of J2131 show signs of possible periodic
modulation. We thus investigated the variations of the light curves by fitting
them with the same sinusoidal function given above, and the emcee was used.
From the fitting, when all parameters were set free, the periods obtained are
larger than those of the optical light curves (see Table 4). Especially for
J1007, the values from the W1 or W2 bands are nearly 10 times or 5 times that
of the optical one, which is because the whole light curves appear upward
curved. The results indicate that the data can not provide good constraints on
periodicities, and no values, similar to those optical ones, can be found.
However, the fitting results do suggest modulation in the light curves. For
example, when we compare the fitting results from the sinusoid model to those
of the constant model, for J0112 that has the smallest $\chi^{2}/{\rm DoF}$
values, $\Delta AIC$ are $\simeq-59$ and $\simeq-74$ in the W1 and W2 band
respectively.

We then tested to fix the periods at the values of the optical ones in our
fitting. The $\chi^{2}/{\rm DoF}$ values from the fitting became slightly
larger (but still smaller than those from a constant model), and the
amplitudes were changed to be smaller (see Table 4). Comparing the results
with those from before, when all parameters for the sinusoidal function were
set free, $\Delta AIC\simeq 20$ and 12, for example, for J0122 in the W1 and
W2 band respectively. This indicates the latter model is preferred. Similar
conclusions can be made for J1007 and J2131. Given the results of the
comparison, no periodicities similar to the optical ones have been found.

Table 4: Results from sinusoidal fitting to the WISE light curves of the three
AGN

Source (Data) | $\chi^{2}/{\rm DoF}$ | DoF |  $P$ (day) | $A$ | $B$ | $C^{\dagger}$
---|---|---|---|---|---|---
J0122+1032 (W1) | 2.4 | 183 | 2773.9 | 0.123 | 1148.3 | 14.752
(W1) | 2.7 | 184 | 2169.7∗ | 0.100 | 2730.8 | 14.751
(W1) | 3.4 | 186 |  |  |  | 14.782
(W2) | 1.3 | 183 | 3200.5 | 0.179 | 2228.8 | 14.259
(W2) | 1.4 | 184 | 2169.7∗ | 0.151 | 2954.1 | 14.279
(W2) | 2.0 | 186 |  |  |  | 14.380
J1007+1248 (W1) | 1.6 | 192 | 19656.3 | 1.631 | 72669.4 | 9.993
(W1) | 5.8 | 193 | 2103.1∗ | 0.067 | 3022.3 | 11.575
(W1) | 11 | 195 |  |  |  | 11.562
(W2) | 1.5 | 192 | 9256.5 | 0.273 | 111314.1 | 10.251
(W2) | 2.4 | 193 | 2103.1∗ | 0.048 | 3066.7 | 10.491
(W2) | 5.0 | 195 |  |  |  | 10.483
J2131$-$1127 (W1) | 10.4 | 187 | 2178.8 | 0.221 | 5710.0 | 13.292
(W1) | 16.7 | 188 | 1462.6∗ | 0.152 | 2425.1 | 13.256
(W1) | 26 | 190 |  |  |  | 13.283
(W2) | 7.1 | 187 | 2404.1 | 0.146 | 9953.7 | 12.231
(W2) | 9.0 | 188 | 1462.6∗ | 0.102 | 3866.0 | 12.216
(W2) | 13 | 190 |  |  |  | 12.247

  * •

∗ Periods are fixed at the best-fit values of the optical light curves.

  * •

† Lines of only $C$ parameter are the fits to a constant.

##  4 Discussion and Summary

By exploring the ZTF data, we have found three good cases of long-term
periodic modulation. These three cases arise from a blazar (J0122), a Seyfert
1 galaxy (J1007), and a candidate AGN (J2131). Combining with the CRTS data,
the periodicities were determined by fitting the CRTS plus ZTF light curves of
the three sources with a sinusoidal function. The periods determined for
J0122, J1007, and J2131 are $\simeq$2170, 2103, and 1463 day, respectively.
While the periods are long, the total lengths of the data for the sources are
$\sim$6000 day, covering relatively well the periodic modulation of each
source. Based on the light-curve fitting results, we consider that the
periodicities have been clearly revealed.

Such sinusoidal modulation seen in AGN has been widely discussed as evidence
for the existence of SMBHBs at the center of AGN. The modulation would reflect
the accretion rate changes according to numerical simulations for SMBHBs,
although how the observed optical modulation is produced is not totally clear
and under discussion (e.g., Graham et al. 2015b). One model that may work
considers the relativistic Doppler effect of emission from a close SMBHB
(D’Orazio et al., 2015), which has been applied to explain AGN periodicities,
quasi-periodicities, and aperiodic variations (e.g., Yan et al. 2018; Charisi
et al. 2018). Taking our cases as the example, the orbital periods
$P^{\prime}_{\rm yr}$ (in units of year) at the host galaxies would be $\sim
6/(1+z)$ or $\sim 4/(1+z)$, and the separation distances between two black
holes would be $\simeq$0.002 pc
$(1+q)^{1/3}(M_{1}/10^{8}\,M_{\odot})^{1/3}(P^{\prime}_{\rm yr})^{2/3}$, where
$M_{1}$ is the mass of the primary black hole and $q=M_{2}/M_{1}$ the mass
ratio of the secondary (mass $M_{2}$) to the primary. Assuming a circular
orbit, the velocity $v_{2}$ of the secondary is $v_{2}/c\simeq
0.1[1/(1+q)](M/10^{9}\,M_{\odot})^{1/3}(P^{\prime}_{\rm yr})^{-1/3}$, where
$c$ is the speed of light and $M=M_{1}+M_{2}$ is the total mass. The observed
flux $F$ would mainly arise from this secondary black hole as it captures most
of the matter coming from the circumbinary disk, and given its relativistic
velocity, $F$ would be orbitally modulated due to the Doppler effect with an
amplitude of $\delta F/F=(3-\alpha)(v_{2}/c)\sin i$, where $\alpha$ is the
power-law index for an emitted spectrum, and $i$ is the inclination angle of
the binary orbit (see details in D’Orazio et al. 2015).

Since only for J1007, there is information of its redshift and estimated black
hole mass, we checked if the model could explain its modulation. We de-
reddened its SDSS spectrum with the Galactic extinction $E(B-V)=0.033$
(Schlegel et al., 1998) and shifted the spectrum back to the rest wavelength
with $z=0.24$. Fitting the $V$, $zg$, and $zr$ (for the latter two, see Bellm
et al. 2019) wavelength ranges of the corrected spectrum, respectively, with a
power-law $f_{\nu}\sim\nu^{\alpha}$, the $\alpha$ values were estimated to be
$\simeq 0.88$, 0.04, and $-0.16$. We then assumed $M=M_{t}$ and inserted
$P\simeq 5.76$ yr (or $P^{\prime}=4.64$ yr), and obtained $\delta
F/F\simeq\\{0.16,0.22,0.23\\}\times[1/(1+q)](M/1.87\times 10^{9})^{1/3}\sin i$
in the three bands. The amplitudes approximately match the observed values
(Table 2), although the $q$ and $i$ terms can significantly decrease the
values if $q\sim 1$ or $i$ is small. We note that the model is rather
simplified, such that the optical continuum could consist of emission not only
from the secondary black hole (e.g., D’Orazio et al. 2015; Farris et al. 2015)
and there are possible large uncertainties in the estimation—for example, the
$\alpha$ values were estimated from a single spectrum which was taken before
the times of the CRTS data and long before those of the ZTF data. For J0122
and J2131, we can derive to find $M=10^{9}\,M_{\odot}(\delta
F/0.1F)^{3}[(1+q)/(3-\alpha)]^{3}(\sin i)^{-3}P^{\prime}_{\rm yr}$. By
inserting their $\delta F/F$ at $V$ band and $P=P^{\prime}(1+z)$ values into
the formula, we obtained $M\sim\\{2\times 10^{11}\,M_{\odot}$, $1.4\times
10^{10}\,M_{\odot}\\}\times[(1+q)/(3-\alpha)]^{3}(\sin i)^{-3}/(1+z)$ for the
first and the latter respectively. In order not to have too large a mass
(e.g., $\geq 10^{11}\,M_{\odot}$) for the SMBHB in J0122, $q$ should be close
to zero, $\alpha\sim 0$, and $i$ should not be small.

For J2131, in addition to the $\sim$4-yr long periodicity, modulation with a
short period of $\simeq 162$ day possibly also existed. This phenomenon, if
verified with continuous observations, would be an intriguing feature to
investigate. We suspect that it could be explained with a hot spot at the
outer edge of the mini-disc around the secondary black hole. The gravitational
radius is $r_{g}=GM_{2}/c^{2}\simeq 1.5\times 10^{13}$ cm for a
10${}^{8}\,M_{\odot}$ black hole. At radius $r\simeq 270r_{g}$ (or $\simeq
0.001$ pc), the Keplerian orbital period is $\simeq 162$ day. As the AGN disc
sizes are found to be mostly in a range of $10^{-3}$–$10^{-2}$ pc (possibly
depending on the black hole mass; e.g., Jha et al. 2022), a hot spot that
orbits the SMBH and is located at the outer edge of the SMBH’s disc could give
rise to the short period signal. Such a hot spot could be produced by the
interaction between the stream of matter from the circumbinary disk with the
outer edge of the mini-disc (e.g., Hayasaki et al. 2007; Roedig et al. 2014).
We realize that there could be alternative explanations for the presence of
two periodicities in J2131, such as the short periodicity being a jet’s
wobbling timescale (Liska et al., 2018). As such discussion would deviate away
from the focus of this work, we refer to Zhang et al. (2022) for the
discussion of alternative scenarios.

Figure 6: Unabsorbed 0.5–10 keV fluxes of J1007 (diamonds: XMM-Newton;
circles: Chandra; square: NuSTAR; triangles: Swift). The flux variations are
possibly consistent with the optical periodicity (blue line).

Finally, Scott et al. (2015) noted that J1007 showed upto 40% X-ray flux
variations, and the variations were likely intrinsic because no significant
variations in the amount of X-ray absorption were seen. We plotted all the
unabsorbed 0.5–10 keV fluxes from Table 5 in Figure 6, and tested to fit the
fluxes with a sinusoidal function, but with the parameters $P$ and $B$ fixed
at the values of the determined optical periodicity. As shown in Figure 6, the
X-ray flux variations are possibly consistent with the optical periodicity,
which may suggest that the X-ray flux variations are actually caused by the
presence of an SMBHB (Serafinelli et al., 2020); the accretion from the
circumbinary disc to the black holes is orbitally modulated. This possibility
predicts very certain X-ray flux variations and can be verified with further
Chandra or XMM-Newton observations.

As a summary, we have discovered three cases of optical modulation in AGN, and
the periodicities are likely caused by the widely-discussed SMBHB scenario.
For J1007 that has redshift and estimated black hole mass, we have found that
the modulation could be explained with the Doppler effect of emission from a
secondary SMBH with a relativistic orbital velocity. We have also analyzed the
MIR light curve variations for the three AGN and found that J2131 was the more
variable one, but whether the variations are related to the optical ones is
not clear. For J2131, its optical light curve showed additional variations
around its long-term modulation, for which a short periodicity was estimated,
and may be explained as the presence of a hot spot at the outer edge of the
mini-disc of the secondary black hole. Also for J1007, which has been observed
multiple times at X-rays, its X-ray flux varied possibly with its optical
periodicity. The verification of this possibility would strengthen its case as
a candidate SMBHB and suggest an accompanying window at X-rays for finding
candidate SMBHBs (Serafinelli et al., 2020). As several large optical
transient survey programs are or will be in operation, more data will
certainly be collected for the three sources. We will be able to keep checking
the persistence of the periodicities, which could potentially establish them
as the good SMBHB candidates.

## Acknowledgements

This work was based on observations obtained with the Samuel Oschin Telescope
48-inch and the 60-inch Telescope at the Palomar Observatory as part of the
Zwicky Transient Facility project. ZTF is supported by the National Science
Foundation under Grant No. AST-2034437 and a collaboration including Caltech,
IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm
University, the University of Maryland, Deutsches Elektronen-Synchrotron and
Humboldt University, the TANGO Consortium of Taiwan, the University of
Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National
Laboratories, and IN2P3, France. Operations are conducted by COO, IPAC, and
UW.

This work made use of data products from the Wide-field Infrared Survey
Explorer, which is a joint project of the University of California, Los
Angeles, and the Jet Propulsion Laboratory/California Insitute of Technology,
funded by the National Aeronautics and Space Administration.

We thank the anonymous referee for very helpful comments. Y. Luo has helped
our understanding of SMBHB evolution scenarios, and M. Gu about the origin of
the AGN optical variability. This research is supported by the Basic Research
Program of Yunnan Province No. 202201AS070005, the National Natural Science
Foundation of China No. 12273033, and the Original Innovation Program of the
Chinese Academy of Sciences (E085021002).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abdollahi et al. (2022) Abdollahi S., Acero F., Baldini L., et al. 2022 ApJS 260, 53
  * Ackermann et al. (2015) Ackermann M., et al., 2015, ApJ, 813, L41
  * Banks & Joyner (2017) Banks H. T., Joyner M. L., 2017, Appl. Math. Lett., 74, 33
  * Bellm et al. (2019) Bellm E. C., et al., 2019, PASP, 131, 018002
  * Charisi et al. (2016) Charisi M., Bartos I., Haiman Z., Price-Whelan A. M., Graham M. J., Bellm E. C., Laher R. R., Márka S., 2016, MNRAS, 463, 2145
  * Charisi et al. (2018) Charisi M., Haiman Z., Schiminovich D., D’Orazio D. J., 2018, MNRAS, 476, 4617
  * Colpi (2014) Colpi M., 2014, Space Sci. Rev., 183, 189
  * Cuadra et al. (2009) Cuadra J., Armitage P. J., Alexander R. D., Begelman M. C., 2009, MNRAS, 393, 1423
  * D’Orazio et al. (2015) D’Orazio D. J., Haiman Z., Schiminovich D., 2015, Nature, 525, 351
  * De Lucia et al. (2006) De Lucia G., Springel V., White S. D. M., Croton D., Kauffmann G., 2006, MNRAS, 366, 499
  * Drake et al. (2009) Drake A. J., et al., 2009, ApJ, 696, 870
  * Evans et al. (2007) Evans P. A., et al., 2007, A&A, 469, 379
  * Evans et al. (2009) Evans P. A., et al., 2009, MNRAS, 397, 1177
  * Evans et al. (2020) Evans P. A., et al., 2020, ApJS, 247, 54
  * Farris et al. (2014) Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2014, ApJ, 783, 134
  * Farris et al. (2015) Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2015, MNRAS, 446, L36
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Förster et al. (2021) Förster F., Cabrera-Vives G., Castillo-Navarrete E., et al. 2021, AJ, 161, 242
  * Graham et al. (2015a) Graham M. J., et al., 2015a, MNRAS, 453, 1562
  * Graham et al. (2015b) Graham M. J., et al., 2015b, Nature, 518, 74
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration et al., 2016, A&A, 594, A116
  * Haiman et al. (2009) Haiman Z., Kocsis B., Menou K., 2009, ApJ, 700, 1952
  * Hayasaki et al. (2007) Hayasaki K., Mineshige S., Sudou H., 2007, PASJ, 59, 427
  * He et al. (2019) He L., et al., 2019, Research in Astronomy and Astrophysics, 19, 098
  * Jha et al. (2022) Jha V. K., Joshi R., Chand H., Wu X.-B., Ho L. C., Rastogi S., Ma Q., 2022, MNRAS, 511, 3005
  * Liska et al. (2018) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S., 2018, MNRAS, 474, L81
  * Liu et al. (2015) Liu T., et al., 2015, ApJ, 803, L16
  * Liu et al. (2019) Liu T., et al., 2019, ApJ, 884, 36
  * Luo et al. (2013) Luo B., et al., 2013, ApJ, 772, 153
  * MacFadyen & Milosavljević (2008) MacFadyen A. I., Milosavljević M., 2008, ApJ, 672, 83
  * Masci et al. (2019) Masci F. J., et al., 2019, PASP, 131, 018003
  * Roedig et al. (2014) Roedig C., Krolik J. H., Miller M. C., 2014, ApJ, 785, 115
  * Rueda et al. (2022) Rueda, H., Glicenstein J., Brun F., 2022, ApJ, 934, 6
  * Sánchez-Sáez et al. (2021) Sánchez-Sáez, P., Reyes, I., Valenzuela, C., et al., 2021, AJ, 161, 141
  * Schlegel et al. (1998) Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525
  * Scott et al. (2015) Scott A. E., Brandt W. N., Miller B. P., Luo B., Gallagher S. C., 2015, ApJ, 806, 210
  * Serafinelli et al. (2020) Serafinelli R., et al., 2020, ApJ, 902, 10
  * Skrutskie et al. (2006) Skrutskie M. F., et al., 2006, AJ, 131, 1163
  * Sobacchi et al. (2017) Sobacchi E., Sormani M. C., Stamerra A., 2017, MNRAS, 465, 161
  * Toba et al. (2014) Toba Y., et al., 2014, ApJ, 788, 45
  * Tu & Wang (2013) Tu X., Wang Z.-X., 2013, Research in Astronomy and Astrophysics, 13, 323
  * Valtonen et al. (2008) Valtonen M. J., et al., 2008, Nature, 452, 851
  * Vanden Berk et al. (2004) Vanden Berk D. E., et al., 2004, ApJ, 601, 692
  * Vestergaard & Peterson (2006) Vestergaard M., Peterson B. M., 2006, ApJ, 641, 689
  * Volonteri et al. (2003) Volonteri M., Haardt F., Madau P., 2003, ApJ, 582, 559
  * Wang et al. (2020) Wang G., Wang Z., Chen L., Zhou J., Xing Y., 2020, PASJ, 72, 9
  * Wills et al. (1999) Wills B. J., Brandt W. N., Laor A., 1999, ApJ, 520, L91
  * Wright et al. (2010) Wright E. L., et al., 2010, AJ, 140, 1868
  * Yan et al. (2018) Yan D., Zhou J., Zhang P., Zhu Q., Wang J., 2018, ApJ, 867, 53
  * Yu (2002) Yu Q., 2002, MNRAS, 331, 935
  * Zhang & Wang (2021) Zhang P., Wang Z., 2021, ApJ, 914, 1
  * Zhang et al. (2022) Zhang P., Wang Z., Gurwell M., Wiita P. J., 2022, ApJ, 925, 207

##  Appendix A X-ray observations of J1007+1248

The two Swift observations of the source field was conducted with XRT in
photon counting mode. We used the online
tools333https://www.swift.ac.uk/user_objects/ for the data analysis and source
detection (Evans et al., 2020). In the first observation, the source was only
found to have a $\geq 3\sigma$ detection significance by comparing its count
rate (4.7${}^{+3.1}_{-2.2}\times 10^{-3}$ cts s-1) with that of the expected
background (Evans et al., 2007, 2009). Its unaborbed flux in 0.5–10 keV was
estimated using the PIMMS tool, where the Galactic $N_{\rm H}$ value
3.5$\times 10^{20}$ cm-2 (HI4PI Collaboration et al., 2016) and photon index
1.0 were assumed. In the second observation, only 6 photons from the source
were detected. A spectrum was created following Evans et al. (2009) and fitted
with a power-law model (in which the same Galactic $N_{\rm H}$ value was
used). The corresponding 0.5–10 keV unabsorbed flux was obtained. The fluxes
from these two observations are given in Table 5.

Table 5: X-ray observations of J1007+1248 and flux measurement results Observation | Date | Obsid | Exposure | $F^{\rm unabs}_{0.5-10}/10^{-13}$
---|---|---|---|---
|  |  | (ks) | (erg cm-2 s-1)
Chandra1 | 2005-01-05 | 5606 | 41 | 6.45$\pm$0.48
Chandra2 | 2014-03-20 | 16034 | 61 | 4.45$\pm$0.37
XMM-Newton1 | 2003-05-04 | 0140550601 | 22 | 4.25$\pm$0.31
XMM-Newton2 | 2013-11-05 | 0728980201 | 66 | 3.51$\pm$0.15
Swift1 | 2012-10-29 | 00080031001 | 2.0 | $2.4^{+1.6}_{-1.1}$
Swift2 | 2014-01-10 | 00033077001 | 1.9 | 1.6${}^{+2.1}_{-1.0}$
NuSTAR | 2012-10-29 | 60001112002 | 33 | $3.26\pm 0.42$