Exploring High-energy Emission from the BL Lacertae Object S5 0716+714 with the Fermi Large Area Telescope

We present the results of an extensive γ-ray data analysis of the emission from the blazar S5 0716+714 with the primary motivation to study its temporal and spectral variability behavior. In this work, we extract a 10 days binned γ-ray light curve from 2008 August 4 to 2016 April 27 in the energy range of 0.1–300 GeV and identify six outburst periods with peak flux of >4 × 10−7 ph cm−2 s−1 from this highly variable source. The brightest flares are identified by zooming in these outburst periods to 1 day binning and using the Bayesian Blocks algorithm. The fastest variability timescale is found to be 1.5 ± 0.3 hr at MJD 57128.01 ± 0.01 with a peak flux above 100 MeV of (26.8 ± 6.9) × 10−7 ph cm−2 s−1. No hint of periodic modulations has been detected for the light curve of S5 0716+714. During the outburst phases, the γ-ray spectrum shows an obvious spectral break with a break energy between 0.93 and 6.90 GeV energies, which may be caused by an intrinsic break in the energy distribution of radiating particles. The five highest-energy photons, with E > 100 GeV, imply that the high-energy emission from this source may originate from a moving emission region in a helical path upstream in the jet. The spectral behavior and temporal characteristics of the individual flares indicate that the location of the emission region lies in the sub-parsec scale (rγ < 0.85 pc).


