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

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.


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 − log 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 -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 ★ E-mail: ruixue@zjnu.edu.cn† E-mail: zerui_wang62@163.com‡ E-mail: hubing.xiao@shnu.edu.cn§ E-mail: fjh@gzhu.edu.cnlines, 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 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 BLR in units of the Eddington luminosity Edd of the central supermassive black hole (SMBH).FSRQs are objects whose BLR 5 × 10 −4 Edd , and BL Lacs are the others.In addition to the FSRQs and BL Lacs classification, Abdo et al. (2010a) divide blazars into high-synchrotronpeaked (HSP) blazars if their synchrotron peak p s > 10 15 Hz, intermediate-synchrotron-peaked (ISP) blazars if their synchrotron peak 10 14 Hz < p s < 10 15 Hz and low-synchrotron-peaked (LSP) blazars if their synchrotron peak p s < 10 14 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 wellknown "blazar sequence".The first physical explanation is proposed by Ghisellini et al. (1998).They suggested that if the intrinsic jet power is connected to 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,a;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. 2016aFan et al. , 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 0 = 69.6 km s −1 Mpc −1 , Ω 0 = 0.29, and Ω Λ = 0.71 are adopted (Bennett et al. 2014).

REVISITING THE BLAZAR SEQUENCE WITH THE 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. 2017Xiong 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 miti-gate 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.

