© The Authors 2026

1. Introduction

Blazars, which comprise flat-spectrum radio quasars (FSRQs) and BL Lac objects, are among the most extreme active galactic nuclei (AGN), and have relativistic jets closely aligned to our line of sight (Blandford et al. 2019). Their broadband emission is nonthermal and strongly Doppler-boosted by the bulk relativistic motion of the jet plasma. The spectral energy distribution (SED) of blazars displays two broad components: a low-energy hump, from radio to optical–UV (and occasionally X-rays), attributed to synchrotron radiation from relativistic electrons; and a high-energy hump, from X-rays to GeV–TeV γ-rays, generally interpreted as inverse Compton scattering (ICS) of either the synchrotron photons themselves, i.e., synchrotron self-Compton (SSC), or photons from external fields, i.e., external Compton (EC) (e.g., Dermer et al. 1992; Sikora et al. 1994; Tavecchio et al. 1998; Boula & Mastichiadis 2022). Variability timescales, combined with relativistic beaming and radiative cooling constraints, provide important clues about the size and location of the emitting regions. Additionally, optical and multiwavelength polarization measurements reveal the magnetic field structure and degree of ordering in the emitting regions, providing complementary constraints on particle acceleration and emission site geometry (Liodakis et al. 2022).

The nonthermal particle populations responsible for these SED components must be continuously replenished by in situ acceleration processes within the jet (for a review, see Matthews et al. 2020). Several such mechanisms have been proposed and investigated in the literature: diffusive shock acceleration at relativistic shocks (Blandford & Ostriker 1978; Kirk & Duffy 1999; Shirin T et al. 2025), magnetic reconnection in highly magnetized plasma (for a review, see Sironi et al. 2025), shear (or shear-driven) acceleration operating in velocity-gradient layers of the flow (for a review, see Rieger 2019), stochastic (second-order Fermi) acceleration in turbulent regions (e.g., Katarzyński et al. 2006), or combinations of them (e.g., Tavecchio 2021). The efficiency and relative importance of these processes depend on local magnetization, the flow Lorentz factor and shear, turbulence amplitude and spectrum, and dissipation physics, and may vary across sources and along the jet (Tavecchio 2021).

A systematic trend in the SED properties – called the blazar sequence by Fossati et al. (1998) and the updated sequence by Ghisellini et al. (2017) – shows that as bolometric luminosity L_bol increases, the synchrotron peak frequency ${\nu }_{\mathrm{pk}}^{\mathrm{syn}}$ decreases, while the Compton dominance (CD; the IC-to-synchrotron luminosity ratio) and γ-ray spectral slope both increase (Finke 2013). Within this framework, the blazar sequence was interpreted as primarily driven by the accretion power onto the central black hole, a single dominant parameter proposed by Böttcher & Dermer (2002) and independently by Cavaliere & D’Elia (2002). The original sequence, based on a sample of 132 objects with only 33 γ-ray detections, has been refined with Fermi-LAT observations, thus expanding the γ-ray–detected blazar sample to over 5000 sources (The Fermi-LAT collaboration 2020) and suggesting that the sequence may ultimately reflect fundamental AGN parameters.

The optically thin nature of blazar GeV emission places the γ-ray production site far from the immediate vicinity of the black hole (BH), possibly at distances up to ∼10⁶R_S (∼10 pc), for M_BH ∼ 10⁸ M_⊙, where R_S is the Schwarzschild radius (Marscher et al. 2010). On these scales, the dusty molecular tori invoked in AGN unification (Antonucci & Miller 1985) may interact significantly with the jet. A compelling picture identifies these structures as magnetohydrodynamic (MHD) accretion-disk winds (Konigl & Kartje 1994; Contopoulos & Lovelace 1994), launched from a few R_S out to the BH sphere of influence (∼10⁶R_S). Observations of warm absorbers, which are blueshifted absorption features successfully modeled as photoionized MHD winds extending to comparable distances (Behar 2009; Fukumura et al. 2010), support this interpretation. Similar winds inferred in stellar-mass black hole systems (Fukumura et al. 2017) suggest a scale-invariant phenomenon.

