# Understanding the phenomenological and intrinsic blazar sequence using a
simple scaling model

Zhu-Jian Wan,1 Rui Xue,1 Ze-Rui Wang,2 Hu-Bing Xiao,3 and Jun-Hui Fan4,5,6
1Department of Physics, Zhejiang Normal University, Jinhua 321004, China
2College of Physics and Electronic Engineering, Qilu Normal University, Jinan
250200, China
3Shanghai Key Lab for Astrophysics, Shanghai Normal University, Shanghai
200234, People’s Republic of China
4Center for Astrophysics, Guangzhou University, Guangzhou 510006, People’s
Republic of China
5Key Laboratory for Astronomical Observation and Technology of Guangzhou,
Guangzhou, 510006, People’s Republic of China
6Astronomy Science and Technology Research Laboratory of Department of
Education of Guangdong Province, Guangzhou, 510006, People’s Republic of China
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The blazar sequence, including negative correlations between radiative
luminosity $L_{\rm rad}$ and synchrotron peak frequency $\nu$, and between
Compton dominance $Y$ and $\nu$, is widely adopted as a phenomenological
description of spectral energy distributions (SEDs) of blazars, although its
underlying cause is hotly debated. In particular, these correlations turn
positive after correcting Doppler boosting effect. In this work, we revisit
the phenomenological and intrinsic blazar sequence with three samples, which
are historical sample (SEDs are built with historical data), quasi-
simultaneous sample (SEDs are built with quasi-simultaneous data) and Doppler
factor corrected sample (a sample with available Doppler factors), selected
from literature. We find that phenomenological blazar sequence holds in
historical sample, but does not exist in quasi-simultaneous sample, and
intrinsic correlation between $L_{\rm rad}$ and $\nu$ becomes positive in
Doppler factor corrected sample. We also analyze if the blazar sequence still
exists in subclasses of blazars, i.e., flat-spectrum radio quasars and BL
Lacertae objects, with different values of $Y$. To interpret these
correlations, we apply a simple scaling model, in which physical parameters of
the dissipation region are connected to the location of the dissipation
region. We find that the model generated results are highly sensitive to the
chosen ranges and distributions of physical parameters. Therefore, we suggest
that even though the simple scaling model can reproduce the blazar sequence
under specific conditions that have been fine-tuned, such results may not have
universal applicability. Further consideration of a more realistic emission
model is expected.

###### keywords:

radiation mechanisms: non-thermal – galaxies: active – galaxies: jets.

††pubyear: 2024††pagerange: Understanding the phenomenological and intrinsic
blazar sequence using a simple scaling model–Understanding the
phenomenological and intrinsic blazar sequence using a simple scaling model

##  1 Introduction

Blazars are the most extreme active galactic nuclei (AGNs) with relativistic
jets aligned to our line of sight (Urry & Padovani, 1995). The spectral energy
distribution (SED) of blazars consists of two broad non-thermal emission
components in the $\log\nu-\log\nu L_{\nu}$ diagram. The low-energy component,
peaking between the infrared to X-ray bands, is believed to originate from the
synchrotron emission of relativistic electrons within the jet. The high-energy
component, peaking in the $\gamma$-ray band, is believed to originate from the
inverse Compton (IC) emission of the same electron population that up-scatter,
either the synchrotron photons emitted by the same population of relativistic
electrons (synchrotron self-Compton, SSC; Konigl, 1981) or external photons
(external Compton, EC; Ghisellini & Maraschi, 1989; Dermer et al., 1992;
Dermer & Schlickeiser, 1993; Sikora et al., 2002) from an accretion disc, a
broad-line region (BLR), or a dusty torus (DT). Based on the equivalent width
(EW) of the emission lines, blazars are divided into flat-spectrum radio
quasars (FSRQs) with strong broad emission lines (EW$>5$ Å), and BL Lacertae
objects (BL Lacs) with absent or weak emission lines (EW$\leqslant 5$ Å;
Marcha et al., 1996; Landt et al., 2004; Xiong & Zhang, 2014). Ghisellini et
al. (2011) introduce a physical distinction between these two subclasses based
on the BLR luminosity $L_{\rm BLR}$ in units of the Eddington luminosity
$L_{\rm Edd}$ of the central supermassive black hole (SMBH). FSRQs are objects
whose $L_{\rm BLR}\gtrsim 5\times 10^{-4}L_{\rm Edd}$, and BL Lacs are the
others. In addition to the FSRQs and BL Lacs classification, Abdo et al.
(2010a) divide blazars into high-synchrotron-peaked (HSP) blazars if their
synchrotron peak $\nu_{\rm s}^{\rm p}>10^{15}~{}\rm Hz$, intermediate-
synchrotron-peaked (ISP) blazars if their synchrotron peak $10^{14}~{}\rm
Hz<\nu_{\rm s}^{\rm p}<10^{15}~{}\rm Hz$ and low-synchrotron-peaked (LSP)
blazars if their synchrotron peak $\nu_{\rm s}^{\rm p}<10^{14}~{}\rm Hz$.

Based on a sample of 126 blazars, Fossati et al. (1998) identified that their
observed radiative luminosity, and the Compton dominance (the ratio of fluxes
or luminosities of high- and low-energy components) are both anti-correlated
to the observed synchrotron peak frequency. These two negative correlations
are referred to as the well-known “blazar sequence”. The first physical
explanation is proposed by Ghisellini et al. (1998). They suggested that if
the intrinsic jet power is connected to $L_{\rm BLR}$, one would expect the
synchrotron peak frequency shifts from high frequency to low frequency by
increasing the intrinsic jet power and increasing the radiative cooling
dominated by the EC scattering. If this interpretation is correct, it provides
a powerful tool for understanding the evolution of blazars’ SEDs, similar to
the evolution of stars on the Hertzsprung–Russel diagram. Theoretical models
suggest that the jet power is generated from accretion and the extraction of
rotational energy or angular momentum from the disc/black hole (Blandford &
Znajek, 1977; Blandford & Payne, 1982), therefore the jet properties and
emissions are linked to the black hole mass and accretion rate. Ghisellini &
Tavecchio (2008) further highlight that the black hole mass and accretion rate
are connected to the jet power and the shape of SED, providing a framework
that could help explain the existence of both low-luminosity ‘blue’ quasars
and low-luminosity ‘red’ quasars. The importance of the accretion rate is also
pointed out by Meyer et al. (2011). They find that AGNs in the blazar sequence
diagram is mainly composed of two groups, the primary difference between these
groups being the accretion rate. The group with a high accretion rate consists
mostly of LSP blazars and FR II radio galaxies, while the group with a low
accretion rate corresponds to the majority of BL Lacs and FR I radio galaxies.
In addition, many other studies endeavor to come up with physical
interpretations of the blazar sequence. For instance, Böttcher & Dermer (2002)
propose a simple analytical model, suggesting that the reduction of gas and
dust in external environment causes the evolutionary sequence from FSRQ to BL
Lacs-LSPs, and then to BL Lacs-HSPs. Björnsson (2010) propose that the blazar
sequence can be well understood in the case of multiple IC scatterings, and
the energy density of relativistic electrons is the main factor influencing
the formation of the blazar sequence. Finke (2013) construct a simple model to
explain the blazar sequence in the 2LAC clean sample, suggesting that the
difference between sources is associated with the magnetic field of the jet’s
dissipation region, the energy density of the external photon field, and the
jet’s viewing angle. Note that all sources in the modeling are assumed to have
the same bulk Lorentz factor. Potter & Cotter (2013b) propose a relativistic
continuous jet model, which suggests that the blazar sequence could be
interpreted if the radius of the transition region and bulk Lorentz factor of
the jet increase with the jet power. According to their interpretation, FSRQs
have large bulk Lorentz factors, and scatter external cosmic microwave
background photons at large distances, resulting in larger Compton dominance
and lower IC peak frequency. On the other hand, BL Lacs, which have lower jet
powers, possess stronger magnetic fields in the bright synchrotron transition
region compared to high-power FSRQs, leading to higher peak frequency
synchrotron and SSC emissions.

Since the discovery of the blazar sequence, its underlying cause and if it is
a result of selection effect have been widely discussed, but even so it is
still frequently adopted as a description of the blazar population. With the
improvement of detector sensitivity and the expansion of samples, the blazar
sequence has been confirmed and developed (Ghisellini & Tavecchio, 2008; Chen
& Bai, 2011; Finke, 2013; Xiong et al., 2015a; Ghisellini, 2016; Mao et al.,
2016; Ghisellini et al., 2017; Savard et al., 2022; Kerby & Falcone, 2023;
Ouyang et al., 2023). Finke (2013) suggests that the redshift selection effect
cannot explain the blazar sequence. In Ghisellini et al. (2017), a revised
version of the blazar sequence, named the “Fermi blazar sequence”, is propose
by using a clean sample of 747 blazars (299 BL Lacs and 448 FSRQs) from the
third catalogue of AGN detected by Fermi-LAT (Ackermann et al., 2015). In this
large blazar sample, the existence of the blazar sequence is confirmed.
However it is also pointed out that when considering BL Lacs and FSRQs
separately, the sequence still holds for BL Lacs but not for FSRQs. On the
other hand, the blazar sequence is questioned by the discovery of “outlier”
blazars (e.g., luminous HSP blazars) and the evidence of selection effects
(Nieppola et al., 2006; Padovani et al., 2012; Giommi et al., 2012b; Giommi et
al., 2012a; Giommi & Padovani, 2015; Raiteri & Capetti, 2016; Cerruti et al.,
2017; Keenan et al., 2021).

In the original physical interpretation (Ghisellini et al., 1998), radiative
cooling is suggested as the main cause of the two negative correlations in the
phenomenological blazar sequence. If so, these two negative correlations
should firstly hold in the comoving frame, therefore it has to be assumed that
all blazars have similar values of the Doppler factor, or the Doppler factor
is not related to other physical parameters (Ghisellini et al., 2017; Prandini
& Ghisellini, 2022). However, some studies argue that the phenomenological
blazar sequence is an artefact of the Doppler boosting, since these two
correlations with Doppler-corrected values turn positive or disappear (Landt &
Bignall, 2008; Nieppola et al., 2008; Wu et al., 2009; Xiong et al., 2015b;
Fan et al., 2016a, 2017; Chen et al., 2021; Yang et al., 2022b). It poses a
challenge to the physical understanding of the blazar sequence.

In this work, we are motivated to revisit the phenomenological and intrinsic
blazar sequence of samples in the literature, and propose a theoretical
interpretation to these correlations. In previous studies, the blazar sequence
is studied with historical SEDs. However, since blazars are highly variable
objects, we also check if the blazar sequence holds for the quasi-simultaneous
SEDs. This paper is organized as follows. In Sect. 2, we revisit the blazar
sequence with three samples. In Sect. 3, we propose a theoretical
consideration to understand the phenomenological and intrinsic blazar
sequence. Finally, we end with the conclusion in Sect. 4. Throughout the
paper, the cosmological parameters $H_{0}=69.6\ \rm km\ s^{-1}Mpc^{-1}$,
$\Omega_{0}=0.29$, and $\Omega_{\Lambda}$= 0.71 are adopted (Bennett et al.,
2014).

##  2 Revisiting the blazar sequence with the $Fermi$ blazar

The original blazar sequence is a phenomenological description of the blazars’
SED, including the negative correlations between the synchrotron peak
frequency and the radiative luminosity, and between the synchrotron peak
frequency and Compton dominance. In this section, we revisit these two
correlations both in the observers’ frame and in the comoving frame with
different samples.

When studying the correlation between the radiative luminosity and the
synchrotron peak frequency in the blazar sequence, radio luminosity was
originally used (e.g., Fossati et al., 1998) since it is less variable than
other bands luminosity (Giommi et al., 2012b) and can serve as a good
indicator of the radiative luminosity (Wang et al., 2017). However, it should
be noted that using radio luminosity may introduces two sources of bias: (i)
radio emission is generally believed to originate from the extended jet
because of the synchrotron self-absorption, while emissions from other bands
come from the inner jet, indicating that they originate from different
locations of the jet and have different physical origins (Ghisellini et al.,
2009; Ghisellini et al., 2010); (ii) blazars are highly variable objects, with
significant variations in both the luminosities and the peak frequencies of
the two humps (e.g., Massaro et al., 2004; Acciari et al., 2011; Xiong et al.,
2017; Xiong et al., 2020). Although radio luminosity may serve as a long-term
proxy for radiative luminosity, there is no quantity available to serve as a
long-term proxy for the synchrotron peak frequency. Therefore, we suggest
using the integrated full-band jet luminosity as a proxy for radiative
luminosity instead of radio luminosity, which can mitigate the first bias. The
bias due to variability is difficult to avoid, but this approach can at least
ensure that the full-band integrated luminosity and the synchrotron peak
frequency originate from the same snapshot of SED, making their physical
connection more intimate. Thanks to the abundant multi-wavelength observations
at present, a more accurate estimate of the radiative luminosity can be
obtained by fitting the full-band SED. To investigate the impact of
variability on correlations of the blazar sequence, it would be intriguing to
test if non-simultaneous and quasi-simultaneous SEDs exhibit similar behavior.

