# Simulations of Stochastic Long-term Variability in Leptonic Models for External-Compton and Synchrotron Self-Compton Dominated Blazars

Hannes Thiersen, Michael Zacharias, Markus Böttcher

► To cite this version:

Hannes Thiersen, Michael Zacharias, Markus Böttcher. Simulations of Stochastic Long-term Variability in Leptonic Models for External-Compton and Synchrotron Self-Compton Dominated Blazars. The Astrophysical Journal, 2022, 925, 10.3847/1538-4357/ac4013. insu-03713231

HAL Id: insu-03713231

https://insu.hal.science/insu-03713231

Submitted on 4 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# Simulations of Stochastic Long-term Variability in Leptonic Models for External-Compton and Synchrotron Self-Compton Dominated Blazars

Hannes Thiersen<fnref n="1" />,*, Michael Zacharias<fnref n="1" />,<fnref n="2" />,†, and Markus Böttcher<fnref n="1" />,‡

## Abstract

In this work we investigate the nature of multiwavelength variability of blazars from a purely numerical approach. We use a time-dependent one-zone leptonic blazar emission model to simulate multiwavelength variability by introducing stochastic parameter variations in the emission region. These stochastic parameter variations are generated by Monte Carlo methods and have a characteristic power-law index of $\alpha = -2$ in their power spectral densities. We include representative blazar test cases for a flat spectrum radio quasar and a high-synchrotron peaked BL Lacertae object for which the high-energy component of the spectral energy distribution is dominated by external-Compton and synchrotron self-Compton emission, respectively. The simulated variability is analyzed in order to characterize the distinctions between the two blazar cases and the physical parameters driving the variability. We show that the variability's power spectrum is closely related to underlying stochastic parameter variations for both cases. Distinct differences between the different progenitor variations are present in the multiwavelength cross-correlation functions.

Unified Astronomy Thesaurus concepts: Blazars (164); Jets (870); Radiative processes (2055)

## 1. Introduction

Blazars are a type of active galactic nuclei (AGN) that possess a relativistic jet of particles moving on an axis aligned close to our line of sight (Urry & Padovani 1995). They have a characteristic nonthermal spectral energy distribution (SED) consisting of a low-energy emission component ranging from radio frequencies up to infrared or X-rays and a high-energy component that spans from X-rays up to γ-rays, in some cases extending into the Very High-Energy (VHE, $E > 100$ GeV) regime. The low-energy emission component is produced through synchrotron radiation of relativistic electrons (and positrons). The dominant mechanism producing the high-energy component on the other hand is still under debate, but is typically explained as inverse-Compton (IC) scattering of ambient photons by the jet electrons within leptonic models. Alternatively, in hadronic models, the high-energy emission is produced through relativistic proton synchrotron and radiation related to proton-photon interactions (Romero et al. 2017). Blazars are subdivided into two subclasses based on optical emission-line features: flat spectrum radio quasars (FSRQs) possessing broad optical emission lines and BL Lac objects for which these features are effectively absent. Another criterion for further classification of BL Lac objects is the value of their SED synchrotron peak frequency for which there are three types: low-synchrotron peaked (LBL; $\nu_{sync,peak} < 10^{14}$ Hz), intermediate-synchrotron peaked (IBL; $10^{14}$ Hz $\le \nu_{sync,peak} < 10^{15}$ Hz), and high-synchrotron peaked (HBL; $\nu_{sync,peak} > 10^{15}$ Hz) BL Lac objects (Abdo et al. 2010a).

Broadband variability on a vast range of timescales is a ubiquitous characteristic of AGN, and in particular, blazars (Wagner & Witzel 1995; Albert et al. 2007; Abdo et al. 2010b;

H.E.S.S. Collaboration et al. 2017). However, in spite of a large number of observational multiwavelength campaigns on many blazars and extensive theoretical studies of blazar variability, the causes of the variability are still unclear. Most research on blazar variability focuses on some significant short-term events, typically flares (e.g., Aharonian et al. 2008; Arlen et al. 2013; Blinov et al. 2018). These flares may be caused either by changes intrinsic to the jet, such as shocks and/or magnetic reconnection (Sironi et al. 2015; Baring et al. 2017; Christie et al. 2019), or external to the jet. Examples of the latter are the ablation of a dust cloud (Zacharias et al. 2017) or interactions of the jet with a stellar wind (Araudo et al. 2010; Barkov et al. 2012). However, such flaring events might not be representative of the underlying long-term variability, which is still present even in the absence of flares.

