Variability and Spectral Behavior of Gamma-Ray Flares of 3C 279

LAT scans the whole sky every three hours in the energy range from 20 MeV to $>$ 300 GeV (Atwood et al., 2009). We selected LAT data from the Fermi Pass 8 database in the time period from 2008 August 4 15:43:36 (UTC) to 2020 Nov 25 23:41:00 (UTC), with energy range in 0.1–300 GeV. Following the recommendations of the LAT team¹¹1https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/, we selected events with zenith angles less than 90 deg to prevent possible contamination from the Earth’s limb. The LAT science tool Fermitools 2.0.8 and instrument response function (IRF) P8R3_SOURCE_V2 were used. For the target 3C 279, a 20${}^{\circ}\times 20^{\circ}$ square region of interest (ROI) centered at its position given in 4FGL-DR2 (Abdollahi et al., 2020b) was selected. The normalization parameters and spectral indices of the sources within 5 deg from 3C 279, as well as sources within the ROI with variable index $\geq$ 72.44 (Acero et al., 2015), were set as free parameters. All other parameters were fixed at their catalog values in 4FGL-DR2. We used the original spectral models in 4FGL-DR2 for the sources in the source model. Gamma-ray spectra and lightcurves binned by hour and minute timescales in subsection 2.3 were derived by an unbinned maximum likelihood analysis with gtlike. Galactic and extragalactic diffuse emission models were added to the source model using the spectral model file gll_iem_v07.fits and iso_P8R3_SOURCE_V2_v1.txt²²2https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html, respectively. The normalizations of the two diffuse emisson components were set as free parameters in the analysis.

We made the long-term light curve of the whole analysis time period above 100 MeV and found three obvious flares in 2018 in addition to the flare in 2015 June. Then we used the parameters in subsection 2.1 and a simple power law for 3C 279 in the source model. We constructed light curves in the daily time bin from MJD 58091 to MJD 58373 by performing a standard binned maximum likelihood analysis. The daily binned light curve is shown in Figure  1, in which when flux data points have the maximum likelihood Test Statistic (TS) values larger than 9 were plotted. The observed gap in Figure  1 is due to a technical issue encountered on 2018 March 16 by Fermi spacecraft. LAT did not carry out any observations between March 16, 2018 and April 08, 2018.

In order to investigate the fast flux variability in detail, we constructed light curves in hour timescale bins (Figure 2-(a)), in which when flux data points have TS values larger than 9 were plotted. We noted that the flare F1 and F3 did not show noteworthy fast component or profile, therefore, we only made minute timescale binned light curve for F2 (Figure 2-(b)). The measured $\gamma$-ray flux in 1-hour bins reached $\sim 4.5\times 10^{-5}\;$${\rm photons}\;{\rm cm}^{-2}\;{\rm s}^{-1}$ on April 19, 2018. Here we defined the flare period of F2 to be between MJD 58227.875 and MJD 58228.500 (April 19, 2018 21:00:00 and April 20, 2018 12:00:00), as indicated in Figure 2-(a). And it concludes 10 bins designated Stage “A” through “J”, respectively. We fitted a constant value to each stage in F2 for both time bins, and calculated a probability ($p$-value) from $\chi^{2}$ for each stage. While stage D, H, I resulted in $p$-values consistent with constant fluxes, we found significant indications of rapid variability on a minute timescale in the 2-min binned light curve. For stage B: ($p$, $\chi^{2}$/dof)$=$($1.789\times 10^{-15}$, 87.11/8), stage C: ($p$, $\chi^{2}$/dof)$=$($3.370\times 10^{-05}$, 47.24/15), stage E: ($p$, $\chi^{2}$/dof)$=$($1.747\times 10^{-05}$, 31.85/6), stage F: ($p$, $\chi^{2}$/dof)$=$($2.145\times 10^{-08}$, 51.45/8), stage G: ($p$, $\chi^{2}$/dof)$=$(0.013, 20.97/9).