Samples
This work focuses on the blazar sequence and uses blazars as samples, since -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: (i) The "historical sample", compiled by Yang et al. (2022aYang et al. ( , 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 them1 , enabling the derivation of Compton dominance , as well as synchrotron peak frequency obs and radiative luminosity obs rad in observers' frame.(ii) Xue et al. (2016) collect quasi-simultaneous multiwavelength 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.
(iii) By cross-matching the Fermi-4FGL blazars collected in Yang et al. (2022aYang et al. ( , 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 D and measured redshifts.Hereafter, this sample is referred to as the " D -corrected sample".Due to the beaming effect, the synchrotron peak frequency and radiative luminosity rad in the comoving frame can be obtained through ≃ obs −1 D and rad ≃ obs rad −4 D , respectively.Please note that this D -corrected sample is highly biased, as there are only three sources with < 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 nonstandard 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 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 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 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 D measured by radio observation to study the intrinsic blazar sequence will inevitably introduce a large uncertainty and potential observational biases.

Correlation Results
We analyze the phenomenological blazar sequence in both the historical and quasi-simultaneous samples, and the intrinsic blazar sequence in D -corrected sample.We do not study the phenomenological blazar sequence in D -corrected sample further, since the D -corrected sample is a subsample of historical sample.The theoretical analysis proposed in Sect. 3 indicates that FSRQs and BL Lacs with ≤ 1 and > 1 may have different correlations and slopes, so we further divide three samples collected in this work using the dividing line = 1.Table 1 shows the statistical correlation results, including slopes of the best linear fitting equations (denoted by 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 > 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( ) and log( obs ) in the observers' frame, we find a strong negative correlation ( = −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 ( = −0.18) is found with = 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 ( = −0.64)and quasi-simultaneous sample ( = −0.53),which are consistent with Finke (2013).Strong negative correlations are also found for BL Lacs with ≤ 1 in the both the historical ( = −0.58)and quasi-simultaneous ( = −0.54)samples.However it should be noted that a large chance probability ( = 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 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 > 1 and the empty symbol represents the source with ≤ 1.If a significant correlation is found ( > 0.1, < 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 > 1, and the dotted red line represents the best-fitting equation for the population of FSRQs with ≤ 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 > 1, and the dotted teal line represents the best-fitting equation for the population of BL Lacs with ≤ 1.
BL Lacs with > 1 in either the historical or quasi-simultaneous sample.The above correlation results are displayed in Figure 1.
The correlations between log( obs rad ) and log( obs ) in the observers' frame for the historical and quasi-simultaneous samples are shown in Figure 2. A moderate negative correlation ( = −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 ( = −0.05) is found in quasi-simultaneous sample.When focusing on FSRQs, weak negative correlations are found both in historical ( = −0.14)and quasi-simultaneous samples ( = −0.2;also a The correlations between log( obs rad ) and log( 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.
large chance probability = 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 ≤ 1.For BL Lacs with > 1, only a weak correlation with a large chance probability ( = −0.17,= 0.07) is found.For BL Lacs in quasi-simultaneous sample, weak correlations are found, but the chance probabilities are significant high ( > 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 -ray emissions, respectively.Given that most FSRQs have > 1, which means the -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 ( = 0.1, = 3.6 × 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 > 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 -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 ob-served 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 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.
We study the intrinsic blazar sequence using the D -corrected sample, and the correlation results are presented in Figure 3.For the correlation between log( ) and log( ), weak correlations are found, but the obtained large chance probabilities ( > 0.1) suggesting that these correlations might not be significant.In contrast, a strong positive correlation ( = 0.54) is found for the correlation between log( rad ) and log( ).Furthermore, when we study FSRQs and BL Lacs (dominated by ones with > 1), respectively, similar correlation results and slopes ( ∼ 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 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.

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 with radius and moving with bulk Lorentz factor Γ = (1 − 2 ) −1/2 , where is the jet speed, along the jet, at a viewing angle obs with respect to observers' line of sight.By assuming obs 1/Γ for blazars, we have D ≈ Γ.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 ∝ , where is the distance between SMBH and the dissipation region.More specifically, = 1 represents the jet has a conical structure, and 0 < < 1 means the jet has a parabolic structure.For the magnetic field , we have ∝ −1 ∝ − 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., D ∝ , where > 0 for an accelerating jet, < 0 for a decelerating jet, and = 0 for a constant-speed jet or D does not related to other parameters.In this way, power-law scaling connections of , , D and 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 Ω( ) at a constant rate given by ( ) = 0 Ω( ), min < < max , where min is the the minimum electron Lorentz factor, max is the maximum electron Lorentz factor, and 0 is the normalization in units of cm −3 s −1 .By giving an electron injection luminosity inj , 0 can be calculated by where is the rest mass of an electron, is the speed of light, and b is the volume of the dissipation region.The radiative luminosity (i.e., the bolometric luminosity) in the comoving frame can be calculated as where ( ) = min{ t dyn t cool , 1} is the cooling efficiency of electrons.More specifically, dyn = / with > 1 is the dynamical timescale 2 , and cool = 3 /[4 T B (1 + )] is the radiative cooling timescale, where T is the Thomson scattering cross-section, B = 2 /(8 ) is the energy density of the magnetic field, and is the Compton dominance.In this work, we aim to study if radia-tive cooling is responsible for the blazar sequence, so it is more appropriate to use the full-band radiative luminosity rad .
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) From Eqs. ( 3) and ( 4), it can be seen that expressions of rad , and (its expression will be given later) consist of several physical parameters coupled together, including , , and D .Considering connections of , , D and in the framework of the simple scaling model, rad becomes and becomes For FSRQs, the Compton dominance can be evaluated as where ext is the energy density of the external photon field in the AGN frame, and ext represents the characteristic distance of external photon fields that is found to scale with the disk luminosity D as ext ∝ 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 ext , ext = D /(4 2 ext ) would be a constant value (Harvey et al. 2020); when the dissipation region is located beyond ext , ext decreases with the index (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 in the comoving frame, and in the observers' frame.If considering ≪ 1 and ≫ 1 3 , above 3 In the following, whenever symbols ' ≪ 1' or ' ≫ 1' appear, it implies that the actual correlation is non-linear in logarithmic space when ∼ 1.
correlations become linear in logarithmic space, i.e., in the comoving frame, and in the observers' frame.It should be noted that due to the actual correlation being non-linear, when ∼ 1, under certain combinations of , , and , the true slope may have significant errors compared to those when ≪ 1 and ≫ 1.From Eq. ( 5), it can be seen that if the first term on the right is dominant, the correlation between and rad depends on if there is a correlation between and inj .
Assuming that the first term on the right is not dominant and values of inj do not have a large dispersion, we have in the comoving frame, and 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 < 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 ∝ (Axford et al. 1977;Bell 1978a,b;Blandford & Ostriker 1978;Kirk & Schneider 1988;Rieger et al. 2007), therefore ∝ 3 2 ∝ − .Then, for the correlation between the synchrotron peak frequency and the Compton dominance, we have in the comoving frame, and in the observers' frame.For the correlation between the synchrotron peak frequency and the radiative luminosity, we have in the comoving frame, and in the observers' frame.
For BL Lacs, the Compton dominance can be evaluated as where syn = syn rad /(4 2 ) is the energy density of synchrotron photons, syn rad represents the radiative synchrotron luminosity where syn ( ) = min{ t dyn t syn , 1} is the cooling efficiency of synchrotron emission.Then can be written as Firstly, we study the correlations of the blazar sequence for = c .If 's the first term on the right is dominant, the correlation between and depends on if there is a correlation between and inj .If the first term on the right is not dominant, we derive the correlation between and , in the comoving frame, and in the observers' frame.Similarly, we derive the correlation between the synchrotron peak frequency and the radiative luminosity, i.e., in the comoving frame, and in the observers' frame.Similar to previous discussion for FSRQs, we also study these two correlations for BL Lacs in the case of < c , i.e., ∝ 3 2 ∝ − .For the correlation between the synchrotron peak frequency and the Compton dominance (Eq.20), we have in the comoving frame, and in the observers' frame.For the correlation between the synchrotron peak frequency and the radiative luminosity, we have in the comoving frame, and in the observers' frame.
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.

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 = c , i.e., cooling is important.As mentioned before, values of D given by radio observation may have a great discrepancy with the actual values of D .If one consider that radio D is significantly influential and unreliable, or if we believe that the distribution of actual D is irregular, i.e., = 0, then it can be seen from Table 2 that the synchrotron peak frequency is both anticorrelated 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 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., > ext and Eq. ( 10): Eq. ( 22): ∝ obs − + Eq. ( 21): ∝ −1 Eq. ( 13): obs rad ∝ obs 4 − + Eq. ( 12): rad ∝ −1 Eq. ( 24): obs rad ∝ obs 4 − + Eq. ( 23): rad ∝ −1 FSRQs, ≫ 1 BL Lacs, ≫ 1 Eq. ( 11): Eq. ( 22): ∝ obs − 3 + Eq. ( 21): Eq. ( 13): obs rad ∝ Eq. ( 24): obs rad ∝ obs 4 −2 +3 Eq. ( 23) Eq. ( 15): Eq. ( 14): Eq. ( 26): ∝ obs − Eq. ( 25): ∝ Eq. ( 17) Eq. ( 16) Eq. ( 28): obs rad ∝  (14,15) are equal to zero, as statistical analysis find that and (or 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 , , and , which are marked by red circles.
2 + 2 − = 0. Consequently, the indexes of Eq. ( 12) between rad and in the comoving frame become -1.In addition, Eq. ( 23) suggests that rad is negatively correlated with for BL Lacs in the comoving frame (even when considering the actual nonlinear correlation).However, moderate and strong positive correlations between rad and 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 < c have the potential to account for the phenomenological and intrinsic blazar sequence simultaneously.
For FSRQs, if we consider < ext , Eqs. (14, 15) suggest that should be equal to − , 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 ≪ 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 = −0.35when ≤ 1 and = −0.18when > 1.Therefore, we need to further discuss the case of > 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 + 2 − = 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 , , and 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 .It implies that the effective parameter space can only be obtained when > 0, i.e., the dissipation region has to be located in the accelerating jet.We emphasize that any combinations of , , and that satisfies 2 + 2 − = 0 are possible.In Figure 4, we present five sets of parameter combinations, which are = 1, ≃ 2.6, ≃ 0.3; = 0.8, ≃ 2.1, ≃ 0.24; = 0.6, ≃ 1.56, ≃ 0.18; = 0.4, ≃ 1.04, ≃ 0.12; and = 0.2, ≃ 0.52, ≃ 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, = 3 has been suggested as a model assumption by Sikora et al. (2009).For the DT, = 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 = 3 or = 4 works for every AGN, our results clearly suggest that the majority of AGNs cannot have an greater than 3. Otherwise, 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. 2010Hao et al. , 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 + 2 − = 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 > 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 > 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 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 > 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 ≤ 1, we do not find enough BL Lacs with measured D in the literature, so we only study the phenomenological blazar sequence.For the phenomenological correlations, no correlation is found between log( obs rad ) vs. log( obs ), indicating the indexes of Eq. (24, 28) are zeros, i.e., 4 − = 0. Since 0 < ≤ 1, it indicates that the dissipation region is in an accelerating jet with 0 < ≤ 0.25 (no matter = c or < c ).If defaulting = 0.25, the obtained strong negative correlation between log( ) and log( obs ) can be explained in the case of = c .Such a cooling scenario is quite possible for BL Lacs with ≤ 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( ) and log( obs ) is lower than that derived by the statistical analysis (∼ −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 = 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 and that satisfies 4 − = 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.

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., = 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 ≤ 1, 500 FSRQs with > 1, 700 BL Lacs with ≤ 1 and 300 BL Lacs with > 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 and of FSRQ and BL Lacs with > 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 < 1, values of 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: 4 Please note that the range of 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 e affects the ratio of inj ( > c )/ rad .As mentioned before Eq. ( 12), if inj ( > c ) dominates, all correlations will depend only on if there is a correlation between and 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 inj ( > c )/ rad ≤ 0.5 to ensure that the final correlation would not be affected too much.Note that, the value of inj ( > c )/ rad is jointly determined by e and c .The removal of blazars with inj ( > c )/ rad > 0.5 does not contradict the condition < c applied in the simulation.
are some differences in the range of uniform distributions of inj , and D .For FSRQs, we set 10 45 erg s −1 L inj 10 46 erg s −1 ; 0.1 G B 1pc 1 G (B 1pc represents the magnetic field at 1 pc); 10 44.5 erg s −1 L D 10 46.5 erg s −1 .For BL Lacs with > 1, we set 10 43 erg s −1 L inj 10 44 erg s −1 ; 0.01 G B 1pc 0.1 G; 10 42.5 erg s −1 L D 10 44.5 erg s −1 .For BL Lacs with < 1, we set 10 43 erg s −1 L inj 10 44 erg s −1 ; 0.05 G B 1pc 0.5 G. 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 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( obs rad ) and log( obs ) for FSRQs with ≤ 1 generated by the model is −0.10 ± 0.04, while the slope obtained from the historical sample is −0.35 ± 0.08.This may be due to the fact that the EC emission of these FSRQs with ≤ 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 ≤ 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( rad ) vs. log( ), all others are weakly correlated or uncorrelated.If the range of parameter distributions continues to expand, e.g., 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 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.

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 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 > 1, we find no correlations between the Compton dominance and the synchrotron peak frequency in either observers' or comoving frame, and no cor- relations 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 ≤ 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.

Figure 1 .
Figure1.The correlations between log( ) and log( 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 > 1 and the empty symbol represents the source with ≤ 1.If a significant correlation is found ( > 0.1, < 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 > 1, and the dotted red line represents the best-fitting equation for the population of FSRQs with ≤ 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 > 1, and the dotted teal line represents the best-fitting equation for the population of BL Lacs with ≤ 1.

Figure 3 .
Figure3.The upper panel shows the correlations between log( rad ) and log( ), and the lower panel shows the correlations between log( ) and log( ) in the D -corrected sample.The meanings of symbols and line styles are the same as in Figure1.Similarly, the best linear fitting equation is only shown when a significant correlation is found.

Figure 4 .
Figure 4. Parameter spaces for FSRQs with different value of when > ext .Solid curves represent that the indexes of Eqs.(14, 15)  are equal to zero, as statistical analysis find that and (or 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 , , and , which are marked by red circles.
0.1 • open 3 • ( open represents the jet half opening angle; Finke 2019); ext r 10 2 pc; the spectral index of electrons energy distribution).In addition, there

Figure 5 .
Figure5.The phenomenological (upper two panels) and intrinsic (lower two panels) blazar sequence generated by the simple scaling model.Similar to Figures1-3, if a significant correlation is found ( > 0.1, < 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 > 1, respectively.Red and teal dotted lines represent the best-fitting equation for the populations of FSRQs and BL Lacs with ≤ 1, respectively.The meanings of symbols are explained in the inset legends.

Table 1 .
Correlation results of historical, quasi-simultaneous and 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.

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

Table 3 .
Correlation results of the phenomenological and intrinsic blazar sequence generated by the simple scaling model.