Extensive multiwavelength campaigns on blazars allowed for investigations into the long-term variability in blazars providing in-depth analyses and characterization of observed variability behavior (e.g., detections of quasi-periodic oscillations and colored noise characteristics of multiwavelength variability (see Bhatta & Zola 2016; Goyal et al. 2017; Goyal 2018; Bhatta & Dhital 2020). However, in this work, we study the underlying long-term blazar variability from a purely theoretical/numerical perspective. We simulate the multiwavelength variability by assuming a stochastic process inducing time-dependent changes to some emission region variables. These parameter changes are fed into a time-dependent leptonic one-zone code in order to obtain the variable light curves in various energy bands. Lastly, we derive the power density spectra for the individual light curves and the cross-correlation functions between the different energy bands. Two representative test cases of blazars are considered: (a) an FSRQ, dominated by IC scattering of external photon fields (external Compton; EC) in its high-energy SED component and (b) an HBL, dominated by synchrotron self-Compton (SSC) emission.

Similar work has been done by Mastichiadis et al. (2013) and Polkas et al. (2021). Mastichiadis et al. (2013) investigated

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Centre for Space Research, North-West University, Potchefstroom, 2520, South Africa; hannesthiersen@gmail.com</footnote>

<footnote n="2">²Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, CNRS, Université de Paris, F-92190 Meudon, France
Received 2021 April 22; revised 2021 November 11; accepted 2021 November 30; published 2022 February 3</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Broadband spectral energy distributions for the representative FSRQ (left) and HBL (right) blazar cases. Spectral components are as indicated in the legend. The shaded areas indicate the frequency ranges for which the optical, X-ray, HE γ-ray, and VHE γ-ray light curves have been derived. We note that there is an overlap in the frequency ranges of the HE and VHE γ-rays.</figcaption></figure>

variability in X-rays and γ-rays for one-zone leptohadronic emission scenarios in which the high-energy SED component is produced by synchrotron radiation of protons or secondary particles from photopair/photopion production, and/or the decay of neutral pions. Random walks were used to generate variations of injection luminosity and maximum particle energy introduced into their emission model to induce variability. Polkas et al. (2021) studied multiwavelength variability in a leptonic framework. Time-dependent variations were generated for magnetic field strength, Doppler factor, external radiation compactness, and electron injection compactness by means of mimicking observed Fermi-LAT light curves. Here, we focus solely on leptonic emission processes, consider wavelengths from optical up to teraelectronvolt γ-rays, and generate independent variations for the electron injection luminosity, the magnetic field strength, the electron injection spectral index, and the maximum electron energy by means of Fourier transforms.

An overview of the blazar emission model, the representative blazar test cases, and the generation of the stochastic variations is given in Section 2. Results are presented and discussed in Section 3.

## 2. Model

### 2.1. Model Setup

This work employs the time-dependent homogenous one-zone leptonic blazar emission model created by Diltz & Böttcher (2014, 2016) and further developed by Zacharias et al. (2017). We refer to these papers for the details of the code. It assumes that the multiwavelength emission is produced by a homogenous region/zone in which a rapid acceleration process (e.g., shocks or magnetic reconnection) injects relativistic electrons with a power-law distribution in energy. The emission region travels along the jet axis at relativistic speed and is pervaded by a tangled magnetic field of uniform strength.

The model uses Fokker-Planck and radiation transfer equations to keep track of the electron and photon populations within the emission region over time. This allows for the extraction of snapshot SEDs and photon-energy-dependent light curves. The relativistic electrons in the emission region produce synchrotron and IC radiation which then, due to the

emission region's relativistic motion, produce Doppler boosted emission in the observer's frame of reference.

Figure 1 shows the SEDs of the two different representative blazar test cases. The parameter sets used to generate those SEDs are defined and listed in Table 1. The low-energy components in both cases are produced by synchrotron radiation (in the FSRQ case from radio up to optical frequencies; in the HBL case from radio up to X-rays). In the FSRQ case X-ray emission (low-energy end of the high-energy component) is dominated by SSC while the γ-rays are dominated by the IC scattering of the broadband region (BLR) photon field. In the HBL case the low-energy synchrotron photons serve as target photons for IC scattering producing the high-energy component through SSC (soft γ-rays up to VHE γ-rays).

The shaded areas in Figure 1 indicate the frequency ranges over which fluxes have been integrated to obtain light curves. The selected frequency ranges are motivated by the ranges of filters and instruments. This also implies the existence of archival light curves from extensive monitoring campaigns on many blazars in these bands. The ranges are the R band in the optical regime, 0.2–10 keV in X-rays (as typically observed with, e.g., Swift-XRT), 20 MeV–300 GeV in high-energy (HE) γ-rays (the sensitivity range of the Fermi-LAT), and 100 GeV–10 TeV in VHE γ-rays (as observed with, e.g., H.E.S.S, MAGIC, and VERITAS). Although detailed light curves also exist in radio wavelengths (approximately gigahertz), we do not consider them in this work since the emission region considered here is optically thick at those frequencies, and most of the radio emission is produced in the larger-scale jet structures.

### 2.2. Stochastic Variations

In order to produce the light curves, we introduce a stochastic process that results in the variation of one of the emission region parameters. This work does not attempt to explain the stochastic process itself. It is merely assumed that a process exists that results in the stochastic variations of individual source parameters considered here. One may speculate that the stochastic variation of the jet parameters is related to the stochastic variations observed in accretion disks (Chatterjee et al. 2018).

The stochastic variations are generated with the algorithm developed by Timmer & Koenig (1995). In the generation
<table><caption>**Table 1**
Initial Model Parameters for the Representative FSRQ and HBL Blazar Cases</caption><thead><tr><th>Definition</th><th>Symbol</th><th>FSRQ</th><th>HBL</th></tr></thead><tbody><tr><td>Magnetic field</td><td>B<sub>0</sub></td><td>1.70 G</td><td>0.40 G</td></tr><tr><td>Blob radius</td><td>R</td><td>3.0 × 10<sup>16</sup> cm</td><td>3.0 × 10<sup>16</sup> cm</td></tr><tr><td>Ratio of the acceleration to escape timescales</td><td>η</td><td>1.00</td><td>1.00</td></tr><tr><td>Escape timescale</td><td>t<sub>esc</sub></td><td>10.0 R/c</td><td>10.0 R/c</td></tr><tr><td>Redshift to the source</td><td>z</td><td>1.0</td><td>1.0</td></tr><tr><td>Minimum Lorentz factor of the electron injection spectrum</td><td>γ<sub>min</sub></td><td>1.0 × 10<sup>2</sup></td><td>1.0 × 10<sup>4</sup></td></tr><tr><td>Maximum Lorentz factor of the electron injection spectrum</td><td>γ<sub>max,0</sub></td><td>1.0 × 10<sup>4</sup></td><td>1.0 × 10<sup>6</sup></td></tr><tr><td>Bulk Lorentz factor</td><td>Γ</td><td>20.0</td><td>20.0</td></tr><tr><td>Observing angle relative to the axis of the BH jet</td><td>θ<sub>obs</sub></td><td>5.0 × 10<sup>-2</sup> rad</td><td>5.0 × 10<sup>-2</sup> rad</td></tr><tr><td>Doppler factor</td><td>δ</td><td>20.0</td><td>20.0</td></tr><tr><td>Electron injection index</td><td>q<sub>0</sub></td><td>2.8</td><td>2.5</td></tr><tr><td>Co-moving injection luminosity of the electron spectrum</td><td>L<sub>inj,0</sub></td><td>5.0 × 10<sup>43</sup> erg s<sup>-1</sup></td><td>1.0 × 10<sup>42</sup> erg s<sup>-1</sup></td></tr><tr><td>Mass of the super massive black hole</td><td>M<sub>BH</sub></td><td>8.5 × 10<sup>8</sup> M<sub>☉</sub></td><td>8.5 × 10<sup>8</sup> M<sub>☉</sub></td></tr><tr><td>Eddington ratio</td><td>l<sub>edd</sub></td><td>1.0 × 10<sup>-1</sup></td><td>1.0 × 10<sup>-4</sup></td></tr><tr><td>Initial location of the blob along jet axis</td><td>d</td><td>6.5 × 10<sup>17</sup> cm</td><td>6.5 × 10<sup>17</sup> cm</td></tr><tr><td>Radius of the BLR</td><td>R<sub>BLR</sub></td><td>6.7 × 10<sup>17</sup> cm</td><td>...</td></tr><tr><td>Effective temperature of the BLR</td><td>T<sub>eff</sub></td><td>5.0 × 10<sup>4</sup> K</td><td>...</td></tr><tr><td>Effective luminosity of the BLR</td><td>L<sub>BLR</sub></td><td>1.0 × 10<sup>45</sup> erg s<sup>-1</sup></td><td>...</td></tr></tbody></table>

Note. These parameters remain fixed unless they represent the varying parameter of the respective simulation realization (see Section 2.2).

procedure an underlying spectrum—specifically, a pure power law in this work—is assumed for the stochastic process’ power spectral density (PSD) in temporal frequency. Using the underlying PSD spectrum as a baseline a unique periodogram<fnref n="3" /> is produced by means of a Monte Carlo scheme in which complex numbers are generated for the periodogram given the power of the underlying PSD at a given temporal frequency. This effectively produces a noisy version of the underlying spectrum, and in the case of an underlying pure power-law spectrum, conserves the power-law index. A complex valued time series is then obtained from the inverse Fourier transform of the generated periodogram whose real values are considered to create parameter variations.

The observed light-curve variability of accreting systems such as AGN and X-ray binaries has shown evidence of colored noise (e.g., Park & Trippe 2014; Nilsson et al. 2018; Goyal 2021). Therefore, in this work all variations were generated for an underlying PSD being a pure power law with an index of $\alpha = -2$ such that $P(f) \propto f^{\alpha}$, where $P(f)$ represents the PSD as function of temporal frequency $f$.

For each generated time series the parameter variations are obtained for the maximum electron Lorentz factor, $\gamma_{max}$, the electron injection luminosity, $L_{inj}$, the magnetic field strength, $B$, and the electron injection spectral index, $q$, and introduced into the emission region in independent realizations of the simulation. These parameters are naturally expected to be subject to changes over time making them prime candidates for causing variability. Induced variability from the underlying accretion flow will most likely impact jet power, and thus, injection luminosity. Magnetic fields are significantly affected by shock compression and/or magnetic reconnection, which in turn will lead to variability. Particle spectral indices and electron Lorentz factors are characteristic of the considered acceleration mechanisms and can even vary for the same mechanisms depending on the state of the environment (e.g.,

shock obliquity in the case of mildly relativistic shocks; see Summerlin & Baring 2012).

The parameter variations are generated by assuming that the time series represents the scale factor of the initial parameter value at each time step, except for the spectral index variations for which the time series represents the constant added to the initial spectral index value instead. The time series is further manipulated by adding a single constant value throughout to produce maximum Lorentz factor, injection luminosity, and magnetic field variations such that the occurrence of non-physical parameter values is prevented (ensuring that $\gamma_{max} > \gamma_{min}$ and $L_{inj}, B > 0$). In the case of spectral index variations the time series is normalized within the range $[-1,1]$ for the same reason. The overall minimum and maximum values of the aforementioned parameters during the variations are as follows: $\gamma_{max}/\gamma_{max,0} \in [0.1, 5.6]$; $L_{inj}/L_{inj,0} \in [0, 5.5]$; and $B/B_0 \in \{0.1, 5.6\}$, where $B_0$ represents the respective parameter’s initial value in Table 1. The minimum and maximum spectral index values are in the range $q \in [q_0 - 1, q_0 + 1]$. Note that the aforementioned modifications to time series do not change the PSD power-law index of the underlying parameter variations. An example of a variation generated with this method is shown in Figure 2.

The time and frequency intervals of the variations produced by the algorithm may be chosen arbitrarily. In the simulations a total of 4000 steps have been used for a fixed time interval of 2 hr in the emission region’s rest frame for both blazar cases, corresponding to 12 minute intervals in the observer’s frame. A total of 100 time series with different random seeds were generated and used to simulate 100 individual realizations for each permutation of blazar case and varied parameter (eight permutations in total).

# 3. Results and Discussions

## 3.1. Light Curves

Figure 3 shows the resulting multiwavelength light curves for one representative realization for each of the eight cases (FSRQ versus HBL; $\gamma_{max}$, $L_{inj}$, $B$, $q$ variation). The figure

<footnote n="3">³ Also typically known as the Fourier transform of a time series that contains power and phase information for each temporal frequency.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Example of a generated stochastic variation. Left: the PSD obtained from the produced periodogram (blue) follows the underlying pure power law (red). Right: the variation obtained from the inverse fast Fourier transform of the periodogram exhibiting red noise. Specifically, this represents the value of $L_{\text{inj}}$ compared to the initial value at each time step.</figcaption></figure>

illustrates that variations in injection luminosity, magnetic field, and spectral index variations produce significant flux variability across all photon frequencies with flux changes of $\pm \sim 50\%$ compared to the average flux. Qualitatively, the variability produced by the different parameter variations is also distinct from one another. Note that the VHE $\gamma$-ray light curve is omitted from the FSRQ results due to its insignificant flux (see Figure 1).

Flux changes due to $\gamma_{\text{max}}$ variations are less volatile with only about $\pm 5\% - 10\%$. However, FSRQ optical and HBL VHE $\gamma$-ray light curves show the most significant drops in flux at the point where $\gamma_{\text{max}}$ is at a minimum. This sensitivity indicates that these respective wavelengths are produced by the highest-energy electrons in the electron population.

Correlated variability for injection luminosity variations can be seen for both the FSRQ and HBL cases evident from the synchronous flux increases and decreases across all wavelengths due to uniform changes in the number of radiating electrons across the whole electron spectrum. Note that the optical and HE $\gamma$-ray light curves in the FSRQ case are practically identical, which is a consequence of being produced by electrons of similar energies.

In the case of magnetic field variations for the FSRQ the optical and X-rays show high sensitivity to such changes, whereas the reaction in HE $\gamma$-rays is weaker. This is expected since the HE $\gamma$-rays are dominated by EC processes that depend only indirectly on the magnetic field through increased synchrotron cooling of the radiating electrons. The HBL case shows comparable variability in all wavelengths; however, in this specific run the SSC radiation is much more sensitive to the magnetic field changes than the synchrotron.

The prominent variability patterns produced by spectral index variations in both blazars are similar to the variability exhibited from variations in injection luminosity since the normalization of the electron distribution depends on both the injection luminosity and the spectral index (see Equation (B2) in Appendix B). However, unlike for the injection luminosity the spectral index variations induce an energy-dependent

variation of the electron distribution resulting in different amplitudes of the variability in the light curves. Only the EC radiation components in the FSRQ show little variability over the course of the parameter variation. This is because the $\gamma$-ray spectrum pivots around a central energy in the high-energy regime and the X-rays are produced by SSC radiation of low-energy electrons. A lowering (hardening) of the spectral index enhances the synchrotron target photon field for SSC emission but depletes low-energy electrons, so that the two effects almost cancel each other out.

## 3.2. PSDs

The average results for the light-curve PSDs of the 100 different realizations for each case are presented in Figure 4 as solid lines while the dotted lines represent the cooling timescales (see Tables 2 and 3) for their respective wavelengths, as indicated by the legend. The individual simulation realizations in comparison to the averages can be found in Figures 6 and 7 in Appendix C, respectively. The average PSDs follow closely the underlying power law of the parameter variation and do not show significant differences in PSD power-law index for the different wavelengths. This result is similar to what Finke & Becker (2014, 2015) found in their analytical analysis of variability in Fourier space with respect to matching PSD power-law indices for different wavelengths. They also predict spectral breaks between $\sim 10^{-6}$ and $\sim 10^{-4}$ Hz that arise from frequencies related to the light-crossing and escape timescales of electrons and photons in the emission region. However, these clear breaks do not appear in our results because these characteristic frequencies, $t_{\text{esc}} \approx 10^7$ s, corresponding to $f_{\text{esc}} \approx 10^{-7}$ Hz, fall outside the bounds of frequencies considered in this work for the chosen model parameters (see Table 1).

Similar PSD power-law indices across the different wavelengths supplement the argument that one-zone models lack the complexity to accurately produce observed variability patterns of blazars, which have been found to exhibit differently colored
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Simulated light curves exhibiting variability for the same stochastic variation applied to different emission region parameters as indicated. Refer to the symbols in Table 1. Note that the optical and HE γ-ray light curves are practically identical for the FSRQ case and injection luminosity variation.</figcaption></figure>

noise between the low-energy synchrotron and high-energy emission components (see Goyal 2020, and references therein). The lack of distinctive differences in the PSDs also makes it unlikely to be an effective diagnostic tool for identifying progenitor and/or cooling mechanisms. However, since they do follow the underlying variation's power law closely they can be a useful tool to probe the stochastic behavior of variations that cause variability.

Notably, several individual light-curve PSDs smooth out at the higher frequencies, which highly contrasts the Monte Carlo generated variations (see Figure 2 with Figures 6 and 7). This is a consequence of the delayed response of the emission processes that cannot keep up with the changes induced by the variations on the shortest timescales. This is especially evident for the emission processes with longer cooling timescales; see, for example, the X-ray PSD of the FSRQ case (Figure 6) and practically all PSDs of the HBL case (Figure 7).

<figure><img src="image_7.png" /></figure>

### 3.3. Cross-correlation Functions

The average multiwavelength cross correlations of the 100 simulation realizations for each case are presented in Figure 5. The individual realizations in comparison to the averages can be found in Figures 8 and 9–11, and 12 in Appendix C, respectively. Furthermore the concentrations of cross-correlation peaks of all realizations are presented in Appendix D with the delay times and correlation strengths of the average peaks listed in Tables 4 and 5. The results show strong distinctive features across the different scenarios.

The maximum Lorentz factor variation shows the least consistent behavior between all simulations in both blazars evident from the weak correlations on average and the large deviations between individual realizations. The only exception in this regard is the HBL HE γ-ray versus optical cross-correlation function, which does present a consistent feature peaking at −30 ± 23 hr with correlation strength 0.8 ± 0.3.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Average light-curve PSDs for the FSRQ (left) and HBL (right) cases for the different parameter variations. Vertical dotted lines indicate average characteristic cooling frequencies (see Tables 2 and 3) corresponding to the waveband indicated by color in the legend.</figcaption></figure>

<table><caption>**Table 2**
FSRQ Cooling Timescales in the Observer's Frame</caption><thead><tr><th>Waveband/Energy</th><th>Timescale (s)</th><th>Temporal Frequency (Hz)</th></tr></thead><tbody><tr><td>R band (658 nm)</td><td>4.42 × 10<sup>2</sup></td><td>2.26 × 10<sup>-3</sup></td></tr><tr><td>X-ray (0.2 keV)</td><td>2.75 × 10<sup>5</sup></td><td>3.64 × 10<sup>-6</sup></td></tr><tr><td>X-ray (10 keV)</td><td>3.89 × 10<sup>4</sup></td><td>2.57 × 10<sup>-5</sup></td></tr><tr><td>HE γ-ray (100 MeV)</td><td>6.35 × 10<sup>3</sup></td><td>1.57 × 10<sup>-4</sup></td></tr><tr><td>HE γ-ray (30 GeV)</td><td>3.67 × 10<sup>2</sup></td><td>2.73 × 10<sup>-3</sup></td></tr></tbody></table>

This is due to the HE γ-rays and optical photons being produced by the low-energy electrons in the HBL that are similarly, weakly affected by changes in $γ_{max}$ leading to a positive correlation. The fact that the optical precedes the HE γ-rays stems from HE γ-rays being SSC radiation, which depends on the optical synchrotron photons, even though the cooling timescales of electrons radiating in these wavebands are comparable. A similar counterpart is absent in the FSRQ case.

<table><caption>**Table 3**
HBL Cooling Timescales in the Observer's Frame</caption><thead><tr><th>Waveband/Energy</th><th>Timescale (s)</th><th>Temporal Frequency (Hz)</th></tr></thead><tbody><tr><td>R band (658 nm)</td><td>5.57 × 10<sup>4</sup></td><td>1.79 × 10<sup>-5</sup></td></tr><tr><td>X-ray (0.2 keV)</td><td>5.41 × 10<sup>3</sup></td><td>1.85 × 10<sup>-4</sup></td></tr><tr><td>X-ray (10 keV)</td><td>7.65 × 10<sup>2</sup></td><td>1.31 × 10<sup>-3</sup></td></tr><tr><td>HE γ-ray (100 MeV)</td><td>1.02 × 10<sup>5</sup></td><td>9.78 × 10<sup>-6</sup></td></tr><tr><td>HE γ-ray (30 GeV)</td><td>5.90 × 10<sup>3</sup></td><td>1.69 × 10<sup>-4</sup></td></tr><tr><td>VHE γ-ray (200 GeV)</td><td>7.41 × 10<sup>3</sup></td><td>1.35 × 10<sup>-4</sup></td></tr><tr><td>VHE γ-ray (10 TeV)</td><td>1.48 × 10<sup>2</sup></td><td>6.74 × 10<sup>-3</sup></td></tr></tbody></table>

The injection luminosity variations produce consistently strong correlations (≈1.0), which deviate little in time (standard deviations ≤2 hr) between all wavebands for both blazars. The time delays of each peak are indicative of the differences in cooling timescales of electrons radiating in the respective wavebands. The FSRQ optical and HE γ-rays produced by
<figure><img src="image_10.png" /><figcaption>Figure 5. Average multiwavelength cross correlations for the EC (left) and SSC (right) cases for the different parameter variations.</figcaption></figure>

<table><caption>Table 4
FSRQ Case Cross-correlation Peaks</caption><thead><tr><th>Variation Parameter</th><th>Wavelengths</th><th>Time Delay (hr)</th><th>Correlation Strength</th></tr></thead><tbody><tr><td rowspan="3">γ<sub>max</sub></td><td>HE γ-rays versus X-rays</td><td>12 ± 76</td><td>-0.2 ± 0.5</td></tr><tr><td>HE γ-rays versus optical</td><td>-27 ± 25</td><td>-0.6 ± 0.4</td></tr><tr><td>Optical versus X-rays</td><td>30 ± 76</td><td>0.5 ± 0.4</td></tr><tr><td rowspan="3">Injection luminosity</td><td>HE γ-rays versus X-rays</td><td>12 ± 2</td><td>1.0 ± 8 × 10<sup>-3</sup></td></tr><tr><td>HE γ-rays versus optical</td><td>-0.8 ± 0.1</td><td>1.0 ± 6 × 10<sup>-5</sup></td></tr><tr><td>Optical versus X-rays</td><td>13 ± 2</td><td>1.0 ± 8 × 10<sup>-3</sup></td></tr><tr><td rowspan="3">Magnetic field</td><td>HE γ-rays versus X-rays</td><td>32 ± 7</td><td>-1.0 ± 1 × 10<sup>-2</sup></td></tr><tr><td>HE γ-rays versus optical</td><td>-5 ± 0.9</td><td>-1.0 ± 2 × 10<sup>-3</sup></td></tr><tr><td>Optical versus X-rays</td><td>38 ± 7</td><td>1.0 ± 2 × 10<sup>-2</sup></td></tr><tr><td rowspan="3">Spectral index</td><td>HE γ-rays versus X-rays</td><td>25 ± 33</td><td>0.8 ± 0.3</td></tr><tr><td>HE γ-rays versus optical</td><td>-3 ± 12</td><td>0.9 ± 8 × 10<sup>-2</sup></td></tr><tr><td>Optical versus X-rays</td><td>30 ± 48</td><td>0.7 ± 0.4</td></tr></tbody></table>

high-energy electrons precede the X-rays produced by low-energy electrons since the cooling timescales of electrons radiating in optical and HE γ-rays are shorter than those of electrons radiating in X-rays (see Table 2). The time delays

presented in the HBL cross correlations further confirm this, as the X-ray band, with the shortest cooling timescale, precedes all other wavebands. Although electrons radiating in VHE γ-rays have similar cooling timescales compared to X-rays, the
<table><caption>**Table 5**
HBL Case Cross-correlation Peaks</caption><thead><tr><th>Variation Parameter</th><th>Wavelength</th><th>Time Delay (hr)</th><th>Correlation Strength</th></tr></thead><tbody><tr><td rowspan="6">γ<sub>max</sub></td><td>HE γ-rays versus X-rays</td><td>27 ± 29</td><td>-0.5 ± 0.5</td></tr><tr><td>HE γ-rays versus VHE γ-rays</td><td>-1 ± 72</td><td>0.3 ± 0.5</td></tr><tr><td>HE γ-rays versus optical</td><td>-30 ± 23</td><td>0.8 ± 0.3</td></tr><tr><td>VHE γ-rays versus X-rays</td><td>14 ± 41</td><td>0.5 ± 0.6</td></tr><tr><td>VHE γ-rays versus Optical</td><td>34 ± 46</td><td>-0.1 ± 0.6</td></tr><tr><td>Optical versus X-rays</td><td>25 ± 51</td><td>-0.4 ± 0.5</td></tr><tr><td rowspan="6">Injection luminosity</td><td>HE γ-rays versus X-rays</td><td>-18 ± 2</td><td>1.0 ± 2 × 10<sup>-2</sup></td></tr><tr><td>HE γ-rays versus VHE γ-rays</td><td>-12 ± 2</td><td>1.0 ± 8 × 10<sup>-3</sup></td></tr><tr><td>HE γ-rays versus optical</td><td>-12 ± 2</td><td>1.0 ± 9 × 10<sup>-3</sup></td></tr><tr><td>VHE γ-rays versus X-rays</td><td>-4 ± 1</td><td>1.0 ± 4 × 10<sup>-3</sup></td></tr><tr><td>VHE γ-rays versus optical</td><td>1 ± 1</td><td>1.0 ± 3 × 10<sup>-3</sup></td></tr><tr><td>Optical versus X-rays</td><td>-6 ± 1</td><td>1.0 ± 3 × 10<sup>-3</sup></td></tr><tr><td rowspan="6">Magnetic field</td><td>HE γ-rays versus X-rays</td><td>-9 ± 3</td><td>-1.0 ± 1 × 10<sup>-2</sup></td></tr><tr><td>HE γ-rays versus VHE γ-rays</td><td>-11 ± 5</td><td>1.0 ± 2 × 10<sup>-2</sup></td></tr><tr><td>HE γ-rays versus optical</td><td>-28 ± 34</td><td>-0.5 ± 0.4</td></tr><tr><td>VHE γ-rays versus X-rays</td><td>1 ± 2</td><td>-1.0 ± 9 × 10<sup>-3</sup></td></tr><tr><td>VHE γ-rays versus optical</td><td>-9 ± 22</td><td>-0.7 ± 0.4</td></tr><tr><td>Optical versus X-rays</td><td>12 ± 27</td><td>0.6 ± 0.4</td></tr><tr><td rowspan="6">Spectral index</td><td>HE γ-rays versus X-rays</td><td>-15 ± 5</td><td>-1.0 ± 2 × 10<sup>-2</sup></td></tr><tr><td>HE γ-rays versus VHE γ-rays</td><td>-28 ± 6</td><td>-0.9 ± 3 × 10<sup>-2</sup></td></tr><tr><td>HE γ-rays versus optical</td><td>-15 ± 4</td><td>1.0 ± 9 × 10<sup>-3</sup></td></tr><tr><td>VHE γ-rays versus X-rays</td><td>10 ± 5</td><td>1.0 ± 3 × 10<sup>-2</sup></td></tr><tr><td>VHE γ-rays versus optical</td><td>13 ± 4</td><td>-1.0 ± 2 × 10<sup>-2</sup></td></tr><tr><td>Optical versus X-rays</td><td>-1 ± 3</td><td>-1.0 ± 9 × 10<sup>-3</sup></td></tr></tbody></table>

production due to SSC depends on the synchrotron production. This explains how the optical variability, with a longer electron cooling timescale, can precede VHE γ-ray variability.

The cross-correlation functions for magnetic field variations also produce strong correlations/anticorrelations (≈±1). However, contrary to injection luminosity variations these correlations/anticorrelations show larger deviations in time delays (standard deviations up to 7 hr) among individual realizations. The exception to these results are HBL cross correlations involving optical. In these cases, a consistent behavior in all simulations is observed, even though the cross-correlation peaks are much weaker (≈0.6) and more scattered (time delay standard deviations up to 34 hr). This can be attributed to the fact that the magnetic field changes have their most prominent effect on the optical emission (see Equation (A4)) combined with the long cooling timescale compared to other wavebands leading to weaker correlations/anticorrelations.

Furthermore, both blazars exhibit correlations between wavebands that share a cooling mechanism and anticorrelations when the cooling mechanisms differ.<fnref n="4" /> The optical versus X-ray correlation in the FSRQ is a special case and seems to contrast the results of the HBL case. However, in the FSRQ case, the X-rays are produced by SSC, and are thus linked to the optical emission—which is synchrotron—resulting in a positive correlation. On the other hand, the HE γ-rays are EC emission, which react only weakly to magnetic field changes providing an anticorrelation with both the optical and the X-ray emission.

For the case of spectral index variations the FSRQ and HBL results differ. Correlation functions produced by the FSRQ are

noticeably less consistent compared to the HBL leading to weaker average correlations (≈0.8) over a much wider range of time delays (standard deviations up to 48 hr), while in the HBL case correlations/anticorrelations are strong (|DCF| > 0.9) and time delays are found in a much narrower range (standard deviations ≤6 hr). This discrepancy is due the specific electron energy range producing both FSRQ X-rays and HE γ-rays. These electrons range from the lowest energies up to intermediate energies and beyond compared to the total electron spectrum. This would then include electron energies higher than the *pivot* electron energy, $\epsilon_{pivot}$, below which the number of electrons does not change significantly in spite of changes in electron spectral index. This explains the small amplitude variability in FSRQ X-rays and HE γ-rays compared to optical (see Figure 3). The HBL results show that correlations between wavebands arise for wavebands produced by equivalent energy electrons and anticorrelations arise otherwise. That is *high-energy* versus *low-energy* emissions (e.g., VHE γ-rays versus optical) exhibit anticorrelations.

## 4. Summary and Conclusions

This work presents a method for simulating long-term multiwavelength variability of blazars by means of varying an emission region parameter in a stochastic manner in the framework of a time-dependent one-zone leptonic blazar emission model. The simulated variability for two representative blazar cases, one dominated in γ-rays by EC and the other by SSC, was investigated.

It was found that the simulated variability PSDs closely resemble the PSD shape of the progenitor variation, in this case a pure power law. Hence, according to the model, the produced variability in blazar light curves likely reflects the variability characteristics of the relevant particle acceleration mechanism,

<footnote n="4">⁴ That is in the HBL case, “SSC” versus “synchrotron” (e.g., HE γ-rays versus optical) ⇒ anticorrelation; in the FSRQ case, “EC” versus “synchrotron” ⇒ anticorrelation.</footnote>
which in this case were represented by the stochastic variations. No significant distinctive differences are present among the different test cases. This indicates that the one-zone emission model might be incapable of producing light-curve variability that posses different PSD indices for different wavebands (Goyal 2020).

The PSDs resulting from our simulations did not present any conclusive spectral breaks at characteristic frequencies corresponding to cooling, light-crossing, or escape timescales. This is likely a consequence of the fact that the probed frequency range in this study was limited, and in most cases, did not include the frequencies where such breaks would be expected (see Finke & Becker 2014, 2015). However, the frequencies corresponding to some electron cooling timescales, were within the probed frequency range. Hence, the results on the absence of corresponding spectral breaks seem unexpected and are still inconclusive.

The cross-correlation results showed distinctive features for most of the different scenarios (FSRQ versus HBL, maximum Lorentz factor, injection luminosity, magnetic field, and spectral index variations). As expected, in the case of injection luminosity and spectral index variations, photons of different wavelengths produced by the same electron populations (low energy versus high energy) are usually strongly correlated. The time delays obtained between the different wavelengths are primarily caused by differences in electron cooling timescales. This illustrates that cross-correlation analysis is a viable strategy to identify the photon frequencies produced by equivalent energy electrons and indicate relative cooling timescales when it is known that the acceleration mechanism that drives the variability modifies the electron spectrum in a similar manner.

Overall, the cross-correlation functions presented similar features in both blazars for the respective variations. Varying the maximum Lorentz factor of the electron spectrum produced inconsistent results across all simulations leading to weak correlations/anticorrelations on average. HBL HE γ-rays versus optical are the only cross-correlation functions that presented a consistent feature due to the almost negligible effect of the maximum Lorentz factor variations on these wavelengths. A similar counterpart is not present in the FSRQ cross-correlation comparisons. Electron injection luminosity variations produced all-round strong correlations across all wavebands in all cases. The variations of magnetic field strength showed strong anticorrelations between synchrotron and IC wavebands and strong correlations when wavebands share a common cooling mechanism. However, comparisons with FSRQ X-rays and HBL optical produce noticeably weaker correlations/anticorrelations. A varying electron spectral index leads to strong correlations between wavebands produced by equivalent energy electrons and strong anticorrelations otherwise except for FSRQ X-rays. These distinctions between the different varied parameters suggest the possibility to identify the underlying changing quantities from multiwavelength variability cross correlations.

This provides a compelling argument to use cross-correlation results as a diagnostic tool to constrain dominant radiation mechanisms as well as changing quantities in the emission region that are the cause of the variability. However, we acknowledge that such an analysis would be difficult to implement in practice since observational data are, more often than not irregularly sampled and true simultaneous multiwavelength data with similar temporal sampling are scarce. Therefore, the theoretical results

presented in this paper are practically challenging to test against currently available observational data. Strategies to test our results against observations are currently under development, and we defer their presentation to a future publication.

We thank the anonymous referee for helpful and constructive feedback. The work of M.B. is supported through the South African Research Chair Initiative of the National Research Foundation<fnref n="5" /> and the Department of Science and Innovation of South Africa, under SARChI Chair grant No. 64789. M.Z. acknowledges postdoctoral financial support from LUTH, Observatoire de Paris.

## Appendix A

### Characteristic Cooling Timescales of Radiation Mechanisms

The total energy loss rate, $P(\gamma)$, of individual electrons with energy $\gamma m_e c^2$—where $\gamma$ is the electron Lorentz factor, $m_e$ the electron mass, and $c$ the speed of light—for synchrotron and IC radiation in the Thomson regime can be evaluated as

$$ P(\gamma)^{\text{total}} = P(\gamma)^{\text{sync}} + P(\gamma)^{\text{C}} \\ = -\frac{4}{3}\sigma_{T}c u_{B}(1 + k_{C})\gamma^{2} \quad (\text{A1}) $$

with the Compton dominance parameter $k_C = \frac{u_{ph}}{u_B} \approx \nu F_{\nu}^{p, IC} / \nu F_{\nu}^{p, sy}$, where $\nu F_{\nu}^{p, sy}$ and $\nu F_{\nu}^{p, IC}$ are the $\nu F_{\nu}$ peak fluxes of the synchrotron and Compton emission components, respectively. Furthermore, $\sigma_T$ represents the Thomson cross section and $u_B$ the magnetic field energy density.

The cooling timescale for the electrons in the emission region rest frame is then

$$ T_{\text{cool}}^{\text{em}}(\gamma) = \frac{\gamma m_e c^2}{|P(\gamma)^{\text{total}}|} \\ = \left[ \frac{4}{3} \sigma_T c \frac{u_B}{m_e c^2} (1 + k_C) \gamma \right]^{-1} \quad (\text{A2}) $$

and is related to the observer's frame cooling timescale through

$$ T_{\text{cool}}^{\text{obs}}(\gamma) = \frac{1+z}{\delta} T_{\text{cool}}^{\text{em}}(\gamma), \quad (\text{A3}) $$

where $\delta = (\Gamma [1 - \beta_{\Gamma} \cos \theta_{\text{obs}}])^{-1}$ is the Doppler factor and $z$ the redshift to the source. The electron Lorentz factor $\gamma$ may be related to an observed photon frequency/energy, depending on the radiation mechanism through which any observed frequency $\nu_{\text{obs}}$ is produced.

For frequencies in the synchrotron SED component, we may set $\nu_{\text{obs}}$ equal to the critical synchrotron frequency of an electron with Lorentz factor $\gamma$, such that $\nu_{\text{obs}} = \frac{\delta}{1+z} \frac{3eB}{4\pi m_e c} \gamma^2$ where $e$ is elementary charge. With the magnetic field energy

<footnote n="5">⁵ Any opinion, finding, and conclusion or recommendation expressed in this material is that of the authors, and the NRF does not accept any liability in this regard.</footnote>
density $u_B = \frac{B^2}{8\pi}$, this then yields

$$T_{\text{cool,sync}}^{\text{obs}}(\nu_{\text{obs}}) = \sqrt{\frac{1+z}{\delta}} \times \left[ \frac{1}{3}\sigma_T(1+k_C)\sqrt{\frac{B^3\nu_{\text{obs}}}{3\pi em_e c}} \right]^{-1}. \quad (A4)$$

We assume that frequencies in X-rays up to HE $\gamma$-rays in the IC SED component are within the Thomson regime and approximate the relation between $\gamma$ and $\nu^{\text{obs}}$ with $\epsilon_s^{\text{em}} = \gamma^2 \epsilon_t^{\text{em}}$ where $\epsilon$ represent the photon energy normalized w.r.t. electron rest energy and leads to $\nu_s^{\text{obs}} = \frac{\delta}{1+z} \nu_t^{\text{em}} \gamma^2$ where the subscripts s and t indicates the scattered and target photons respectively. The cooling timescale is then

$$T_{\text{cool,T}}^{\text{obs}}(\nu_{\text{obs}}) = \sqrt{\frac{1+z}{\delta}} \left[ \frac{1}{6}\sigma_T \frac{B^2}{\pi m_e c} (1+k_C) \sqrt{\frac{\nu_{\text{obs}}}{\nu_t^{\text{em}}}} \right]^{-1}. \quad (A5)$$

For frequencies in VHE $\gamma$-rays we approximate IC radiation in the Klein–Nishina regime, $\epsilon = \gamma$ leading to the relation $\nu^{\text{obs}} = \frac{\delta m_e c^2}{1+z h} \gamma$, where h is Planck's constant, giving the following result for the cooling timescale:

$$T_{\text{cool,KN}}^{\text{obs}}(\nu_{\text{obs}}) = \left[ \frac{1}{6} \sigma_T \frac{B^2}{\pi m_e^2 c^3} (1 + k_C) h \nu_{\text{obs}} \right]^{-1}. \quad (A6)$$

<figure><img src="image_11.png" /><figcaption>Figure 6. Average and individual PSD comparisons for FSRQ simulation realizations.</figcaption></figure>

## Appendix B
Injection Luminosity within the One-zone Model

In the model the emission region is continuously injected with electrons at each time step, $t$, where the rate of electrons injected with Lorentz factor $\gamma$ follow a power law:

$$Q(\gamma, t) = Q_0(t)\gamma^{-q(t)}H(\gamma; \gamma_{\text{min}}, \gamma_{\text{max}}) \quad (B1)$$

with the normalization factor $Q_0(t)$, and the Heaviside function $H(\gamma; \gamma_{\text{min}}, \gamma_{\text{max}})$, which is 1 for $\gamma_{\text{min}} \le \gamma \le \gamma_{\text{max}}$ and 0 otherwise, and the electron injection spectral index, $q(t)$. Normalization with respect to the electron injection luminosity, $L_{\text{inj}}$, and the size of the emission region, $V_b$, yields

$$Q_0(t) = \begin{cases} \frac{L_{\text{inj}}(t)}{V_b m_e c^2 \gamma_{\text{max}}^{2-q(t)} - \gamma_{\text{min}}^{2-q(t)}} & \text{if } q(t) \neq 2 \\ \frac{L_{\text{inj}}(t)}{V_b m_e c^2 \ln \frac{\gamma_{\text{max}}}{\gamma_{\text{min}}}} & \text{if } q(t) = 2 \end{cases}. \quad (B2)$$

This shows how changing the injection luminosity and/or the electron injection index can have a large affect on the number of electrons injected into the emission region. Increasing the electron injection luminosity increases the number of electrons of all energies equally, while a change in the spectral index results in an energy-dependent variation of the electron density.

## Appendix C
Individual Simulation Realization Results

A collection of all individual simulation realization PSD (Figures 6 and 7) and cross-correlation (Figures 8 to 12) results plotted against the calculated averages for the respective cases.
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Average and individual PSD comparisons for HBL simulation realizations.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 8. Average and individual cross-correlation comparisons for FSRQ simulation realizations.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 9. Average and individual cross-correlation comparisons for HBL simulation realizations with maximum Lorentz factor variations.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 10. Average and individual cross-correlation comparisons for HBL simulation realizations with injection luminosity variations.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 11. Average and individual cross-correlation comparisons for HBL simulation realizations with magnetic field variations.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Average and individual cross-correlation comparisons for HBL simulation realizations with spectral index variations.</figcaption></figure>
## Appendix D
### Cross-correlation Peak Analyses

Figure 13 shows the peak in each cross-correlation function within a time delay range of [−100, 100] hr for each of the

simulations from which average cross-correlation peaks (crosses) are calculated.

<figure><img src="image_18.png" /><figcaption>Figure 13. Individual cross-correlation peaks of all simulations as well as the average cross-correlation peaks.</figcaption></figure>
ORCID iDs

Hannes Thiersen https://orcid.org/0000-0002-4024-3280

Michael Zacharias https://orcid.org/0000-0001-5801-3945

Markus Böttcher https://orcid.org/0000-0002-8434-5692

References

Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010a, ApJ, 716, 30

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010b, ApJ, 722, 520

Aharonian, F., Akhperjanian, A., Bazer-Bachi, A., et al. 2008, ApJL, 664, L71

Albert, J., Aliu, E., Anderhub, H., et al. 2007, ApJ, 669, 862

Araudo, A. T., Bosch-Ramon, V., & Romero, G. E. 2010, A&A, 522, A97

Arlen, T., Aune, T., Beilicke, M., et al. 2013, ApJ, 762, 92

Baring, M. G., Böttcher, M., & Summerlin, E. J. 2017, MNRAS, 464, 4875

Barkov, M. V., Aharonian, F. A., Bogovalov, S. V., Kelner, S. R., & Khangulyan, D. 2012, ApJ, 749, 119

Bhatta, G., & Dhital, N. 2020, ApJ, 891, 120

Bhatta, G., Zola, S., Stawarz, Ł., et al. 2016, ApJ, 832, 47

Blinov, D., Pavlidou, V., Papadakis, I., et al. 2018, MNRAS, 474, 1296

Chatterjee, R., Roychowdhury, A., Chandra, S., & Sinha, A. 2018, ApJL, 859, L21

Christie, I. M., Petropoulou, M., Sironi, L., & Giannios, D. 2019, MNRAS, 482, 65

Diltz, C., & Böttcher, M. 2014, JHEAp, 1, 63

Diltz, C., & Böttcher, M. 2016, ApJ, 826, 54

Finke, J. D., & Becker, P. A. 2014, ApJ, 791, 21

Finke, J. D., & Becker, P. A. 2015, ApJ, 809, 85

Goyal, A. 2020, MNRAS, 494, 3432

Goyal, A. 2021, ApJ, 909, 39

Goyal, A., Stawarz, Ł., Ostrowski, M., et al. 2017, ApJ, 837, 127

Goyal, A., Stawarz, Ł., Zola, S., et al. 2018, ApJ, 863, 175

H.E.S.S. Collaboration, Abdalla, H., Abramowski, A., et al. 2017, A&A, 598, A39

Mastichiadis, A., Petropoulou, M., & Dimitrakoudis, S. 2013, MNRAS, 434, 2684

Nilsson, K., Lindfors, E., Takalo, L. O., et al. 2018, A&A, 620, A185

Park, J.-H., & Trippe, S. 2014, ApJ, 785, 76

Polkas, M., Petropoulou, M., & Vasilopoulos, G. 2021, MNRAS, 505, 6103

Romero, G. E., Boettcher, M., Markoff, S., & Tavecchio, F. 2017, SSRv, 207, 5

Sironi, L., Petropoulou, M., & Giannios, D. 2015, MNRAS, 450, 183

Summerlin, E. J., & Baring, M. G. 2012, ApJ, 745, 63

Timmer, J., & Koenig, M. 1995, A&A, 300, 707

Urry, C. M., & Padovani, P. 1995, PASP, 107, 803

Wagner, S. J., & Witzel, A. 1995, ARA&A, 33, 163

Zacharias, M., Böttcher, M., Jankowsky, F., et al. 2017, ApJ, 851, 72