On the Origin of Gamma-ray Flares from Bright Fermi Blazars

The origin of gamma-ray flares observed from blazars is one of the major mysteries in jet physics. We have attempted to address this problem following a novel spectral energy distribution (SED) fitting technique that explored the flaring patterns identified in the broadband SEDs of two gamma-ray bright blazars 3C 279 (z=0.54) and 3C 454.3 (z=0.86), using near-simultaneous radio-to-gamma-ray observations. For both sources, the gamma-ray flux strongly correlates with the separation of the SED peaks and the Compton dominance. We propose that spectral hardening of the radiating electron population and/or enhancement of the Doppler factor can naturally explain these observations. In both cases, magnetic reconnection may play a pivotal role in powering the luminous gamma-ray flares.


INTRODUCTION
The γ-ray emission from flat-spectrum radio quasars (FS-RQs) offers a powerful diagnostic tool to probe the innermost regions of jets and their surroundings. FSRQs are a special class of active galactic nuclei (AGN) having relativistic jets oriented close to the line of sight to the observer and characterized by broad optical emission lines. The typical spectral energy distribution (SED) of an FSRQ is dominated by the non-thermal jet emission and is characterized by two broad bumps. In the widely accepted leptonic emission scenario, the low-energy bump produced by the synchrotron process peaks in the sub-millimeter (mm) or the infrared (IR)-optical energy regime. The high energy bump, on the other hand, produced by inverse Compton scattering, peaks in the MeV range, leading to a steep γ-ray spectrum in the GeV energy band (e.g., Ajello et al. 2020). These intriguing objects exhibit flux variations at all accessible wavelengths with the most extreme variability, on the timescale of hours or less, observed at GeV energies (see, e.g., Foschini et al. 2011). During high-activity states, many FSRQs show a rising γ-ray spectrum in νF ν , indicating a shift of the inverse Compton peak from MeV to GeV energies (cf. Paliya et al. 2016). In a few cases, the associated synchrotron peak is also observed to shift from far-IR to optical-UV energies (e.g., Angioni et al. 2019). On the other hand, there are observations indicating that only the inverse Compton peak shifted to higher frequencies without any significant change in the synchrotron peak (Paliya et al. 2016). Such contrasting behaviors indicate the complex interplay between particle acceleration and radiation mechanisms in relativistic jets and, at the same time, provide tantalizing clues to explore the underlying causes responsible for the flux variability. The key to understanding the origin of γ-ray flares lies in unraveling the spectral patterns shown by FSRQs across the electromagnetic spectrum.
We have attempted to address this outstanding problem by carrying out a systematic study of the broadband SEDs of two well-known γ-ray detected blazars, namely 3C 279 (z = 0.54) and 3C 454.3 (z = 0.86). The motivation is to understand the evolution of the low-and high-energy SED peaks as a function of the γ-ray activity following a novel SED fitting approach which treated the SED peak transitions separately. We further developed a theoretical model to decipher the observed patterns revealing the origin of γ-ray flares. We selected 3C 279 and 3C 454.3 because these are among the brightest γ-ray blazars with hundreds of existing multi-  Figure 1. The 3-day binned γ-ray light curves of 3C 279 (left) and 3C 454.3 (right). Vertical dotted lines refer to the epochs when there is at least one simultaneous observation available each in the radio, NIR-to-UV, X-and γ-ray bands.
wavelength observations supplemented with more-than-adecade of continuous monitoring with the Fermi-Large Area Telescope (Fermi-LAT). The dense sampling has allowed us to study their multi-epoch SED evolution using contemporaneous data covering all accessible electromagnetic bands. Moreover, near-IR (NIR) to optical-ultraviolet (UV) spectra of both objects are strongly dominated by non-thermal synchrotron emission from the jet, with signatures of direct accretion disk emission only emerging in very low activity states (e.g., Pian et al. 1999) and at frequencies far beyond the synchrotron peak. Thus, we are able to reliably estimate the location of the synchrotron peak based on the observed NIR -UV spectra. The following cosmology parameters were adopted: H 0 = 67.8 km s −1 Mpc −1 , Ω m = 0.308, and Ω Λ = 0.692 (Planck Collaboration et al. 2016) 2. DATA REDUCTION AND ANALYSIS Gamma-rays: We analyzed 0.1−300 GeV data taken with the Fermi-LAT between 2008 August 4 and 2021 May 5 (MJD 54683-59339) to generate 3-day binned γ-ray light curves and spectra. We followed the standard data reduction procedure 1 fully described in various works (cf. Paliya et al. 2019) and we briefly elaborate it here. A 15 • region of interest (ROI) was chosen to select SOURCE class events which were subjected to the relational filter 'DATA QUAL>0 and LAT CONFIG==1' to determine good time intervals. Additionally, we adopted a zenith angle cut of z max < 90 • to avoid contamination from Earth's albedo. A binned likelihood fitting technique was used to determine the spectral parameters of sources of interest. For this purpose, we considered all γ-ray sources located within 20 • from the center of the ROI and included in the fourth catalog of the Fermi-LAT detected sources (4FGL-DR2; Abdollahi et al. 2020;Ballet et al. 2020). The standard Galactic and isotropic background templates 2 were also used. The spectral parameters of all sources lying within the ROI were kept free in the likelihood fit, whereas, those lying outside ROI were fixed to their 4FGL-DR2 values. The significance of the γ-ray detec-1 http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/ 2 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html tion was quantified using maximum likelihood test statistic TS=2∆ log(L) where L is the ratio of the likelihood values for models with and without a γ-ray point object. We first performed a likelihood fit for the whole duration covered in this work and then froze the spectra parameters of all sources whose TS values were found to be <25. This updated sky model was then used to generate 3-day binned light curve and spectra. In this study, we considered only those epochs during which the target quasars were significantly detected ( 5σ) in the γ-ray band. The generated light curves are shown in Figure 1.
X-rays: We used the "Swift-XRT data products generator 3 " (Evans et al. 2009) to extract 0.3−10 keV spectra from the observations taken with the Neil Gehrels Swift X-ray telescope in the time interval MJD 54683 -59339. Based on the current count rate, this tool automatically takes into account the possible pile-up effects and choose the source and background regions of appropriate sizes. An absorbed power-law model was adopted to derive the source flux for each observation using the Galactic neutral hydrogen column density from Kalberla et al. (2005). The spectral fitting was carried out in XSPEC (Arnaud 1996).
NIR-UV: Data from the Swift Ultraviolet and Optical telescope (UVOT -V, B, U, W 1, M 2 and W 2 filters), Steward observatory (V and R bands; Smith et al. 2009) and Small and Medium Aperture Research Telescope System (SMARTS -B, V, R, J and K filters; Bonning et al. 2012) were used to extract the source magnitudes in the NIR, optical, and UV bands. They were corrected for the Galactic extinction following Schlafly & Finkbeiner (2011) and converted to flux units using the standard zero points (Bessell et al. 1998;Breeveld et al. 2011).
Radio: We considered the multi-frequency radio fluxes made available by the University of Michigan Radio Astronomy Observatory (UMRAO -4.8, 8, and 14.5 GHz;Aller et al. 1985), the Submillimeter Array 4 (230 and 350 GHz; Gurwell et al. 2007), and the Atacama Large Milime-  In the appendix, we provide the output of the multiwavelength data reduction in Table 1 and 2 for 3C 279 and 3C 454.3, respectively.