We constructed light curves in 2-min, 3-min, and 5-min bins for the whole flare period of F2, and eventually found minute-scale flare profile only in stage B, which was shown in Figure 2-(b). Then we used an exponential equation to fit the 2-min binned light curve profile, which is given by

where $F_{\rm c}$ and $F_{0}$ are the constant flux and height of a peak, respectively, $t_{0}$ is the flux peak time, $T_{r}$ and $T_{d}$ are used to measure the rise and decay time separately. The best fitting result gives $\chi^{2}$/dof$=$1.71/5 with the flare peak time $t_{0}$ on 26.16 (2018 Apr 19 22:16:10(UTC)). Flare peak flux with the best fitting at $t_{0}$ is $F_{\rm p}=4.32\times 10^{-5}\;$${\rm photons}\;{\rm cm}^{-2}\;{\rm s}^{-1}$ . We calculate the average flux value of the four data points before flare as the flux in quiescence ($F_{\rm q}=1.90\times 10^{-5\;}$ ${\rm photons}\;{\rm cm}^{-2}\;{\rm s}^{-1}$ ), which is shown in Figure 2-(b). The time when the flux equals to $F_{\rm q}$ before and after the flare is 15.66 and 50.99, respectively. Therefore the flux decay time (from $F_{\rm p}$ to $F_{\rm q}$) is 24.83 min. We calculate the upper limit of $F_{\rm c}$ from the best fitting result of the exponential equation and take it as $F_{\rm q}$, which equals to $2.25\times 10^{-5}\;$ ${\rm photons}\;{\rm cm}^{-2}\;{\rm s}^{-1}$ . Then we can get that the lower limit of flux decay time is 15.04 min.

We fit the $\gamma$-ray spectra of 3C 279 described with a log-parabola model ($dN/dE\propto(E/E_{0})^{-\alpha-\beta\log(E/E_{0})}$, where $\alpha$ and $\beta$ are spectral parameters, here $E_{0}$ = 466 MeV). We plotted the spectral data points with TS$>$4 and upper limits when TS$<$4 in Figure 3. Figure 3-(a) shows the spectra in the ‘pre-outburst’ and ‘post-flare’ periods as defined in Figure 2 for comparison. Figure 3-(b) shows the $\gamma$-ray spectra as measured by Fermi-LAT for each stage. Stages B and C, and stages D, E, and F, were combined because they showed similar fluxes and spectral fitting results. The spectral energy flux peaks in each stage correspond with the light curve during the flare F2. The flux and spectral fitting results above 100 MeV for Each stage (A−-J) during flare F2 are shown in Table 1.

The $\gamma$-ray flux (0.1–300 GeV) vs. the $\gamma$-ray power law photon index $\Gamma$ of 3C 279 with daily time bin based on the LAT data during the flare F2 is given in Figure  4. A correlation between the $\gamma$-ray flux and photon index was visible to suggest a ‘harder-when-brighter (BWB)’ trend. The tendency of spectral hardening during the flaring state was seen in blazar flares by many authors (e.g., Britto et al. 2016; Shukla et al. 2018; Shah et al. 2019). A similar trend pattern was also previously seen for 3C 279 by Hayashida et al. (2015), Paliya et al. (2015b) and Prince (2020). The BWB could be explained in many ways, such as two-component (one stable component and a variable one) model (Fiorucci et al., 2004), fresh electrons injection (Kirk et al., 1998; Mastichiadis & Kirk, 2002), etc.