###  2.1 Samples

This work focuses on the blazar sequence and uses $Fermi$ blazars as samples,
since $\gamma$-ray luminosity is a good indicator of radiative luminosity
(Wang et al., 2017) and is crucial for calculating Compton dominance. We study
three samples selected from the literature using linear correlations:

  1. 1.

The “historical sample”, compiled by Yang et al. (2022a, 2023) from the Fourth
Fermi-LAT 12-year Source catalog (Abdollahi et al., 2022), contains 750 FSRQs
and 844 BL Lacs with measured redshifts. The historical SEDs of all blazars
are constructed, and a parabolic equation is used to fit them111The low-energy
component is fitted in Yang et al. (2022a) and the high-energy component is
fitted in Yang et al. (2023)., enabling the derivation of Compton dominance
$Y$, as well as synchrotron peak frequency $\nu^{\rm obs}$ and radiative
luminosity $L_{\rm rad}^{\rm obs}$ in observers’ frame.

  2. 2.

Xue et al. (2016) collect quasi-simultaneous multi-wavelength data of 279
blazars from the second Fermi-LAT AGN catalog (Ackermann et al., 2011). In
their sample, 81 FSRQs and 28 BL Lacs have the full-waveband quasi-
simultaneous SEDs, which are also fitted using a parabolic equation. This
sample is referred to as the “quasi-simultaneous sample” hereafter.

  3. 3.

By cross-matching the Fermi-4FGL blazars collected in Yang et al. (2022a,
2023) with the sample of Liodakis et al. (2018), Yang et al. (2022b) find 180
blazars, including 129 FSRQs and 51 BL Lacs, with available Doppler factors
$\delta_{\rm D}$ and measured redshifts. Hereafter, this sample is referred to
as the “$\delta_{\rm D}$-corrected sample”. Due to the beaming effect, the
synchrotron peak frequency $\nu$ and radiative luminosity $L_{\rm rad}$ in the
comoving frame can be obtained through $\nu\simeq\nu^{\rm obs}\delta_{\rm
D}^{-1}$ and $L_{\rm rad}\simeq L_{\rm rad}^{\rm obs}\delta_{\rm D}^{-4}$,
respectively. Please note that this $\delta_{\rm D}$-corrected sample is
highly biased, as there are only three sources with $Y<1$. This implies that
there are clear selection biases when using this sample for correlation
studies. If a more complete sample can be collected in the future, the
corresponding correlation results could be more reliable and meaningful.

In the above samples, the synchrotron peak frequency, the radiative
luminosity, and Compton dominance are all obtained by fitting SEDs with
parabolic equations. However, one may argue that fitting a non-standard
parabola SED with the parabolic equation may introduce bias in correlation
studies. In Xue et al. (2016), they respectively fit apparent asymmetrical
SEDs of some specific blazars using the parabolic equation and the cubic
polynomial, and find that the difference in obtained luminosities is
insignificant. In Yang et al. (2022a), they compare the synchrotron peak
frequency derived from parabolic equation with those derived from the cubic
polynomial of Abdo et al. (2010b). They find that the difference of the
average peak frequency obtained by different fitting methods is minor.
Therefore, we suggest that it is reasonable to use the parameters obtained by
the parabolic equation to revisit the blazar sequence.

In the following, we use historical and quasi-simultaneous samples to revisit
the phenomenological blazar sequence and compare their correlation results.
The $\delta_{\rm D}$-corrected sample enables us to study the intrinsic blazar
sequence. It should be noted that the obtained Doppler factor is based on the
variability brightness temperature of radio observations (Liodakis et al.,
2018). Although Finke (2019) finds that Doppler factors in parsec-scale jets
obtained from the conical jet model (Blandford & Königl, 1979) are basically
in agreement with those obtained from radio observations (Hovatta et al.,
2009), the agreement is within quite significant errors. It should be noted
that obtaining accurate $\delta_{\rm D}$ for the dominant dissipation region
is incredibly difficult, and even the values measured by radio observations at
different times for the same source can vary greatly. For example, for the
same sample of blazars, values of $\delta_{\rm D}$ obtained by Liodakis et al.
(2018) has a large discrepancy with those obtained by Hovatta et al. (2009)
and Liodakis et al. (2017). Therefore, using $\delta_{\rm D}$ measured by
radio observation to study the intrinsic blazar sequence will inevitably
introduce a large uncertainty and potential observational biases.

Figure 1: The correlations between $\log(Y)$ and $\log(\nu^{\rm obs})$ in the
historical sample (upper panel), and in the quasi-simultaneous sample (lower
panel). As explained in the inset legends, the red symbol represents the FSRQs
and the teal symbol represents the BL Lacs, the full symbol represents the
source with $Y>1$ and the empty symbol represents the source with $Y\leq 1$.
If a significant correlation is found ($\tau>0.1,~{}p<0.01$) for a specific
(sub-) sample statistically, the best linear fitting equation is also shown.
The solid black line represents the best-fitting equation for the whole
sample. The solid red line represents the best-fitting equation for the whole
population of FSRQs, the dashed red line represents the best-fitting equation
for the population of FSRQs with $Y>1$, and the dotted red line represents the
best-fitting equation for the population of FSRQs with $Y\leq 1$. The solid
teal line represents the best-fitting equation for the whole population of BL
Lacs, the dashed teal line represents the best-fitting equation for the
population of BL Lacs with $Y>1$, and the dotted teal line represents the
best-fitting equation for the population of BL Lacs with $Y\leq 1$.

###  2.2 Correlation Results

We analyze the phenomenological blazar sequence in both the historical and
quasi-simultaneous samples, and the intrinsic blazar sequence in $\delta_{\rm
D}$-corrected sample. We do not study the phenomenological blazar sequence in
$\delta_{\rm D}$-corrected sample further, since the $\delta_{\rm
D}$-corrected sample is a subsample of historical sample. The theoretical
analysis proposed in Sect. 3 indicates that FSRQs and BL Lacs with $Y\leq 1$
and $Y>1$ may have different correlations and slopes, so we further divide
three samples collected in this work using the dividing line $Y=1$. Table 1
shows the statistical correlation results, including slopes of the best linear
fitting equations (denoted by $s$ hereafter). Based on the Spearman rank
correlation test, our results are described as follows: correlation
coefficients between $0.10$ and $0.29$ indicate a weak correlation,
coefficients between $0.30$ and $0.49$ represent a moderate correlation, and
coefficients between $0.50$ and $1$ represent a strong correlation. If the
chance probabilities are $p>0.01$, it is suggested that the correlation is not
established statistically, which implies that we cannot rule out the
possibility that the correlation result is due to chance factors (Cohen et
al., 2013).

For the correlation between $\log(Y)$ and $\log(\nu^{\rm obs})$ in the
observers’ frame, we find a strong negative correlation ($\tau=-0.59$) for all
blazars in historical sample, consistent with previous studies (Fossati et
al., 1998; Chen & Bai, 2011). While, for the quasi-simultaneous sample, only a
weak negative correlation ($\tau=-0.18$) is found with $p=0.06$, which might
be attributed to the smaller sample size. When considering FSRQs separately,
the correlation disappears in historical sample, while the correlation result
for the quasi-simultaneous sample remains basically unchanged, since this
sample is dominated by FSRQs. This result is consistent with that of Finke
(2013) using data from the second catalogue of AGN detected by Fermi-LAT. For
BL Lacs, we find strong negative correlations in both the historical
($\tau=-0.64$) and quasi-simultaneous sample ($\tau=-0.53$), which are
consistent with Finke (2013). Strong negative correlations are also found for
BL Lacs with $Y\leq 1$ in the both the historical ($\tau=-0.58$) and quasi-
simultaneous ($\tau=-0.54$) samples. However it should be noted that a large
chance probability ($p=0.09$) is obtained for the quasi-simultaneous sample,
which may caused by the small sample size. Also it is interesting that no
significant correlation is found for BL Lacs with $Y>1$ in either the
historical or quasi-simultaneous sample. The above correlation results are
displayed in Figure 1.

Figure 2: The correlations between $\log(L_{\rm rad}^{\rm obs})$ and
$\log(\nu^{\rm obs})$ in the historical sample (upper panel), and in the
quasi-simultaneous sample (lower panel). The meanings of symbols and line
styles are the same as in Figure 1. Similarly, the best linear fitting
equation is only shown when a significant correlation is found. Figure 3: The
upper panel shows the correlations between $\log(L_{\rm rad})$ and
$\log(\nu)$, and the lower panel shows the correlations between $\log(Y)$ and
$\log(\nu)$ in the $\delta_{\rm D}$-corrected sample. The meanings of symbols
and line styles are the same as in Figure 1. Similarly, the best linear
fitting equation is only shown when a significant correlation is found.

The correlations between $\log(L^{\rm obs}_{\rm rad})$ and $\log(\nu^{\rm
obs})$ in the observers’ frame for the historical and quasi-simultaneous
samples are shown in Figure 2. A moderate negative correlation ($\tau=-0.48$)
is found in historical sample, which is consistent with many previous studies
(Fossati et al., 1998; Chen & Bai, 2011; Finke, 2013; Xiong et al., 2015a; Fan
et al., 2016b; Ghisellini et al., 2017), but no correlation ($\tau=-0.05$) is
found in quasi-simultaneous sample. When focusing on FSRQs, weak negative
correlations are found both in historical ($\tau=-0.14$) and quasi-
simultaneous samples ($\tau=-0.2$; also a large chance probability $p=0.08$ is
obtained for quasi-simultaneous sample). For BL Lacs in historical sample, no
correlation is found for the whole sample and for subsample with $Y\leq 1$.
For BL Lacs with $Y>1$, only a weak correlation with a large chance
probability ($\tau=-0.17,~{}p=0.07$) is found. For BL Lacs in quasi-
simultaneous sample, weak correlations are found, but the chance probabilities
are significant high ($p>0.3$), suggesting that these correlations are not
significant. Among the above results of the whole sample in the historical
sample, our main focus is on the correlation between the full-band radiative
luminosity and the peak frequency. In contrast, Yang et al. (2022a) focuses on
the correlations between the peak frequency and radio, optical, X-ray, and
$\gamma$-ray emissions, respectively. Given that most FSRQs have $Y>1$, which
means the $\gamma$-ray luminosity dominates the full-band radiative
luminosity, our results are basically consistent with Yang et al. (2022a).
However, for BL Lacs, the band that dominates the full-band radiative
luminosity is uncertain. Therefore, it is difficult to directly compare our
results with Yang et al. (2022a), which focuses on the correlation of a
specific band. Our correlation study suggests that no correlation is found for
the whole sample of BL Lacs. As can be seen from the upper panel of Figure 2,
no correlation found for the whole sample of BL Lacs may be due to the
discovery of many high-luminosity HSPs and low-luminosity LSPs in the Fourth
Fermi-LAT 12-year Source catalog. It is worth noting that our findings differ
from those obtained using data from the third catalogue of AGN detected by
Fermi-LAT (3LAC), where correlations are found in BL Lacs instead of FSRQs
(Ghisellini et al., 2017). In addition, Fan et al. (2016b) find marginal
correlations in BL Lacs-HSPs ($\tau=0.1,~{}p=3.6\times 10^{-3}$) instead of
FSRQs or BL Lacs-LSPs, and Mao et al. (2016) find strong correlations without
particular separation between FSRQs and BL Lacs.

Based on the results above, there are similarities and differences in
correlation studies of the phenomenological blazar sequence between the
historical and quasi-simultaneous samples. Some of these differences may be
attributed to the fact that the quasi-simultaneous sample contains different
states of SEDs (Xue et al., 2016), while others may be due to the small size
of the quasi-simultaneous sample, resulting in $p>0.01$. Moreover, for the
quasi-simultaneous SED given by (Xue et al., 2016), the observation time for
radio, optical, and X-ray is within one week, while the $\gamma$-ray data has
been integrated over several months. However, variabilities on the scales of
days (e.g., Bonning et al., 2012), hours (e.g., Hayashida et al., 2015), and
even minutes (e.g., Begelman et al., 2008) have been widely observed in
different energy bands of blazars. Therefore, using the quasi-simultaneous
data collected by (Xue et al., 2016) will introduce potential bias. It is
interesting that if the correlation holds both for the historical and quasi-
simultaneous samples, their $s$ values are also quite close. In the subsequent
theoretical interpretation (see Sect. 3), we choose to rely on the results of
the historical sample. A larger quasi-simultaneous sample with high-quality
and same state (high or low) SEDs is needed to check these results in detail.