EVOLUTION OF THE SED PEAKS
Our primary objective is to study the evolution of the SED peaks as a function of the γ-ray activity of 3C 279 and 3C 454.3. Therefore, from the γ-ray light curves, only those epochs were chosen that have at least one observation each in the radio and NIR-UV (covering the synchrotron peak), and X-ray and γ-ray (constraining the inverse Compton peak) bands. This criterion of selecting near-simultaneous data across the electromagnetic spectrum has led to the identification of 155 and 149 epochs for 3C 279 and 3C 454.3, respectively ( Figure 1). To determine the peak frequencies and luminosities, we fitted a third-order polynomial following the Markov Chain Monte Carlo technique implemented in the software emcee (Foreman-Mackey et al. 2013). Both SED peaks were fitted independently and we derived the uncertainties based on the 16th, 50th, and 84th percentiles of the sample in the marginalized distribution. We show the results of this fitting procedure in Figure 2. Figure 3 shows the variations of the synchrotron peak frequency (ν syn pk ), inverse Compton peak frequency (ν IC pk ), and their separation measured during various γ-ray flux states, normalized to those derived from their average activity state SEDs by Paliya et al. (2021). The strength of the correlation of the plotted quantities was determined by computing Spearman's rank correlation coefficient (ρ) taking uncertainties into account. Both ν syn pk and ν IC pk have been found to correlate positively with the γ-ray activity for 3C 279. For 3C 454.3, ν IC pk follows the same trend though ν syn pk does not show any significant correlation with the γ-ray flux. Interestingly, the separation of both SED peaks strongly correlates with the γ-ray flux for both sources. This can be understood by comparing the ν syn pk and ν IC pk panels in Figure 3 which reveals that the extent of variation is larger for the latter. In other words, ν IC pk shifts considerably more than ν syn pk during γ-ray flares, thus leading to an overall increase in the SED peak separation.
We have also derived the ratio of the inverse Compton to synchrotron peak luminosities, also known as Compton dominance, and show it as a function of the γ-ray flux in Figure 4. In both cases, a strong correlation is found.
One possible caveat of this work could be the fact that we have considered near-simultaneous radio data to constrain the synchrotron peak. One can argue that the observed radio emission may originate from more extended regions, further down the jet, than the compact one radiating in NIR-toγ-ray bands due to the synchrotron self-absorption process. Though our analysis method does not rely on any such assumptions, we performed the following test to reveal its impact on the results. Since the radio emission produced by the plasma blob radiating NIR-to-γ-rays is expected to be ≤ the observed flux, we recomputed ν syn pk and determine the strength of correlations described above by dividing the observed radio flux by 10. No assumption was made on the spectral shape of the self-absorbed radio flux. In all cases, the results followed the same trends and the derived correlations remain significant. The changes were noticed only in the derived ratios plotted in Figure 3 and 4. For example, decreasing the radio flux reduced the corresponding synchrotron peak luminosity thereby increasing the Compton dominance which still showed a strong correlation with the γ-ray flux, ρ = 0.58 ± 0.01 and 0.43 ± 0.01 for 3C 279 and 3C 454.3, respectively.