The physical properties of these MHD winds depend primarily on a small set of parameters and, most notably, on the mass accretion rate. Previous work showed that varying only this parameter can, in principle, reproduce the main trends of the blazar sequence (Boula et al. 2019a, 2019b). To explain the observed multiwavelength spectrum, we adopted a two-zone framework. In the compact acceleration region considered here, particle energization can occur, for instance, through diffusive shock acceleration, after which electrons escape into a larger-scale cooling zone where radiative losses dominate. Building on this idea, we developed a two-zone emission model in which electrons are accelerated near the central engine and subsequently escape into a cooling zone, where interactions with the ambient photon field associated with the MHD wind can dominate the high-energy emission.

This Letter is structured as follows. In Sect. 2 we describe the emission model and its scalings with mass accretion rate. In Sect. 3 we present the two-zone emission framework. In Sect. 4 we summarize our conclusions.

2. Scaling of particle emission and acceleration with mass accretion rate

In our previous work, Boula et al. (2019a), we modeled blazar emission with relativistic electrons radiating via synchrotron and Compton processes, the latter involving seed photons originating from the MHD wind. By varying only the mass accretion rate, the model reproduced the observed correlations between the synchrotron peak, the gamma-ray index, and the Compton dominance. A natural next step is to include explicit particle acceleration physics, which can arise from shocks, turbulence, or shear flows (e.g., Kirk et al. 1998; Rieger & Mannheim 2002).

To capture these effects, we introduced a characteristic acceleration timescale t_acc, parameterized to depend on $\dot{m}$. For first-order Fermi acceleration at relativistic shocks, the timescale can be written $t_{\mathrm{acc},\mathrm{FI}}\gtrsim {\left(\frac{c}{u_s}\right)}^2\frac{\lambda }{c}\simeq \frac{r_{g}c}{u_s^2}$, with u_s the shock speed in the comoving frame, λ the particle mean free path, and r_g = γmc²/(eB) the gyroradius. The spectral shape of blazars depends on the ratio of magnetic to photon energy densities. Following Boula et al. (2019a), we assume partial equipartition between magnetic and accretion power, with accretion luminosity $P_{\mathrm{acc}}=\dot{m}\mathcal{M}{L}_{\mathrm{Edd}}$, where L_Edd = 1.28 × 10³⁸ erg s⁻¹, $\dot{m}$ is the Eddington-normalized accretion rate, and $\mathcal{M}=M_{\mathrm{BH}}/M_{\odot }$. The magnetic energy density at distance z then scales as $U_{\mathrm{B}}\propto {\eta }_{\mathrm{b}}\dot{m}{\mathcal{M}}^{-1}$, with η_b a proportionality constant. We assume that the external photon field originates from accretion-disk photons scattered by the MHD wind. The scattering region extends from a radial distance R₁ to R₂ along the wind, with electron density n(r)=n₀(r/R_s)⁻¹ following our MHD wind prescription (Contopoulos & Lovelace 1994). The energy density of scattered photons in the jet comoving frame is U_ext = Γ²U_sc, where U_sc ≃ L_diskτ_T/(4πR₂²c), and the wind Thomson depth is τ_T = n₀σ_TR_sln(R₂/R₁), where $n_0\propto \dot{m}$. Here $L_{\mathrm{disk}}\propto ϵ{\dot{m}}^{\alpha }\mathcal{M}{L}_{\mathrm{Edd}}$, with α = 1 for radiatively efficient disks (FSRQs) and α = 2 for ADAF-like disks (BL Lac obects), and ϵ is the radiative efficiency. Therefore, the photon field energy density scales as $U_{\mathrm{ext}}\propto {\mathrm{\Gamma }}^2ϵ{\dot{m}}^{\alpha +1}{\mathcal{M}}^{-1}$.