Table 1: Correlation results of historical, quasi-simultaneous and $\delta_{\rm D}$-corrected samples. Columns from left to right: (1): the sample studied in correlation analysis; (2): the number of blazars in the sample; (3) and (6) are the Spearman test correlation coefficients; (4) and (7) are the chance probabilities; (5) and (8) are the slopes of the best linear fitting equations. in the observers’ frame | N | $\log(Y)~{}{\rm vs.}~{}\log(\nu^{{\rm obs}})$ | $\log(L_{{\rm rad}}^{{\rm obs}})~{}{\rm vs.}~{}\log(\nu^{{\rm{obs}}})$
---|---|---|---
$\tau$ | $p$ | $s$ | $\tau$ | $p$ | $s$
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8)
ALL, historical sample | 1594 | -0.59 | $5.5\times 10^{-152}$ | $-0.31\pm 0.009$ | -0.48 | $2.0\times 10^{-91}$ | $-0.35\pm 0.02$
ALL, quasi-simultaneous sample | 109 | -0.18 | 0.06 | $-0.34\pm 0.07$ | -0.05 | 0.63 | $-0.38\pm 0.12$
FSRQs, historical sample |  |  |  |  |  |  |
ALL | 750 | -0.002 | 0.95 | $-0.04\pm 0.04$ | -0.14 | $7.8\times 10^{-5}$ | $-0.29\pm 0.06$
$Y\leq 1$ | 288 | -0.005 | 0.93 | $-0.02\pm 0.03$ | -0.20 | $5.0\times 10^{-4}$ | $-0.35\pm 0.08$
$Y>1$ | 462 | $5.28\times 10^{-4}$ | 0.98 | $8.05\times 10^{-4}\pm 0.04$ | -0.14 | 0.002 | $-0.18\pm 0.07$
FSRQs, quasi-simultaneous sample |  |  |  |  |  |  |
ALL | 81 | -0.22 | 0.05 | $-0.25\pm 0.13$ | -0.20 | 0.08 | $-0.29\pm 0.20$
BL Lacs, historical sample |  |  |  |  |  |  |
ALL | 844 | -0.64 | $9.6\times 10^{-97}$ | $-0.28\pm 0.01$ | -0.04 | 0.19 | $-0.02\pm 0.03$
$Y\leq 1$ | 726 | –0.58 | $3.3\times 10^{-65}$ | $-0.22\pm 0.01$ | 0.07 | 0.07 | $0.08\pm 0.04$
$Y>1$ | 118 | -0.19 | 0.04 | $-0.06\pm 0.02$ | -0.17 | 0.07 | $-0.21\pm 0.14$
BL Lacs, quasi-simultaneous sample |  |  |  |  |  |  |
ALL | 28 | -0.53 | 0.004 | $-0.32\pm 0.08$ | -0.18 | 0.37 | $-0.22\pm 0.19$
$Y\leq 1$ | 11 | -0.54 | 0.09 | $-0.20\pm 0.10$ | 0.21 | 0.54 | $0.13\pm 0.27$
$Y>1$ | 17 | -0.03 | 0.90 | $-0.14\pm 0.14$ | 0.15 | 0.56 | $0.20\pm 0.31$
in the comoving frame |  | $\log(Y)~{}{\rm vs.}~{}\log(\nu)$ | $\log(L_{{\rm rad}})~{}{\rm vs.}~{}\log(\nu)$
ALL, $\delta_{\rm D}$-corrected sample | 180 | -0.10 | 0.17 | $-0.03\pm 0.05$ | 0.54 | $8.8\times 10^{-15}$ | $0.91\pm 0.09$
FSRQs, $Y>1$ | 126 | -0.23 | 0.12 | -$0.24\pm 0.08$ | 0.40 | $4.4\times 10^{-6}$ | $0.89\pm 0.17$
BL Lacs, $Y>1$ | 51 | 0.16 | 0.15 | $0.06\pm 0.05$ | 0.76 | $8.8\times 10^{-15}$ | $0.95\pm 0.12$

We study the intrinsic blazar sequence using the $\delta_{\rm D}$-corrected
sample, and the correlation results are presented in Figure 3. For the
correlation between $\log(Y)$ and $\log(\nu)$, weak correlations are found,
but the obtained large chance probabilities ($p>0.1$) suggesting that these
correlations might not be significant. In contrast, a strong positive
correlation ($\tau=0.54$) is found for the correlation between $\log(L_{\rm
rad})$ and $\log(\nu)$. Furthermore, when we study FSRQs and BL Lacs
(dominated by ones with $Y>1$), respectively, similar correlation results and
slopes ($s\sim 1$) are derived as well. The obtained positive correlations are
consistent with those obtained by some previous works (Nieppola et al., 2008;
Wu et al., 2009; Fan et al., 2016a; Chen et al., 2021), and the obtained
slopes are consistent with our previous studies (Fan et al., 2017; Yang et
al., 2022b). It is noteworthy that $\delta_{\rm D}$ used by Chen et al. (2021)
are derived from the conventional one-zone leptonic model (Chen, 2018).
However, it is unfortunate that the slope of the best linear fitting equation
is omitted in their study.

##  3 Theoretical analysis

###  3.1 Simple scaling model

In the study of blazars, many phenomenological models have been proposed.
However, most of these models contain numerous free parameters, even the
simplest conventional one-zone model has at least seven free parameters that
are coupled to each other. This may not be inappropriate for studying the
global properties of blazars, such as the blazar sequence. In radio and X-ray
observations, bright knots have been observed along the jet of blazars and
radio galaxies (e.g., Marscher & Jorstad, 2011; Mertens et al., 2016). In
other words, significant dissipation regions can potentially appear at any
location along the jet (Worrall et al., 2007; Ghisellini & Tavecchio, 2009;
Liu et al., 2023). This kind of framework that considers the dominant
dissipation region appearing at different positions has been applied by Potter
(2018) and Wang et al. (2022) to explain the orphan flares of blazars.
Following this phenomenological framework, we assume that one single spherical
dissipation region appearing at different position of the jet dominates the
whole jet emission, which is certainly an over-simplified assumption (for a
multi-zone view, please see Liu et al., 2023). This dominant dissipation
region is composed of a plasma of charged particles in a uniform magnetic
field $B$ with radius $R$ and moving with bulk Lorentz factor
$\Gamma=(1-\beta^{2})^{-1/2}$, where $\beta c$ is the jet speed, along the
jet, at a viewing angle $\theta^{\rm obs}$ with respect to observers’ line of
sight. By assuming $\theta^{\rm obs}\lesssim 1/\Gamma$ for blazars, we have
$\delta_{\rm D}\approx\Gamma$. Let us envision that the dominant dissipation
regions of all blazars can be placed within one jet. In this way, these
dominant dissipation regions are distributed at different positions within
this jet. Based on some reasonable assumptions that usually applied in the
continuous jet model (e.g., Blandford & Königl, 1979), physical parameters of
these dissipation regions are connected with each other. Radio observations
have found that AGNs jets have an approximate conical (Kovalev et al., 2007;
Sokolovsky et al., 2011) or parabolic structure (Nakamura & Asada, 2013;
Algaba et al., 2017; Giovannini et al., 2018; Hodgson et al., 2018). If
assuming that the dissipation region occupies the entire jet cross section, we
have $R\propto r^{\alpha}$, where $r$ is the distance between SMBH and the
dissipation region. More specifically, $\alpha=1$ represents the jet has a
conical structure, and $0<\alpha<1$ means the jet has a parabolic structure.
For the magnetic field $B$, we have $B\propto R^{-1}\propto r^{-\alpha}$ by
assuming the magnetic power is conserved in the jet (Blandford & Königl,
1979), which is a crucial assumption for reproducing the flat radio spectrum
(Potter & Cotter, 2012). In addition, if considering the jet’s continuous
acceleration or deceleration as suggested by phenomenological model,
magnetohydrodynamic simulation, and observation (Blandford & Königl, 1979;
Hawley & Krolik, 2006; Potter & Cotter, 2013a; Mościbrodzka et al., 2014;
Roychowdhury et al., 2024), we assume a simplified acceleration/deceleration
profile, i.e., $\delta_{\rm D}\propto r^{x}$, where $x>0$ for an accelerating
jet, $x<0$ for a decelerating jet, and $x=0$ for a constant-speed jet or
$\delta_{\rm D}$ does not related to other parameters. In this way, power-law
scaling connections of $R$, $B$, $\delta_{\rm D}$ and $r$ can be established.
With this simple scaling model, the potential physical origins of the
statistical correlation of the blazar sequence can be investigated. Of course,
one should bear in mind that the jet each dissipation region located in is
actually different, therefore the normalization of each power-law correlations
is still a free parameter. In this section, parameters with superscript “obs”
are measured in observers’ frame, whereas the parameters without the
superscript are measured in the comoving frame, unless specified otherwise.

Assuming that accelerated relativistic electrons are injected in the
dissipation region with a specific spectral shape $\Omega(\gamma)$ at a
constant rate given by $\dot{Q}(\gamma)=\dot{Q}_{0}\Omega(\gamma),\gamma_{\rm
min}<\gamma<\gamma_{\rm max},$ where $\gamma_{\rm min}$ is the the minimum
electron Lorentz factor, $\gamma_{\rm max}$ is the maximum electron Lorentz
factor, and $\dot{Q}_{0}$ is the normalization in units of $\rm
cm^{-3}~{}s^{-1}$. By giving an electron injection luminosity $L_{\rm inj}$,
$\dot{Q}_{0}$ can be calculated by

$\dot{Q}_{0}=\frac{L_{\rm inj}}{V_{\rm b}mc^{2}\int_{\gamma_{\rm
min}}^{\gamma_{\rm max}}\gamma\Omega(\gamma)d\gamma},$ (1)

where $m$ is the rest mass of an electron, $c$ is the speed of light, and
$V_{\rm b}$ is the volume of the dissipation region. The radiative luminosity
(i.e., the bolometric luminosity) in the comoving frame can be calculated as

$L_{\rm rad}=V_{\rm b}mc^{2}\int_{\gamma_{\rm min}}^{\gamma_{\rm
max}}f(\gamma)\gamma\dot{Q}(\gamma)d\gamma,$ (2)

where $f(\gamma)=\rm min\\{\frac{\it{t}_{\rm dyn}}{\it{t}_{\rm cool}},1\\}$ is
the cooling efficiency of electrons. More specifically, $t_{\rm dyn}=\eta R/c$
with $\eta>1$ is the dynamical timescale222Undoubtedly, this is clearly an
oversimplification, since the modeling of the particle escape requires a
detailed specification of the geometry, the boundary conditions, the magnetic-
field configuration. The particle’s dynamical timescale is likely to be energy
dependent of the energy-dependent diffusion. In phenomenological blazar model,
the form applied here is to represent the average dynamical timescale (e.g.,
Ghisellini, 2013), when considering the escaping electrons at different
position of the spherical blob., and $t_{\rm
cool}=3mc/\left[4\gamma\sigma_{\rm T}U_{\rm B}(1+Y)\right]$ is the radiative
cooling timescale, where $\sigma_{\rm T}$ is the Thomson scattering cross-
section, $U_{\rm B}=B^{2}/(8\pi)$ is the energy density of the magnetic field,
and $Y$ is the Compton dominance. In this work, we aim to study if radiative
cooling is responsible for the blazar sequence, so it is more appropriate to
use the full-band radiative luminosity $L_{\rm rad}$.

Substituting Eq. (1) into Eq. (2), we have

$\begin{split}L_{\rm rad}&=\frac{L_{\rm inj}}{\int_{\gamma_{\rm
min}}^{\gamma_{\rm max}}\gamma\Omega(\gamma)d\gamma}\int_{\gamma_{\rm
min}}^{\gamma_{\rm max}}f(\gamma)\gamma\Omega(\gamma)d\gamma\\\ &=\frac{L_{\rm
inj}}{\int_{\gamma_{\rm min}}^{\gamma_{\rm
max}}\gamma\Omega(\gamma)d\gamma}\bigg{\\{}\int_{\gamma_{\rm c}}^{\gamma_{\rm
max}}\gamma\Omega(\gamma)d\gamma+\frac{f(\gamma)}{\gamma}\int_{\gamma_{\rm
min}}^{\gamma_{\rm c}}\gamma^{2}\Omega(\gamma)d\gamma\bigg{\\}}\\\ &\propto
L_{\rm inj}(\gamma>\gamma_{\rm c})+L_{\rm inj}(\gamma<\gamma_{\rm
c})B^{2}R(1+Y)\end{split}$ (3)

where $\gamma_{\rm c}=\frac{3mc^{2}}{4\sigma_{\rm T}RU_{\rm B}(1+Y)}$
represents the cooling break at which $t_{\rm dyn}=t_{\rm cool}$.

Following the physical interpretation of Ghisellini et al. (1998), if we
believe that the peak of the low-energy component is caused by radiative
cooling, the peak frequency in the comoving frame can be estimated using the
monochromatic approximation (Tavecchio et al., 1998)

$\begin{split}\nu&\approx 3.7\times 10^{6}\gamma_{\rm c}^{2}B\\\ &\propto
B^{-3}R^{-2}(1+Y)^{-2}\end{split}.$ (4)

From Eqs. (3) and (4), it can be seen that expressions of $L_{\rm rad}$, $\nu$
and $Y$ (its expression will be given later) consist of several physical
parameters coupled together, including $R$, $B$, and $\delta_{\rm D}$.
Considering connections of $R$, $B$, $\delta_{\rm D}$ and $r$ in the framework
of the simple scaling model, $L_{\rm rad}$ becomes