PHYSICAL INTERPRETATION
The multiwavelength variability and correlation studies of both 3C 279 and 3C 454.3 have been extensively studied, especially during major flares (e.g., Giuliani et al. 2009;Vercellone et al. 2009;Raiteri et al. 2011;Hayashida et al. 2012;Jorstad et al. 2013;Paliya et al. 2015;Bottacini et al. 2016;Larionov et al. 2020;Zhou et al. 2021).   In most works discussing the leptonic modeling of blazar SEDs, the γ-ray emission from FSRQs is dominated by inverse-Compton scattering of a target radiation field external to the jet, such as line-dominated optical -UV radiation from the Broad Line Region or infrared emission from a dusty torus from a relativistic non-thermal electron population assumed to be isotropic in the co-moving frame of an emission region moving with bulk Lorentz factor Γ along the jet (e.g., Ghisellini et al. 2010;Böttcher et al. 2013). This results in relativistic Doppler boosting of the radiation with respect to the co-moving frame of the emission region by a factor of δ = (Γ [1 − β Γ cos θ obs ]) −1 in frequency and a fac-tor of δ 4 in νF ν peak flux. In such models, therefore, we find that the synchrotron peak is determined by the source parameters as ν syn pk ∝ B γ 2 p δ (1) and νF syn ν,pk ∝ N e (γ p ) γ 3 p B 2 δ 4 where B is the magnetic field and γ p is the peak of the electron spectrum in a γ 3 N e (γ) representation. In the case of FSRQs with an inverse-Compton peak in the MeV -GeV regime, inverse-Compton scattering to photon energies around the peak is likely dominated by Compton scattering in the Thomson regime. Therefore, the inverse-Compton peak location is determined by where u ext and ν ext are the energy density and peak frequency of the external radiation field in the stationary frame of the AGN.
The peak separation ratio therefore depends on source parameters as and the Compton dominance as As evident from Eq. (4), γ-ray flaring may be caused by several factors, which we will discuss in the following: (1.) an increase in the number of relativistic electrons N e (γ); (2.) a hardening of the radiating electron distribution, resulting in an increase in γ p ; (3.) an increase in the energy density of the external radiation field, u ext , and/or (4.) an increase in the bulk Lorentz factor Γ and Doppler factor δ.