The emission is assumed to arise from a stationary spherical zone of radius R_b located at some distance z along the jet. The blob moves with bulk Lorentz factor Γ at angle θ to the line of sight, giving Doppler factor δ = [Γ(1 − βcosθ]⁻¹. Particles are accelerated on a timescale $t_{\mathrm{acc}}\propto \gamma A_{\mathrm{acc}}(\dot{m})$. In the framework of diffusive (Fermi I) shock acceleration, the acceleration timescale depends on the particle gyroradius and the shock speed, and can be indirectly linked to the accretion rate through the magnetic field strength, since $B\propto {\dot{m}}^{1/2}$ in standard jet–disk coupling scenarios. We initially tested this physical dependence, but found that it could not reproduce the observed trend of softer spectral peaks in high accretion rate sources. To investigate this effect without introducing additional free parameters, we therefore adopted a phenomenological scaling in which the acceleration timescale increases linearly with the accretion rate, $t_{\mathrm{acc}}\propto \dot{m}$. This prescription effectively captures the possibility that more powerful (high $\dot{m}$) systems may experience slower particle acceleration, for example due to reduced turbulence or more complex shock geometries.

The electron distribution in the acceleration region is governed by the kinetic equation

$$\begin{array}{c}\hfill \frac{\partial n_e(\gamma ,t)}{\partial t}+\frac{n_e(\gamma ,t)}{t_{\mathrm{esc}}}=Q_{\mathrm{inj}}+{\mathcal{L}}_{\mathrm{syn}}(\gamma ,t)+{\mathcal{L}}_{\mathrm{ICS}}(\gamma ,t)+\frac{\partial }{\partial \gamma }\left[\frac{\gamma n_e(\gamma ,t)}{t_{\mathrm{acc}}}\right],\end{array}$$(1)

where n_e(γ, t) is the differential distribution, t_esc the escape timescale, Q_inj the injection term, and ℒ_syn and ℒ_ICS the synchrotron and inverse Compton loss terms (Mastichiadis & Kirk 1995). Solving the coupled electron–photon equations for $\dot{m}\le 0.1$, we find that reproducing the blazar sequence requires t_acc to scale nearly linearly with $\dot{m}$ (Fig. 1, inset).

Fig. 1. Bottom: Calculated spectral energy distributions of BL Lac objects for various values of the normalized accretion rate $\dot{m}$. Inset (top center): Dependence of the acceleration timescale tacc on $\dot{m}$, showing the nearly linear scaling required to explain BL Lac sources. The colored points correspond to the SED curves, indicating that higher $\dot{m}$ produces more luminous sources.

This section establishes a one-zone leptonic framework where both the acceleration timescale and the external photon field depend on the mass accretion rate. In the next section we extend this approach by embedding the acceleration region in a larger cooling zone, moving from a one-zone description to a two-zone emission model.

3. Particle escape: A two-zone emission model

In order to take into account that the high-energy electrons escaping the acceleration zone still undergo ICS losses from the scattered disk photons on the wind, we extend the one-zone model presented in the previous section onto a two-zone scenario. Thus, we consider that particles are accelerated in a compact region near the central engine (acceleration zone, zone I), then escape into a larger volume farther along the jet (cooling zone, zone II), where they continue to radiate via synchrotron and inverse Compton processes (Fig. 2).

Fig. 2. Sketch of the two-zone model. Particles accelerate and radiate in the acceleration zone (dark blue). A fraction escapes into the cooling zone (light blue), where it continues to radiate.

The kinetic equation for electrons in the acceleration zone is described by Eq. (1). In the cooling zone, the kinetic equation is

$$\begin{array}{c}\hfill \frac{\partial n_{e_{\mathit{II}}}(\gamma ,t)}{\partial t}+\frac{n_{e_{\mathit{II}}}(\gamma ,t)}{t_{{\mathrm{esc}}_{\mathit{II}}}}=Q_{\mathrm{inj}}+{\mathcal{L}}_{\mathrm{loss}}(\gamma ,t),\end{array}$$(2)

with t_escII = R_II/c set by the light-crossing time of the cooling zone, and $Q_{\mathrm{inj}}=\frac{n_{e_I}}{t_{{\mathrm{esc}}_I}}$ representing electrons escaping from zone I and injected into zone II. The acceleration and cooling regions in our model differ in size by about one order of magnitude (R_II ≃ 10R_I). Regarding the transition between these zones, we explicitly compared the relevant energy loss timescales. While the volume expansion implies adiabatic cooling, this channel is negligible for the high-energy particle population. The external photon field is assumed constant along the jet, originating from isotropic scattering of accretion-disk photons by the MHD wind (Section 2). The magnetic field decreases with distance as B ∝ 1/z, consistent with magnetic flux conservation in a conical jet.

The two-zone framework is depicted in Fig. 3. BL Lac and FSRQ objects exhibit distinct emission patterns, reflecting their different mass accretion rates and the resulting variations in key physical parameters such as magnetic field, external photon energy density, and acceleration efficiency (see Table A.1). In BL Lac objects the magnetic energy density dominates over the external photon field in both the acceleration and cooling zones. In the acceleration zone, synchrotron emission is therefore the primary radiative process. In contrast, in the cooling zone, although the magnetic field decreases with distance, it still satisfies U_BII ≫ U_ext. As a result, the cooling zone contributes little to the total flux, and the overall spectrum is primarily determined by the acceleration zone (Fig. 3). Instead, FSRQs display a more complex behavior. Synchrotron emission dominates in the acceleration zone, where U_BI >  U_ext, but in the cooling zone the magnetic field decreases. In contrast, the external photon field remains approximately constant, Boula et al. (2019a), leading to U_ext ≫ U_BII. Consequently, external Compton scattering dominates the energy losses in the cooling zone, giving rise to the characteristic two-component SED observed in FSRQs, with low-energy synchrotron emission from the acceleration zone and high-energy external Compton emission from the cooling zone.

Fig. 3. Results for FSRQ and BL Lac objects according to the two-zone model. Straight lines represent the emission from the acceleration zone, and dotted lines represent the emission from the cooling zone. The shaded region depicts the Fermi γ-ray energy band.

Fig. 4. Theoretical blazar sequence according to the superposition of the two-zone emission by varying only the mass accretion rate (see Table A.1 for the values of the input parameters).

In Fig. 4 we present our results for a theoretical blazar sequence in the two-zone model (see Table A.1 for the input parameters and Table A.2 for the resulting spectral peak properties). In all classes, particles escape from the acceleration zone. FSRQs exhibit a strong γ-ray signature because electrons injected into the cooling zone interact with a dense external photon field, whereas in BL Lacs the escaping particles interact only weakly with the external field, and the cooling zone contributes less to the total flux. This two-zone approach is fully consistent with the statistical trends observed by Fermi-LAT, including the Compton dominance, the γ-ray spectral indices, and the positions of the synchrotron and inverse Compton peaks. The values in Tables A.1 and A.2 show that these trends can be reproduced by varying only the mass accretion rate, which emerges as the key physical driver of the sequence, as demonstrated in Boula et al. (2019a). The same study shows that the black hole mass cannot account for the observational spread of the data. Although the black hole mass and the location of the emission region can both affect the details, their influence is subdominant compared to the role of $\dot{m}$. Furthermore, by allowing variations in the acceleration timescale, the electron injection properties, and the magnetic field efficiency, the two-zone framework can reproduce not only the systematic sequence, but also the observed spread in the data, consistent with our previous one-zone study (Boula et al. 2019a, 2019b). This demonstrates that the combination of mass accretion rate scaling and a two-zone treatment of particle escape provides a robust framework for capturing both the systematic and statistical features of the blazar population.

Figure 5 illustrates the application of the two-zone model to the FSRQ 3C273, focusing on the X-ray and γ-ray bands. SSC primarily produces X-ray emission in the compact acceleration zone, where relativistic electrons upscatter synchrotron photons. The cooling zone dominates the γ-ray emission through external Compton scattering, where electrons interact with ambient photons from the MHD wind or the broad-line region. This two-zone scenario naturally explains the distinct emission components and spectral features observed in 3C273. For the parameters adopted (e.g., B_II ≈ 0.1 G, R_II ≈ 3 × 10¹⁷ cm, and u_exp ≈ 0.01c), the synchrotron cooling timescale t_syn is significantly shorter than the adiabatic expansion timescale t_ad ≈ R/u_exp for electrons with γ ≳ 200. Consequently, the radiating electrons lose their energy via synchrotron emission before they undergo significant adiabatic decompression.

Fig. 5. Application of the two-zone model to the FSRQ 3C273, focusing on the high-energy spectral range. The acceleration zone is located at z = 0.01 pc, with a magnetic field strength of B = 1 G and radius RI = 5 × 1015 cm, the cooling zone is one order of magnitude more distant and larger than the acceleration zone. The external photon field arises from isotropic scattering of disk photons on wind particles between radii R1 = 9 × 1014 cm and R2 = 3 × 1018 cm, with energy density Uext = 2.5 × 10−3 erg cm−3. The bulk Lorentz factor is Γ = 30, the Doppler factor is δ = 15, and the characteristic disk temperature is Tdisk = 3 × 104 K. The observational data are from Giommi et al. (2012).

4. Conclusion

For this work we successfully reproduced the theoretical blazar sequence by extending the model of Boula et al. (2019a), using the mass accretion rate as the primary physical parameter to account for the broad diversity observed in blazar phenomenology. This demonstrates that a simple physically grounded approach where changes in accretion rate alone drive the differences between flat-spectrum radio quasars (FSRQs) and BL Lacertae objects (BL Lacs) is consistent with the blazar unification paradigm.

Our model solves the coupled kinetic equations for electrons and photons self-consistently, including Fermi acceleration, synchrotron, and inverse Compton losses, within a two-zone emission framework. Electrons accelerate in a compact region and partially escape into a larger cooling zone, naturally producing the distinct spectral components characteristic of blazars. The competition between magnetic and external photon field energy densities, modulated by the accretion rate, explains the dominance of synchrotron and SSC emission in BL Lacs versus strong external Compton emission in FSRQs.

Importantly, this approach reproduces not only the main trends of the blazar sequence, but also the observed statistical distributions reported by Fermi-LAT, including Compton dominance, γ-ray spectral indices, and the locations of the synchrotron and inverse Compton peaks. Application of the model to the well-studied FSRQ 3C273 captures key features of its broadband SED, confirming that variations in accretion power are sufficient to account for its high-energy emission.

An interesting implication of the two-zone framework is that, in a flaring state, FSRQs may exhibit a time delay in the γ-ray emission relative to the X-rays. Initially, γ-rays are produced simultaneously with X-rays in the acceleration zone, but additional γ-ray emission arises from the cooling zone as electrons propagate and radiate. The duration of the γ-ray pulse is therefore set primarily by the electron cooling timescale since freshly accelerated ones do not immediately replace cooled electrons. Detailed modeling of this effect requires a time-dependent approach. A similar argument applies to the radio emission, which is expected to appear later due to the longer cooling and propagation timescales of the lower-energy electrons.

Overall, our results underscore the central role of the mass accretion rate in shaping blazar spectra within a unified theoretical framework grounded in particle acceleration physics and radiative transfer. The inclusion of particle escape and spatially structured emission zones provides a more realistic interpretation of blazar spectra. Future refinements, including jet stratification, time variability, and polarization studies, will further enhance the predictive power of the model. However, the present work establishes a robust foundation that links accretion physics, jet emission, and blazar phenomenology.

Acknowledgments

We thank the anonymous referee for their useful comments. We thank Yannis Liodakis for the helpful feedback on the manuscript and Fabrizio Tavecchio for valuable discussions. We have used the following Python libraries: Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), Scipy (Virtanen et al. 2020).