$L_{\rm rad}\propto L_{\rm inj}(\gamma>\gamma_{\rm c})+L_{\rm
inj}(\gamma<\gamma_{\rm c})(1+Y)r^{-\alpha},$ (5)

and $\nu$ becomes

$\nu\propto r^{\alpha}(1+Y)^{-2}.$ (6)

For FSRQs, the Compton dominance $Y$ can be evaluated as

$\begin{split}Y&\simeq\frac{\delta_{\rm D}^{2}U_{\rm ext}}{U_{\rm B}}\\\
&\propto\left\\{\begin{array}[]{rl}r^{2\alpha+2x},r<r_{\rm ext}\\\
r^{2\alpha+2x-n},r>r_{\rm ext}\end{array}\right.,\end{split}$ (7)

where $U_{\rm ext}$ is the energy density of the external photon field in the
AGN frame, and $r_{\rm ext}$ represents the characteristic distance of
external photon fields that is found to scale with the disk luminosity $L_{\rm
D}$ as $r_{\rm ext}\propto{L_{\rm D}}^{1/2}$ through reverberation mapping
(Suganuma et al., 2006; Kishimoto et al., 2011; Bentz et al., 2013; Pozo Nuñez
et al., 2014). When the dissipation region is located within $r_{\rm ext}$,
$U_{\rm ext}=L_{\rm D}/(4\pi r_{\rm ext}^{2}c)$ would be a constant value
(Harvey et al., 2020); when the dissipation region is located beyond $r_{\rm
ext}$, $U_{\rm ext}$ decreases with the index $n$ (Sikora et al., 2009;
Hayashida et al., 2012). From Eqs. (6) and (7), the non-linear correlation
between the synchrotron peak frequency and the Compton dominance in
logarithmic space can be obtained

$\begin{split}\nu\propto\left\\{\begin{array}[]{rl}Y^{\frac{\alpha}{2\alpha+2x}}(1+Y)^{-2},r<r_{\rm
ext}\\\ Y^{\frac{\alpha}{2\alpha+2x-n}}(1+Y)^{-2},r>r_{\rm ext}\\\
\end{array}\right.\end{split}$ (8)

in the comoving frame, and

$\begin{split}\nu^{\rm
obs}\propto\left\\{\begin{array}[]{rl}Y^{\frac{1}{2}}(1+Y)^{-2},r<r_{\rm
ext}\\\ Y^{\frac{\alpha+x}{2\alpha+2x-n}}(1+Y)^{-2},r>r_{\rm ext}\\\
\end{array}\right.\end{split}$ (9)

in the observers’ frame. If considering $Y\ll 1$ and $Y\gg 1$333In the
following, whenever symbols ‘$Y\ll 1$’ or ‘$Y\gg 1$’ appear, it implies that
the actual correlation is non-linear in logarithmic space when $Y\sim 1$.,
above correlations become linear in logarithmic space, i.e.,

$\begin{split}Y\propto\left\\{\begin{array}[]{rl}{\nu}^{\frac{2\alpha+2x}{\alpha}},r<r_{\rm
ext}\\\ {\nu}^{\frac{2\alpha+2x-n}{\alpha}},r>r_{\rm
ext}\end{array}\right.,Y\ll 1\\\
Y\propto\left\\{\begin{array}[]{rl}{\nu}^{\frac{2\alpha+2x}{-3\alpha-4x}},r<r_{\rm
ext}\\\ {\nu}^{\frac{2\alpha+2x-n}{2n-3\alpha-4x}},r>r_{\rm
ext}\end{array}\right.,Y\gg 1\end{split}$ (10)

in the comoving frame, and

$\begin{split}Y\propto\left\\{\begin{array}[]{rl}{\nu^{\rm obs}}^{2},r<r_{\rm
ext}\\\ {\nu^{\rm obs}}^{\frac{2\alpha+2x-n}{\alpha+x}},r>r_{\rm
ext}\end{array}\right.,Y\ll 1\\\ Y\propto\left\\{\begin{array}[]{rl}{\nu^{\rm
obs}}^{-2/3},r<r_{\rm ext}\\\ {\nu^{\rm
obs}}^{\frac{2\alpha+2x-n}{2n-3\alpha-3x}},r>r_{\rm
ext}\end{array}\right.,Y\gg 1\end{split}$ (11)

in the observers’ frame. It should be noted that due to the actual correlation
being non-linear, when $Y\sim 1$, under certain combinations of $\alpha$, $n$,
and $x$, the true slope may have significant errors compared to those when
$Y\ll 1$ and $Y\gg 1$. From Eq. (5), it can be seen that if the first term on
the right is dominant, the correlation between $\nu$ and $L_{\rm rad}$ depends
on if there is a correlation between $\nu$ and $L_{\rm inj}$. Assuming that
the first term on the right is not dominant and values of $L_{\rm inj}$ do not
have a large dispersion, we have

$\begin{split}L_{\rm
rad}\propto\left\\{\begin{array}[]{rl}&{\nu}^{-1},\it{Y}\ll 1\\\
&{\nu}^{\frac{2x+\alpha}{-3\alpha-4x}},\it{Y}\gg 1,r<r_{\rm ext}\\\
&{\nu}^{\frac{2x+\alpha-n}{2n-3\alpha-4x}},\it{Y}\gg 1,r>r_{\rm
ext}\end{array}\right.\end{split}$ (12)

in the comoving frame, and

$\begin{split}L_{\rm rad}^{\rm
obs}\propto\left\\{\begin{array}[]{rl}&{\nu^{\rm
obs}}^{\frac{4x-\alpha}{x+\alpha}},\it{Y}\ll 1\\\ &{\nu^{\rm
obs}}^{\frac{6x+\alpha}{-3\alpha-3x}},\it{Y}\gg 1,r<r_{\rm ext}\\\ &{\nu^{\rm
obs}}^{\frac{6x+\alpha-n}{2n-3\alpha-3x}},\it{Y}\gg 1,r>r_{\rm
ext}\end{array}\right.\end{split}$ (13)

in the observers’ frame. Note that the above analysis is based on the
assumption that radiative cooling induces the spectrum break. However cooling
induced change of spectral index is rarely discovered in observations (Baheeja
et al., 2022). Therefore, it is necessary to re-study these two correlations
for $\nu<\nu_{\rm c}$. In this scenario, electrons are cooled in slow cooling
regime, so the break of the electron energy distribution corresponding to the
peak frequency of the low-energy component might be ascribed to multiple
acceleration processes (Bell, 1978b; Tan et al., 2023). In the diffusive shock
acceleration, the electron Lorentz factor can be evaluated as $\gamma\propto
BR$ (Axford et al., 1977; Bell, 1978a, b; Blandford & Ostriker, 1978; Kirk &
Schneider, 1988; Rieger et al., 2007), therefore $\nu\propto B^{3}r^{2}\propto
r^{-\alpha}$. Then, for the correlation between the synchrotron peak frequency
and the Compton dominance, we have

$Y\propto\left\\{\begin{array}[]{rl}{\nu}^{-\frac{2\alpha+2x}{\alpha}},r<r_{\rm
ext}\\\ {\nu}^{-\frac{2\alpha+2x-n}{\alpha}},r>r_{\rm ext}\end{array}\right.$
(14)

in the comoving frame, and

$Y\propto\left\\{\begin{array}[]{rl}{\nu^{\rm
obs}}^{\frac{2\alpha+2x}{x-\alpha}},r<r_{\rm ext}\\\ {\nu^{\rm
obs}}^{\frac{2\alpha+2x-n}{x-\alpha}},r>r_{\rm ext}\end{array}\right.$ (15)

in the observers’ frame. For the correlation between the synchrotron peak
frequency and the radiative luminosity, we have

$\begin{split}L_{\rm rad}\propto\left\\{\begin{array}[]{rl}&\nu,\it{Y}\ll 1\\\
&{\nu}^{-\frac{2x+\alpha}{\alpha}},\it{Y}\gg 1,r<r_{\rm ext}\\\
&{\nu}^{-\frac{2x+\alpha-n}{\alpha}},\it{Y}\gg 1,r>r_{\rm
ext}\end{array}\right.\end{split}$ (16)

in the comoving frame, and

$\begin{split}L_{\rm rad}^{\rm
obs}\propto\left\\{\begin{array}[]{rl}&{\nu^{\rm
obs}}^{\frac{4x-\alpha}{x-\alpha}},\it{Y}\ll 1\\\ &{\nu^{\rm
obs}}^{\frac{6x+\alpha}{x-\alpha}},\it{Y}\gg 1,r<r_{\rm ext}\\\ &{\nu^{\rm
obs}}^{\frac{6x+\alpha-n}{x-\alpha}},\it{Y}\gg 1,r>r_{\rm
ext}\end{array}\right.\end{split}$ (17)

in the observers’ frame.

For BL Lacs, the Compton dominance $Y$ can be evaluated as

$Y\simeq\frac{U_{\rm syn}}{U_{B}},$ (18)

where $U_{\rm syn}=L_{\rm rad}^{\rm syn}/(4\pi R^{2}c)$ is the energy density
of synchrotron photons, $L_{\rm rad}^{\rm syn}$ represents the radiative
synchrotron luminosity

$L_{\rm rad}^{\rm syn}=V_{\rm b}mc^{2}\int_{\gamma_{\rm min}}^{\gamma_{\rm
max}}f_{\rm syn}(\gamma)\gamma\dot{Q}(\gamma)d\gamma,$ (19)

where $f_{\rm syn}(\gamma)=\rm min\\{\frac{\it{t}_{\rm dyn}}{\it{t}_{\rm
syn}},1\\}$ is the cooling efficiency of synchrotron emission. Then $Y$ can be
written as

$Y\propto L_{\rm inj}(\gamma>\gamma_{\rm c})+L_{\rm inj}(\gamma<\gamma_{\rm
c})r^{-\alpha}.$ (20)

Firstly, we study the correlations of the blazar sequence for $\nu=\nu_{\rm
c}$. If $Y$’s the first term on the right is dominant, the correlation between
$\nu$ and $Y$ depends on if there is a correlation between $\nu$ and $L_{\rm
inj}$. If the first term on the right is not dominant, we derive the
correlation between $\nu$ and $Y$,

$\begin{split}Y\propto\left\\{\begin{array}[]{rl}&{\nu}^{-1},Y\ll 1\\\
&{\nu}^{-\frac{1}{3}},\it{Y}\gg 1\end{array}\right.\end{split}$ (21)

in the comoving frame, and

$\begin{split}Y\propto\left\\{\begin{array}[]{rl}&{\nu^{\rm
obs}}^{-\frac{\alpha}{\alpha+x}},Y\ll 1\\\ &{\nu^{\rm
obs}}^{-\frac{\alpha}{3\alpha+x}},\it{Y}\gg 1\end{array}\right.\end{split}$
(22)

in the observers’ frame. Similarly, we derive the correlation between the
synchrotron peak frequency and the radiative luminosity, i.e.,

$\begin{split}L_{\rm
rad}\propto\left\\{\begin{array}[]{rl}&{\nu}^{-1},\it{Y}\ll 1\\\
&{\nu}^{-2/3},\it{Y}\gg 1\end{array}\right.\end{split}$ (23)

in the comoving frame, and

$\begin{split}L_{\rm rad}^{\rm
obs}\propto\left\\{\begin{array}[]{rl}&{\nu^{\rm
obs}}^{\frac{4x-\alpha}{x+\alpha}},\it{Y}\ll 1\\\ &{\nu^{\rm
obs}}^{\frac{4x-2\alpha}{x+3\alpha}},\it{Y}\gg 1\end{array}\right.\end{split}$
(24)

in the observers’ frame.

Similar to previous discussion for FSRQs, we also study these two correlations
for BL Lacs in the case of $\nu<\nu_{\rm c}$, i.e., $\nu\propto
B^{3}r^{2}\propto r^{-\alpha}$. For the correlation between the synchrotron
peak frequency and the Compton dominance (Eq. 20), we have

$Y\propto{\nu}$ (25)

in the comoving frame, and

$Y\propto{\nu^{\rm obs}}^{\frac{\alpha}{\alpha-x}}$ (26)

in the observers’ frame. For the correlation between the synchrotron peak
frequency and the radiative luminosity, we have

$\begin{split}L_{\rm rad}\propto\left\\{\begin{array}[]{rl}&\nu,\it{Y}\ll 1\\\
&\nu^{2},\it{Y}\gg 1\end{array}\right.\end{split}$ (27)

in the comoving frame, and

$\begin{split}L_{\rm rad}^{\rm
obs}\propto\left\\{\begin{array}[]{rl}&{\nu^{\rm
obs}}^{\frac{4x-\alpha}{x-\alpha}},\it{Y}\ll 1\\\ &{\nu^{\rm
obs}}^{\frac{4x-2\alpha}{x-\alpha}},\it{Y}\gg 1\end{array}\right.\end{split}$
(28)

in the observers’ frame.

Table 2: Deduced equations of the correlations in the blazar sequence.

$\nu=\nu_{\rm c}$ | in the observers’ frame | in the comoving frame |  | in the observers’ frame | in the comoving frame
---|---|---|---|---|---
FSRQs, $Y\ll 1$ |  |  |  BL Lacs, $Y\ll 1$ |  |
|  Eq. (11): $Y\propto\left\\{\begin{array}[]{rl}{\nu^{\rm obs}}^{2},r<r_{\rm ext}\\\ {\nu^{\rm obs}}^{\frac{2\alpha+2x-n}{\alpha+x}},r>r_{\rm ext}\end{array}\right.$ |  Eq. (10): $Y\propto\left\\{\begin{array}[]{rl}{\nu}^{\frac{2\alpha+2x}{\alpha}},r<r_{\rm ext}\\\ {\nu}^{\frac{2\alpha+2x-n}{\alpha}},r>r_{\rm ext}\end{array}\right.$ |  |  Eq. (22): $Y\propto{\nu^{\rm obs}}^{-\frac{\alpha}{\alpha+x}}$ |  Eq. (21): $Y\propto{\nu}^{-1}$
|  Eq. (13): $L_{\rm rad}^{\rm obs}\propto{\nu^{\rm obs}}^{\frac{4x-\alpha}{x+\alpha}}$ |  Eq. (12): $L_{\rm rad}\propto{\nu}^{-1}$ |  |  Eq. (24): $L_{\rm rad}^{\rm obs}\propto{\nu^{\rm obs}}^{\frac{4x-\alpha}{x+\alpha}}$ |  Eq. (23): $L_{\rm rad}\propto{\nu}^{-1}$
FSRQs, $Y\gg 1$ |  |  |  BL Lacs, $Y\gg 1$ |  |
|  Eq. (11): $Y\propto\left\\{\begin{array}[]{rl}{\nu^{\rm obs}}^{-2/3},r<r_{\rm ext}\\\ {\nu^{\rm obs}}^{\frac{2\alpha+2x-n}{2n-3\alpha-3x}},r>r_{\rm ext}\end{array}\right.$ |  Eq. (10): $Y\propto\left\\{\begin{array}[]{rl}{\nu}^{\frac{2\alpha+2x}{-3\alpha-4x}},r<r_{\rm ext}\\\ {\nu}^{\frac{2\alpha+2x-n}{2n-3\alpha-4x}},r>r_{\rm ext}\end{array}\right.$ |  |  Eq. (22): $Y\propto{\nu^{\rm obs}}^{-\frac{\alpha}{3\alpha+x}}$ |  Eq. (21): $Y\propto{\nu}^{-\frac{1}{3}}$
|  Eq. (13): $L_{\rm rad}^{\rm obs}\propto\left\\{\begin{array}[]{rl}&{\nu^{\rm obs}}^{\frac{6x+\alpha}{-3\alpha-3x}},r<r_{\rm ext}\\\ &{\nu^{\rm obs}}^{\frac{6x+\alpha-n}{2n-3\alpha-3x}},r>r_{\rm ext}\end{array}\right.$ |  Eq. (12): $L_{\rm rad}\propto\left\\{\begin{array}[]{rl}&{\nu}^{\frac{2x+\alpha}{-3\alpha-4x}},r<r_{\rm ext}\\\ &{\nu}^{\frac{2x+\alpha-n}{2n-3\alpha-4x}},r>r_{\rm ext}\end{array}\right.$ |  |  Eq. (24): $L_{\rm rad}^{\rm obs}\propto{\nu^{\rm obs}}^{\frac{4x-2\alpha}{x+3\alpha}}$ |  Eq. (23): $L_{\rm rad}\propto{\nu}^{-2/3}$
$\nu<\nu_{\rm c}$ |  |  |  |  |
FSRQs |  |  | BL Lacs |  |
|  Eq. (15): $Y\propto\left\\{\begin{array}[]{rl}{\nu^{\rm obs}}^{\frac{2\alpha+2x}{x-\alpha}},r<r_{\rm ext}\\\ {\nu^{\rm obs}}^{\frac{2\alpha+2x-n}{x-\alpha}},r>r_{\rm ext}\end{array}\right.$ |  Eq. (14): $Y\propto\left\\{\begin{array}[]{rl}{\nu}^{-\frac{2\alpha+2x}{\alpha}},r<r_{\rm ext}\\\ {\nu}^{-\frac{2\alpha+2x-n}{\alpha}},r>r_{\rm ext}\end{array}\right.$ |  |  Eq. (26): $Y\propto{\nu^{\rm obs}}^{\frac{\alpha}{\alpha-x}}$ |  Eq. (25): $Y\propto{\nu}$
|  Eq. (17): $L_{\rm rad}^{\rm obs}\propto\left\\{\begin{array}[]{rl}&{\nu^{\rm obs}}^{\frac{4x-\alpha}{x-\alpha},\it{Y}\ll 1}\\\ &{\nu^{\rm obs}}^{\frac{6x+\alpha}{x-\alpha},\it{Y}\gg 1},r<r_{\rm ext}\\\ &{\nu^{\rm obs}}^{\frac{6x+\alpha-n}{x-\alpha},\it{Y}\gg 1},r>r_{\rm ext}\end{array}\right.$ |  Eq. (16): $L_{\rm rad}\propto\left\\{\begin{array}[]{rl}&\nu,\it{Y}\ll 1\\\ &{\nu}^{-\frac{2x+\alpha}{\alpha}},\it{Y}\gg 1,r<r_{\rm ext}\\\ &{\nu}^{-\frac{2x+\alpha-n}{\alpha}},\it{Y}\gg 1,r>r_{\rm ext}\end{array}\right.$ |  |  Eq. (28): $L_{\rm rad}^{\rm obs}\propto\left\\{\begin{array}[]{rl}&{\nu^{\rm obs}}^{\frac{4x-\alpha}{x-\alpha}},\it{Y}\ll 1\\\ &{\nu^{\rm obs}}^{\frac{4x-2\alpha}{x-\alpha}},\it{Y}\gg 1\end{array}\right.$ |  Eq. (27): $L_{\rm rad}\propto\left\\{\begin{array}[]{rl}&\nu,\it{Y}\ll 1\\\ &\nu^{2},\it{Y}\gg 1\end{array}\right.$

As deduced above, we obtain equations of the blazar sequence for FSRQs and for
BL Lacs in the comoving and observers’ frames, respectively. All equations are
summarized in Table 2. By taking the logarithm of the parameters in the
correlations, we can get the corresponding slopes, which can be compared with
the statistical linear regression results.

The simple scaling model proposed in this subsection is quite similar to the
conventional one-zone model, where all the jet emission is represented by one
single dissipation region. Compared to the conventional one-zone model, the
advantage of this simple scaling model is that power-law relationships have
been established among several physical parameters based on reasonable
assumptions. However, normalizations of power-law relationships for each
blazar remains unknown and varies. This will greatly increase the dispersion
of the correlation intercept, even if the correlation slope predicted by the
model is the same. Moreover, the simple scaling model, compared to the
conventional one-zone model, takes many shortcuts. For example, to establish
the power-law correlations among multiple physical parameters, a crude
assumption is introduced, i.e., the electron injected luminosities of all
blazars do not have a large dispersion. However, as suggested by Ghisellini &
Tavecchio (2008), the black hole mass and the accretion rate would be crucial
for interpreting the blazar sequence, since the jet power is always linked
with the accretion rate. Therefore, such an assumption is likely to generate
potential bias. Also, blazars have a complex environment with various external
photon fields, such as the BLR, the DT, etc, each with its own unique
influence distance. Gathering information on each blazar’s external photon
field is tough, so we make an simplification that there is just one main
external photon field. When using the conventional one-zone model to interpret
blazars’ SEDs, the evolution of the relativistic electron energy distribution
is a factor that cannot be ignored. However, in our attempt to understand the
blazar sequence using the simple scaling model, we do not specify the
particular shape for the electron energy distribution, nor do we study its
evolution. This simplification is feasible because the study of the blazar
sequence only requires information about the synchrotron peak frequency and
integrated luminosity, both of which can be obtained through analytical
calculations.

###  3.2 Theoretical Implications

In this subsection, we apply the simple scaling model described in Sect. 3.1
to interpret the statistical correlation results of the phenomenological and
intrinsic blazar sequence derived in Sect. 2.

First of all, let us investigate if statistical results can be interpreted
under the conditions that $\nu=\nu_{\rm c}$, i.e., cooling is important. As
mentioned before, values of $\delta_{\rm D}$ given by radio observation may
have a great discrepancy with the actual values of $\delta_{\rm D}$. If one
consider that radio $\delta_{\rm D}$ is significantly influential and
unreliable, or if we believe that the distribution of actual $\delta_{\rm D}$
is irregular, i.e., $x=0$, then it can be seen from Table 2 that the
synchrotron peak frequency is both anti-correlated with the Compton dominance
and the radiative luminosity in the observers’ and comoving frames. This is
consistent with the original physical interpretation proposed by Ghisellini et
al. (1998). If here we trust the radio $\delta_{\rm D}$ and believe in the
correlation results of the intrinsic blazar sequence, then the physical
explanation of cooling might faces difficulties. From the statistical results
in Table 1, it can be seen that no correlation is found between the Compton
dominance and the synchrotron peak frequency for all (sub-)samples of FSRQs,
whether in the observers’ frame or the comoving frame. It indicates that
indexes of the deduced Eqs. (10, 11) are zero, i.e., $r>r_{\rm ext}$ and
$2\alpha+2x-n=0$. Consequently, the indexes of Eq. (12) between $L_{\rm rad}$
and $\nu$ in the comoving frame become -1. In addition, Eq. (23) suggests that
$L_{\rm rad}$ is negatively correlated with $\nu$ for BL Lacs in the comoving
frame (even when considering the actual nonlinear correlation). However,
moderate and strong positive correlations between $L_{\rm rad}$ and $\nu$ are
found for FSRQs and BL Lacs, respectively, implying a different physical
origin of these correlations. In the following, we will discuss if deduced
equations under the condition that $\nu<\nu_{\rm c}$ have the potential to
account for the phenomenological and intrinsic blazar sequence simultaneously.

Figure 4: Parameter spaces for FSRQs with different value of $\alpha$ when
$r>r_{\rm ext}$. Solid curves represent that the indexes of Eqs. (14, 15) are
equal to zero, as statistical analysis find that $Y$ and $\nu$ (or $\nu^{\rm
obs}$) are not correlated. Dashed, dotted, and dash-dotted curves are obtained
by setting the indexes of Eqs. (16, 17) equal to slopes obtained from
statistical analysis, respectively. Since errors of slopes are also taken into
account, the intersection area between two curves of each type represents the
effective parameter space. After considering the intersection of all line
types, we find the corresponding combinations of $n$, $x$, and $\alpha$, which
are marked by red circles.

For FSRQs, if we consider $r<r_{\rm ext}$, Eqs. (14, 15) suggest that $x$
should be equal to $-\alpha$, since no correlation is found between the
Compton dominance and the synchrotron peak frequency both in the observers’
frame and the comoving frame. Consequently, the index of Eq. (16) becomes
unity (including the case $Y\ll 1$), which is consistent with the statistical
results. However, the index of Eq. (17) will be 5/2, which is inconsistent
with the statistical results that found $s=-0.35$ when $Y\leq 1$ and $s=-0.18$
when $Y>1$. Therefore, we need to further discuss the case of $r>r_{\rm ext}$.
In this case, since statistical studies find that the Compton dominance and
the synchrotron peak frequency are not correlated, Eqs. (14, 15) suggest that
$2\alpha+2x-n=0$. Moreover, by making the indexes of Eqs. (16, 17) equal to
the slopes obtained in the statistical studies (including errors), we can get
combinations of $x$, $n$, and $\alpha$ that satisfies all conditions, which
are shown in Figure 4. It can be seen that the size of effective intersection
area is proportional to the value of $\alpha$. It implies that the effective
parameter space can only be obtained when $x>0$, i.e., the dissipation region
has to be located in the accelerating jet. We emphasize that any combinations
of $x$, $n$, and $\alpha$ that satisfies $2\alpha+2x-n=0$ are possible. In
Figure 4, we present five sets of parameter combinations, which are
$\alpha=1,n\simeq 2.6,x\simeq 0.3$; $\alpha=0.8,n\simeq 2.1,x\simeq 0.24$;
$\alpha=0.6,n\simeq 1.56,x\simeq 0.18$; $\alpha=0.4,n\simeq 1.04,x\simeq
0.12$; and $\alpha=0.2,n\simeq 0.52,x\simeq 0.06$. Please note that our
discussion here is based on the assumption of a single external photon field,
which is a simplification. In the actual AGNs’ environment, it is generally
believed that there are two external photon fields, i.e., the BLR and the DT,
both of which could have significant implications for jet emission. For the
BLR, $n=3$ has been suggested as a model assumption by Sikora et al. (2009).
For the DT, $n=4$ has been found by Hayashida et al. (2012) for a specific
observation of an FSRQ 3C 279. While it may not be the case that $n=3$ or
$n=4$ works for every AGN, our results clearly suggest that the majority of
AGNs cannot have an $n$ greater than 3. Otherwise, $\alpha$ would be larger
than 1, implying that the jet profile becomes hyperbolic. This is evidently in
contradiction with the conical or parabolic structure found by radio
observations (see Blandford et al., 2019, for a review). On the other hand,
the covering factors, which represent the fractions of the disk luminosity
reprocessed into the BLR and DT radiation, evidently vary among different
AGNs. For instance, values such as 0.1, 0.2, and 0.5 have been suggested
(Maiolino et al., 2007; Ghisellini & Tavecchio, 2009; Hao et al., 2010, 2011).
This degree of variability introduces complexity into the model interpretation
attempted here, and may even preclude the possibility of reaching definitive
conclusions. Therefore, the assumption of a single external photon field in
this work is a simplifying approximation. In general, our model suggests that
the phenomenological and intrinsic blazar sequence of FSRQs can be explained
only when the condition $2\alpha+2x-n=0$ is satisfied, implying that the
dissipation region is in an accelerating jet located beyond the external
photon field in the slow-cooling regime.

For BL Lacs, the physical interpretation of phenomenological and intrinsic
blazar sequence appears to be complicated. For BL Lacs with $Y>1$, the derived
correlation results and slopes (if correlations exist) are similar to those of
FSRQs. On the other hand, we find that BL Lacs with $Y>1$ in different samples
are dominated by LSPs (61/118 = 52 % for historical sample, 14/17 = 82 % for
quasi-simultaneous sample, and 23/51 = 45 % for $\delta_{\rm D}$-corrected
sample), which may suggest that their high-energy components are mainly from
the EC emission rather than the SSC emission (Böttcher, 2007; Böttcher et al.,
2013; Wang & Xue, 2021; Deng & Jiang, 2023). Blazars are classified as FSRQs
and BL Lacs based on whether broad emission lines are detected. However, there
are many blazars with comparable jet and broad emission lines intensities that
are classified into either subclass depending on the jet activity during
observation. For instance, if the jet emission is in a low state during
observation, the blazar will be observed with broad emission lines and
classified as FSRQs. On the other hand, if the jet emission is in a high
state, the emission lines will be masked, then the blazar will be classified
as BL Lacs. Such blazars are known as “changing-look” blazars (Matt et al.,
2003; Bianchi et al., 2005). Additionally, some blazars with broad emission
lines are classified as BL Lacs because their broad emission lines are
outshone by the jet emission. These blazars are suggested as “masquerading” BL
Lacs (Giommi et al., 2013). In our sample, eighteen blazars (including 4FGL
J0238.6+1637, 4FGL J0334.2-4008, 4FGL J0407.5+0741, 4FGL J0428.6-3756, 4FGL
J0438.9-4521, 4FGL J0516.7-6207, 4FGL J0538.8-4405, 4FGL J0629.3-1959, 4FGL
J0710.9+4733, 4FGL J0831.8+0429, J1001.1+2911, 4FGL J1058.4+0133, 4FGL
J1147.0-3812, 4FGL J1751.5+0938, 4FGL J1800.6+7828, 4FGL J1954.6-1122, 4FGL
J2134.2-0154, and 4FGL J2152.5+1737) are suggested as changing-look blazars
(Xiao et al., 2022), which can be seen as the direct evidence of presence of
external photon fields. Therefore, we suggest that the phenomenological and
intrinsic blazar sequence of BL Lacs with $Y>1$ can also be explained by the
dissipation region in accelerating jet (located beyond the external photon
field) emitting in the slow-cooling regime, as in FSRQs. For BL Lacs with
$Y\leq 1$, we do not find enough BL Lacs with measured $\delta_{\rm D}$ in the
literature, so we only study the phenomenological blazar sequence. For the
phenomenological correlations, no correlation is found between $\log(L_{{\rm
rad}}^{{\rm obs}})~{}{\rm vs.}~{}\log(\nu^{{\rm{obs}}})$, indicating the
indexes of Eq. (24, 28) are zeros, i.e., $4x-\alpha=0$. Since $0<\alpha\leq
1$, it indicates that the dissipation region is in an accelerating jet with
$0<x\leq 0.25$ (no matter $\nu=\nu_{\rm c}$ or $\nu<\nu_{\rm c}$). If
defaulting $x=0.25$, the obtained strong negative correlation between
$\log(Y)$ and $\log(\nu^{{\rm{obs}}})$ can be explained in the case of
$\nu=\nu_{\rm c}$. Such a cooling scenario is quite possible for BL Lacs with
$Y\leq 1$ since HSPs are dominant in the historical (472/726=65%) and quasi-
simultaneous (5/11=45%) samples. On the other hand, it is necessary to note
that the model predicted slope ($-0.8$) for the correlation between $\log(Y)$
and $\log(\nu^{{\rm{obs}}})$ is lower than that derived by the statistical
analysis ($\sim-0.2$). We suppose that this might be due to the BL Lacs in two
samples having mixed cooling regimes, since the model predicted correlations
under two cooling scenarios are opposite as shown in Table 2. Here we
emphasize that the above discussion assumes $x=0.25$ merely for the
convenience of showing that the statistical correlations could be reproduced
based on some specific conditions. In fact, any combination of $x$ and
$\alpha$ that satisfies $4x-\alpha=0$ is plausible. Furthermore, the attemp of
introducing mixed cooling scenarios arises from the single cooling scenario’s
inadequacy to account for the slope derived from statistical analysis.

###  3.3 Reproducing the blazar sequence

In previous, our theoretical analysis of the blazar sequence focuses on the
power-law relation between various parameters (summarized in Table 2), as the
obtained indexes can be compared with slopes given by correlation studies.
However, specific values of normalization in power-law relations are ignored,
which will inevitably increase the dispersion of the correlation. In addition,
we do not provide specific values for each physical parameter, so one may
worry that special and unreasonable physical parameters will be introduced. On
the other hand, as shown in Sect. 3.1, the two correlations of the blazar
sequence under some conditions are nonlinear, which will have a certain impact
on the slope and dispersion of the correlation. In this subsection, we would
like to apply the simple scaling model to reproduce the blazar sequence shown
in Figures 1–3, and provide the distribution of important physical parameters.
In the following, we default that the jet has a conical structure (i.e.,
$\alpha=1$) as the benchmark case, as found in observations (Kovalev et al.,
2007; Sokolovsky et al., 2011) and assumed in theoretical models (Kaiser,
2006; Potter & Cotter, 2012).

Applying the simple scaling model, 1800 blazars, including 300 FSRQs with
$Y\leq 1$, 500 FSRQs with $Y>1$, 700 BL Lacs with $Y\leq 1$ and 300 BL Lacs
with $Y>1$, are generated that are similar in composition to the historical
sample. The physical parameters of these FSRQs and BL Lacs are assigned by
generating random numbers from a certain range of values, following either a
normal or uniform distribution. It is expected that model generated slopes are
consistent with those in Table 1, so we assume that values of $n$ and $x$ of
FSRQ and BL Lacs with $Y>1$ respectively conform to normal distributions with
a mean of 2.6 and a standard deviation of 0.1, and a mean of 0.3 and a
standard deviation of 0.05. For BL Lacs with $Y<1$, values of $x$ conform to a
normal distribution with a mean of 2.5 and a standard deviation of 0.1. Values
of the other physical parameters conform to uniform distributions within
physically plausible ranges. For all blazars, uniform distributions of some
parameters have the same range: $0.1^{\circ}\leqslant\theta_{\rm
open}\leqslant 3^{\circ}$ ($\theta_{\rm open}$ represents the jet half opening
angle; Finke, 2019); $r_{\rm ext}\leqslant\it{r}\leqslant\rm 10^{2}~{}pc$;
$1\leqslant\eta\leqslant 10$ (Inoue & Takahara, 1996; Mücke & Protheroe, 2001;
Gao et al., 2017); $1\leqslant\delta_{\rm D,0}\leqslant 3$ ($\delta_{\rm D,0}$
represents the initial $\delta_{\rm D}$); $2\leqslant s_{\rm e}\leqslant
4$444Please note that the range of $s_{\rm e}$ here is for the convenience of
simulation, and it at odds with the simplest diffusion shock acceleration
(e.g., Bednarz & Ostrowski, 1998). The value of $s_{\rm e}$ affects the ratio
of $L_{\rm inj}(\gamma>\gamma_{\rm c})/L_{\rm rad}$. As mentioned before Eq.
(12), if $L_{\rm inj}(\gamma>\gamma_{\rm c})$ dominates, all correlations will
depend only on if there is a correlation between $\nu$ and $L_{\rm inj}$,
which will inevitably change the correlation or increase the dispersion of the
correlation largely. When applying the simple scaling model to reproduce the
blazar sequence, we only retain blazars with $L_{\rm inj}(\gamma>\gamma_{\rm
c})/L_{\rm rad}\leq 0.5$ to ensure that the final correlation would not be
affected too much. Note that, the value of $L_{\rm inj}(\gamma>\gamma_{\rm
c})/L_{\rm rad}$ is jointly determined by $s_{\rm e}$ and $\gamma_{\rm c}$.
The removal of blazars with $L_{\rm inj}(\gamma>\gamma_{\rm c})/L_{\rm
rad}>0.5$ does not contradict the condition $\nu<\nu_{\rm c}$ applied in the
simulation. ($s_{\rm e}$ represents the spectral index of electrons energy
distribution). In addition, there are some differences in the range of uniform
distributions of $L_{\rm inj}$, $B$ and $L_{\rm D}$. For FSRQs, we set
$10^{45}~{}\rm erg~{}s^{-1}\leqslant\it{L}_{\rm inj}\leqslant\rm
10^{46}~{}erg~{}s^{-1}$; $0.1~{}\rm G\leqslant\it{B}_{\rm 1pc}\leqslant\rm
1~{}G$ ($\it{B}_{\rm 1pc}$ represents the magnetic field at 1 pc);
$10^{44.5}~{}\rm erg~{}s^{-1}\leqslant\it{L}_{\rm D}\leqslant\rm
10^{46.5}~{}erg~{}s^{-1}$. For BL Lacs with $Y>1$, we set $10^{43}~{}\rm
erg~{}s^{-1}\leqslant\it{L}_{\rm inj}\leqslant\rm 10^{44}~{}erg~{}s^{-1}$;
$0.01~{}\rm G\leqslant\it{B}_{\rm 1pc}\leqslant\rm 0.1~{}G$; $10^{42.5}~{}\rm
erg~{}s^{-1}\leqslant\it{L}_{\rm D}\leqslant\rm 10^{44.5}~{}erg~{}s^{-1}$. For
BL Lacs with $Y<1$, we set $10^{43}~{}\rm erg~{}s^{-1}\leqslant\it{L}_{\rm
inj}\leqslant\rm 10^{44}~{}erg~{}s^{-1}$; $0.05~{}\rm G\leqslant\it{B}_{\rm
1pc}\leqslant\rm 0.5~{}G$.

(a)  (b)

Figure 5: The phenomenological (upper two panels) and intrinsic (lower two
panels) blazar sequence generated by the simple scaling model. Similar to
Figures 1–3, if a significant correlation is found ($\tau>0.1$, $p<0.01$) for
a specific subsample statistically, the best linear fitting equation is also
shown. Red and teal dashed lines represent the best-fitting equation for the
populations of FSRQs and BL Lacs with $Y>1$, respectively. Red and teal dotted
lines represent the best-fitting equation for the populations of FSRQs and BL
Lacs with $Y\leq 1$, respectively. The meanings of symbols are explained in
the inset legends.

Based on the chosen ranges and distributions of parameters described above,
the phenomenological and intrinsic blazar sequence generated by the simple
scaling model are shown in Figure 5, and the corresponding correlation results
are given in Table 3. It can be seen that although the model generated
correlation results and slopes are basically consistent with those of
historical sample and $\delta_{\rm D}$-corrected sample, there are
morphological differences between the data distribution obtained from the
model and the observed data when viewed with the naked eye. For the
correlation results, one major difference is that the slope of the correlation
between $\log(L_{{\rm rad}}^{{\rm obs}})$ and $\log(\nu^{{\rm{obs}}})$ for
FSRQs with $Y\leq 1$ generated by the model is $-0.10\pm 0.04$, while the
slope obtained from the historical sample is $-0.35\pm 0.08$. This may be due
to the fact that the EC emission of these FSRQs with $Y\leq 1$ no longer
dominates the high-energy peak, which has changed the slope to some extent. On
the other hand, our model, under the current settings, gives the predicted
correlation results of FSRQs and BL Lacs with $Y\leq 1$, respectively, which
can be checked in the future with larger and more complete samples.
Distributions of physical parameters are shown in Figure 6. It can be seen
that parameters distributions are generally reasonable and consistent with
those suggested in the conventional leptonic model (Chen, 2018; Cerruti, 2020;
Tan et al., 2020). However, it should be noted that the parameters
distributions obtained here are only used to show that the phenomenological
and intrinsic blazar sequence can be reproduced by the simple scaling model
within a reasonable parameter range, and may not necessarily represent the
actual physical properties of blazars. In fact, the blazar sequence currently
reproduced using this simple scaling model may not necessarily have universal
applicability. Although power-law correlations are built between physical
parameters with reasonable assumptions, the normalizations are all random
numbers, which greatly increases the dispersion and weakens the correlation.
For example, as shown in Table 3, except for $\log(L_{{\rm rad}})~{}{\rm
vs.}~{}\log(\nu)$, all others are weakly correlated or uncorrelated. If the
range of parameter distributions continues to expand, e.g., $L_{\rm inj}$, all
the relevant correlations will disappear or the slope will change largely.
Actually, Ghisellini & Tavecchio (2008) suggest that the accretion rate, which
is directly related to $L_{\rm inj}$, is important for the explanation of the
blazar sequence. Therefore, the model generated results are highly sensitive
to the chosen ranges and distributions of parameters. Further study is needed
to identify the parameters that truly dominate. In addition, the actual blazar
radiation mechanism must be more complex than the model we currently used. For
example, the AGN environment contains multiple external photon fields, rather
than the current single external photon field assumption; the jet acceleration
profile is more complex and not even continuous; observations indicate that
some jets have a parabolic structure (Nakamura & Asada, 2013; Algaba et al.,
2017; Giovannini et al., 2018). As a result, correlations of the historical
sample are mostly weakly correlated or uncorrelated as well, and the
correlation obtained from different samples also varies, as discussed in Sect.
2.2.

Table 3: Correlation results of the phenomenological and intrinsic blazar sequence generated by the simple scaling model. | $r$ | $p$ | $s$
---|---|---|---
FSRQ, $Y>1$ |  |  |
$\log(L_{{\rm rad}}^{{\rm obs}})~{}{\rm vs.}~{}\log(\nu^{{\rm{obs}}})$ | -0.21 | $2.5\times 10^{-5}$ | $-0.16\pm 0.04$
$\log(Y)~{}{\rm vs.}~{}\log(\nu^{{\rm obs}})$ | -0.08 | 0.08 | $-0.05\pm 0.03$
$\log(L_{{\rm rad}})~{}{\rm vs.}~{}\log(\nu)$ | 0.64 | $2.2\times 10^{-53}$ | $0.74\pm 0.04$
$\log(Y)~{}{\rm vs.}~{}\log(\nu)$ | -0.09 | 0.05 | $-0.05\pm 0.02$
FSRQ, $Y\leq 1$ |  |  |
$\log(L_{{\rm rad}}^{{\rm obs}})~{}{\rm vs.}~{}\log(\nu^{{\rm{obs}}})$ | -0.15 | $6.7\times 10^{-3}$ | $-0.10\pm 0.04$
$\log(Y)~{}{\rm vs.}~{}\log(\nu^{{\rm obs}})$ | -0.04 | 0.51 | $-0.02\pm 0.03$
$\log(L_{{\rm rad}})~{}{\rm vs.}~{}\log(\nu)$ | 0.68 | $6.6\times 10^{-43}$ | $0.77\pm 0.05$
$\log(Y)~{}{\rm vs.}~{}\log(\nu)$ | -0.07 | 0.22 | $-0.03\pm 0.02$
BL Lac, $Y>1$ |  |  |
$\log(L_{{\rm rad}}^{{\rm obs}})~{}{\rm vs.}~{}\log(\nu^{{\rm{obs}}})$ | -0.14 | 0.02 | $-0.14\pm 0.05$
$\log(Y)~{}{\rm vs.}~{}\log(\nu^{{\rm obs}})$ | -0.08 | 0.16 | $-0.05\pm 0.04$
$\log(L_{{\rm rad}})~{}{\rm vs.}~{}\log(\nu)$ | 0.64 | $3.5\times 10^{-35}$ | $0.76\pm 0.05$
$\log(Y)~{}{\rm vs.}~{}\log(\nu)$ | -0.12 | 0.05 | $-0.07\pm 0.03$
BL Lac, $Y\leq 1$ |  |  |
$\log(L_{{\rm rad}}^{{\rm obs}})~{}{\rm vs.}~{}\log(\nu^{{\rm{obs}}})$ | 0.07 | 0.05 | $0.08\pm 0.02$
$\log(Y)~{}{\rm vs.}~{}\log(\nu^{{\rm obs}})$ | -0.34 | $4.8\times 10^{-24}$ | $-0.28\pm 0.02$
$\log(L_{{\rm rad}})~{}{\rm vs.}~{}\log(\nu)$ | -0.22 | $4.6\times 10^{-10}$ | $-0.24\pm 0.03$
$\log(Y)~{}{\rm vs.}~{}\log(\nu)$ | -0.22 | $4.6\times 10^{-10}$ | $-0.20\pm 0.02$

(a)  (b)  (c)  (d)

(e)  (f)  (g)

Figure 6: Distributions of the electron injection luminosity $\log L_{\rm
inj}$, the distance between SMBH and the dissipation region $\log r$, the
radius of the dissipation region $\log R$, the electron Lorentz factor
corresponding to the synchrotron peak frequency $\log\gamma_{\rm b}$, the
magnetic field in the dissipation region $\log B$, the Doppler factor of the
dissipation region $\log\delta_{\rm D}$, and the disk luminosity $\log L_{\rm
D}$. Mean values for FSRQs are, $\bar{L}_{\rm inj}=3.2\times
10^{45}\rm~{}erg~{}s^{-1}$, $\bar{r}=3.5~{}\rm pc$, $\bar{R}=1.9\times
10^{17}\rm~{}cm$, $\bar{\gamma}_{\rm b}=4.8\times 10^{3}$,
$\bar{B}=0.3\rm~{}G$, $\bar{\delta}_{\rm D}=16.8$, $\bar{L}_{\rm D}=3.6\times
10^{45}\rm~{}erg~{}s^{-1}$, respectively. Mean values for BL Lacs are,
$\bar{L}_{\rm inj}=3.3\times 10^{43}\rm~{}erg~{}s^{-1}$,
$\bar{r}=0.14\rm~{}pc$, $\bar{R}=7.8\times 10^{15}\rm~{}cm$,
$\bar{\gamma}_{\rm b}=2.3\times 10^{3}$, $\bar{B}=1.7\rm~{}G$,
$\bar{\delta}_{\rm D}=21.2$, $\bar{L}_{\rm D}=8.7\times
10^{43}\rm~{}erg~{}s^{-1}$, respectively. The meanings of histograms with
different colors are explained in the inset legends.

##  4 Conclusion

In this work, we revisit correlations in the phenomenological and intrinsic
blazar sequence across three samples, which are the historical sample, the
quasi-simultaneous sample and the $\delta_{\rm D}$-corrected sample, selected
from literature. In attempting to interpret these correlations, we propose the
simple scaling model, in which physical parameters of the dissipation region
are connected to the location of the dissipation region. Our conclusions are
as follows:

Statistical correlation results: When considering all types of blazars as a
whole, the phenomenological blazar sequence holds in the historical sample. In
which, a strong negative correlation is found between the Compton dominance
and the synchrotron peak frequency, as well as a moderate negative
correlations is found between the radiative luminosity and the synchrotron
peak frequency. However, the phenomenological blazar sequence does not exist
in the quasi-simultaneous sample. It might be caused by the fact that the
influence of variability that could cause massive changes in the SED is not
considered in Xue et al. (2016). Their main purpose is to make use of the
maximum availability of simultaneous data coverage, so SEDs in both low and
high states are included. For the intrinsic blazar sequence, correlations
between the Compton dominance and the synchrotron peak frequency, and between
the radiative luminosity and the synchrotron peak frequency display no
correlation and strong positive correlation, respectively. For FSRQs, we find
no correlations between the Compton dominance and the synchrotron peak
frequency in either observers’ or comoving frame. In addition, we find weak
negative correlations (excluding the quasi-simultaneous sample) in the
observers’ frame and a moderate positive correlation in the comoving frame
between the radiative luminosity and the synchrotron peak frequency. For BL
Lacs with $Y>1$, we find no correlations between the Compton dominance and the
synchrotron peak frequency in either observers’ or comoving frame, and no
correlations in the observers’ frame and strong positive correlations in the
comoving frame between the radiative luminosity and the synchrotron peak
frequency. The derived correlation results and slopes (when correlations
exist) are similar to those of FSRQs. For BL Lacs with $Y\leq 1$, we find
strong negative correlations (excluding the quasi-simultaneous sample) between
the Compton dominance and the synchrotron peak frequency, and no correlations
between the radiative luminosity and the synchrotron peak frequency in the
observers’ frame.

Theoretical implications: In this work, in attempting to reproduce the
correlations of the blazar sequence, we propose a simple scaling model. In
this model, a dominant dissipation region, which takes into account radiative
cooling, is considered to occur along the jet. Consequently, under reasonable
assumptions, the physical parameters of the dissipation region are linked to
the location of the dissipation region itself. In the modeling, we find that
the correlations in the blazar sequence cannot be reproduced satisfactorily
unless considering some specific conditions that have been fine-tuned. This
implies that radiative cooling alone may not be the primary cause of the
blazar sequence. It further indicates that additional physical processes not
considered in the simple scaling model is needed to interpret the blazar
sequence within a more physically plausible framework. Based on a sensible
range of physical parameters, we employ our simple scaling model to generate a
population of blazars. The objective is to ascertain whether the observed
correlations in the blazar sequence can be accurately replicated with the
generated blazar population. Whilst this method is promising to test different
hypotheses of the underlying physical mechanism of the blazar sequence, we
find that the model generated results are so sensitive to the chosen ranges
and distributions of parameters that it may not be able to accurately
reproduce the broad properties of the observed blazar population. This
demonstrates that this simple scaling model lacks the complexity required to
interpret the blazar sequence. Further consideration of black hole mass,
accretion rate and a more realistic emission calculation may be required to
explain the underlying physics of the blazar sequence.

## Acknowledgements

We thank the anonymous referees for insightful comments and constructive
suggestions. We thank Prof. J. H. Yang for his kindly help and sharing data.
This work is supported by the National Natural Science Foundation of China
(NSFC) under the Grants No. 12203043. Z.R.W. acknowledges the support by the
NSFC under the the Grants No. 12203024. H.B.X. acknowledges the support by the
NSFC under the the Grants No. 12203034 and by Shanghai Science and Technology
Fund under the Grants No. 22YF1431500. J.H.F. acknowledges the support by the
NSFC under the the Grants No. U2031201.

## Data Availability

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

## References

  * Abdo et al. (2010a) Abdo A. A., et al., 2010a, ApJ, 710, 1271
  * Abdo et al. (2010b) Abdo A. A., et al., 2010b, ApJ, 716, 30
  * Abdollahi et al. (2022) Abdollahi S., et al., 2022, ApJS, 260, 53
  * Acciari et al. (2011) Acciari V. A., et al., 2011, ApJ, 729, 2
  * Ackermann et al. (2011) Ackermann M., et al., 2011, ApJ, 743, 171
  * Ackermann et al. (2015) Ackermann M., et al., 2015, ApJ, 810, 14
  * Algaba et al. (2017) Algaba J. C., Nakamura M., Asada K., Lee S. S., 2017, ApJ, 834, 65
  * Axford et al. (1977) Axford W. I., Leer E., Skadron G., 1977, in International Cosmic Ray Conference. p. 132
  * Baheeja et al. (2022) Baheeja C., Sahayanathan S., Rieger F. M., Jagan S. K., Ravikumar C. D., 2022, MNRAS, 514, 3074
  * Bednarz & Ostrowski (1998) Bednarz J., Ostrowski M., 1998, Phys. Rev. Lett., 80, 3911
  * Begelman et al. (2008) Begelman M. C., Fabian A. C., Rees M. J., 2008, MNRAS, 384, L19
  * Bell (1978a) Bell A. R., 1978a, MNRAS, 182, 147
  * Bell (1978b) Bell A. R., 1978b, MNRAS, 182, 443
  * Bennett et al. (2014) Bennett C. L., Larson D., Weiland J. L., Hinshaw G., 2014, ApJ, 794, 135
  * Bentz et al. (2013) Bentz M. C., et al., 2013, ApJ, 767, 149
  * Bianchi et al. (2005) Bianchi S., Guainazzi M., Matt G., Chiaberge M., Iwasawa K., Fiore F., Maiolino R., 2005, A&A, 442, 185
  * Björnsson (2010) Björnsson C. I., 2010, ApJ, 723, 417
  * Blandford & Königl (1979) Blandford R. D., Königl A., 1979, ApJ, 232, 34
  * Blandford & Ostriker (1978) Blandford R. D., Ostriker J. P., 1978, ApJ, 221, L29
  * Blandford & Payne (1982) Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433
  * Blandford et al. (2019) Blandford R., Meier D., Readhead A., 2019, ARA&A, 57, 467
  * Bonning et al. (2012) Bonning E., et al., 2012, ApJ, 756, 13
  * Böttcher (2007) Böttcher M., 2007, Ap&SS, 309, 95
  * Böttcher & Dermer (2002) Böttcher M., Dermer C. D., 2002, ApJ, 564, 86
  * Böttcher et al. (2013) Böttcher M., Reimer A., Sweeney K., Prakash A., 2013, ApJ, 768, 54
  * Cerruti (2020) Cerruti M., 2020, Galaxies, 8, 72
  * Cerruti et al. (2017) Cerruti M., Benbow W., Chen X., Dumm J. P., Fortson L. F., Shahinyan K., 2017, A&A, 606, A68
  * Chen (2018) Chen L., 2018, ApJS, 235, 39
  * Chen & Bai (2011) Chen L., Bai J. M., 2011, ApJ, 735, 108
  * Chen et al. (2021) Chen Y., et al., 2021, ApJ, 906, 108
  * Cohen et al. (2013) Cohen J., Cohen P., West S. G., Aiken L. S., 2013, Applied multiple regression/correlation analysis for the behavioral sciences. Routledge
  * Deng & Jiang (2023) Deng J., Jiang Y., 2023, MNRAS, 521, 6210
  * Dermer & Schlickeiser (1993) Dermer C. D., Schlickeiser R., 1993, ApJ, 416, 458
  * Dermer et al. (1992) Dermer C. D., Schlickeiser R., Mastichiadis A., 1992, A&A, 256, L27
  * Fan et al. (2016a) Fan X.-L., Bai J.-M., Mao J.-R., 2016a, Research in Astronomy and Astrophysics, 16, 173
  * Fan et al. (2016b) Fan J. H., et al., 2016b, ApJS, 226, 20
  * Fan et al. (2017) Fan J. H., et al., 2017, ApJ, 835, L38
  * Finke (2013) Finke J. D., 2013, ApJ, 763, 134
  * Finke (2019) Finke J. D., 2019, ApJ, 870, 28
  * Fossati et al. (1998) Fossati G., Maraschi L., Celotti A., Comastri A., Ghisellini G., 1998, MNRAS, 299, 433
  * Gao et al. (2017) Gao S., Pohl M., Winter W., 2017, ApJ, 843, 109
  * Ghisellini (2013) Ghisellini G., 2013, Radiative Processes in High Energy Astrophysics. Vol. 873, doi:10.1007/978-3-319-00612-3,
  * Ghisellini (2016) Ghisellini G., 2016, Galaxies, 4, 36
  * Ghisellini & Maraschi (1989) Ghisellini G., Maraschi L., 1989, ApJ, 340, 181
  * Ghisellini & Tavecchio (2008) Ghisellini G., Tavecchio F., 2008, MNRAS, 387, 1669
  * Ghisellini & Tavecchio (2009) Ghisellini G., Tavecchio F., 2009, MNRAS, 397, 985
  * Ghisellini et al. (1998) Ghisellini G., Celotti A., Fossati G., Maraschi L., Comastri A., 1998, MNRAS, 301, 451
  * Ghisellini et al. (2009) Ghisellini G., Tavecchio F., Ghirlanda G., 2009, MNRAS, 399, 2041
  * Ghisellini et al. (2010) Ghisellini G., Tavecchio F., Foschini L., Ghirlanda G., Maraschi L., Celotti A., 2010, MNRAS, 402, 497
  * Ghisellini et al. (2011) Ghisellini G., Tavecchio F., Foschini L., Ghirlanda G., 2011, MNRAS, 414, 2674
  * Ghisellini et al. (2017) Ghisellini G., Righi C., Costamante L., Tavecchio F., 2017, MNRAS, 469, 255
  * Giommi & Padovani (2015) Giommi P., Padovani P., 2015, MNRAS, 450, 2404
  * Giommi et al. (2012a) Giommi P., Padovani P., Polenta G., Turriziani S., D’Elia V., Piranomonte S., 2012a, MNRAS, 420, 2899
  * Giommi et al. (2012b) Giommi P., et al., 2012b, A&A, 541, A160
  * Giommi et al. (2013) Giommi P., Padovani P., Polenta G., 2013, MNRAS, 431, 1914
  * Giovannini et al. (2018) Giovannini G., et al., 2018, Nature Astronomy, 2, 472
  * Hao et al. (2010) Hao H., et al., 2010, ApJ, 724, L59
  * Hao et al. (2011) Hao H., Elvis M., Civano F., Lawrence A., 2011, ApJ, 733, 108
  * Harvey et al. (2020) Harvey A. L. W., Georganopoulos M., Meyer E. T., 2020, Nature Communications, 11, 5475
  * Hawley & Krolik (2006) Hawley J. F., Krolik J. H., 2006, ApJ, 641, 103
  * Hayashida et al. (2012) Hayashida M., et al., 2012, ApJ, 754, 114
  * Hayashida et al. (2015) Hayashida M., et al., 2015, ApJ, 807, 79
  * Hodgson et al. (2018) Hodgson J. A., et al., 2018, MNRAS, 475, 368
  * Hovatta et al. (2009) Hovatta T., Valtaoja E., Tornikoski M., Lähteenmäki A., 2009, A&A, 494, 527
  * Inoue & Takahara (1996) Inoue S., Takahara F., 1996, ApJ, 463, 555
  * Kaiser (2006) Kaiser C. R., 2006, MNRAS, 367, 1083
  * Keenan et al. (2021) Keenan M., Meyer E. T., Georganopoulos M., Reddy K., French O. J., 2021, MNRAS, 505, 4726
  * Kerby & Falcone (2023) Kerby S., Falcone A. D., 2023, ApJ, 951, 133
  * Kirk & Schneider (1988) Kirk J. G., Schneider P., 1988, A&A, 201, 177
  * Kishimoto et al. (2011) Kishimoto M., Hönig S. F., Antonucci R., Millour F., Tristram K. R. W., Weigelt G., 2011, A&A, 536, A78
  * Konigl (1981) Konigl A., 1981, ApJ, 243, 700
  * Kovalev et al. (2007) Kovalev Y. Y., Lister M. L., Homan D. C., Kellermann K. I., 2007, ApJ, 668, L27
  * Landt & Bignall (2008) Landt H., Bignall H. E., 2008, MNRAS, 391, 967
  * Landt et al. (2004) Landt H., Padovani P., Perlman E. S., Giommi P., 2004, MNRAS, 351, 83
  * Liodakis et al. (2017) Liodakis I., et al., 2017, MNRAS, 466, 4625
  * Liodakis et al. (2018) Liodakis I., Hovatta T., Huppenkothen D., Kiehlmann S., Max-Moerbeck W., Readhead A. C. S., 2018, ApJ, 866, 137
  * Liu et al. (2023) Liu R.-Y., Xue R., Wang Z.-R., Tan H.-B., Böttcher M., 2023, MNRAS, 526, 5054
  * Maiolino et al. (2007) Maiolino R., Shemmer O., Imanishi M., Netzer H., Oliva E., Lutz D., Sturm E., 2007, A&A, 468, 979
  * Mao et al. (2016) Mao P., Urry C. M., Massaro F., Paggi A., Cauteruccio J., Künzel S. R., 2016, ApJS, 224, 26
  * Marcha et al. (1996) Marcha M. J. M., Browne I. W. A., Impey C. D., Smith P. S., 1996, MNRAS, 281, 425
  * Marscher & Jorstad (2011) Marscher A. P., Jorstad S. G., 2011, ApJ, 729, 26
  * Massaro et al. (2004) Massaro E., Perri M., Giommi P., Nesci R., 2004, A&A, 413, 489
  * Matt et al. (2003) Matt G., Guainazzi M., Maiolino R., 2003, MNRAS, 342, 422
  * Mertens et al. (2016) Mertens F., Lobanov A. P., Walker R. C., Hardee P. E., 2016, A&A, 595, A54
  * Meyer et al. (2011) Meyer E. T., Fossati G., Georganopoulos M., Lister M. L., 2011, ApJ, 740, 98
  * Mościbrodzka et al. (2014) Mościbrodzka M., Falcke H., Shiokawa H., Gammie C. F., 2014, A&A, 570, A7
  * Mücke & Protheroe (2001) Mücke A., Protheroe R. J., 2001, Astroparticle Physics, 15, 121
  * Nakamura & Asada (2013) Nakamura M., Asada K., 2013, ApJ, 775, 118
  * Nieppola et al. (2006) Nieppola E., Tornikoski M., Valtaoja E., 2006, A&A, 445, 441
  * Nieppola et al. (2008) Nieppola E., Valtaoja E., Tornikoski M., Hovatta T., Kotiranta M., 2008, A&A, 488, 867
  * Ouyang et al. (2023) Ouyang Z., et al., 2023, ApJ, 949, 52
  * Padovani et al. (2012) Padovani P., Giommi P., Rau A., 2012, MNRAS, 422, L48
  * Potter (2018) Potter W. J., 2018, MNRAS, 473, 4107
  * Potter & Cotter (2012) Potter W. J., Cotter G., 2012, MNRAS, 423, 756
  * Potter & Cotter (2013a) Potter W. J., Cotter G., 2013a, MNRAS, 429, 1189
  * Potter & Cotter (2013b) Potter W. J., Cotter G., 2013b, MNRAS, 436, 304
  * Pozo Nuñez et al. (2014) Pozo Nuñez F., et al., 2014, A&A, 561, L8
  * Prandini & Ghisellini (2022) Prandini E., Ghisellini G., 2022, Galaxies, 10, 35
  * Raiteri & Capetti (2016) Raiteri C. M., Capetti A., 2016, A&A, 587, A8
  * Rieger et al. (2007) Rieger F. M., Bosch-Ramon V., Duffy P., 2007, Ap&SS, 309, 119
  * Roychowdhury et al. (2024) Roychowdhury A., Meyer E. T., Georganopoulos M., Kollmann K., 2024, MNRAS, 527, 10262
  * Savard et al. (2022) Savard K., Ruan J. J., Haggard D., 2022, MNRAS, 509, 4620
  * Sikora et al. (2002) Sikora M., Błażejowski M., Moderski R., Madejski G. M., 2002, ApJ, 577, 78
  * Sikora et al. (2009) Sikora M., Stawarz Ł., Moderski R., Nalewajko K., Madejski G. M., 2009, ApJ, 704, 38
  * Sokolovsky et al. (2011) Sokolovsky K. V., Kovalev Y. Y., Pushkarev A. B., Lobanov A. P., 2011, A&A, 532, A38
  * Suganuma et al. (2006) Suganuma M., et al., 2006, ApJ, 639, 46
  * Tan et al. (2020) Tan C., Xue R., Du L.-M., Xi S.-Q., Wang Z.-R., Xie Z.-H., 2020, ApJS, 248, 27
  * Tan et al. (2023) Tan H.-B., Liu R.-Y., Böttcher M., 2023, arXiv e-prints, p. arXiv:2311.13873
  * Tavecchio et al. (1998) Tavecchio F., Maraschi L., Ghisellini G., 1998, ApJ, 509, 608
  * Urry & Padovani (1995) Urry C. M., Padovani P., 1995, PASP, 107, 803
  * Wang & Xue (2021) Wang Z.-R., Xue R., 2021, Research in Astronomy and Astrophysics, 21, 305
  * Wang et al. (2017) Wang Z., Xue R., Du L., Xie Z., Xiong D., Yi T., Xu Y., Liu W., 2017, Ap&SS, 362, 191
  * Wang et al. (2022) Wang Z.-R., Liu R.-Y., Petropoulou M., Oikonomou F., Xue R., Wang X.-Y., 2022, Phys. Rev. D, 105, 023005
  * Worrall et al. (2007) Worrall D. M., Birkinshaw M., Laing R. A., Cotton W. D., Bridle A. H., 2007, MNRAS, 380, 2
  * Wu et al. (2009) Wu Z.-Z., Gu M.-F., Jiang D.-R., 2009, Research in Astronomy and Astrophysics, 9, 168
  * Xiao et al. (2022) Xiao H., Fan J., Ouyang Z., Hu L., Chen G., Fu L., Zhang S., 2022, ApJ, 936, 146
  * Xiong & Zhang (2014) Xiong D. R., Zhang X., 2014, MNRAS, 441, 3375
  * Xiong et al. (2015a) Xiong D., Zhang X., Bai J., Zhang H., 2015a, MNRAS, 450, 3568
  * Xiong et al. (2015b) Xiong D., Zhang X., Bai J., Zhang H., 2015b, MNRAS, 451, 2750
  * Xiong et al. (2017) Xiong D., Bai J., Zhang H., Fan J., Gu M., Yi T., Zhang X., 2017, ApJS, 229, 21
  * Xiong et al. (2020) Xiong D., et al., 2020, ApJS, 247, 49
  * Xue et al. (2016) Xue R., et al., 2016, MNRAS, 463, 3038
  * Yang et al. (2022a) Yang J. H., et al., 2022a, ApJS, 262, 18
  * Yang et al. (2022b) Yang W. X., et al., 2022b, ApJ, 925, 120
  * Yang et al. (2023) Yang J., et al., 2023, Science China Physics, Mechanics, and Astronomy, 66, 249511