1.
A pure increase of the electron number density (without significant spectral changes or other parameter changes) will lead to an equal increase of both radiation components (synchrotron and inverse-Compton). It will not shift the peak frequency of either component and will not alter the Compton dominance. Hence, such a scenario is ruled out as the dominant cause of γ-ray variability by the significant correlation of both the peak separation and the Compton dominance with flux in both sources.

2.
A spectral hardening of the electron distribution, as expected in most particle acceleration scenarios, such as internal shocks or magnetic reconnection, shifting γ p to higher values, without any other parameter changes, will lead an equal increase of the synchrotron and Compton peak frequencies and peak fluxes. In order to increase the peak separation, other parameter changes would need to be invoked. It has been found in several modeling studies of blazar variability (see, e.g., Zhang et al. 2015;Böttcher & Baring 2019) that a reduction in the magnetic field accompanied by more efficient particle acceleration can successfully be invoked to reproduce the observed increase of the Compton dominance and peak separation. In particular, it may explain the lack of a correlation of the synchrotron peak frequency with the γ-ray flux, as we found for 3C454.3. Such a reduction of the magnetic field accompanying γ-ray flaring activity could indicate magnetic energy conversion, as in magnetic reconnection, transferring magnetic-field energy into relativistic particle energy.
An interesting alternative explanation has recently been suggested by Sobacchi & Lyubarsky (2019) and Sobacchi et al. (2021): Those authors propose that stochastic particle acceleration mediated by magneto-hydrodynamic (MHD) turbulence will cause the highest-energy particles to maintain small pitch angles with respect to the magnetic field. This will naturally reduce both the synchrotron peak frequency and the synchrotron flux. It may result in flaring due to enhanced particle acceleration, shifting γ p to higher values, not leading to a significant increase of the synchrotron peak, while the inverse-Compton peak frequency and flux may increase significantly.
3. Blazar variability models invoking a temporary increase of the energy density of the external target photon field for inverse-Compton scattering have been proposed in the form of a synchrotron mirror model, in which part of the jet synchrotron emission is reflected by a cloud near the jet trajectory (Böttcher & Dermer 1999;Tavani et al. 2015) or by a stationary jet feature, as in the "ring of fire model" (MacDonald et al. 2017).
If the particle acceleration characteristics in the emission region are not changing in such a scenario, as the emission region passes the mirror, the enhanced radiative cooling is naturally expected to decrease γ p , thus, in fact, leading to a decrease of both the synchrotron and inverse-Compton peak frequencies as well as a reduction of the synchrotron peak flux. If these were the only changes in such a scenario, one would expect to see the observed increase in Compton dominance correlated with increasing flux, but not the increasing peak separation and the positive correlation of the peak frequencies with flux. In order to reproduce those latter correlations, a reduction of the magnetic field, along with increased particle acceleration efficiency (increasing γ p ) would have to occur along with the passing of the mirror. While additional particle injection and acceleration may be the result of a shock caused by the mechanical interaction of the emission region with the obstacle (see, e.g., Barkov et al. 2012;Zacharias et al. 2017), one would expect a compres-sion and, thus, increase of the magnetic field to result as well, leading to a decrease of the peak separation.
Thus, it seems difficult to imagine a plausible scenario of an increasing external radiation field to reproduce the flux correlations found in our study.