Introduction
Blazars, including BL Lac objects and flat-spectrum radio quasars, whose relativistic jets are pointed close to our line of sight, represent a small subclass of active galactic nuclei (AGNs), which are extremely variable objects in the sky (e.g., Blandford & Königl 1979;Urry & Padovani 1995). The highly luminous, rapid, and variable broadband flux of blazars is observed throughout the entire electromagnetic spectrum. The spectral energy distributions (SEDs) of blazars are dominated by nonthermal emission and consist of two distinct components: a low-energy component with a peak between submillimeter and UV/X-rays, and a high-energy component with a peak at MeV-GeV energies. It is widely believed that the low-energy component of the SED is due to synchrotron emission from relativistic electrons, whereas the origin of the high-energy component is still a matter of debate. Therefore, two classes of models, leptonic models (see, e.g., Finke et al. 2008;Ghisellini et al. 2010;Ding et al. 2019;van den Berg et al. 2019) and hadronic models (see, e.g., Böttcher et al. 2013;Cao & Wang 2014;Dermer et al. 2014;Cao et al. 2020), are widely used to explain the high-energy emission.
The blazar S50716+714, discovered in 1979 by Perley et al. (1980), is one of the brightest BL Lac objects observed in the γ-ray band. No emission or absorption lines were identified in this source, and it is highly variable from radio to γ-ray energies (Wagner et al. 1996). Its redshift, z=0.31±0.08, was estimated by using its host galaxy as a standard candle (Nilsson et al. 2008). Danforth et al. (2013) set a statistical upper bound of z < 0.32 with a 95% confidence for this source by using Hubble Space Telescope observations. This source is classified as an intermediate-synchrotron-peaked blazar (ISP:10 14 n < < Hz 1 0 peak s 15 Hz) in the 1LAC, 2LAC, and 4LAC catalogs (Abdo et al. 2010b(Abdo et al. , 2010cAckermann et al. 2011), but a low-synchrotron-peaked blazar in the 3LAC catalog (LSP: n < 10 Hz peak s 14 ). 9 It is a famous intraday variability (IDV) source in the radio and optical bands with different timescales from hours, months, to even years (see, e.g., Xie et al. 2002;Gupta et al. 2012;Dai et al. 2015). Its IDV in radio bands is attributed to interstellar scintillation (Bignall et al. 2003). High optical polarization of ∼29%-30% and periodic variability on different timescales have also been observed (see, e.g., Raiteri et al. 2003;Rani et al. 2010Rani et al. , 2011Rani et al. , 2013a. Rani et al. (2013a) found that the optical spectrum of S50716+714 tends to be bluer with increasing brightness.
The first γ-ray detection and variability measurement of S50716+714 came from the Energetic Gamma Ray Experiment Telescope on board the Compton Gamma Ray Observatory (Lin et al. 1995). It is also well detected by other γ-ray observatories, such as AstroRivelatore Gamma a Immagini Leggero (AGILE; Chen et al. 2008), the Major Atmospheric Gamma Imaging Cherenkov Telescopes (MAGIC; Anderhub et al. 2009), and the Fermi Gamma-ray Space Telescope (see, e.g., Abdo et al. 2009aAbdo et al. , 2009bAbdo et al. , 2009cAbdo et al. , 2010c. Its high-energy spectra (X-ray and γ-ray) can be well described by a logparabola (LP) or a broken power law (BPL; Ackermann et al. 2011;Rani et al. 2013b;Wierzcholska & Siejkowski 2016). The combined GeV-TeV spectrum of S50716+714 might exhibit absorption-like features in the energy range of 10-100GeV. (Sentürk et al. 2011(Sentürk et al. , 2013 fitted the GeV-TeV spectrum using a power law (PL) with absorption from line emission of H I + He II. However, a double-absorption scenario does not provide an improved fit over the He II single-line absorption. In a full broad-line region (BLR) absorption model, they found that absorption from the He II complex seems to be dominant.
The high-energy γ-ray flaring behaviors from 0716+714 and other blazars are complex. Some studies that attempted to understand broadband flaring activity of S50716+714 have been done (see, e.g., Rani et al. 2014Rani et al. , 2015Chandra et al. 2015). However, we still do not have a clear understanding of the exact origin of blazar variability, especially the location and origin of the γ-ray flares. Short-timescale γ-ray variability is one important step in understanding the physical γ-ray mechanisms of blazars (Aharonian et al. 2017;Böttcher 2019). The γ-ray properties of some blazars have been studied with minute-scale and hour/sub-hour scale γ-ray variability (see, e.g., Aleksić et al. 2011;Ackermann et al. 2016;Prince et al. 2017;Shukla et al. 2018;Ding et al. 2019). These studies implied that the short-timescale γ-ray variability may be produced in a compact high-energy emission region located far away from the black hole, at the edge or outside of the BLR. Moreover, the observed rapid variability may suggest the dissipation of a magnetic island or collimated proton beams from turbulent plasma at the end of the magnetic nozzle (Shukla et al. 2018). Therefore, it is valuable to investigate short time variability and spectral properties of S50716+714 in the γ-ray band, which will allow us to constrain the location of the γ-ray emission region and explore the origin of highenergy emission.
In the present work, we report a long-term study of the source using GeV observations made by the Fermi Large Area Telescope (LAT). Specifically, we focus on the individual outburst phases observed between 2008 August and 2016 April. Although some outburst phases of S50716+714 have been studied before, a comprehensive analysis that includes all its outburst phases, short-timescale characteristics, and spectral evolution has not been presented in the literature. We identify six major outburst phases in the energy band >100 MeV during 2008-2016, and we analyze the short-timescale characteristics and spectral evolution from four of these outburst phases in detail. Hour-scale timescale and a sharp break energy with high significance are first found, and their possible physical implications in the γ-ray band are discussed.
The framework of this paper is as follows. In Section 2, we provide a brief description of the observations and the data reduction. The results on the flux variability of S50716+714 are shown in Section 3. The γ-ray spectra are presented in Section 4, and a discussion of our results is provided in Section 5. Finally, a summary is given in Section 6. A flat Λ cold dark matter (ΛCDM) cosmology with H 0 =69.6 km s −1 Mpc −1 , Ω M =0.29, Ω Λ =0.71, Ω κ =0 is used in this paper (Planck Collaboration et al. 2014).

Fermi-LAT Observations and Data Reduction
The data used in this study were obtained from the Fermi-LAT public data server. 10 The LAT is the primary instrument on the Fermi Gamma-ray Space Telescope. It is a pairconversion γ-ray telescope with an energy range from 20MeV to more than 300GeV, and it can scan the entire sky in approximately 3hr with a field of view of ∼2.4sr. This configuration allows for the study of the high-energy properties of AGN and the short time evolution of γ-ray sources. For a detailed description of the LAT detector, see Atwood et al. (2009).
In this paper, we analyze γ-ray data collected between 2008 August 4 and 2016 April 27 (MJD=54682-57505) in the 0.1-300 GeV energy range. We use the standard procedures described in the Fermi-LAT documentation 11 to analyze the data with the LAT analysis software ScienceTools v10r0p5 with P8R2_SOURCE_V6 instrument response functions. Photons are selected in a circular region of interest (ROI) of 10°radius centered at the location of S50716+714. Only photon-like events classified as evclass=128, evtype= 3 are selected. A zenith angle cut of 90°is used to remove the contamination of the background γ-rays from the Earth limb. We include all the sources from the Fermi-LAT Third Source Catalog (3FGL; ) within 20°of S50716+714 in our model file. The parameters of point sources within 10°of S50716+714, either the normalization only, or other spectral parameters, are left to vary freely during the likelihood fitting. The parameters of all sources outside the ROI are fixed to the values published in the 3FGL catalog. We use an unbinned maximum likelihood technique to investigate the photon flux and index variations. There are four sources with a detection significance (maximum likelihood test statistic (TS) 12 ) of TS > 25 in the ROI, namely, 3FGLJ0841.4+7053, 3FGLJ0855.4+7142, 3FGLJ0814.7+6428, and 3FGLJ0805.4 +7534. The normalized parameters for these sources are freed in the fits with 7 days and 1 days binning. When some data points exhibit non-convergence in the process of gtlike, we fix photon indexes and normalized parameters for these sources with TS > 25 lying within 10°except for the target source. For sources with TS values less than 25, their photon indexes and normalized parameters are fixed regardless of their location. The isotropic background model iso_P8R2_SOURCE_V6_v06. txt and the Galactic diffuse-emission model gll_iem_v06. fits are used 13 ). In addition, we reanalyze the data, including all sources from the latest Fermi-LAT Fourth Source Catalog (4FGL catalog; Ajello et al. 2020) to check if any of the new 4FGL sources in the region around S5 0716+714 are significant enough to affect the analysis. We find no significant difference in the two cases.
We perform different spectral analyses of several epochs of flaring activity by fitting the γ-ray spectra with the PL, LP (E 0 : the peak energy fixed at 440 MeV), or BPL function forms over the 0.1-300 GeV range by using the Unbinned Likelihood analysis 10 https://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi 11 https://fermi.gsfc.nasa.gov/ssc/data/analysis/ 12 = -  TS 2 log log 0 ( () ), where  and  0 ( ) are the likelihood with and without the source. 13 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html package. These functions are described on the Fermi Science Support Center website. 14 We use the TS curve value, which is calculated as , to evaluate the significance of spectral curvature (Nolan et al. 2012).

Long-term Light Curve
In Figure 1(a), we plot a 10 days binned light curve of S50716+714 above 100 MeV using the procedures described in Section 2. The light curve exhibits significant flux variability throughout the whole 2008-2016 (∼8yr) period. Six highactivity periods with peak fluxes above 4×10 −7 ph cm −2 s −1 are found (see Table 1), which allow us to investigate the evolution of outbursts with a shorter temporal resolution, down to ∼3hr (i.e., the single sky survey time for LAT).
We calculate the probabilities that the detected photons with energies E > 10 GeV are associated with S50716+714 by using gtsrcprob. Figure 1(b) shows those photons that have probabilities greater than 3σ. The high-energy photons during the outburst phases are listed in Table 1. The five highestenergy photons, with measured energies greater than 100 GeV, are displayed and numbered in Figure 1(b). The highest-energy photon was detected with the energy of 207 GeV, implying that the electrons in the jet are accelerated to higher energies, assuming a leptonic scenario for production.

Identifying the Flares of S50716+714
For the six outburst periods, when we zoom-in to 1 day binning, some substructures, and various phases are clearly seen in Figure 2, including the quiescent, pre-flare, inter-flare, flare, and post-flare. The interval before the flaring is defined as "pre-flare," the interval between the two flares is defined as "inter-flare," and the interval after the last flare is defined as the "post-flare." 1 < E γ < 300 GeV during 2008 August 4∼2016 April 27 (MJD 54682-57505) using the LP model with 10 days bins. Bottom: the distribution of photons with E > 10 GeV with different significance levels of 3σ (green circles) and 4σ (red triangles). The numbers 1-5 represent photons with energy of E > 100 GeV, implying S50716+714 has Very High Energy (VHE; E > 100 GeV) emission. Shaded regions indicate the six active phases studied in this work. The gray dashed line presents E > 100 GeV. true disp ). The strong flares exceeding the threshold value are identified from the block representation, which are shown with the blue-violet-shaded regions. The green-shaded regions represent monotonically decreasing sets of adjacent block, which are a part of identified flares. Here, considering Flare 1 of Phases IV and V exhibit a significant outburst behavior with the duration of ∼7 days, they are not a part of identified Flare 2 of Phases IV and V. The horizontal orange dashed lines represent the average flux. The blue lines represent the BB with the false alarm rate parameter p 0 =0.05. Points with TS < 4 are not shown. Nalewajko (2013) defined the flare as a continuous period of time with the flux greater than half of a given flux peak. However, this definition requires that every two flares must be separated and thus cannot treat overlapping flares. Here, we adopt the Bayesian Blocks (BB) algorithm to identify the flares for the source S50716+714. The BB method proposed by Scargle (1998) can divide a data set comprising N p photon elementary cells into N b longer blocks. A simple nonparametric analysis model with an improved and generalized version of the BB algorithm has been proposed to find the optimal segmentation of the data in the whole observation interval (Scargle et al. 2013). This algorithm is also included in the astropy.stats 15 and astroML packages. The BB algorithm with some improvements has been applied to all the photon event data in high-energy astrophysics, including timetagged events, binned counts, and time-to-spill data, with no lower limit on the timescale (see, e.g., Ahlgren et al. 2019;Kerr 2019;Meyer et al. 2019). In this work, we adopted directly the BB algorithm from the astropy.stats to implement identification of flares, with the false alarm rate parameter p 0 =0.05 for 1 day binning light curves.
In order to obtain the time range of the flare, we consider monotonically decreasing sets of adjacent blocks. We adopt the "burst_def" function from Ivan Kramarenkoʼs algorithm to select time ranges of the flares. 16 In this algorithm, an iterative approach is used to split the points into two sets, low-state points and everything else (anti-set) as implemented in the "burst_def" function. At the start of the algorithm, every point is in the low-state set. In each iteration, the time averaged flux of the low-state set F is calculated, as well as the absolute deviation of each point in the low-state set from this average (abs(F i -F )). The maximum of this deviation is called F disp and the associated point is moved from the low-state set into the anti-set. This is iterated until F disp /F anti < 0.25, where F anti is the mean flux of the anti-set. The flare group is also obtained by calculating the dispersion of the low-state group. If the flux levels of the block from the flare groups are less than ( +´-F F 2 true disp (dispersion of low state)) for the 1 day binned light curves, these flares are too weak and we remove them from the flares group. Finally, we can distinguish the flare groups from the low-state groups in all the BBs.
We also extract 12 hr and 6 hr binned light curves of all flares and find that Phases III and VI do not exhibit a sufficiently significant outburst in 12 hr and 6 hr binning. Therefore, we do not study short-timescale variability of Phases III and VI. Because the statistical fluctuations may tend to dominate in the short-timescale variability, we adopt the false alarm rate parameter both = p 0.32 0 and p 0 =0.05 to identify peaks in 12 hr, 6 hr, and 3 hr binning. Here, we adopted p 0 =0.32 to select peaks. If the flux levels of the block from the flare groups are less than ( +´-F F 1.5 true disp ) for the 12 hr and 6 hr binned light curves, these flares are removed from the flares group. Moreover, data points with TS < 4 are rejected from the light curve analysis. (Note the factor 1.5 Table 2 Best-fitting Parameters of the 12 hr, 6 hr, and 3 hr Binned Light Curves of Flares 3 and 4 of Phase V, as Marked in Figure 2 Note. The standard deviations of T r and T f are obtained by light curve fitting with Equation (2). As the timescales are strictly positive (i.e., T r and T f 0), the maximum value of lower errors is equal to the values of T r and T f . We adopted asymmetric error bars, namely, the upper errors and the maximum lower error, and corresponding doubling times are also considered. The error in lower frequencies is estimated from the rms scatter of the frequency points in the logarithmic bins, and for the higher frequencies the error is based on the rms scatter of the 6 PSDs used to compute the averaged PSD. As the number of Fourier frequencies in each logarithmic bin decreases toward the lower frequencies the errors become correspondingly more uncertain. The lowest frequency bins lack error bars since they contain only one Fourier frequency. compared to 2.0 for the daily binning. Empirically, we find that larger values, like 2.0 or 3.0 missed some flares, while a smaller value identifies too many weak features.) Since Flares 1 and 2 during Phase IV and Flares 1 and 2 and Flares 3 and 4 during Phase V are located in a continuous time interval, the flares are combined for these phases with 12 hr and 6 hr binning to conveniently study the evolution characteristics of the shorttimescale variability and the spectra from the flares.

Flux Variability
To probe the evolution of outbursts with a finer time resolution that depends on the γ-ray brightness of the source, we use the following function to fit the time profile of a single flare (Abdo et al. 2010d): In this function, F 0 is the flux at time T 0 representing the approximate flare amplitude. T 0 approximately describes the time of the peak (i.e., it corresponds to the actual maximum only for symmetric flares), while T r and T f represent the characteristic timescale of the rise and decay of the flares, respectively. In fact, it is hard to handle the statistical fluctuations, especially for the very short timescales. Some multiple peaks are also found in the light curves of S50716+714. In order to fit the multiple peaks, we first identify each individual peak by the BB algorithm. Then we fit each individual peak component by the function provided in Equation (1) and obtain the best-fit parameters of the function. Finally, the total function for the multiple peaks is obtained by including n individual peak components and one constant background component. Therefore, the form of the multiple peaks can be defined as The constant flux (F C ) is obtained by fitting a constant. Although the minimum flux observed over time also could be taken as the constant flux, it is a rather crude measure. The subsequent fitting of light curves was performed by fixing the peak position.
Following the fit method above, we fix the location of each individual peak component to fit the 12 hr, 6 hr, and 3 hr binned light curves with Equation (2). In general, the values of the χ 2 divided by the number of degrees of freedom (ndf) are between ∼1 and ∼2 (see Table 8 in AppendixA). The 12 hr and 6 hr binned light curves cannot be well fitted by Equation (2) if the values of the χ 2 /ndf are too large. In this case, we add some substructures peaks to improve the fit. In fact, the profiles of γ-ray flares with the short-timescale variability can be described better by adding some substructures, which allows us to better describe the flareʼs properties and explore the γ-ray emission of S50716+714.  Figure 1(a), which were fitted by the PL (orange lines), LP (violet dashed-dotted lines), and BPL (magenta dashed lines) spectral models. Their best-fitting parameters are given in Table 3. The time (T p ) of maximum of a flare is calculated by the following equations : where T p is equal to T 0 when T r =T f . The total duration of a flare can be estimated by . The symmetry of a flare can be described as The parameter ξ can be used to define three different profiles: (1) if x < 0.3 | | , the flares have a symmetric temporal profile; (2) for moderately asymmetric flares, and (3) for markedly asymmetric flares, 3.3.1. Variability of Phase I Figure 2(a) presents the light curve of Phase I with 1 day binning. Three flaring periods (Flare 1, Flare 2, and Flare 3) can be clearly identified by the BB algorithm, and have peak fluxes of greater than 8×10 −7 ph cm −2 s −1 , which allow us to analyze the evolution of outbursts with fine time resolution. A pre-flare phase is observed with increasing flux from MJD 55592.66-55621.15. Inter-flares 1 and 2 with some substructures are observed from MJD 55632.67-55745.97 and 55761.00-55849.98, during which the flux is less than 5×10 −7 ph cm −2 s −1 . A post-flare is observed from MJD 55859.66-55882.66, during which the flux decreased gradually.
We also extract light curves of all flaring periods with 12 hr and 6 hr bins, which are shown in the top panel and middle panel of Figure 8 in Appendix A. We find that Flares 1 and 2 have three distinctive peaks: P1, P2 and P3, where the flux exceeds 4.7×10 −7 ph cm −2 s −1 . Flare 3 consists of two peaks, P1 and P2, as shown in the bottom panel of Figure 8 in Appendix A. The highest fluxes are (14.1±4.2) and (12.1±3.9)×10 −7 ph cm −2 s −1 at MJD 55854.11 and 55854.35 for the 12 hr and 6 hr bins, respectively. The shortest variability timescale during Phase I is (2.3±1.7 hr) with a symmetric temporal profile. All the values of the fitted parameters are shown in Table 5 in Appendix A.

Variability of Phase II
Similar to Phase I, we also show the 1 day binned light curves of Phase II in Figure 2(b). Four flaring periods (marked as Flares 1-4) that have duration of a few days to two weeks, can clearly be seen in Figure 2(b). Pre-flare and post-flare Table 3 Results of the Spectra Using Model Fitting for Phases I -VI Note. D ln L ( ) is half of TS curve . We follow Wilks' theorem to estimate the significance, i.e., twice the difference between the log(likelihood) values for the two spectral models is formally distributed as χ 2 with Dn degrees of freedom, where Δn is the difference in the number of degrees of freedom between the two models (Rani et al. 2013c). phases were also observed between MJD 56092.66-56109.66 and 56185.66-56212.66, during which substructures could be clearly seen. Inter-flare 1 was observed from MJD 56122.66-56137.66, during which the flux was close to a constant (average flux: 2.4×10 −7 ph cm −2 s −1 ).
The 12 hr and 6 hr binned light curves of the flares are shown in the top panels of Figures 9 and 10 in Appendix A. Three peaks, P1, P2, and P3, are observed during the Flare 1 period with 12 hr and 6 hr binning (see the top panel of Figure 9 from Appendix A). These peaks are identified by the BB algorithm, where the maximum observed fluxes are (6.7±1.7) and (7.3±2.3)× 10 −7 ph cm −2 s −1 at MJD 56112.70 and 56116.38 with 12 hr and 6 hr binning, respectively.

Variability of Phase IV
The 1 day binned light curve of Phase IV is shown in Figure 2(c). A pre-flare phase, where the flux is close to the quiescent state (average flux: 1.5×10 −7 ph cm −2 s −1 ), extends from MJD 56692-56749, after which the source entered into an outburst phase. Two flares (Flare 1 and Flare 2) can be clearly seen during MJD 56750.21-56759.09 and 56759.09-56772.66.

Variability of Phase V
The 1 day binned light curve of Phase V displays seven variability patterns (a pre-flare, Flare 1, Flare 2, an inter-flare, Flare 3, Flare 4, and a post-flare) in Figure 2(d). The pre-flare and post-flare phases are observed with some substructures between MJD 57030.66-57049.66 and 57113.66-57133.66, respectively. The pre-flare phase includes two states, namely, an active state and a quiescent state (average flux of 1.5×10 −7 ph cm −2 s −1 ).
We also extract 12 hr and 6 hr binned light curves of Flares 1 and 2 (see the top panels of Figure 3). Three peaks, P1, P2, and P3, are observed in the 12 hr and 6 hr binned light curves, where the maximum observed fluxes are (8.1±2.4) and (8.9±2.9)×10 −7 ph cm −2 s −1 at MJD 57048.65 and 57048.72 with 12 hr and 6 hr binning, respectively. Their light curves exhibit symmetric temporal profiles except for P2 of 12 hr binning.
The 12 hr, 6 hr, and 3 hr binned light curves of Flares 3 and 4 are also extracted in Figure 3 (bottom panel). Two peaks, P1 and P2, are clearly seen for the light curves of three different time bins. The peak fluxes of P2 are larger than the flux values of P1. The highest peak flux of (26.8±6.6)× 10 −7 ph cm −2 s −1 is found at MJD 57128.01 with 3 hr binning. The flux-doubling times of P2 are (1.6±2.8) hr, (1.8±0.8) hrm and (1.5±0.3) hr with 12 hr, 6 hr and 3 hr binning, respectively. All the peaks in the 6 hr and 3 hr binned light curves exhibit a symmetric temporal profile. The modeling parameters are given in Table 2. Table 8 in Appendix A shows additional information about the fitting of the flares in the various phases. The χ 2 /ndf indicates that the fitting is reasonably good in all cases.

Power Spectrum Density
To explore the nature of the variability and to search for periodicity, a power spectral density (PSD) analysis method is used (Vaughan 2005;Chidiac et al. 2016). The PSD of the full 10 days binned light curve is shown together with the average PSD of the six phases in the 1 day binned light curves ( Figure 4). The shift between the two PSDs in the frequency range where they overlap corresponds to a 1.7 times higher fractional variance (variance/(mean square)) in the flaring state (blue line: the 1 day binned data segments) relative to that of the total light curve (red line: the 10 days binned light curve). The power-law fit (magenta line) to the high-frequency part of the outburst PSD has a slope of α=−1.5±0.2 (uncertainty based on standard deviations from the bin averaging). Although the combined PSD seems to indicate some feature or break below 1/ -10 2 days, this feature may be an artifact of the selection of flare light curves due to the increasing statistical and systematic uncertainties toward lower frequencies. The PSD analysis does not indicate the presence of periodic variations in the source.

Time-resolved SEDs
We fit SEDs of S50716+714 during different active states with three different models (PL, LP, and BPL). The SEDs of Phases I to VI as marked in Figure 1(a) are shown in Figures 5(a)-(f). The different colored lines (orange lines, violet dashed-dotted lines, magenta dashed lines) represent the different fitting models (PL, LP, and BPL), respectively. Their corresponding fitted parameters are given in Table 3. They show spectral breaks, where the range of break energies is between 0.93 GeV and 6.90 GeV. Their TS curve values with respect to the PL model are between 16.4 and 29.2 (corresponding to significances that exceed 4 σ) except for Phases III & IV. We therefore conclude that both the BPL and LP models better describe the γ-ray spectral shapes than the PL model.
The same spectral models are also applied to shorter time intervals of Phases I, II, IV, and V. The corresponding γ-ray SEDs and the fitted parameters are shown in Figures 13-16 and in Tables 10-12 in Appendix C, respectively. These spectra also show a hint of the curved or break shapes at the different states, although the corresponding significance levels are less than 3 σ. A possible reason is the poor photon statistics at shorter timescales, but we cannot also rule out other possible causes, including a combination of different factors such as jet dynamics or the geometry of substructures in the emission region (Tanihata et al. 2001;Kushwaha et al. 2014).
A significant spectral hardening with higher flux is also seen during Phases I, II, IV, and V. Using the PL fit for these short intervals, the photon index (Γ) changes are 2.17±0.05 to 1.93±0.04 (Figure 6(a)), 2.07±0.04 to 1.92±0.04 (Figure 6(b)), 2.10±0.08 to 1.92±0.04 (Figure 6(c)) and 2.09±0.06 to 1.91±0.03 (Figure 6(d)), respectively. Similar behavior was also found for some other bright Fermi blazars, like 3C279 (Paliya 2015a), S50836-71 (Paliya 2015b), 3C454.3 (Britto et al. 2016), PKS1510-089 and CTA102 (Prince et al. 2017(Prince et al. , 2018. The variations of the difference of spectral slopes (ΔΓ) with respect to the photon flux, the break energy (E break ), or the highest photon energy for the six strong outburst phases are shown in Figure 7. We do not find any significant variation in ΔΓ with respect to the flux, E break , or the highest photon energy during Phases I to VI. For the different activity periods during Phases I, II, IV, and V, we also do not find any significant variation in ΔΓ with respect to the flux and to E break , neither in E break with respect to the photon flux (see Figure 17 of Appendix C). This is similar to what has been found for other Fermi blazars (Ackermann et al. 2010;Rani et al. 2013c).

Discussion
We have used 8 years of observations with the Fermi LAT to explore the high-energy emission properties of the BL Lac S50716+714. The motivations for this study are: (i) to investigate both the short and long-term variations in the source, (ii) to explore its spectral variations, and (iii) to pinpoint the location and origin of the observed γ-ray flares in the source.
We do not find any hint of periodic modulation in the light curves. The PSD analysis (Section 3.3.5) of the long-term (10 days binned) and short (1 day binned) flux variations is found to be consistent with a red-noise slope/index of −1.5. The fractional variance for the 1 day binned data is comparatively high, which suggests higher variability power during the episodes of rapid flares.
Another feature of the light curves is the temporal profile of flares. For S50716+714, all the flares with the 12 hr and 6 hr binned light curves are symmetric within the error bars (see Table 2 and Tables 5-7 of Appendix A), which is consistent with the results of long-term outbursts (Chatterjee et al. 2012;Roy et al. 2019) and is slightly different from results of shorttimescale variability (Roy et al. 2019). This symmetric temporal profile suggests that the rise and decay timescales are dominated by the crossing time of radiation or a disturbance through the γ-ray emission region. An asymmetric temporal profile would be expected if caused by a fast injection of accelerated particles and a slower radiative cooling or an escape from the emission region (Sikora et al. 2001). Roy et al. (2019) found that the majority of the short-term (∼1 day) flares in their study exhibit an asymmetric temporal profile, indicating that the asymmetric temporal profile may be due to the gradual acceleration of the particles to the GeV energy band or the change of the width of the emission region and the bulk Lorentz factor of the plasma. The symmetric flares of S50716 +714 suggest that the relevant timescale is the light-crossing time of the emission region, which can be explained by the superposition of several episodes of short duration. It has also been suggested that different amounts of Doppler boosting for different shells in the emission region may also be responsible for the asymmetry of γ-ray flares (Nalewajko et al. 2014). Moreover, the flare with the slower rise timescales than the decay timescales (negative ξ) may be attributed to the gradual acceleration of particles to the GeV energy band, indicating that the radiation cooling timescale is shorter than the acceleration timescale. The increased cooling time may cause a faster rise timescale than the decay timescale (positive ξ), which may imply a change in the bulk Lorentz factor of the plasma or the width of the shocked region (Roy et al. 2019). Therefore, the observed γ-ray flares in the S50716+714 may originate from a combination of different physical processes.

γ-Ray Doppler Factor
We can constrain the Doppler factor by using the highestenergy detected γ-ray photon. The high-energy photons can interact with lower-energy photons by pair production, if the bulk of the high-energy emission (γ-rays and X-rays) is produced in the same emission region. The maximum cross section of this process is ∼σ T /5 (Svensson 1987), where σ T is the Thomson cross section. If the optical depth t n gg ( ) of the photon field is t n < ( ) is the dimensionless energy of a γ-ray photon with highest energy E max when the optical depth of the emitting region is t n = gg 1 ( ) . The luminosity distance d L that corresponds to z=0.31±0.08 is d L =1.55 Gpc (Jorstad et al. 2017), t var is the shortest observed doubling/halving timescale, and it is approximately equal to T r ×ln2 (see, e.g., Rani et al. 2013c). We obtain δ γ 4.6-5.4 (see Table 4), which is quantitatively consistent with the result (δ γ 5.9) from Dondi & Ghisellini (1995). Rani et al. (2013a) discussed the variability properties of S50716+714 from radio to γ-ray wavelengths and derived a Doppler factor δ γ 9.1 above 200 GeV and δ γ 9.8 above 400 GeV. The estimated δ γ value is comparable to δ VLBI =6-21 observed by Rani et al. (2014Rani et al. ( , 2015.

Size of the Emission Region
The fastest flux-doubling timescale of the source can be used to constrain the size of the emission region. The size of the γ-ray emitting region can be estimated by Using δ=15.6±4.0 estimated by radio VLBA observations (Jorstad et al. 2017), we derive an upper limit on the size of the γ-ray emission region radius R∼2.0-3.5×10 15 cm during the different outburst phases. The upper limit on the angular size [mas] of the emission region can be calculated by Rani et al. (2013aRani et al. ( , 2013c as: We obtain θ≅0.4-0.7 μas during the different outburst phases, which is much smaller than the range of the core sizes 0.07-0.09 mas at 15 GHz and 0.01 mas at 230 GHz observed by Lee et al. (2017). This result implies that the origin of emission is likely different for the γ-ray and radio bands.

Origin of the Spectral Break/Curvature in S50716+714
The γ-ray spectra of the source significantly deviate from a PL (see Figure 5) and show GeV spectral breaks. The origin of spectral breaks is still a matter of debate. Many theoretical models are proposed in the literature to explain the observed spectral breaks (Ding et al. 2019).
The attenuation of γ-rays via photon-photon pair production on He II Lyman recombination lines within the BLR may be responsible for the observed breaks (Poutanen & Stern 2010). Sentürk et al. (2011) proposed that the γ-γ absorption of the full BLR also can produce the observed spectral breaks if the γ-ray dissipation region lies inside the BLR. The GeV spectral breaks also observed in many bright Fermi blazars, like 3C 454.3, 3C279, 4C+21.35 etc. were also interpreted by the attenuation of γ-rays via photon-photon pair production and γ-γ absorption of the full BLR models. On the other hand, MAGIC Collaboration et al. (2018) reported that the VHE emission from S50716+714 originated in the entrance and exit of a superluminal knot in and out of a recollimation shock in the inner jet, which is attributed to a shock in the helical jet downstream of the core. This TeV result suggests that the origin of the γ-ray emission within the BLR may be ruled out.
Alternatively, the GeV spectral breaks could also be explained as a transition of inverse Compton scattering from the accretion disk (in the Thomson regime) to disk emission reprocessed in the BLR (taking place in the Klein-Nishina regime) (Finke & Dermer 2010). This model focuses on the powerful FSRQs with luminous BLRs. However, for S50716 +714 the BLR is probably weak due to an inefficient accretion process in BL Lac objects. Abdo et al. (2009c) proposed that the spectral break could be attributed to radiative cooling. A key feature of the radiative cooling break is that the change of spectral slope (ΔΓ) is close to 0.5. For the six major outburst phases and the different activity states, only a few values are close to 0.5. Therefore, it is difficult to reconcile the constancy of the break energy with respect to the photon flux variations within the cooling break. From this we conclude that the observed spectral breaks in 0716+714 are unlikely to have an intrinsic origin associated with radiative cooling.
A spectral break can be produced if there is a cutoff/ curvature in the energy distribution of particles. Using an equipartition approach, the GeV break may arise from the onset of Klein-Nishina effects on the Compton scattering of BLR photons, and with the continuously curving electron energy distribution given by a log-parabola function, this continuously curving electron energy distribution derives from stochastic acceleration processes with radiation and escape (Cerruti et al. 2013;Dermer et al. 2014). Therefore, the GeV spectral breaks of S50716+714 are likely due to a break in the energy distribution of particles, based on the measured γ-ray spectral shapes.

Location of the Gamma-Ray Emission Region
We can constrain the location of the emission region by using the observed fastest variability timescale, although this approach may not apply to cases with very rapid variability (e.g., Aleksić et al. 2011). The location of the emission region can be estimated by < G + The energy dependence of a falling timescale T f of a flare can also be used to estimate the location of the γ-ray emission region (Dotson et al. 2012). By considering that the maximum decay time difference is between a high-energy E HE and a lowenergy E LE , the condition that Δt max T f (LE)-T f (HE) can be used to derive an upper limit for the location r γ of the γ-ray region [cm]: For S50716+714, we find that the 12 hr binned light curves do not show a significant flare profile, in particular for the 1-300 GeV energy band. Therefore, we fit the 12 hr data from Flare 3 of Phase I, Flares 2-4 of Phase II, Flares 1 and 2 of Phase V and Flares 3 and 4 of Phase V, because these light curves show a significant flare profile for the 0.1-1 GeV and 1-300 GeV energy bands. Only Flares 3 and 4 of Phase V among these light curves show a significant flare profile in the 6 hr binned time for the 0.1-1 GeV and 1-300 GeV energy bands, so we only fitted 6 hr binned data from the Flares 3 and 4 as marked in the 1 day binned light curve of Figure 2(d). The light curves fitted by Equation (2) and the best-fitting results are listed in Figures 11,  12 and Table 9 of Appendix B. Therefore, the allowed range of r γ from the outburst Phases I, II, IV, and V by using the above methods is 0.016-0.23, 0.019-0.85, 0.028-0.38, and 0.016-0.61 pc (see Table 11 of Appendix C), respectively.
Our results therefore suggest that the location of the γ-ray flaring activity observed in S50716+714 lies within 0.016-0.85 pc of the central engine. Based on the upper limit to the disk luminosity of S50716+714 of 1.8×10 44 erg s −1  pc. Therefore, the location of the emission region lies in the sub-parsec scale near the outer edge or well outside the BLR, which is consistent with results from Rani et al. (2014Rani et al. ( , 2015 (where the distance of the γ-ray emission region cannot be larger than 6.5 pc). Rani et al. (2013a) found that the SSC+external Compton (EC) model can well describe the broadband SED of the source if the external radiation field is dominated by Ly-α emission from a putative BLR with an external radiation field energy density of 10 −6 -10 −5 erg cm −2 s −1 , which is factor of ∼1000 lower than what we expect for a typical quasar. This low BLR energy density may explain the γ-ray spectral breaks in the S50716 +714 and be consistent with the non-detection of emission lines. However, this model does not provide a good fit to the radio data and the VHE data.

Summary
In this paper, we presented a detailed investigation of the γ-ray flux and spectral variations of the blazar S50716+714 for 8 years of Fermi-LAT observations, from 2008 August 4 to 2016 April 27. The source displays significant flaring activity after 2011, with six major outburst phases and many substructures found in the 10 days binned light curve. Each individual outburst phase is further studied with the 1 day, 12 hr, and 6 hr binning. The shortest variability timescale is 1.5±0.3 hr, with peak flux of (26.9± 6.1)×10 −7 ph cm −2 s −1 , which is close to the lightcrossing time, indicating a very compact and anisotropic emission region in the inner jet. The shortest variability timescale is also comparable to the fast γ-ray flares observed in other Fermi-LAT blazars, which can put a constraint on the size of the emission region of R∼2.0-3.5×10 15 cm and a lower limit of the γ-ray Doppler factor of δ γ 4.6-5.4. Our results also indicate that the γ-ray emission region is located at r γ ∼0.016-0.23 to 0.019-0.85 pc, suggesting the emission region lies near or beyond the outer edge of the BLR.
The short-timescale γ-ray flares show symmetric properties, suggesting the variability may be dominated by the crossing time of radiation or a disturbance through the activity region rather than by the acceleration or energy-loss timescales of the radiating electrons (Chatterjee et al. 2012). Despite the presence of flaring activity, we do not detect any periodic behavior in the source light curve. The variations are well characterized by a red-noise slope of −1.5. The shorttimescale variability may be triggered by the turbulence/ magnetic reconnection, shock-in-jet or interactions of the jet with external media (see, e.g., Aharonian et al. 2017;Böttcher 2019;Ding et al. 2019). Multiband synergistic observational data would be needed to identify such rapid flaring activity.
The energy of the highest-energy photon observed from the source in this study was 207 GeV. The five highest-energy photons with energies greater than 100 GeV have theoretical implication about the nature of the high-energy γ-ray emission region, indicating that (1) the γ-ray emission may be associated to a component entering and exiting the core region (MAGIC Collaboration et al. 2018), (2) the evolution of the spectra of flares lasting approximately months may follow the shock-injet model, (3) the γ-ray emission region has a magnetic field of a few mG (Rani et al. 2013a;Lee et al. 2017).
The SEDs are fitted with three different functional forms: PL, LP, and BPL. Both the BPL and LP models can better describe the γ-ray spectral shapes than the PL model. GeV spectral breaks are observed in the source. The break energies varied between 0.93 and 6.90 GeV, and show evidence for curvature with >4σ significance for Phases I, II, V, and VI. This curvature could be attributed to a break in the energy distribution of particles. The source also follows "harder-whenbrighter" behavior during the different activity states.
Our results seem to support the idea that the γ-ray emission region of S50716+714 varies from one flare to another, suggesting the presence of multiple γ-ray emitting sites in the source. The phenomenon of emission of many HE photons may be explained in terms an emission region moving on a helical path of the jet (Rani et al. 2014(Rani et al. , 2015Raiteri et al. 2003Raiteri et al. , 2017. The γ-ray flaring activity may be triggered by the interaction of moving blobs of plasma and shock (Agudo et al. 2011). Further combined multiband contemporaneous observations are needed to identify the origin of the γ-ray flaring activity clearly and put a stronger constraint on the location of the emitting region. The coming era of multimessenger astronomy will provide more high-quality multiband synergistic observational data (see, e.g., Burns et al. 2019;Rani et al. 2019).    (2). The royal blue lines represent the total fit. The color lines correspond to the contribution of single components in the total fit. The orange shaded regions represent the 68% confidence bands. The gray lines represent the constant flux.

Table 5
Best-fitting Parameters of the 12 hr and 6 hr Light Curves of Flares 1-3 of Phase I, as Marked in Figure 2(a) Peak Note. The standard deviation of T r and T f are same as Table 2.

Table 6
Best-fitting Parameters of the 12 hr and 6 hr Light Curves of Flares 1 to 4 of Phase II, as Marked in Figure 2 Table 7 Best-fitting Parameters of the 12 hr and 6 hr Light Curves of Flares 1 and 2 of Phase IV, as Marked in Figure 2  Sections 3.3.1, 3.3.2 and 3.3.3. We also identify peaks in the short-timescale light curves with the BB algorithm, and with the false alarm rate parameter p 0 =0.05 and 0.32. These peaks are fitted by Equation (2). All plotted light curves are shown in Figures 8, 9, and 10. All fitted parameters are described in Tables 5, 6, and 7. The constant flux and χ 2 /ndf are provided in Table 8, which imply that the fitting is reasonably good in all light curves with 12 hr and 6 hr binning. Here, the light curves cannot be well fitted by the Equation (2) if the values of the χ 2 /ndf and residuals are too large. In this case, we add some substructures or split a peaks into two peaks to improve the fit by visual inspection of the light curves and residuals. So, for the 12 hr and 6 hr binned light curves of Flare 2 from Phase I, we add some substructures to improve the fit. For the 12 hr and 6 hr binned light curves of Flare 3 from Phase I and Flares 1 and 4 from Phase I, we split a peaks into two peaks to improve the fit.
Appendix B Light Curves of Phases I, II, IV, and V for the 0.1-1 GeV and 1-300 GeV Energy Bands We present the light curves of Phase I, II, IV, and V for the 0.1-1 GeV and 1-300 GeV energy bands. Because of the poor statistics in the 1-300 GeV band with the short time bins, many flares do not exhibit a significant flare profile with 12 hr and 6 hr binning. The 12 hr binned light curves of Flare 3 during Phase I, Flares 2-4 during Phase II and Flares 1-2 during Phase IV are presented. The 12 hr and 6 hr binned light curves of Flares 3-4 during Phase V are also presented. These light curves are fitted by Equation (2). All plotted light curves are shown in Figures 11 and 12. All fitted parameters are described in Table 9. The constant flux and χ 2 /ndf are also provided, which imply that the fitting is reasonably good in all light curves with 12 hr and 6 hr binning.     Figure 14. Gamma-ray SEDs of S50716+714 for different activity states during Phase II, as defined in Figure 2(b). These were fitted by the PL (orange lines), LP (violet dashed-dotted lines), and BPL (magenta dashed lines) spectral models. Their best-fitting parameters are given in Table 11. Figure 15. Gamma-ray SEDs of S50716+714 for different activity states during Phase IV, as defined in Figure 2(c), fitted by the PL (orange lines), LP (violet dashed-dotted lines), and BPL (magenta dashed lines) spectral models. Their best-fitting parameters are given in Table 11.