The question of the location of the blazar $\gamma$-ray dissipation region ($d_{\rm diss}$) is intriguing to researchers. To directly determine the $\gamma$-ray location, we need an extremely large $\gamma$-ray telescope to fulfill the required resolution at its wavelength. Some indirect methods have been proposed to solve this issue, one of the popular ways is to measure the variability timescale and calculate the size of the dissipation region ($R_{\rm diss}$) by assuming the variability timescale accounts for the light travel across the dissipation region. Timescale from hours to minutes have been reported for blazars (Abdo et al., 2011; Tanaka et al., 2011; Brown, 2013; Shukla et al., 2018) that constrain a very compact dissipation region. Then, the size of the dissipation region is used to calculate distance to the central black hole through $d_{\rm diss}=R_{\rm diss}\psi$ with an assumption that the dissipation region diameter is identical to the width of the jet at this distance, where $\psi$ is the semi-aperture opening angle of the jet (Dermer et al., 2009; Ghisellini & Tavecchio, 2009). However, the distance is usually constrained to the base of the jet, where locates on the inner side of BLR, according to this method. The BLR is full of optical/UV photons, which will cause absorption to the high energy $\gamma$-ray photons through pair production ($\gamma\gamma\to e^{\pm}$) (Poutanen & Stern, 2010) that arises a spectral cut-off at GeV energies. Acharyya et al. (2021) studied the $\gamma$-ray spectra for 18 flares and found that 11 of them with very high energy emission are incompatible with a BLR origin.

In the present work, we have studied the $\gamma$-ray spectra for the flare on April 19, 2018, see in Figure 3. Our results in Table 1 indicate that a log-parabola model is enough to fit the spectra and no significant cut-off at GeV energies is observed during the flaring state, suggesting the dissipation region is outside the BLR and the energy dissipation may be due to the IC process of scattering the DT infrared photons. This result is consistent with the one suggested by Tolamatti et al. (2022).

We calculate the broadband SED for ‘pre-outburst’, ‘flare’, and ‘post-flare’ states of F2. We employ a log-parabolic-power-law (LPPL) function, which is a combination of a power-law function in the lower energy range and a log-parabola function at a higher energy range (Tramacere et al., 2009, 2011; Tramacere, 2020), for the electron energy distribution (EED) in the dissipation blob (Massaro et al., 2006). The LPPL function is expressed as follows:

where $N$ is the normalization parameter, and $\gamma_{\rm 0}$ is the turn-over energy of the electron spectrum, $s$ is the spectral index, and $r$ is the spectral curvature.

In the framework of leptonic scenario (Ghisellini et al., 1996; Tavecchio et al., 1998), blazar broadband SED is formed through coupling the EED with the radiation mechanism (synchrotron radiation and inverse Compton scattering), and is usually displayed in the diagram of ${\rm log}\,\nu f_{\nu}-{\rm log}\,\nu$. In the present work, we introduce the Jets SED modeler and fitting Tool (JetSet) ³³3https://jetset.readthedocs.io/en/latest/index.html (Tramacere et al., 2009, 2011; Tramacere, 2020) to generate SEDs.

There are nine basic parameters to describe the SED model. Magnetic field ($B$), dissipation region size ($R_{\rm diss}$) and Doppler beaming factor ($\delta$) to describe the dissipation region environment and the jet geometry. The rest of the six parameters, the maximum energy of electrons ($\gamma_{\rm max}$), the minimum energy of electrons ($\gamma_{\rm min}$), $\gamma_{\rm 0}$, $N$, $s$ and $r$ to describe the EED. In the case of FSRQs, the entire $\gamma$-ray emission is a combination of self-Compton scattering (SSC), in which the soft photons are provided by the same population of relativistic electrons in synchrotron process, and external Compton scattering (EC), in which the soft photons can have multi-suppliers e.g., accretion disk (Dermer & Schlickeiser, 1993), the broad emission line region (BLR) (Sikora et al., 1994; Fan et al., 2006), and the dusty torus (DT) (Błażejowski et al., 2000; Arbeiter et al., 2002; Sokolov & Marscher, 2005).

We consider the EC model with soft photons origin from the DT, and there are 5 additional parameters to define the accretion disk and the DT, the disk luminosity ($L_{\rm disk}$), the disk temperature ($T_{\rm disk}$), the distance of the DT from the central engine ($R_{\rm DT}$), the temperature of the DT ($T_{\rm DT}$), and the opacity of the DT ($\tau_{\rm DT}$).