4.
A change in the bulk Lorentz factor and Doppler factor will lead to enhanced Doppler boosting of both the peak frequencies and fluxes of both radiation components. If such a change were effected by a change in the viewing angle, not changing Γ, both radiation components would be affected in the same way, in contrast to the observed peak separation vs. flux and Compton dominance vs. flux correlations. Thus, more plausible would be an increase of the Doppler factor, which, for initial viewing angles θ obs < 1/Γ, will also lead to an increase of δ. Positive correlations of all quantities (both peak frequencies, peak separation, and Compton dominance) are expected in this case, as observed for 3C279 (see also Larionov et al. 2020) In order to explain the lack of a correlation between ν sy pk and γ-ray flux for 3C 454.3, a reduction of the magnetic field, as in scenarios 1. and 2., would need to be invoked.

SUMMARY
In this work, we have studied the evolution of the SED peaks as a function of the γ-ray activity of two bright blazars with good multi-wavelength coverage. The separation of the SED peaks and the Compton dominance are found to strongly correlate with the γ-ray flux activity for both sources, indicating a similar mechanism to be responsible for the observed γ-ray flares. We have interpreted these findings based on the spectral hardening of the radiating electron population and/or enhancement of the Doppler factor. In both cases, magnetic reconnection may play a major role in powering the γ-ray flares. Similar results have been reported in various studies focused on the major flaring events detected from 3C 279 and 3C 454.3 (cf. Bonnoli et al. 2011;Shah et al. 2019). Furthermore, these objects have occasionally shown 'unusual' SEDs during flaring episodes that might not have been included in our analysis due to lack of nearsimultaneous data e.g., at radio or X-ray wavelengths (e.g., 2013 December flare of 3C 279; Paliya et al. 2016). Such pe-culiar flaring events are crucial in understanding the origin of γ-ray flares and associated diverse radiative mechanisms by applying various physically-motivated models rather than the polynomial function used in this study. Therefore, the results presented here should be considered as representing the overall broadband spectral behavior of bright Fermi blazars as a function of their γ-ray activities. The proposed hypotheses can be tested for a larger sample of blazars with the accumulation of contemporaneous radio-to-γ-ray data made available by the current-and next-generation of telescopes.

ACKNOWLEDGMENTS
We are grateful to the journal referee for constructive feedback. The work of M. Böttcher is supported by the South African Research Chairs Initiative (grant no. 64789) of the Department of Science and Innovation and the National Research Foundation 5 of South Africa. V.S.P. is grateful to H. Aller and C. Raiteri for providing radio data taken with UMRAO and published in Larionov et al. (2020), respectively. This work is partly based on data taken and assembled by the WEBT collaboration and stored in the WEBT archive at the Osservatorio Astrofisico di Torino -INAF (http://www.oato.inaf.it/blazars/webt/). This paper has made use of up-to-date SMARTS optical/near-infrared light curves that are available at www.astro.yale.edu/smarts/glast/home. php. Data from the Steward Observatory spectropolarimetric monitoring project were used. This program is supported by Fermi Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX15AU81G. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center. This research has made use of NASA's Astrophysics Data System Bibliographic Services.

APPENDIX
The multiwavelength data used in this paper is provided in Table 1 and 2 for 3C 279 and 3C 454.3, respectively. NOTE-All information are same as in Table 1.
(This table is available in its entirety in a machine-readable form in the online journal. A portion is shown here for guidance regarding its form and content.)