References

Appendix A: Parameters and results of the theoretical blazar sequence

In Tables A.1 and A.2 we present the input parameters of the theoretical blazar sequence and the main characteristics of the produced spectral energy distributions that are in agreement with the observational trends of fermi blazars (The Fermi-LAT collaboration 2020).

All Tables

All Figures

Fig. 1. Bottom: Calculated spectral energy distributions of BL Lac objects for various values of the normalized accretion rate $\dot{m}$. Inset (top center): Dependence of the acceleration timescale tacc on $\dot{m}$, showing the nearly linear scaling required to explain BL Lac sources. The colored points correspond to the SED curves, indicating that higher $\dot{m}$ produces more luminous sources.

In the text

	Fig. 2. Sketch of the two-zone model. Particles accelerate and radiate in the acceleration zone (dark blue). A fraction escapes into the cooling zone (light blue), where it continues to radiate.
In the text

	Fig. 3. Results for FSRQ and BL Lac objects according to the two-zone model. Straight lines represent the emission from the acceleration zone, and dotted lines represent the emission from the cooling zone. The shaded region depicts the Fermi γ-ray energy band.
In the text

	Fig. 4. Theoretical blazar sequence according to the superposition of the two-zone emission by varying only the mass accretion rate (see Table A.1 for the values of the input parameters).
In the text

Fig. 5. Application of the two-zone model to the FSRQ 3C273, focusing on the high-energy spectral range. The acceleration zone is located at z = 0.01 pc, with a magnetic field strength of B = 1 G and radius RI = 5 × 1015 cm, the cooling zone is one order of magnitude more distant and larger than the acceleration zone. The external photon field arises from isotropic scattering of disk photons on wind particles between radii R1 = 9 × 1014 cm and R2 = 3 × 1018 cm, with energy density Uext = 2.5 × 10−3 erg cm−3. The bulk Lorentz factor is Γ = 30, the Doppler factor is δ = 15, and the characteristic disk temperature is Tdisk = 3 × 104 K. The observational data are from Giommi et al. (2012).

In the text