The simultaneously observed data that we used for SED calculation are collected from Tolamatti et al. (2022), in which the multi-wavelength data of duration MJD 58218–58243 is concluded. This duration spans 25 days and has a full coverage of flare F2. Thus, we collect radio, optical/UV, and X-ray data from their work to combine with our GeV $\gamma$-ray data of the three states. During the SED fitting, we have six parameters fixed as tabulated in Table 2 and the rest seven parameters to vary in Table 3, and the quantities related to jet power and SED structure for all the three states are listed in Table 4 .

The blazar, 3C 279 at the redshift 0.536, is one of the extreme variable celestial objects in multi-wavelengths. In 2018, it demonstrated three outbursts, namely F1, F2, and F3, with GeV $\gamma$-ray fluxes exceeding the previous flare in 2015. Especially, an unprecedented high flux was observed in flare F2.

We, in particular, analyzed the spectral behavior of flare F2, which pronounced ‘pre-outburst’, ‘flare’ and ‘post-flare’, seen in Figure 3. The spectral fitting results above 100 MeV are shown in Table 1. It showed clearly that a ‘harder-when-brighter’ trend was observed for F2, which suggests the detection of a high-energy photon during the high state.

During the broadband SED modeling, we fixed parameters and used those values as in Tolamatti et al. (2022). In their work, those parameters were optimized during the fitting procedure and consistent with those values in literature (Hayashida et al., 2015; Paliya et al., 2015b). In addition, we also assumed a fraction $\tau_{\rm DT}=0.15$ of the disk luminosity was intercepted and re-radiated from dust as IR emission, as a comparison, this value was used as 0.5 for Aleksić et al. (2011b) and in Paliya et al. (2015b) and was used as 0.1 in Hayashida et al. (2012, 2015).

One should be able to calculate the jet power with these given fixed parameters and those parameters obtained through SED modeling. The jet power $P_{\rm jet}$, which is a critical parameter to understand the formation and composition of jets, is usually assumed to be carried by relativistic electrons ($P_{\rm e}$), magnetic field ($P_{\rm B}$), radiation ($P_{\rm r}$) and cold proton kinetic ($P_{\rm p}$) (Celotti & Ghisellini, 2008; Tan et al., 2020) as:

where $U_{i}$ represents the energy density of the magnetic field ($i=B$), the relativistic electron ($i=e$), cold proton ($i=p$) and bolometric radiation ($i=r$) (Ghisellini & Tavecchio, 2010; Ding et al., 2019) in the comoving frame, and can be estimated in Tan et al. (2020):

where we applied a standard assumption of one cold proton every ten electrons;

where $L_{\rm obs}$ is the observed non-thermal bolometric luminosity, which is estimated by integrating over the broadband SED. The $P_{\rm e}$, $P_{\rm p}$, $P_{\rm B}$, $P_{\rm r}$ are calculated and tabulated in Table 3.

Our results suggest the magnetic field stays stable, but we find a slightly higher magnetic field for low activity states and a lower magnetic field for the high emission state, which has also been reported in previous works (Shah et al., 2019; Prince, 2020). Besides, the Doppler factor varies and the relativistic electron distribution changed significantly during the ‘pre-outburst’, ‘flare’, and ‘post-flare’ states. The Doppler factor gives the largest value for the ‘flare’ state, while giving the smallest value for the ‘post-flare’ state among the three states; This is consistent with the fact found in Shah et al. (2019), in which they studied the flare of 3C 279 in 2018 January and found that the bulk Lorentz factor increased when the flux increased. During the ‘flare’ state, the electron energy distribution is significantly changed, $\gamma_{\rm min}$, $\gamma_{\rm 0}$, and $s$ are becoming smaller, compared to those of the ‘pre-outburst’ state and the ‘post-flare’ state.

Different explanations of this flare have been purposed by previous studies. Shukla & Mannheim (2020) analyzed the light curve of F2 and found a $\sim$8 minute-timescale flare, which was superimposed on a longer duration envelope, and a peak-in-peak structure in this flare. They suggested a magnetic reconnection-derived particle acceleration to form the fast flare as well as the observed peak-in-peak light curve structure. However, Tolamatti et al. (2022) suggested that the flare should not be caused by a magnetic reconnection process, because they found the jet is particle dominated ($U_{\rm e}/U_{\rm B}\gg 1$) and the magnetization ratio ($\sigma_{\rm B}=P_{\rm B}/(P_{\rm e}+P_{\rm p})$) is too low to accelerate particles via magnetic reconnection, base on what they suggested the shock-in-jet may still responsible for the particle acceleration in the jet.

