Blazars exhibit multiwavelength variability a phenomenon whose underlying mechanisms remain elusive. This study investigates the origin of such variability through leptonic blazar emission simulations focusing on stochastic fluctuations in environmental parameters. By analyzing the spectral indices of the power spectral densities of the variability we assess their relationship with the underlying fluctuations. Our findings reveal that the variability spectral indices remain almost independent of the variations responsible for their emergence. This suggests a complex interplay of factors contributing to the observed multiwavelength variability in blazars.
['Blazars', 'Variable radiation sources', '164', '1759', 'Astrophysics - High Energy Astrophysical Phenomena']
The Relation between Simulated Multiwavelength Blazar Variability and Stochastic Fluctuations
{'The relation between Simulated Multiwavelength Blazar Variability and Stochastic Fluctuations': 'Hannes Thiersen , 1 Michael Zacharias , 2, 1 and Markus Bottcher 1  \n1  \nCentre for Space Research, North-West University, Potchefstroom, 2520, South Africa 2 Landessternwarte, Universitat Heidelberg, Konigstuhl 12, 69117 Heidelberg, Germany  \n(Received ...; Revised ...; Accepted ...)  \nSubmitted to ApJ', 'ABSTRACT': 'Blazars exhibit multiwavelength variability, a phenomenon whose underlying mechanisms remain elusive. This study investigates the origin of such variability through leptonic blazar emission simulations, focusing on stochastic fluctuations in environmental parameters. By analyzing the spectral indices of the power spectral densities of the variability, we assess their relationship with the underlying fluctuations. Our findings reveal that the variability spectral indices remain almost independent of the variations responsible for their emergence. This suggests a complex interplay of factors contributing to the observed multiwavelength variability in blazars.  \nKeywords: -', '1. INTRODUCTION': "Blazars are a subclass of active galactic nuclei known to be some of the brightest persistent multi-wavelength sources in the sky (Urry & Padovani 1995; Abdollahi et al. 2020). Blazars possess jets of relativistic particles propagating along their central region's rotation axis which is closely aligned to our line of sight. The relativistic bulk motion of the emission region leads to highly Doppler-boosted observed emission. Radio up to UV/X-ray radiation is produced by leptonic synchrotron emission, while the mechanisms producing the high-energy spectral component (X-ray through γ -rays, in some cases up to TeV energies) are still up for debate. Some candidate mechanisms include inverse Compton scattering (Begelman & Sikora 1987; Maraschi et al. 1992; Dermer & Schlickeiser 1993; Sikora  \net al. 1994; Ghisellini et al. 2010; Bottcher et al. 2013), proton synchrotron radiation, and photopion production and subsequent decay of pions and muons, initiating electromagnetic cascades (Mannheim & Biermann 1992; Mannheim 1993; Mucke & Protheroe 2000; Bottcher et al. 2013).  \nBlazar emission exhibits unpredictable variability across the entire electromagnetic spectrum. This variability is observed on both short (minutes up to days) and long (weeks up to years) time scales (e.g., Albert et al. 2007; Aharonian et al. 2007; Ackermann et al. 2016; Arlen et al. 2013; Chatterjee et al. 2021). Multiwavelength blazar observation campaigns (e.g., Nilsson et al. 2018; Abdollahi et al. 2020; Pe˜nil et al. 2024; MAGIC Collaboration et al. 2024) elucidated the phenomenology of blazar variability in some detail. However, the exact cause(s)  \nof this phenomenon still elude(s) the scientific community.  \nOne of the characteristics routinely found in long-term observations of blazar variability as well as Galactic X-ray binaries is either simple or broken power-law spectra in the power spectral densities (PSDs) of their light curves (H. E. S. S. Collaboration et al. 2017; Goyal 2020; Bhatta & Dhital 2020; Goyal et al. 2022). Therefore it is plausible that some type of stochastic process is at work producing the multiwavelength variability such that emission power and temporal frequency are anticorrelated, i.e. P ( f ) ∝ f -α , also commonly referred to as coloured noise. Furthermore, some works quantifying the power law indices have found that these indices do not necessarily agree between different wavebands (e.g. Goyal et al. 2022).  \nVarious models have been developed that exhibit colored noise variability, but are often constrained by the duration of particular events, rendering them effective only temporarily. Very few models attempt to explain continuous sustained variability over the large range of observed variability time scales (e.g., Finke & Becker 2015). Recent papers, including Tavecchio et al. (2020); Adams et al. (2022); Brill (2022), employ stochastic differential equations to model high-energy blazar variability with some success, adding to the growing evidence for the role of stochastic processes in these phenomena.  \nThis work attempts to understand how variability can arise assuming some self-sustaining process manipulates key parameters in the emitting region in a stochastic manner. Details of the model and methodology are explained in Section 2. Section 3 presents the main results of the paper followed by a discussion and conclusions in Section 4.", '2.1. Leptonic one-zone Emission Model': "We use a leptonic one-zone blazar emission model to simulate multi-wavelength blazar variability (Diltz & Bottcher 2014; Zacharias et al. 2017, 2022; Thiersen et al. 2022; Zacharias 2023). The model assumes that all emission is produced by a single homogeneous volume/zone located close to the central region, containing isotropically distributed relativistic leptons (electrons and possibly positrons). Due to the relativistic nature of the leptons they produce non-thermal emission via synchrotron (radio up to UV and/or X-ray radiation) and inverseCompton radiation (X-ray up to very highenergy [VHE] γ -rays). Furthermore, the blob travels at relativistic speed with bulk Lorentz factor Γ with respect to the central region. The observer's line of sight is closely aligned to the jet axis ( θ obs < 10 · ). This results in strongly Doppler boosted non-thermal emission, characterized by the Doppler factor δ = (Γ [1 -β Γ cos θ obs ]) -1 , where β Γ is the velocity (normalized to the speed of light) corresponding to Γ. The νF ν peak fluxes are boosted by a factor δ 4 into the observer's frame. Additionally, radiation frequencies are boosted by a factor δ and timescales in the co-moving frame are shortened in the observer's frame by a factor δ -1 .  \nOur code solves the time dependent FokkerPlanck and radiation-transfer equations for the lepton and photon populations in the emission region. The calculations account for radiative cooling from synchrotron emission and inverse Compton scattering of the synchrotron spectrum itself (synchrotron-self Compton; SSC; Konigl (1981); Marscher & Gear (1985); Ghisellini & Maraschi (1989)) as well as external photon fields, i.e. accretion disc, broad-line region and/or dusty torus photons (Ghisellini & Madau 1996; Bottcher et al. 1997; Bglyph[suppress]la˙zejowski et al. 2000; Ghisellini & Tavecchio 2008). Additionally γ -γ pair production and external pho-t  \nbsorption effects are accounted for in the model. Accelerated relativistic leptons are injected at each time step with a power law energy spectrum, representing the typical non-thermal electron spectra resulting in most known astrophysical acceleration mechanisms.", '2.2. Blazar Parameters': 'Table 1 shows the parameters for two synthetic blazar types used in this work, similar to the ones used in Thiersen et al. (2022): a flat spectrum radio quasar (FSRQ) and a high synchrotron-peaked BL Lac object (HBL). The parameters for these blazars were chosen to be as closely matched as possible while maintaining the characteristics of their classification. This is done to eliminate unintentional side effects due to differences in blazar parameters as far as possible. Although observed HBLs generally exhibit redshifts z ≲ 0 . 5, we use the same redshift for both source types in order to allow for a direct comparison. Placing the simulated HBL at z = 1 . 0 implies decreased fluxes and a slight shift of the SED towards lower frquencies, which do not meaningfully impact the variability results.  \nThe differences in parameters are mainly the different lepton spectra of the two blazars. Compared to the leptons in the FSRQ model, the leptons in HBL models possess higher Lorentz factors to produce the necessary significant SSC emission in X-rays up to VHE γ -rays. However, the synchrotron emission in radio up to soft X-rays also produced by HBL-leptons would be overestimated with magnetic fields strengths similar to the FSRQ case, which therefore warrants the lower magnetic field strength for the HBL. The higher Lorentz factors of the leptons in the HBL case increase the peak frequency of the synchrotron emission such that the spectrum can be classified as HBL.  \nFurthermore, the differences in electron injection index and co-moving injection luminosity of the electron spectrum are adjusted to bal-  \nance the spectral energy distributions (SEDs) to be realistic. Figure 1 shows the steady-state SEDs of the two blazars, illustrating also which cooling mechanisms dominate in each case.  \nWe derive light curves in the optical ( R band), X-ray (0.2 - 10 keV), HE (100 MeV - 300 GeV), and VHE (20 GeV - 300 TeV) γ -ray domains as indicated in Fig. 1. These domains are chosen to be representative of current and future instruments such as Swift -XRT, Fermi -LAT and CTAO.', '2.3. Stochastic Variations': "We assume that a self-sustaining process exists that can change a parameter or parameters in the emission region following our generated stochastic variations. For simplicity and to isolate the signatures of variations of individual parameters, the stochastic variations are limited to vary only a single parameter in each realization of the model.  \nFour different parameters were used as varying parameters: the maximum electron Lorentz factor, γ max , the electron injection luminosity, L inj , the magnetic field strength, B , and the electron spectral index, q . All these parameters are plausible candidates to change in different acceleration processes (Summerlin & Baring 2012; Mo'scibrodzka & Falcke 2013; Zech & Lemoine 2021; Zhang et al. 2022).  \nThe algorithm for producing synthetic light curves by Timmer & Koenig (1995) is employed here to produce the variations. The algorithm takes a power spectral density (PSD) as input to produce a periodogram - the Fourier transform of a time series - that maps to a noisy version of the input PSD. This noise is generated by means of mapping normally distributed random variables that depend on the input PSD. Applying an inverse Fourier transform on the periodogram produces a time series which represents the synthetic light curve.  \nTable 1. Initial model parameters for the representative FSRQ and HBL blazar cases. These parameters remain fixed unless they represent the varying parameter of the respective simulation realization. The subscript 0 indicates the parameter's initial value.  \nFigure 1. Steady-state spectral energy distributions (SEDs) for the FSRQ ( left ) and HBL ( right ). The gray bands depict the energy ranges over which the flux is integrated to generate the optical, X-ray, and HE and VHE γ -ray light curves.  \n<!-- image -->  \nThe time series is manipulated appropriately to represent a scale factor or an offset of the varying parameters and to fall within physical limits both numerically and physically imposed. The respective parameter scale factors span the following ranges: γ max /γ max,0 ∈ [0 . 1 , 4], L inj /L inj,0 ∈ [0 , 4], B/B 0 ∈ [10 -2 , 3], q ∈ [ q 0 -1 , q 0 +1]. 1 These manipulations are such that the underlying indices in the PSDs are not affected. Five different values of the PSD index α of the variations were explored: α = 1 . 0 , 1 . 5 , 2 . 0 , 2 . 5, and 3 . 0.", '2.4. Variation Permutations': "Variability is generated for each permutation of blazar parameters, variation PSD indices and varying parameters - 40 permutations in total. For each permutation a total of 100 simulations were run to obtain good statistics. Furthermore, in order to test the robustness of our results and ensure that there are no artefacts due to the finite time resolution of the simulations, we run every simulation twice, once with a short time step and once with a long time step, with the respective time steps differing by a factor of ∼ 20. In addition to testing the statistical robustness of our simulations, this also allows us to explore an extended range of temporal frequencies, referred to as short- and long-timescale frequencies for the short- and long-timesteps, respectively.  \nThe short-time-step case used a co-moving time step of ∆ t = 2 hours which translates to ∆ t obs = 720 seconds in the observer's frame. Cooling time scales for the radiation mechanisms in the observer's frame given the blazar parameters in Table 1 range between ∼ 10 2 seconds up to ∼ 10 5 seconds depending on the energy of the leptons (Thiersen et al. 2022). Due to Doppler and redshift effects this time step duration translates to the effective temporal fre-  \nin the observer's frame falling within the range [7 . 18 × 10 -7 , 2 . 9 × 10 -3 ] Hz.  \nSimilarly, in the long-time-step case ∆ t obs = 1 . 5 × 10 4 s (∆ t = 1 . 5 × 10 5 seconds in the co-moving frame) which is much longer than most of the cooling time scales of the emission processes involved. This case probes temporal frequencies in the range [3 . 48 × 10 -9 , 1 . 41 × 10 -5 ] Hz in the observer's frame. The overlap, [7 . 18 × 10 -7 , 1 . 41 × 10 -5 ] Hz, allows for testing the consistency of the model for different time step sizes.  \nLight curves were generated for a total of 8000 simulations. For each of the permutations a waveband-specific average fractional variability (Vaughan et al. 2003; Schleicher et al. 2019), and PSD power law index is calculated from the light curves. The average fractional variability for a permutation is calculated from the fractional variability of each of the light curves of its 100 simulations. The standard deviation of the fractional variability data is used as error bar. The waveband-specific average PSD power law index of a permutation is derived by first calculating an average PSD. This is done by calculating the average and standard deviation for the power in each temporal frequency over the 100 simulations of a permutation. A power law is then fitted to the average PSD to determine the average PSD index.", '3. RESULTS': 'Figures 2 and 3 show the averaged PSDs of simulated variability for each permutation of blazar parameters and varying parameter with input PSD power law index of 2 . 0. The variability indices within the overlapping frequency range are compatible within error margins. The power differences between the short- and longtime-scale simulations are due to keeping the magnitude of input parameter variations the same for both time step cases. This leads to a proportional shift in the power of each temporal frequency. Most notably, the variability indices  \nFigure 2. Variability PSDs for the FSRQ case with variation power law index 2 . 0. Blue : long-time-scale frequencies. Red : short-time-scale frequencies. Varying parameters as indicated in the titles and wavelengths as indicated in the panels.  \n<!-- image -->  \nfor all cases are larger than the variation index of 2 . 0.  \nWe do not show PSDs for the other variation indices, as they look rather similar to Figs. 2 and 3. The main results of these plots are the variability indices in the overlap region of the shortand long-time-scale simulations. The comparison is shown in Figs. 4 and 5 for all parameter variations and variation indices. The results show that the short- and long-time-scale indices are mostly consistent with each other. With few exceptions, the comparison between the variation and variability PSD indices shows no strong evidence of proportionality as most results are  \ncompatible with a constant index of the variability within error margins of 2 σ .  \nThe exceptions from the aforementioned compatibility with constant variability index are for the FSRQ the X-ray variability indices for varying (a.1) magnetic field and (a.2) electron spectral index; and for the HBL the (b.1) X-ray and (b.2) VHE γ -ray variability for varying maximum Lorentz factor, (c) optical variability for magnetic field variations and (d) HE γ -rays for electron spectral index variations.  \nIn some cases the variability produced in different wavebands in the model is of so small amplitude that it can be considered practically unobservable. Figures 6 and 7 show the aver-  \nFigure 3. Same as Fig. 2, but for the HBL case.  \n<!-- image -->  \naged fractional variability, F var , as a function of frequency for the different permutation cases. The fractional variability of the simulated light curves is low compared to typical observations, in some cases with values as low as 10 -3 . This is likely due to relatively low amplitudes of the input parameter variations, but is not expected to have a significant impact on the PSD characteristics of the light curves. As expected, the fractional variability increases with steeper variation temporal indices, as this leads to increased power at low frequencies of the underlying parameter variations.', '4. DISCUSSION AND CONCLUSION': 'This work shows that a time-dependent onezone leptonic blazar model with stochastic pa-  \nrameter variations produces variability patterns that are almost independent of the PSD powerlaw index of the input parameter variations and of the time step used for the simulations. However, it is unphysical for such a power law trend to extend indefinitely to lower temporal frequencies. Therefore, some cut-off is expected, but no clear indications of such a cut-off were found in our results.  \nThe consistency of a constant temporal spectral index for the variability suggests that it is unlikely that the PSDs of blazar variability are dominated by the temporal characteristics of the underlying physical conditions in blazars, specifically the changes in the injected electron spectrum (i.e. maximum Lorentz factor, injection luminosity, and electron spectral in-  \nFigure 4. PSD index comparisons for the FSRQ simulations. Blue : long-time-scale frequencies; Red : short-time-scale frequencies. Varying parameters as indicated in the titles and wavebands as indicated in each panel.  \n<!-- image -->  \ndex) and magnetic field variations tested in this work. As a caveat, we note that the magneticfield variations simulated in our work do not reflect the drastic changes expected in case of magnetic reconnection.  \nThis leads to the conclusion that the variability power law index is not intrinsically determined by the temporal spectrum of underlying single parameter fluctuations in the emission region. However, the temporal variations of physical parameters are still expected to be a plausible cause for small amplitude flux variations which continuously occur even in quiescent states in blazars.  \nAfurther caveat is that the fractional variability produced by our simulations is significantly smaller than typically observed in blazar light curves. While this could be remedied by largeramplitude parameter variations in our simulations, we do not expect this difference to have a significant impact on the PSD characteristics of the simulated light curves.  \nIn our previous work (Thiersen et al. 2022) we studied only a single variation index (2 . 0) and concluded that the resulting PSD index was similar to the variation index. However, a fit of the PSD was not done, and thus the results in both studies are compatible. However, we concluded in Thiersen et al. (2022) that there is  \nFigure 5. Same as Fig. 4, but for the HBL case.  \n<!-- image -->  \na high likelihood that the variability PSD reflected the PSD of the varying parameter. This is refuted in this work given that the variability index as a function of the variation index is compatible with a constant value throughout the various cases. As in our previous work, we find that the model does not produce wavebanddependent variability spectral indices as found in observations. This is an indication that in reality blazar emission does not originate from a single emission zone; see, e.g., Aharonian et al. (2023) for an example of observational evidence for this. Therefore, a multi-zone emission model with emission zones quasi-causally connected would be better suited for simulation of variability that is more representative of observed variability patterns. Furthermore, one can also ar-  \ngue that single parameter variations are highly unlikely and not physically plausible which supports the hypothesis of using multi-parameter variations. This is left for future investigations.  \nAsimilar study was conducted by Polkas et al. (2021), but they employed observed light curves from Fermi-LAT in order to generate single parameter variations. Similar to our conclusions, they found that single parameter variations are unable to reproduce the observed variability.  \nFigure 6. Fractional variability for the FSRQ case with colours indicating the short-/long-time-step simulations and PSDs fit by power-laws with indices provided in the legend.  \n<!-- image -->  \nThe authors acknowledge the following: the Centre for High Performance Computing (CHPC) in Cape Town, South Africa, for providing access to computational resources and technical support that facilitated the simulations presented in this paper. MZ acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) project number 460248186 (PUNCH4NFDI).', 'REFERENCES': 'Abdollahi, S., Acero, F., Ackermann, M., et al.  \n2020, ApJS, 247, 33,  \ndoi: 10.3847/1538-4365/ab6bcb  \nAckermann, M., Anantua, R., Asano, K., et al.  \n2016, ApJL, 824, L20,  \ndoi: 10.3847/2041-8205/824/2/L20  \nFigure 7. Same as Fig. 6, but for the HBL case.  \n<!-- image -->  \nAdams, C. B., Batshoun, J., Benbow, W., et al. 2022, ApJ, 924, 95,  \ndoi: 10.3847/1538-4357/ac32bd  \nAharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, ApJL, 664, L71,  \ndoi: 10.1086/520635  \nAharonian, F., Benkhali, F. A., Aschersleben, J., et al. 2023, ApJL, 952, L38,  \ndoi: 10.3847/2041-8213/ace3c0  \nAlbert, J., Aliu, E., Anderhub, H., et al. 2007, ApJ, 669, 862, doi: 10.1086/521382  \nArlen, T., Aune, T., Beilicke, M., et al. 2013, ApJ, 762, 92, doi: 10.1088/0004-637X/762/2/92  \nBegelman, M. C., & Sikora, M. 1987, ApJ, 322, 650, doi: 10.1086/165760  \nBhatta, G., & Dhital, N. 2020, ApJ, 891, 120, doi: 10.3847/1538-4357/ab7455  \nBglyph[suppress]la ˙zejowski, M., Sikora, M., Moderski, R., & Madejski, G. M. 2000, ApJ, 545, 107, doi: 10.1086/317791  \nB¨ottcher, M., Mause, H., & Schlickeiser, R. 1997, A&A, 324, 395, doi: 10.48550/arXiv.astro-ph/9604003 B¨ottcher, M., Reimer, A., Sweeney, K., & Prakash, A. 2013, ApJ, 768, 54, doi: 10.1088/0004-637X/768/1/54 Brill, A. 2022, ApJ, 936, 147, doi: 10.3847/1538-4357/ac887a Chatterjee, R., Das, S., Khasnovis, A., et al. 2021, Journal of Astrophysics and Astronomy, 42, 80, doi: 10.1007/s12036-021-09709-3 Dermer, C. D., & Schlickeiser, R. 1993, ApJ, 416, 458, doi: 10.1086/173251 Diltz, C., & B¨ottcher, M. 2014, Journal of High Energy Astrophysics, 1, 63, doi: 10.1016/j.jheap.2014.04.001 Finke, J. D., & Becker, P. A. 2015, ApJ, 809, 85,  \ndoi: 10.1088/0004-637X/809/1/85  \nGhisellini, G., & Madau, P. 1996, MNRAS, 280, 67, doi: 10.1093/mnras/280.1.67'}