In analogy, we calculate the $U_{\rm e}/U_{\rm B}$ and $\sigma_{\rm B}$ and report them in Table 5, our results suggest that the particle energy density dominates over magnetic energy density. And the low values of $\sigma_{\rm B}\sim 0.01$ indicate that magnetic reconnection is less likely to account for particle acceleration, but it can accelerate electrons via relativistic shocks as suggested by Baring et al. (2017), in which electrons would be efficiently accelerated by relativistic shocks in blazar jets with $\sigma_{\rm B}$ ranges from $\sim 10^{-4}$ to 0.06.

Our results of Table 5 meet the conclusion from Tolamatti et al. (2022), but the contradiction between the shock acceleration in the emission blob and the minute-scale variability still exist that mentioned by Shukla & Mannheim (2020). Besides, the EED with a spectral index $\sim 2$ predicted by the non-relativistic shock acceleration and $\sim 2.2$ for the case of a classic relativistic shock (Baring et al., 1999; Kirk et al., 2000; Ellison & Double, 2004; Hu et al., 2020), but we obtained a much harder electron spectrum with $s=1.65$ for the flaring state.

Instead of both the shock acceleration and the magnetic reconnection, we consider this flare is possibly caused by an injection of higher-energy electrons from outside of the blob, these electrons are accelerated to higher energy outside and injected into the blob in a very short time. This injection of higher-energy electrons changes the distribution of electrons to a harder spectrum and the increases amount of electrons with energy above $\gamma_{\rm 0}$. These higher-energy electrons could be accelerated by a more comprehensive mechanism that can explain the unique light curve behaviour (i.e., ‘peak-in-peak’ structure), for instance, they are accelerated by the magnetic reconnection outside the blob and also obtain energy from the shock inside the blob. Otherwise, a multi-zone model should be involved to explain all these phenomenon.

In the leptonic scenario, the jet power is mostly dominated by the kinetic power of particles, including the relativistic electrons (or positrons) and cold protons. Energy dissipation of relativistic electrons is more efficient through the IC process than the synchrotron radiation resulting in a $\gamma$-ray dominant emission from the jet. The observed spectra (see in Figure 3 and Table 1) of the flare state favor a model of EC_DT, and the broadband SED has been constructed in the previous subsection. The cooling time of electrons is able to be estimated in the observer frame through

where $m_{\rm e}$ is the electron rest mass, $\gamma$ is the electron energy in unit of Lorentz factor, and the $U^{{}^{\prime}}$ is the energy density in the jet frame. The Lorentz factor $\gamma$ can be derived by the IC peak frequency: $\nu_{\rm IC}\approx\nu_{\rm ext}\gamma^{2}\Gamma\delta/(1+z)$, then

the external photon energy density $U^{{}^{\prime}}=3\times 10^{-4}\,\Gamma^{2}\ {\rm erg\cdot cm^{-3}}$ and the bulk Lorentz factor ($\Gamma$) is believed to be close to the Doppler factor ($\delta$) for blazars due to a small viewing angle (Ghisellini & Tavecchio, 2008; Tavecchio et al., 2010; Foschini et al., 2011). We calculate the cooling time scale for our EC_DT model and get $t^{\rm ob}_{\rm cool}=30.1\,{\rm min},\,19.6\,{\rm min},\,31.8\,{\rm min}$ for the three states, respectively. An average cooling time scale $t^{\rm ob}_{\rm cool}=27.2\pm 6.6\,{\rm min}$, which meets the flux decay time scale that we obtained in the previous section (see also in Figure 2). The result suggests that the flux decay is dominated by the EC cooling of scattering DT photons and proves the dissipation region is located outside the BLR.