Characterizing the Gamma-Ray Variability of the Brightest Flat Spectrum Radio Quasars Observed with the Fermi LAT

Our goal is to characterize the long-term γ-ray behavior of the selected FSRQs, and the γ-ray behavior of the brightest flares. To this end, we select γ-rays that were measured with the Fermi LAT between 2008 August 4 and 2018 January 30, yielding data over an interval of 114 months, or almost 9.5 yr. The Fermi LAT is a pair conversion telescope designed to measure γ-rays with energies from 20 MeV to above 300 GeV (Atwood et al. 2009).

We follow the standard data selection recommendations⁴ and restrict ourselves to γ-rays in the energy range between 100 MeV and 316 GeV.⁵ Below 100 MeV, the effective area of the LAT quickly decreases, and the point-spread function increases to above ∼6°,⁶ making a point-source analysis challenging. Since FSRQs usually have soft γ-ray spectra, we do not expect significant detection of these sources above our chosen maximum energy.

To mitigate contamination of γ-rays originating from the Earth limb, we further limit the sample to events that have arrived at a zenith angle less than 90°, and we excise periods of bright γ-ray bursts and solar flares that have been detected with a test statistic (TS) >100. The TS is defined as $\mathrm{TS}\,=-2\mathrm{ln}({{ \mathcal L }}_{1}/{{ \mathcal L }}_{0})$, i.e., the log-likelihood ratio between the maximized likelihoods ${{ \mathcal L }}_{1}$ and ${{ \mathcal L }}_{0}$ for the hypotheses with and without an additional source, respectively (Mattox et al. 1996). We use the latest Pass 8 instrumental response functions and Monte Carlo simulations (Atwood et al. 2013) and select γ-rays that pass the P8R2 SOURCE event selection. For each source, we analyze a 10° × 10° region of interest (ROI) centered on the position of each source as provided in the 3FGL (Acero et al. 2015). We choose a spatial binning of 01 pixel^–1 and eight energy bins per decade.

Our analysis proceeds iteratively, starting from the full time range and zooming in on bright flares and shorter timescales (see Section 2.3). In a first step, we optimize the global γ-ray model of each ROI using the Fermi Science Tools version 11-05-03⁷ and fermipy version 0.16.0+188⁸  (Wood et al. 2017). The initial model consists of all γ-ray point sources within 15° of the ROI center included in the 3FGL, as well as the standard templates for isotropic and Galactic diffuse emission.⁹ After an initial optimization, we free the spectral normalization of all sources in the model. The spectral shape parameters, such as power-law indices, curvature, or exponential cutoff energies, are free to vary for sources within 5° of the ROI center. We freeze all spectral parameters for sources with TS < 1 or a predicted number of γ-rays after the initial optimization less than 10⁻³ counts. The normalizations of the diffuse backgrounds are left free during the fit,¹⁰ together with the spectral index of the Galactic diffuse background template. After the fit has converged successfully, we relocalize the central γ-ray source and refit all spectral model parameters. The relocalized source positions are provided in Table 1. After this step, we generate a TS map to search for additional point sources. For each pixel in the ROI, we add a putative point source with a power-law spectrum with index Γ = 2 and calculate its TS. If $\sqrt{\mathrm{TS}}\geqslant 5$, i.e., a detection with a significance of just over ∼4σ (Acero et al. 2015), we permanently add the source at the position of the highest TS value and reoptimize the spectral parameters for the whole ROI. This step is repeated until no further sources are found.

With the best-fit model for each ROI, we compute the γ-ray light curves for the FSRQs with an initial binning of 7 days. In each light-curve bin, we leave the spectral parameters free during the fit for sources within 3° of the ROI center and, additionally, the normalizations of the Galactic and isotropic emission. If any of these sources have TS < 1 or the number of predicted photons is less than 10⁻³, their parameters are fixed to the average values.

The 9.5 yr light curves for all considered FSRQs are shown in Figure 1. If the source is detected with TS < 9 within one time bin, or the flux in one bin is equal to or smaller than its statistical uncertainty F_i ≤ σ_i, we show upper limits¹¹ at the 2σ confidence level instead.

Zoom In Zoom Out Reset image size

The average source fluxes with their 1σ statistical uncertainties, $\overline{F}\pm {\sigma }_{\overline{F}}$, derived from the likelihood maximization over the full 9.5 yr are shown as gray bands. These flux measurements and uncertainties determine the optimal step-function representations of the light curves using the Bayesian block (BB) algorithm (Scargle et al. 2013) maximizing the overall fitness function appropriate to point measurements. From all possible partitions of the data into blocks, this algorithm finds the unique one maximizing the total fitness of the resulting step-function model.

These blocks provide an objective way to detect significant local variations in the light curve. Several strong flares exceeding the average flux level are easily identified from this block representation.

There is no generally accepted consensus on the best way to determine which data points belong to a flaring state and which characterize the quiescent level. Nalewajko (2013) suggested the simple definition that a flare is a continuous time interval associated with a flux peak in which the flux is larger than half the peak flux value. This definition is intuitive, however, and it is unclear how to treat overlapping flares and identify flux peaks in an objective way. Here we use a simple two-step procedure tailored to the block representation: (1) identify a block that is higher than both the previous and subsequent blocks as a peak, and (2) proceed downward from the peak in both directions as long as the blocks are successively lower. The two resulting monotonically decreasing sets of adjacent blocks are analogous to the watershed concept of topological data analysis. In fact, this approach was suggested by the HOP algorithm (Eisenstein & Hut 1998), which is based on a bottom-up hill-climbing concept of great use in higher dimensions (e.g., Way et al. 2011).

The time-series segments shown in Figure 1 are the result of feeding the block representations of the light curves to the HOP algorithm. The combination of BB and the HOP algorithm provides an objective way to split a light curve into groups of quiescent and flaring episodes; we will refer to one connected flare episode as a HOP group of consecutive BBs.

We iteratively zoom in on time ranges with bright γ-ray activity by identifying HOP groups where the peak BB fulfills the condition ${F}_{\mathrm{BB}}\geqslant {F}_{\max }=5\times \overline{F}$ and include adjacent blocks within the group with ${F}_{\mathrm{BB}}\geqslant {F}_{\min }=\overline{F}$. This relatively conservative definition gives reliable group shape information at the cost of slightly underestimating the extent of the groups and overestimates that of the quiescent intervals. Furthermore, we prefer a criterion based on peak flux rather than, for instance, integrated flare luminosity. This is because our approach promises the most straightforward way to find those time ranges with sufficient photon statistics to search for short-timescale flux variations and exponential cutoffs due to γ-ray absorption. If multiple adjacent HOP groups fullfill our criteria, they are combined into one time interval. The final time ranges are then extended by one time bin on either side. For the identified time span, we reoptimize the spectral model of the ROI in the same way as described in Section 2.2 but without relocalizing the central FSRQ or adding new point sources. The results of the best-fit spectra for the different time ranges are provided in the Appendix. Subsequently, we calculate a light curve with finer binning and again select the time ranges of the highest γ-ray activity. We repeat this procedure twice, down to a binning equal to the good time intervals (GTIs) of the Fermi satellite, which correspond to one passage of the source through the field of view of the satellite during one ∼95 minute orbit. The choices of time binnings and values for F_max and F_min are summarized in Table 2, together with the number of identified high γ-ray activity states (which might consist of several flares, as indicated by the HOP groups).

Table 2.  Thresholds for BB Fluxes in One HOP Group to Select Time Ranges of High γ-Ray Activity, Together with Selected Time Binning and Number of Selected Time Ranges

Binning	Fmin	Fmax	${N}_{\mathrm{time}\mathrm{ranges}}$
7 days	$\overline{F}$	$5\times \overline{F}$	20
1 day	$\overline{F}$	$\max ({10}^{-5}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1},1.5\times \overline{F})$ a	21
GTI	$\overline{F}$	$2\times \overline{F},\,\mathrm{TS}\geqslant 150$ b	7

Notes. The criteria are applied to all sources. Furthermore, if no interval fulfills the F_BB ≥ F_max criterion for the weekly or daily light curves, we include the HOP group of the maximum flux if that flux exceeds 5 × 10⁻⁶ cm⁻² s⁻¹; i.e., we change F_max to the maximum value F_BB.

^aWe choose here the absolute flux (rather than the flux relative to the average) as a threshold in order to be consistent with our initial source selection. However, because of the high average flux of 3C 454.3, we also include the max argument. If ${F}_{\max }=1.5\times \overline{F}$, we set F_min = 10⁻⁵ cm⁻² s⁻¹. ^bMotivated by the high TS found for the flare of 3C 279  (Ackermann et al. 2016b), we also demand that at least one GTI of each HOP group be detected with TS ≥ 150 in order to ensure enough statistics to search for variability on timescales of minutes.

Download table as:  ASCIITypeset image

The values of the threshold fluxes F_max and F_min are somewhat arbitrary and are a compromise between including as many flares as possible and keeping the overall number of flares manageable. Note that the sole purpose of this exercise is to select the brightest flares for further analysis; consideration of a more complete statistical sample is postponed to the future.

The time ranges of the identified flares for the weekly light curves are plotted as colored horizontal lines in Figure 1. The intermediate daily light curves are provided in Figure 2. Figure 5 shows the GTI (equivalent to orbital) light curves with exponential fits of flare profiles that we will discuss in Section 4. The source exposure can vary significantly between two adjacent orbits, as the satellite rocks between the celestial north and south poles between orbits. This explains the large error bars on some of the time bins of the GTI light curves.

Zoom In Zoom Out Reset image size

In a last step, we derive light curves on sub-GTI timescales. The time bin size is calculated from the adaptive binning method of Lott et al. (2012), where we choose bins of constant flux uncertainty of ∼20%. In this step, we use the spacecraft information in time steps of 1 s instead of 30 s. Additionally, we compute the effective area in five bins of the azimuthal spacecraft coordinates because, on such short timescales, the exposure dependence on the azimuth should not be averaged over.¹² We discuss these light curves in Section 4.2.

Light curves and fit results for the different time intervals and binnings are provided online.¹³

As mentioned in the Introduction, there is no rigorous way to distinguish the flux in flares from the flux in a quiescent background (QB). It is not even guaranteed that there is a corresponding astrophysical distinction. Nevertheless, some progress can be made, especially given the assumption that the QB is truly constant.

The minimum flux observed over time might be taken as an estimate of the true QB, but it is a rather crude lower limit, biased downward because of the scatter due to observational errors. On the other hand, long overlapping tails of flares could contribute an approximately constant flux level yielding an upward bias to some QB estimates. For the present, we assume that the overlap of flare tails is negligible and propose a statistical procedure addressing the observational error bias.

This algorithm is based on an approximate separation of the distribution function of the observed fluxes into two components, (a) a low end and (b) a high end, dominated by the QB and the flares, respectively. While not relying on (b) being devoid of any QB flux, we do assume that (a) contains almost entirely QB flux. In the picture outlined above, this amounts to the assumption that the flares are isolated (no overlap) and the intervening intervals are essentially pure QB. The algorithm implements this separation using only flux values with no regard to their time sequence. It is based on finding the flux interval that maximizes the symmetry of the resulting distribution.

In the following pseudo-code, the term "distribution" refers to the distribution of low values (a) in a putative flux interval defining the QB.

• 1.  
  Sort the N flux values (with TS ≥ 9) in increasing order: ${F}_{1}\leqslant {F}_{2}\leqslant ...{F}_{N-1}\leqslant {F}_{N}$.
• 2.  
  Define a search grid of flux values ϕ_m, m = 1, ..., M to serve as candidates for the peak of the distribution. We take these values on an evenly spaced grid between the minimum and mean flux $\bar{F}$: ${\phi }_{m}\in [{F}_{1}+\epsilon (\bar{F}-{F}_{1});\bar{F}]$. By definition, the QB is very unlikely to be less than the minimum flux or exceed the mean. The factor is a small number to avoid edge effects near the generally sparse lower end of the flux distribution.
• 3.  
  For each m, define two intervals of equal length straddling ϕ_m:

  • (i)  
    LOW: from F₁ to ϕ_m; and
  • (ii)  
    HIGH: from ϕ_m to ϕ_m + (ϕ_m − F₁).
• 4.  
  Construct fine grids of flux values spanning these two intervals.
• 5.  
  Compute the cumulative distributions (CDFs) of the corresponding flux values.
• 6.  
  Reverse the CDF for the HIGH interval.
• 7.  
  Normalize both CDFs to unity at the peak ϕ_m.
• 8.  
  Compute the ratio of the posterior probabilities that the CDFs come from different distributions or the same distribution, using the algorithm of Wolpert (1996).
• 9.  
  Find the value of ϕ_m that minimizes this ratio, which measures the asymmetry of the total CDF of HIGH and LOW.
• 10.  
  Report the median flux, F_QB, in this optimally symmetric distribution as the QB estimate.

In the limit of the moderately finely gridded bins used in step 5, the CDF estimate is effectively unbinned (little or no dependence on the binning).

We show our estimate for the QB level in Figure 1 as black dashed lines and report the values in Table 3. In general, these F_QB values match the visual flux baselines extremely well. In the case of 3C 454.3, it is slightly higher than the minimum flux level observed for this source around MJD 55,800–56,200. The reason is our applied TS cut and a contamination of F_QB from the tails of the flares. We have tested the latter point with simulations drawing random numbers from a uniform distribution and Cauchy distributions to emulate flares. For the Cauchy distributions, we randomize the height, width, and position. Applying our algorithm to these simulated histograms, we find that F_QB slightly overestimates the true uniform background. Therefore, we conclude that F_QB can be seen as a firm upper limit on a truly constant QB level. We also note that using the peak or mean of the distribution in step 10, instead of the median, only changes the results marginally. We have furthermore tested different metrics instead of the algorithm of Wolpert (1996), namely, the Kolmogorov–Smirnov (K-S) test and the least squares between the CDFs. We again find comparable results. However, the K-S test estimates underpredict the true QB level in simulations, while the least squares give estimates that are too high (higher than the Wolpert estimate). We also note that the F_QB values are either consistent with or lower than the break flux of the BPL fit, F_br. This is expected, as F_QB marks the median of the QB, while F_br probably indicates the transition from the QB to the flaring states.

We further characterize the global γ-ray light curves in terms of their PSD, which usually can be described with simple power laws in frequency, PSD ∝ 1/ν^β. An analysis by Abdo et al. (2010b) of the first 11 months of LAT data from 106 blazars revealed that these objects have β values between 1 and 2, the intermediate regime between flicker noise (β = 1) and Brownian motion (β = 2). In addition to the noise behavior, we will use the derived PSDs in Section 5.3 to simulate γ-ray light curves in order to calculate the significance of a correlation between radio and γ-ray emission.

The best-fit PSDs are estimated from the periodograms and simulated light curves following the method described in detail in Emmanoulopoulos et al. (2013) and Max-Moerbeck et al. (2014b) and summarized briefly below. The observed periodograms P(ν) as a function of frequency (inverse time) ν are calculated from the absolute square of the Fourier transformation of the light curve (Equation (3) in Max-Moerbeck et al. 2014b). We include all data points detected with TS ≥ 9 and perform a linear interpolation between gaps in the light curve to guarantee an even sampling. Since we are using weekly binned light curves and bright FSRQs, the gaps are small and at most six consecutive data points long (42 days) in the case of PKS B1222+216. The number of nondetected bins is less than ∼13% for all sources. In contrast to Max-Moerbeck et al. (2014b), we do not apply a window function (see the discussion below).

We compare the observed periodogram with simulated light curves, which we produce with the method of Timmer & Koenig (1995) for power-law PSDs with values 0 ≤ β ≤ 3 in steps of Δβ = 0.05. For each β value, 100 light curves are generated, each one 100 times longer than the actual observation to account for possible red-noise leakage. Splitting the simulated light curves (without overlap) leaves us with 10⁴ realizations. The light curves are initially produced with a regular time binning equal to 0.7 days and rebinned into 7 day light curves through averaging. The same observational gaps and interpolation as in the observed light curves are applied. The periodograms are then calculated for the simulated light curve in the same way as for the observed one.

To fix the normalization of the PSD model, Max-Moerbeck et al. (2014b) suggested variance matching; i.e., they rescaled the simulated flux data points with a factor A⁻¹, where ${A}^{2}={\sigma }_{\mathrm{sim}}^{2}/({\sigma }_{\mathrm{obs}}^{2}-\bar{{\sigma }_{i}^{2}})$, with ${\sigma }_{\mathrm{sim}}^{2}$ (${\sigma }_{\mathrm{data}}^{2}$) the variance of the simulated (observed) light curve and $\bar{{\sigma }_{i}^{2}}$ the variance of the observational noise. For the γ-ray light curves, we choose to follow Emmanoulopoulos et al. (2013) instead and iteratively match the probability distribution of the simulated fluxes to the observed ones, given by the dN/dF distributions shown in Figure 3. The reason is that the algorithm of Timmer & Koenig (1995), which is used by Max-Moerbeck et al. (2014b), produces light curves with Gaussian-distributed fluxes, which is clearly not the case at γ-ray energies.¹⁵ In a final step, we add uncertainties to the light curves by randomly drawing with replacements from the observed uncertainties σ_i and adding a Gaussian random number ${ \mathcal N }(0,{\sigma }_{i})$ to the simulated flux values.

The peridograms of the observed and simulated light curves, P_obs and P_sim, are averaged in 31 logarithmic bins¹⁶ (Papadakis & Lawrence 1993) between ν_min = 1/T (where T is the full duration of the light curve with N measurements) and ν_max = N/(2T) (the Nyquist frequency) and compared by means of a χ² test (Max-Moerbeck et al. 2014b),

$$\begin{eqnarray}&&{\chi }^{2}(\beta )=\displaystyle \sum _{{\nu }_{\min }}^{{\nu }_{\max }}\displaystyle \frac{{\left({P}_{\mathrm{obs}}(\nu )-{\overline{P}}_{\mathrm{sim}}(\nu ,\beta )\right)}^{2}}{{\rm{\Delta }}{\overline{P}}_{\mathrm{sim}}{\left(\nu ,\beta \right)}^{2}}.\end{eqnarray} \tag{ 2 }$$

Here ${\overline{P}}_{\mathrm{sim}}(\nu )$ and ${\rm{\Delta }}{\overline{P}}_{\mathrm{sim}}{(\nu )}^{2}$ are the mean and variance of the periodograms obtained from the simulated light curves. The averaged periodograms of the simulated light curves and the observed ones are shown in Figure 4, and the best-fit average periodogram is shown as a thick solid line. The quality of the best-fit value $\hat{\beta }$ with a corresponding minimum χ² value ${\hat{\chi }}^{2}\equiv {\chi }^{2}(\hat{\beta })$ is evaluated from the light curves simulated with $\beta =\hat{\beta }$ in the following way. We form the distributions of simulated χ² values, ${\chi }_{\mathrm{sim}}^{2}$, by replacing P_obs(ν) with P_sim(ν, β) in Equation (2),

$$\begin{eqnarray}&&{\chi }_{\mathrm{sim}}^{2}(\beta ,\beta ^{\prime} )=\displaystyle \sum _{{\nu }_{\min }}^{{\nu }_{\max }}\displaystyle \frac{{\left({P}_{\mathrm{sim}}(\nu ,\beta )-{\overline{P}}_{\mathrm{sim}}(\nu ,\beta ^{\prime} )\right)}^{2}}{{\rm{\Delta }}{\overline{P}}_{\mathrm{sim}}{\left(\nu ,\beta ^{\prime} \right)}^{2}},\end{eqnarray} \tag{ 3 }$$

and calculate the p-value, p_β, as the fraction of simulations that result in ${\chi }_{\mathrm{sim}}^{2}(\hat{\beta },\hat{\beta })\gt {\hat{\chi }}^{2}$.¹⁷ The confidence interval for $\hat{\beta }$ is derived by determining the ${\rm{\Delta }}{\chi }_{\mathrm{sim}}^{2}(\hat{\beta },\beta )$ value from simulations such that 95% of the time, the simulated (true) β value is contained within Δχ².¹⁸ The same Δχ² value is then applied to the observed χ² curve. The results of our PSD analysis are summarized in Table 3, where we also report the value of β obtained from a linear regression in log–log space. In general, the periodograms are well fit by our method, as indicated by the p_β-values and observed in Figure 4. The only exception is 3C 279, where only two of the 10⁴ simulated light curves result in a ${\chi }_{\mathrm{sim}}^{2}(\hat{\beta },\hat{\beta })\gt {\hat{\chi }}^{2}$. The steep χ² curve for this source also explains the small error bars on the reconstructed value of β. The reason might be the variation of the underlying PSD with time, as found by Ackermann et al. (2016b). Another possibility is the specific 7 day binning we have chosen here.

Zoom In Zoom Out Reset image size

Our results are compatible with the PSD slopes found by Nakagawa & Mori (2013) at the 1σ–2σ level but less so for the PSD power laws obtained by Sobolewska et al. (2014). The reason for the discrepancy with the latter analysis might be the different binning schemes and time ranges (Sobolewska et al. 2014 used an adaptive binning for 4 yr of LAT data instead of the constant binning and 9.5 yr used here).

We again use BBs and HOP groups to identify different states in the orbital light curves (see Figure 5). To assess the time profile of the flares, we fit each HOP group i with a sum of exponential profiles, F_flare,i, using a χ² minimization, and

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{flare},i}(t) & = & \displaystyle \sum _{j=1}^{{N}_{i}}{F}_{0,{ij}}\\ & & \times \,{\left[\exp \left(\displaystyle \frac{t-{t}_{0,{ij}}}{{\tau }_{\mathrm{rise},{ij}}}\right)+\exp \left(\displaystyle \frac{{t}_{0,{ij}}-t}{{\tau }_{\mathrm{decay},{ij}}}\right)\right]}^{-1},\end{array}\end{eqnarray} \tag{ 4 }$$

where t_0,ij are the times when the flare flux is equal to F_0,ij/2, and τ_rise,ij and τ_decay,ij are the flare rise and decay times, respectively. All light-curve points are included that fulfill TS ≥ 9 and F_i ≥ 3σ_i/2. The number of flare profiles per HOP group, N_i, is either one or two and determined during the fit using the Bayesian information criterion (BIC), defined as $\mathrm{BIC}={n}_{\mathrm{par}}\mathrm{ln}(n)+{\chi }^{2}$, where n_par is the number of fit parameters (n_par = 4 for N_i = 1), and n is the number of data points within one HOP group i. Two flare profiles are selected if the difference between the two BIC values is ΔBIC = BIC(N_i = 2) − BIC(N_i = 1) < 0. The reason for allowing N_i > 1 is that the flare profile in Equation (4) does not capture the long-lasting plateaus of a flare (see, e.g., all flares of 3C 279 or the last panel with the flare of 3C 454.3 starting at 55,551.65 MJD in Figure 5).

After each HOP group is fitted individually and N_i is determined, we refit the entire light curve, which consists of N_HOP groups, with the function

$$\begin{eqnarray}&&{F}_{\mathrm{flare}}(t)=\displaystyle \sum _{i=1}^{{N}_{\mathrm{HOP}}}{F}_{\mathrm{flare},i}(t)+{F}_{\mathrm{bkg}}(t),\end{eqnarray} \tag{ 5 }$$

where F_bkg(t) is an order 2 polynomial to describe a slow varying background. The fit results are shown as black solid lines in Figure 5. In general, the χ² values divided by the degrees of freedom (dof) are between 1 and 2 (see the legends in Figure 5). Given the large values of dof, the fit qualities are generally poor. This is not unexpected, as we only allow up to two flare profiles per HOP group and no arbitrary functions. Already with this choice, there are probably some spurious flares identified, see, e.g., the second flare profile in the first PKS 1510–089 flare (starting at 54,908.65 MJD). Nevertheless, the overall light-curve evolution is well captured, which allows us to describe the local flare properties from the ensembles of flare profiles.

We show the distribution of rise and decay times in Figure 6 for flares with a time-integrated flux >10⁻⁷ cm⁻². Remarkably, all sources show values of τ_rise and τ_decay that are close to or below the horizon-crossing timescale of the central supermassive black hole, t_g = r_g/c, where r_g = 2GM_•/c² is the gravitational radius, G is the gravitational constant, M_• is black hole mass (taken from Table 1), and c is the speed of light. This results in values between ∼0.5 and ∼2.5 hr using the black hole masses listed in Table 1.

Zoom In Zoom Out Reset image size

From the rise and decay times, we can calculate the flare asymmetry as

$$\begin{eqnarray}&&A=\displaystyle \frac{{\tau }_{\mathrm{rise}}-{\tau }_{\mathrm{decay}}}{{\tau }_{\mathrm{rise}}+{\tau }_{\mathrm{decay}}},\end{eqnarray} \tag{ 6 }$$

so that A < 0 for fast-rise exponential-decay (FRED) type flares, as expected from an injection of energetic particles on timescales faster than the subsequent cooling through radiative processes such as inverse Compton (IC) scattering or synchrotron emission. The asymmetry is shown versus integrated flux, peak flux, and flare duration T₉₀, defined as the time around the flare peak that contains 90% of the integrated flux, in Figure 7. The peak flux for each flare of each HOP group is derived from the maximum of Equation (4) with respect to time (suppressing indices ij):

$$\begin{eqnarray}&&{F}_{\mathrm{peak}}=\displaystyle \frac{{F}_{0}{\tau }_{\mathrm{rise}}}{{\tau }_{\mathrm{rise}}+{\tau }_{\mathrm{decay}}}{\left(\displaystyle \frac{{\tau }_{\mathrm{decay}}}{{\tau }_{\mathrm{rise}}}\right)}^{\tfrac{{\tau }_{\mathrm{decay}}}{{\tau }_{\mathrm{rise}}+{\tau }_{\mathrm{decay}}}}.\end{eqnarray} \tag{ 7 }$$

The error bars on the peak flux and asymmetry are derived from standard Gaussian error propagation from the fit uncertainties.

Zoom In Zoom Out Reset image size

The median of the asymmetry is found to be −0.195; i.e., FRED-type flares are more common than the opposite. In general, the flares show versatile behavior, and no clear trends are seen from Figure 7. This is also reflected in the fact that we do not find any significant correlation between the asymmetry, integrated and peak flux, rise and decay times, and flare duration using Kendall's τ.

We also investigate whether subsequent flares in each panel of Figure 5 show a trend with time in peak flux, asymmetry, or duration. For consecutive flares, we calculate the difference between, e.g., the peak fluxes, and we calculate the p-value of a binomial distribution assuming an equal probability of finding negative and positive differences. For 32 values of differences, the p-values for the peak flux, asymmetry, and flare duration are close to 0.1 (14, 13, and 13 positive values for 32 trials, respectively), indicating no particular evolution of these quantities with time. As a systematic check, we have repeated the entire profile fit for time-reversed versions of the light curves. In general, we find good agreement between the fitted profiles. Correcting for the sign reversal, the median of the asymmetry changes from −0.195 to −0.198, suggesting that the systematic error is of the order of 2%.

We also find complex behavior of the spectral evolution during the flares. Evidence for a "harder-when-brighter" trend is found for some sources and flares, which is, however, not significant. We therefore cannot draw any firm conclusions from the spectral evolution, which we show in Figure 8.

Zoom In Zoom Out Reset image size

We search for suborbital variability in a subset of orbital light curves, namely, those where at least one orbital bin is detected with TS ≥ 150. In this way, we ensure high photon statistics and reduce the number of trials when searching for minute-scale variability (for comparison, the orbital light-curve bin for which Ackermann et al. 2016b measured the minute-scale variability in 3C 279 is detected in our analysis with TS ∼ 400). The selected time regions are indicated with solid horizontal lines in Figure 5, whereas the dashed lines show the time intervals selected with the criteria in Table 2 that do not pass the additional TS cut.

The resulting light curves, binned such that the uncertainty in each bin is of the order of ∼20% (using the adaptive binning introduced by Lott et al. 2012), are shown in Figure 9. In order to make an objective selection of GTIs that we test against the hypothesis of a constant flux, we consider only those where the BBs indicate a significant flux change within the GTI. Naively, one could expect that a BB change within one GTI would correspond to a significance of 95% for a nonconstant flux, since this is the threshold we have selected in the BB algorithm (Scargle et al. 2013). However, the BB algorithm also takes data before and after the particular GTI into account and only provides qualitative evidence for minute-scale variability. Therefore, we test each bin selected in this way against the hypothesis of constant flux using a simple χ² test. The best-fit constant flux is given by $\hat{F}={(\sum ({F}_{i}/{\sigma }_{i}^{2}))(\sum {F}_{i}^{-2})}^{-1}$. For the GTIs where the constant fit results in a pretrial p-value of less than 0.1, we show the sub-GTI light curves in Figure 10 and report the pre- and posttrial fit probabilities in Table 4. We count each tested GTI as one trial. We also provide the minimum values for the variability times for pairs of fluxes F_i and F_j measured at t_i and t_j, respectively, given by Zhang et al. (1999),

$$\begin{eqnarray}&&{t}_{\mathrm{var},{ij}}=\displaystyle \frac{{F}_{i}+{F}_{j}}{2}\left|\displaystyle \frac{{t}_{i}-{t}_{j}}{{F}_{i}-{F}_{j}}\right|.\end{eqnarray} \tag{ 8 }$$

The pretrial p-values for rejecting the constant flux hypothesis are all around 2σ. The trial correction leaves only one GTI for 3C 279 and two GTIs for CTA 102 close to or above the 2σ evidence for a variable flux. However, inspecting the light curve of the second GTI of CTA 102 (starting at MJD 57,758.86; see Figure 10), the high χ² value might be the result of our chosen binning; the first bin actually spans a long time range for which the exposure is mostly zero. For the other two GTIs, the suggested variability timescales are between 3 and 4 minutes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In comparison to previous results for 3C 279 and CTA 102, our method results in lower-significance detections of minute-scale variability. Furthermore, Shukla et al. (2018) found evidence for minute-scale variability in a different orbit compared to our results. This is due to trial correction but probably also to differences in the analysis. For example, we use a finer binning for the exposure in the azimuthal direction to take this dependence for such short observations into account. The change in exposure is, however, below 10%. More importantly, we use a different binning within one GTI, which can also change the significance. Taking the χ² test and the BBs together, we conclude that we find evidence that minute-scale variability is a common phenomenon during bright FSRQ flares.

5.1.1. γ-Ray Attenuation

The attenuation due to pair production on a radiation field of soft photons should manifest itself as a cutoff feature in the γ-ray spectrum. The cutoff energy depends on the distance of the γ-ray emitting region to the central black hole, r, and the photon density of the considered photon field. For FSRQs, the photon densities of the external radiation fields of the accretion disk, BLR, and extended dust torus usually dominate those of internal synchrotron emission (see, e.g., Dotson et al. 2012). The most relevant external photon field for the γ-ray energies that can be probed with Fermi LAT is the BLR. Pair production on photons from the dust torus or the accretion disk only becomes important at energies beyond 1 TeV, even when the γ-rays are produced close to the central black hole (Finke 2016).

We can search for BLR absorption features by fitting the observed γ-ray spectra during the brightest flares with functions of the form

$$\begin{eqnarray}\begin{array}{rcl}f(E,{\boldsymbol{\pi }},r,z) & = & {f}_{\mathrm{int}}(E,{\boldsymbol{\pi }})\\ & & \times \,\exp \left[-\left({\tau }_{\gamma \gamma }^{\mathrm{BLR}}(E,r)+{\tau }_{\gamma \gamma }^{\mathrm{EBL}}(E,z)\right)\right],\end{array}\end{eqnarray} \tag{ 9 }$$

where ${f}_{\mathrm{int}}(E,{\boldsymbol{\pi }})$ describes the intrinsic spectrum at observed γ-ray energy E emitted by the source, which depends on spectral source parameters, ${\boldsymbol{\pi }}$, such as, e.g., flux normalization, power-law index, and spectral curvature (see also the Appendix for the definitions of the spectral models), and ${\tau }_{\gamma \gamma }^{\mathrm{BLR}/\mathrm{EBL}}$ is the optical depth due to interactions of γ-rays with photons of the BLR and extragalactic background light (EBL), respectively. For the EBL optical depth, which depends on E and the source redshift, z, we use the EBL model of Domínguez et al. (2011).¹⁹ The BLR optical depth is described by the stratified BLR model introduced by Finke (2016), who modeled the BLR as a collection of either shells or rings perpendicular to the jet axis in order to emulate a flattened BLR. Each shell or ring is assumed to have infinitesimal thickness and emit a monochromatic UV or optical emission line. The radii R_li of the shells and rings, as well as the line luminosities L_li = ξ_liL_disk, are taken from templates of average spectra obtained in reverberation mapping campaigns and provide values relative to the radius and luminosity of the Hβ line (see Finke 2016 for further details). With the Hβ luminosities listed in Table 1, we fix the absolute luminosities (or, conversely, ξ_Hβ) and radii of all lines included in the model. Together with the masses of the supermassive black holes, we can then calculate ${\tau }_{\gamma \gamma }^{\mathrm{BLR}}$ for both geometries as a function of r and observed γ-ray energy E. In the BLR model, the absorption is dominated by pair production with Lyα photons at a rest-frame energy of _Lyα ∼ 10.2 eV emitted at radii between ∼8 × 10¹⁶ and 2 × 10¹⁷ cm. In the ring geometry, the corresponding energy density, which is assumed to be isotropic in the stationary frame of the galaxy, becomes (Finke 2016)

$$\begin{eqnarray}&&{u}_{\mathrm{BLR}}\approx {u}_{\mathrm{Ly}\alpha }=\displaystyle \frac{{\xi }_{\mathrm{Ly}\alpha }{L}_{\mathrm{disk}}}{4\pi c({R}_{\mathrm{Ly}\alpha }^{2}+{r}^{2})},\end{eqnarray} \tag{ 10 }$$

and takes values of ∼5 × 10⁻² erg cm⁻³ regardless of the source for r = R_Lyα. These numbers can be compared against typical values for the BLR radius, R_BLR ∼ 10¹⁷ cm (L_disk/10⁴⁵ erg s⁻¹)^1/2 (e.g., Kaspi et al. 2007; Bentz et al. 2009), and energy density u_BLR ∼ 10⁻² erg cm⁻³ (again in the stationary galaxy frame) assuming L_BLR = ξ_BLRL_disk with ξ_BLR ∼ 0.1. The chosen BLR model gives values broadly consistent with typical values within a factor of a few.

Typically, FSRQ spectra show intrinsic curvature, even below energies at which BLR absorption becomes important (see, e.g., the 3FGL). Therefore, we chose a log-parabola for the intrinsic spectral function and also test a power law with a superexponential cutoff (see Equations (19) and (20) for the definition of the models). In the fit, we only include energy bins above 1 GeV, as we expect the BLR cutoff at energies ≳10 GeV. In this way, we avoid having the best fit determined mainly by the high photon statistics at lower energies. Additionally, we select narrow time intervals around the brightest flares (see Section 2.3 and Figure 5). This is a compromise between sufficient photon statistics to probe energies above 10 GeV and avoiding the mixing of different activity states with potentially different spectral states. From Figure 8, we see that for the time bins with the highest fluxes, spectral variability is only marginally present, which should render our results robust against potential variations of the intrinsic spectra. Also, in the fit, we only include energy bins detected with TS > 0 and skip flares where the absorption is below 80% in the highest energy bin with TS > 0 and for the smallest BLR distance tested (r = 10⁻²R_Lyα). This excludes all flaring periods from 3C 273, for which we cannot obtain any limits from the γ-ray spectra.

We derive the best-fit values for the spectral parameters ${\boldsymbol{\pi }}$ and the distance r with a likelihood maximization of the bin-by-bin likelihood curves, which we extract with fermipy²⁰ and are shown as gray shaded bands in Figure 11. The flux points in the figure coincide with the maximum likelihood. Also shown are the best-fit spectra and BLR attenuation for different values of r (colored curves).

Zoom In Zoom Out Reset image size

For both tested BLR geometries, the best-fit value of r is always close to or coincides with the maximum tested value, $r=10{R}_{\mathrm{Ly}\alpha }$, and hence no significant absorption is found (dashed black lines in Figure 11). Consequently, we use Minos to derive the profile likelihood as a function of r from which we determine the 95% lower limit on r, r_lim (dashed–dotted black lines). The limit values are reported for each flare in Figure 12 and summarized in Table 5 for the ring BLR geometry and log-parabola spectrum. Assuming instead a power law with a superexponential cutoff yields consistent results. For the BLR shell geometry, the lower limits are a factor of ∼2–3 higher because this geometry predicts stronger absorption (Finke 2016). The ring geometry is therefore the conservative choice.

Zoom In Zoom Out Reset image size

Table 5.  Results from BLR Absorption Fits to γ-Ray Spectra

t0	Δt	rlim	rlim	rlim	EHEP	${E}_{{\tau }_{\gamma \gamma }=1}$	${t}_{\mathrm{cool},\mathrm{BLR}}$	${t}_{\mathrm{cool},\mathrm{dt}}$	τdecay
[MJD]	[days]	[1017 cm]	[RLyα]	[rg]	[GeV]	[GeV]	[minutes]	[hr]	[hr]
PKS B1222+216
55,364.68	3.42	1.33	1.40	609	75.39	69.69	8.2	2.3	47.4 ± 8.3
3C 279
57,188.07	1.87	0.49	0.64	867	56.03	42.91	2.7	19.0	0.5 ± 0.9
58,133.34	5.32	1.45	1.91	2580	92.56	107.91	9.0	19.0	8.2 ± 6.3
PKS 1510–089
57,114.16	1.42	0.51	0.66	1088	66.54	54.99	0.6	4.5	0.4 ± 0.3
57,243.84	4.53	0.74	0.97	1591	75.93	65.39	0.8	4.5	44.4 ± 9.4
CTA 102
57,737.41	1.67	1.41	0.86	562	36.25	21.23	1.0	1.5	0.3 ± 0.5
57,749.10	4.99	3.20	1.95	1275	73.80	37.94	2.8	1.5	8.7 ± 1.2
57,757.55	4.66	2.76	1.67	1096	39.19	32.38	2.2	1.5	24.6 ± 2.3
57,861.71	2.42	1.95	1.18	776	34.73	24.94	1.4	1.5	1.2 ± 0.7
3C 454.3
55,516.55	8.93	3.19	1.36	1598	41.19	28.73	4.2	16.8	2.6 ± 1.0

Note. The HEPs are given for source probabilities >0.99. The decay times are given for the flare component with the highest peak flux as determined in the fit to the orbital light curves in Section 4. For the cooling times, an observed γ-ray energy of 10^8.5 eV ≈ 316 MeV is assumed.

Download table as:  ASCIITypeset image

As can be seen from Table 5 and Figure 12, the limits are of the order of r_lim ∼ 10¹⁷ cm, which translates to a distance close to or even beyond the Lyα-emitting ring and, consequently, the BLR itself. In terms of gravitational radii, the emission regions are located at distances of at least ∼10³r_g. Table 5 also reports the energy of the highest-energy photon (HEP) associated with the FSRQ with at least 99% probability. For all but one source, this energy is larger than the energy where the optical depth due to absorption in the BLR exceeds ${\tau }_{\gamma \gamma }^{\mathrm{BLR}}\gt 1$ (assuming r = r_lim). Our limits generally agree with the results of Costamante et al. (2018), who could limit the maximum value of ${\tau }_{\gamma \gamma }^{\mathrm{BLR}}$ to ∼1 for 3C 454.3 and PKS B1222+216 and ∼0.2 for CTA 102.

The limits sensitively depend on the detection of the source at energies as high as possible. To demonstrate that the detections are not spurious, we also report the detection significance and the number of detected γ-rays (associated with the source with a probability >85%) for each energy bin below and above the flux points in Figure 11, respectively. The highest-energy bins only contain a handful of source photons (one to four), which underlines the necessity to use the full Poisson likelihood information. Doing so, the energy bins are indeed detected with significances of $\sim \sqrt{\mathrm{TS}}\gtrsim 4\sigma$. The reason is that in the considered energy interval and short time spans (see Table 5), the number of expected background events is small.

We also compare the limits from fits to γ-ray spectra to considerations from variability arguments in Figure 12. If the emission region ${R}_{\mathrm{blob}}^{{\prime} }$ (the prime denotes the comoving frame) is causally connected during the flare, the shortest variability time t_var sets a lower limit on its size (e.g., Begelman et al. 2008),

$$\begin{eqnarray}&&{R}_{\mathrm{blob}}^{{\prime} }\leqslant \displaystyle \frac{{{ct}}_{\mathrm{var}}{\delta }_{{\rm{D}}}}{1+z},\end{eqnarray} \tag{ 11 }$$

where ${\delta }_{{\rm{D}}}={{\rm{\Gamma }}}_{{\rm{L}}}^{-1}{(1-\beta \cos {\theta }_{\mathrm{obs}})}^{-1}$ is the Doppler boost factor with the bulk Lorentz factor of the flow, Γ_L; $\beta =\sqrt{1-{{\rm{\Gamma }}}_{{\rm{L}}}^{-2}}$ is the associated velocity; and θ_obs is the angle between the line of sight and the jet axis. Clausen-Brown et al. (2013) found a correlation, ${\theta }_{{\rm{j}}}\sim \rho {{\rm{\Gamma }}}_{{\rm{L}}}^{-1}$, with ρ ∼ 0.2 from the data of the MOJAVE very long baseline interferometry (VLBI) blazar monitoring program. Similar values are also found from Very Long Baseline Array (VLBA) monitoring observations (Jorstad et al. 2017). Using this correlation and under the assumption that the plasma blob occupies half of the jet's cross section, we find ${\theta }_{{\rm{j}}}\sim 0.1{{\rm{\Gamma }}}_{{\rm{L}}}^{-1}$ and obtain an upper limit on the distance to the black hole, $r\sim 10{R}_{\mathrm{blob}}^{{\prime} }{{\rm{\Gamma }}}_{{\rm{L}}}\equiv {r}_{{\rm{j}}}$. The values are plotted in Figure 12 for the minimum of the rise and decay times of the brightest flares found in Figure 5. The average values for δ_D and Γ_L obtained from VLBA observations are used (Jorstad et al. 2017; see also Table 1), and the total uncertainty is obtained by summing the uncertainties on δ_D and Γ_L and the fit uncertainty of t_var in quadrature. It should be noted that the underlying assumption is that the γ-rays are produced cospatially with the radio emission for which the Doppler factors are measured. In general, we find that r_lim ≳ r_j, indicating that the emission regions are at larger distances to the black hole than predicted from the conical jet scenario.

In our BLR model, we use a simplified BLR geometry and do not include the hydrogen or He ii recombination continua or emission lines of the He ii Ly series (as done in, e.g., Poutanen & Stern 2010; Stern & Poutanen 2014). Comparing the optical depths of our model to the results of the sophisticated modeling of Abolmasov & Poutanen (2017; see, in particular, their Figure 11), who described the BLR as a collection of ionized gas clouds irradiated by the accretion disk, we find that our values of ${\tau }_{\gamma \gamma }^{\mathrm{BLR}}$ reach unity at energies a factor of ∼1.5 higher, but around 100 GeV, the optical depth in the two models is similar. Also, the BLR model of Finke (2016) reproduces the trend observed in the sophisticated model that the absorption sets in at higher energies and is overall weaker for larger values of r. As we do not observe significant cutoffs in the spectra, we therefore conclude that the ring geometry adopted here provides a conservative limit on r.

Furthermore, evidence exists that the the luminosity of BLR emission lines is variable in 3C 454.3, PKS 1510–089, and PKS B1222+216 (León-Tavares et al. 2013; Isler et al. 2015) and correlates with the γ-ray emission (León-Tavares et al. 2013, 2015), which could indicate that γ-rays are produced through IC scattering with BLR photons (see also Section 5.2). Additionally, León-Tavares et al. (2013) found that the BLR brightening coincides with the passage of a superluminal jet component through the radio core, which could indicate that BLR clouds are located at larger distances, ≳1 pc, than assumed here. A brighter BLR emission during γ-ray flares and BLR material located at larger distances would mean stronger γ-ray attenuation, which would shift our limits to even larger distances.

5.1.2. Jet Shielding by a Plasma Sheath

5.2.1. Broad Emission Line Radiation

With the limits on r, it is possible to derive constraints on the energy density of external photon fields in the comoving frame, which could be responsible for the γ-ray emission due to IC scattering with relativistic electrons in the emission region. Because the cooling time depends on the energy density and, in turn, on r, a comparison between the predicted IC cooling times and the observed decay times can provide further information on where the γ-rays are emitted. We first focus on the BLR photon field, but the discussion also applies for IC scattering with photons of the dust torus, which we will discuss at the end of this section.

In the galaxy frame, the energy density of the BLR in the ring geometry is approximately given by Equation (10); hence, the photon number density is n_BLR ≈ u_BLR/_Lyα ∼ 10⁹ cm⁻³. Assuming that the BLR photon field is isotropic and in the limit Γ_L ≫ 1, the energy density in the comoving frame becomes ${u}_{\mathrm{BLR}}^{{\prime} }=(4/3){{\rm{\Gamma }}}_{{\rm{L}}}^{2}{u}_{\mathrm{BLR}}$ (Dermer & Schlickeiser 1994, 2002). We calculate the energy loss of the electrons, ${\dot{\gamma }}_{\mathrm{BLR}}^{{\prime} }$, due to IC scattering in the comoving frame numerically, in order to incorporate Klein–Nishina effects following Blumenthal & Gould (1970). The observed cooling time is then given by

$$\begin{eqnarray}&&{t}_{\mathrm{cool},\mathrm{BLR}}=\displaystyle \frac{1+z}{{\delta }_{{\rm{D}}}}\displaystyle \frac{\gamma ^{\prime} }{{\dot{\gamma }}_{\mathrm{BLR}}^{{\prime} }}.\end{eqnarray} \tag{ 12 }$$

In the Thomson regime, this becomes

$$\begin{eqnarray}&&{t}_{\mathrm{cool},\mathrm{BLR}}=\displaystyle \frac{1+z}{{\delta }_{{\rm{D}}}}\displaystyle \frac{3{m}_{e}{c}^{2}}{4c{\sigma }_{{\rm{T}}}{u}_{\mathrm{BLR}}^{{\prime} }{\gamma }_{\mathrm{BLR}}^{{\prime} }},\end{eqnarray} \tag{ 13 }$$

where m_e is the electron mass and σ_T is the Thomson cross section. In what follows, we approximate the electron Lorentz factor with (e.g., Dermer & Menon 2009; Finke 2016)

$$\begin{eqnarray}&&{\gamma }_{\mathrm{BLR}}^{{\prime} }=\displaystyle \frac{1}{{\delta }_{{\rm{D}}}}\sqrt{\displaystyle \frac{{E}_{\mathrm{obs},\mathrm{BLR}}(1+z)}{2{\epsilon }_{\mathrm{Ly}\alpha }}},\end{eqnarray} \tag{ 14 }$$

where ${E}_{\mathrm{obs},\mathrm{BLR}}$ is the observed γ-ray energy of IC-scattered BLR photons. From Equations (10), (13), and (14), it becomes clear that the cooling time scales as ${t}_{\mathrm{cool},\mathrm{BLR}}\propto {r}^{2}{{\rm{\Gamma }}}^{-2}$, i.e., the cooling becomes less efficient for large distances. Furthermore, if the γ-ray emission is produced very far away from the BLR, the external photons will appear as a point source illuminating the emission region from behind, so that ${u}_{\mathrm{BLR}}^{{\prime} }=(1/4){{\rm{\Gamma }}}^{-2}{u}_{\mathrm{BLR}}$ (Dermer & Schlickeiser 1994), leading to an additional decrease of the cooling time. The BLR cooling time for one flare of CTA 102 is shown for r = r_lim and r = 1 pc as a function of γ' and ${E}_{\mathrm{obs},\mathrm{BLR}}$ in Figure 13, assuming again the average values for δ_D ∼ 31 and Γ_L ∼ 22. Klein–Nishina effects become important for ${\gamma }_{\mathrm{BLR}}^{{\prime} }\gtrsim 2\times {10}^{3}$ or ${E}_{\mathrm{obs},\mathrm{BLR}}\gtrsim 1\,\mathrm{GeV}$ and a clear departure from the Thomson regime in which ${t}_{\mathrm{cool}}\propto {({\gamma }_{\mathrm{BLR}}^{\prime} )}^{-1/2}$ becomes visible. The decay times derived from the fit of the two bright flares of CTA 102 around MJD 57,738 are also shown in Figure 13 (see the first solid horizontal line in the panel in the fourth row and second column in Figure 5), where one decay time is shown as an upper limit due to its large uncertainty. If the decay is indeed caused by IC scattering with BLR photons, a distance of r ∼ 1 pc is still compatible with the observed decay times. Similar conclusions can be drawn for the other sources. We provide the cooling times at 300 MeV and the shortest decay times in Table 5.

Zoom In Zoom Out Reset image size

Figure 13 also shows the variability time t_var derived for the suborbital period in Section 4.2, which is, however, in the rising part of the flare (see Figure 9). If equally short decay times were observed, this would suggest a distance r ∼ r_lim.

As noted by Dotson et al. (2012), the energy dependence of the cooling times could further reveal the dominant photon field responsible for IC scattering: while Klein–Nishina effects become important already at 1 GeV for scattering with BLR photons, the Thomson regime should be valid to higher energies for IC scattering with photons of the dusty torus (dt). Again following Finke (2016), we assume that the torus also has a ring geometry and emits monochromatic photons with energy _dt = 2.7k_BT_dt, with k_B the Boltzmann constant. Generic values for the dust temperature T_dt are around 1000 K, and the dust luminosity is taken to be ξ_dt L_disk, with a dust scattering fraction ξ_dt = 0.1. We adopt these values for all sources except PKS B1222+216, 3C 273, and CTA 102. For PKS B1222+216 and CTA 102, Malmrose et al. (2011) found T_dt = 1200 K and a dust luminosity of 7.9 × 10⁴⁵ and 7 × 10⁴⁵ erg s⁻¹, respectively. Hao et al. (2005) observed silicate emission in 3C 273 from which they deduced a silicate temperature of 140 K and luminosity of ∼10⁴⁵ erg s⁻¹ but were unable to derive a temperature of the dust (see also the discussion in Malmrose et al. 2011). Soldi et al. (2008) instead found a dust temperature of 1200 K. We adopt the latter value and again set ξ_dt = 0.1. The sublimation radius of the dust torus, R_dt = 3.5 × 10¹⁸ cm(L_disk/10⁴⁵ erg s⁻¹)^1/2(T_dt/10³ K)^−2.6, is used as the ring radius. Making the appropriate substitutions in Equations (10) and (12)–(14), we plot the cooling time ${t}_{\mathrm{cool},\mathrm{dt}}$ for r = r_lim as a gray line in Figure 13 (note the additional x-axis, since ${E}_{\mathrm{obs},\mathrm{BLR}}\ne {E}_{\mathrm{obs},\mathrm{dt}}$). Indeed, Klein–Nishina effects only become relevant at ${E}_{\mathrm{obs},\mathrm{dt}}\gtrsim {10}^{2}\,\mathrm{GeV}$ or γ' ≳ 10⁵. The cooling times at ${E}_{\mathrm{obs},\mathrm{dt}}=316\,\mathrm{MeV}$ are also provided in Table 5.

To further investigate the energy dependence of the decay times, we split the energy range of our analysis into three energy bins, 100–300 MeV, 300 MeV–1 GeV, and 1–100 GeV, and recompute the orbital light curves. The energy bins are chosen as a compromise between the number of bins and sufficient photon statistics in each bin. The light curves for which at least two BBs are identified in each energy bin are shown in Figure 14. We repeat the fits of the exponential profiles to the energy-dependent light curves but allow only one flare profile per HOP group. The resulting decay times are also plotted in Figure 13. Only for the 300 MeV–1 GeV energy bin is the double peak of the flare resolved, and in general, the fit qualities are rather poor, with χ² dof^–1 between 1.67 and 2.26. The steep χ² curves also explain the rather small error bars on τ_decay. From the fit values, we cannot draw a conclusion as to whether the decay times evolve with (E_obs)^−1/2, as expected in the Thomson regime. Due to the lack of high photon statistics at energies beyond 10 GeV, where the differences in cooling times between the dust torus and BLR become more pronounced, we are not able to use the method suggested by Dotson et al. (2012) to determine the photon field dominating the IC scattering. This conclusion also holds for the other sources.

Zoom In Zoom Out Reset image size

Interestingly, the BBs for the energy-dependent light curves of the flares of 3C 279 and PKS 1510–089 and the last flare of CTA 102 seem to show time lags between the energy bands, with the high-energy emission leading the low-energy γ-rays. For the 3C 279 flare around MJD 57,188, Paliya (2015) could not find any time lags between the energy bins 0.1–1 GeV and above 1 GeV using the Z-transformed discrete correlation function (DCF; Alexander 1997, 2013). Using the same methodology, we show the DCFs for our energy-dependent light curves with at least two BBs per energy bin in Figure 15. We mark the time lags τ with horizontal lines at the maximum DCF values if $\max (\mathrm{DCF})\gt 2\sqrt{\mathrm{Var}(\mathrm{DCF})}$. In contrast to Paliya (2015), we find evidence that the emission above 1 GeV leads the emission at lower energies with ∼0.1 days. However, from the fits to the light curves, the decay time at higher energies is actually longer (0.45 ± 0.16 days above 1 GeV versus 0.22 ± 0.03 between 0.1 and 0.3 GeV). Therefore, the lag might not be associated with cooling but rather with a changing particle injection. From the spectral variation (Figure 8), it seems that the time bin before the peak of the flare has a harder spectrum; however, the uncertainties are too large to draw firm conclusions. For the CTA 102 flare around MJD 57,758, we also find that the high-energy emission is leading, whereas for the flare at MJD 57,749, the picture is reversed. For 3C 454.3, the DCF also indicates that the low-energy emission is leading the high-energy emission, again suggesting that these lags are connected to the injection of particles rather than radiative cooling.

Zoom In Zoom Out Reset image size

5.2.2. Synchrotron γ-Rays

Cross-correlating γ-ray light curves with radio light curves provides an alternative method to locate the γ-ray emitting region (e.g., Fuhrmann et al. 2014). Under the assumption that the flares are produced in a common compact emission region moving down the jet (e.g., Max-Moerbeck et al. 2014a), the distance ${d}_{\gamma ,r\nu }$ between the γ-ray sphere, where the γ-ray opacity due to, e.g., absorption in the BLR, becomes less than unity (Blandford & Levinson 1995), and the radio core, where synchrotron self-absorption becomes negligible (Königl 1981), can be estimated from the time lag between the light curves,

$$\begin{eqnarray}&&{d}_{\gamma ,r\nu }=\displaystyle \frac{{\rm{\Gamma }}{\delta }_{{\rm{D}}}\beta c{\tau }_{\mathrm{peak},\gamma ,r\nu }}{1+z},\end{eqnarray} \tag{ 15 }$$

where ${\tau }_{\mathrm{peak},\gamma ,r\nu }$ is the time lag corresponding to a peak in the cross-correlation function between the γ-ray and radio light curves obtained at frequency ν. Under the assumption that the radio emission lags the γ-rays, the distance of the γ-ray emission region to the central black hole is thus ${d}_{\gamma }={d}_{\mathrm{core},r\nu }-{d}_{\gamma ,r\nu }$, where ${d}_{\mathrm{core},r\nu }$ is the position of the radio core at frequency ν. The core position itself is frequency-dependent (the core shift effect; see, e.g., Lobanov 1998),

$$\begin{eqnarray}&&{d}_{\mathrm{core},r\nu }=\displaystyle \frac{{{\rm{\Omega }}}_{r\nu }}{{\nu }^{1/{k}_{r}}\sin {\theta }_{\mathrm{obs}}},\end{eqnarray} \tag{ 16 }$$

where k_r depends on the electron energy spectrum and the magnetic field in the emitting region (Königl 1981) and

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{r\nu }=4.85\times {10}^{-9}\displaystyle \frac{{\rm{\Delta }}{r}_{\mathrm{mas}}{d}_{{\rm{L}}}}{1+z}\displaystyle \frac{{\nu }^{1/{k}_{r}}{\nu }_{0}^{1/{k}_{r}}}{{\nu }_{0}^{1/{k}_{r}}-{\nu }^{1/{k}_{r}}},\end{eqnarray} \tag{ 17 }$$

where d_L is the luminosity distance and Δr_mas is the offset between the radio cores in milliarcseconds at frequencies ν and ν₀. The offest is related to the time lag between two radio light curves ${\tau }_{\mathrm{peak},r\nu ,r{\nu }_{0}}$ through ${\rm{\Delta }}{r}_{\mathrm{mas}}=\mu {\tau }_{\mathrm{peak},r\nu r{\nu }_{0}}$, where μ is the jet proper motion.

The proper motion and core position at 15 GHz have been determined from the MOJAVE VLBI blazar monitoring program (Pushkarev et al. 2012; Lister et al. 2016). The following distances ${d}_{\mathrm{core},15\mathrm{GHz}}$ were determined under the assumption that k_r = 1: for PKS B1222+216, ${d}_{\mathrm{core},15\mathrm{GHz}}=23.41\,\mathrm{pc}$; for 3C 279, ${d}_{\mathrm{core},15\mathrm{GHz}}\lt 7.88\,\mathrm{pc}$; for PKS 1510–089, ${d}_{\mathrm{core},15\mathrm{GHz}}\,=17.71\,\mathrm{pc}$; for CTA 102, ${d}_{\mathrm{core},15\mathrm{GHz}}=46.7\,\mathrm{pc}$; and for 3C 454.3, ${d}_{\mathrm{core},15\mathrm{GHz}}=20.36\,\mathrm{pc}$. Dedicated analyses have also been carried out and found for 3C 454.3, k_r = 0.6–0.8 and ${d}_{\mathrm{core},15\mathrm{GHz}}\sim 38\,\mathrm{pc}$, and, since ${d}_{\mathrm{core},r\nu }\propto {\nu }^{-1/{k}_{r}}$, ${d}_{\mathrm{core},43\mathrm{GHz}}\,\sim 9\,\mathrm{pc}$ (Kutkin et al. 2014). For 3C 273, Vol'vach et al. (2013) found k_r = 1.4 and ${d}_{\mathrm{core},r\nu }=134{\nu }^{-1/1.4}$ using radio observations at frequencies between 4.8 and 362 GHz. Lastly, Fromm et al. (2015) conducted VLBA observations of CTA 102 ranging from 5 to 86 GHz and found k_r = 1.0 as a best-fit value and ${d}_{\mathrm{core},86\mathrm{GHz}}\sim 7\,\mathrm{pc}$. Provided that we can estimate ${\tau }_{\mathrm{peak},15\mathrm{GHz},r\nu }$, it is possible with the above results to estimate the core position at an arbitrary radio frequency using Equations (16) and (17). In order to arrive at an estimate for d_γ, the only remaining task is to perform a cross-correlation study between γ-ray and radio light curves.

We search for time lags between the Fermi LAT light curves and radio light curves obtained with the Owens Valley Radio Observatory (OVRO) at 15 GHz, the Atacama Large submillimeter/millimeter Array (ALMA) between 84 and 116 GHz (Band 3, 3.6–2.6 mm), and the Submillimeter Array (SMA) at 230 GHz (1.3 mm). All of the studied FSRQs are included in an ongoing blazar monitoring program at OVRO (Richards et al. 2011) and SMA (Gurwell et al. 2007) and also serve as calibrators for SMA and ALMA (Bonato et al. 2018) at millimeter wavelengths.²¹ We show all radio and γ-ray light curves in Figure 16. It is evident that the OVRO light curves show variations on longer timescales and with less flicker noise behavior. For 3C 454.3, there appears to be a correlation between the radio and γ-ray flux, at least for the giant flare in 2010 (around MJD 55,500). Due to scarce and uneven sampling, we do not use the ALMA and SMA light curves of PKS B1222+216 and PKS 1510–089. We also do not include the SMA light curve of CTA 102 in the following analysis for the same reason.

Zoom In Zoom Out Reset image size

To quantify the correlations, we again closely follow the methodology laid out by Max-Moerbeck et al. (2014b). For two light curves with fluxes a_i and b_j measured at times t_ai and t_bj, we compute the local cross-correlation function (LCCF)

$$\begin{eqnarray}&&\mathrm{LCCF}(\tau )=\displaystyle \frac{1}{M}\displaystyle \frac{\sum ({a}_{i}-{\bar{a}}_{\tau })({b}_{j}-{\bar{b}}_{\tau })}{{\sigma }_{a\tau }{\sigma }_{b\tau }},\end{eqnarray} \tag{ 18 }$$

where the sum runs over the M pairs for which τ ≤ t_ai − t_bj < τ + Δt for some chosen time step Δt, and ${\bar{a}}_{\tau }$ and ${\bar{b}}_{\tau }$ and ${\sigma }_{a\tau }$ and ${\sigma }_{b\tau }$ are the flux averages and standard deviations over the M pairs, respectively (Welsh 1999). The LCCF is bound between −1 and 1 and has much larger efficiency in recovering linear correlations between light curves compared to the DCF (Max-Moerbeck et al. 2014b). For the binning of the time lags τ, we choose half the maximum of the median of the time separations between consecutive data points in the two light curves. The minimum and maximum values of τ are chosen to be ±0.5 times the length of the shortest light curve (Max-Moerbeck et al. 2014a).

We determine the significance of a peak in the LCCF by cross-correlating pairs of simulated light curves. For the γ-ray light curves, the simulation proceeds in the same way as described in Section 3, where we use the best-fit values $\hat{\beta }$ for the assumed PSDs. For the radio light curves, we proceed in a similar way. First, we determine the best-fit PSDs similar to the γ-ray light curves. In order to achieve a good fit between the observed and simulated PSDs, we change the methodology used for the γ-ray light curves in the following ways. Instead of matching the flux probability distribution of the light curves, as suggested by Emmanoulopoulos et al. (2013), we use variance matching (Max-Moerbeck et al. 2014b).²² Furthermore, we do not apply uncertainties to the simulated light curves. Doing so generally leads to a strong flattening of the periodograms at high frequencies when they become dominated by white noise introduced by the uncertainties. This is not observed in the periodogram derived from observations. The reason might be a correlation between uncertainties and flux, which is not taken into account by the adopted simulation scheme and can lead to an overestimation of the simulated uncertainties. Lastly, the radio light curves are unevenly sampled and can show large observational gaps. Simply applying the interpolation scheme used for the the γ-ray light curves would mean that most data points that enter the calculation of the periodogram are actually interpolated. To mitigate this problem, we split the light curves where they show large gaps. We found that a split at gaps that are 20 (4.5) times larger than the median separation between consecutive measurements for the OVRO and SMA (ALMA) light curves provides a good compromise between minimizing the number of splits and too few data points within a light-curve segment. Furthermore, for the interpolated light curves, we use a time step equal to the 80% quantile of the observed separation (the median would correspond to the 50% quantile). In this way, we lose sensitivity to the highest frequencies but end up with interpolated light curves with roughly the same number of data points as the observed ones. We do not average the interpolated flux points as in the γ-ray case, since radio observations are usually short in duration and report flux densities instead of integrated fluxes. The periodograms of the individual light-curve segments are finally log-averaged following Papadakis & Lawrence (1993).

We report the best-fit slopes of the assumed power-law PSDs $\hat{\beta }$, their confidence interval, and the p_β-value of the fits in Table 6. The confidence interval and p_β-value are determined in the same way as described in Section 3. All fits show a high fit quality. For the SMA and OVRO light curves for 3C 454.3, as well as for the 3C 273 light curve obtained with ALMA, we are only able to provide upper bounds on β. In general, we confirm the trend that the OVRO light curves show a softer PSD $\hat{\beta }\gtrsim 2$ for all sources. Moving to higher frequencies with ALMA and SMA, the PSD hardens and becomes more flicker noise–like. Comparing our results for the OVRO light curves to previous analyses of Max-Moerbeck et al. (2014a), who used 4 yr of data, we find them to be consistent within the uncertainties.

Table 6.  Results from PSD Analysis of Radio Light Curves, as well as γ-Ray and Radio LCCF Results

Source	$\hat{\beta }$	pβ	τpeak [days]	pτ	${d}_{\gamma ,r}$ [pc]
OVRO
PKS B1222+216	${1.92}_{-0.59}^{+0.39}$	0.59	⋯	⋯	⋯
3C 273	${2.38}_{-0.97}^{+0.30}$	0.94	$-{416.5}_{-140.0}^{+217.0}$	0.0068	$10.96\,[5.2,14.6]\,\pm 4.4$
3C 279	${2.29}_{-0.94}^{+0.32}$	0.71	⋯	⋯	⋯
PKS 1510–089	${1.89}_{-0.84}^{+0.45}$	0.34	⋯	⋯	⋯
CTA 102	${2.23}_{-0.92}^{+0.26}$	0.84	⋯	⋯	⋯
3C 454.3	${2.20}_{-2.20}^{+0.36}$	0.40	$-{101.5}_{-112.0}^{+49.0}$	0.0156	$15.39\,[8.0,32.4]\,\pm 2.8$
ALMA Band 3
3C 273	${2.12}_{-2.12}^{+0.40}$	0.73	⋯	⋯	⋯
3C 279	${1.82}_{-0.45}^{+0.38}$	0.89	⋯	⋯	⋯
CTA 102	${1.94}_{-1.33}^{+0.42}$	0.45	$-{216.0}_{-11.0}^{+209.0}$	0.0092	$58.85\,[1.9,61.8]\,\pm 7.3$
3C 454.3	${1.73}_{-0.30}^{+0.36}$	0.25	$-{27.0}_{-30.0}^{+30.0}$	0.0164	$4.09\,[-0.5,8.6]\,\pm 0.7$
SMA 1.3 mm
3C 273	${1.48}_{-0.33}^{+0.40}$	0.17	$-{122.5}_{-7.0}^{+84.0}$	0.0088	$3.22\,[1.0,3.4]\,\pm 1.3$
3C 279	${1.61}_{-0.28}^{+0.16}$	0.97	⋯	⋯	⋯
3C 454.3	${1.64}_{-1.64}^{+0.31}$	0.21	${10.5}_{-28.0}^{+21.0}$	0.0002	$-1.59\,[-4.8,2.7]\,\pm 0.3$

Note. For the LCCF analysis, we only report time lags <0 and with a significance p_τ < 0.05. The range of ${d}_{\gamma ,r}$ values reported in square brackets is due to the uncertainty on τ_peak; the remaining uncertainties are the propagated errors on Γ and δ_D.

Download table as:  ASCIITypeset image

Having determined the best-fit PSDs, we use them to create artificial light curves in the same way as for fitting the periodogram itself. We then calculate the LCCF between 5000 pairs of uncorrelated simulated light curves and derive confidence bands on the LCCF. We use the confidence bands to determine the p_τ-value, which gives the probability of finding an LCCF value at a given τ greater or equal to the observed value under the assumption that the light curves are uncorrelated.

For observed LCCFs where we find time lags with a significance 1 − p_τ > 0.95, we estimate the uncertainty on the peak time using flux randomization and random subsample selection (drawing 1000 samples) following Peterson et al. (1998), as suggested by Max-Moerbeck et al. (2014b).

We show the observed LCCF between γ-ray and radio light curves in Figure 17 for the cases where we find a significant peak in the LCCF. This is the case for 3C 273 (OVRO and SMA), CTA 102 (ALMA), and 3C 454.3 (all radio and millimeter light curves correlate with the γ-ray light curve). The correlation between the SMA and γ-ray light curve of 3C 454.3 is the most significant, with p_τ = 2 × 10⁻⁴ (3.72σ). The peak and uncertainties are marked by a dotted line and shaded region in Figure 17. They are also summarized, together with the p_τ values, in Table 6. In the table, we only consider peaks with τ < 0 (the γ-ray light curve leads the radio light curve); however, for 3C 273 and 3C 454.3, peaks with τ > 0 are also visible. For 3C 454.3, these peaks occur at large values of τ that might be due to the several small flares observed at γ-ray energies. Since the overlap between the light curves is smaller for larger values of $| \tau |$, we deem these peaks less credible. For 3C 273, the situation is less clear. The SMA light curve shows two prominent flares around the γ-ray flare, which gives rise to the two peaks in the LCCF. The low state of the source in recent years and the less dense sampling of the SMA light curve render it difficult to draw firm conclusions. As we see below, even the peak at τ < 0 does not lead to constraints on the position of the γ-ray emitting region.

Zoom In Zoom Out Reset image size

Using Equation (15) and the peaks identified in the LCCF, we can now derive the distance between the γ-ray and radio emitting regions. For 3C 454.3, the time lags between the γ-ray and millimeter light curves of ALMA and SMA are consistent with zero; hence, ${d}_{\gamma ,r}$ is consistent with zero as well. For the correlation with the OVRO light curve, a longer time lag of ${\tau }_{\mathrm{peak},\gamma ,15\mathrm{GHz}}=-102$ days is found, placing ${d}_{\gamma ,r}$ between ∼5 and 35 pc. This time lag is consistent with the recent DCF analysis carried out by Liodakis et al. (2018), who found a time lag of (115 ± 6) days at 2.5σ significance using 8 yr of OVRO data and the Fermi LAT weekly monitored light curve.

For CTA 102, we only find a significant lag between the ALMA and γ-ray light curves. The time lag translates into ${d}_{\gamma ,r}\sim [-5;69]\,\mathrm{pc}$, taking uncertainties on ${\tau }_{\mathrm{peak},\gamma ,100\mathrm{GHz}}$, as well as Γ and δ_D, into account. Hence, the distance is also consistent with zero.

We also find a significant correlation between the γ-ray light curve of 3C 273 and both light curves of SMA and OVRO. For the 15 GHz OVRO light curve, ${d}_{\gamma ,r}$ is found between 0.8 and 19 pc, whereas for 230 GHz, the distance falls between −0.7 and 4.7 pc and is consistent with zero.

In order to determine the core positions at the ∼100 GHz (ALMA) and 230 GHz (SMA) cores, we carry out a cross-correlation between the radio light curves to derive ${\tau }_{\mathrm{peak},15\mathrm{GHz},r\nu }$. The results are reported in Table 7. Using Equations (16) and (17), we arrive at the new core position, which we also show in Table 7, where we use the average values of the measured jet proper motion (see Table 4 in Lister et al. 2016), the observation angle θ_obs reported in Jorstad et al. (2017), and the values of k_r found in the dedicated analyses discussed above. For 3C 273 and 3C 454.3, we find that our core positions are consistent with values calculated from the ${d}_{\mathrm{core}}\propto {\nu }^{-1/{k}_{r}}$ relation obtained by Vol'vach et al. (2013) and Kutkin et al. (2014), respectively. This is a nontrivial result, since we have combined our time lags with jet proper motions from MOJAVE and jet properties from VLBA observations.

Table 7.  Results for Time Lags and Core Positions from a Radio/Radio LCCF Analysis

Source	${\tau }_{\mathrm{peak},r{\nu }_{1},r{\nu }_{2}}$ [days]	pτ	${d}_{\mathrm{core},r{\nu }_{1}}$ [pc]
ALMA Band 3 and OVRO
3C 273	$-{161}_{-36}^{+72}$	0.0222	$1.5\,[0.8,1.8]\,\pm 0.6$
3C 279	$-{622}_{-168}^{+154}$	0.0036	$8.1\,[6.1,10.3]\,\pm 2.6$
CTA 102	⋯	⋯	⋯
3C 454.3	$-{667}_{-30}^{+100}$	0.0782	$4.6\,[3.9,4.8]\,\pm 2.6$
SMA 1.3 mm and OVRO
3C 273	$-{427}_{-87}^{+293}$	0.0256	$4.0\,[1.2,4.8]\,\pm 1.5$
3C 279	$-{165}_{-125}^{+12}$	0.0152	$2.2\,[2.0,3.8]\,\pm 0.7$
3C 454.3	$-{110}_{-0}^{+14}$	0.0924	$0.8\,[0.7,0.8]\,\pm 0.4$
SMA 1.3 mm and ALMA Band 3
3C 273	${1}_{-36}^{+9}$	0.0000	⋯
3C 279	$-{188}_{-119}^{+175}$	0.0002	⋯
3C 454.3	${3}_{-20}^{+10}$	0.0004	⋯

Note. For the LCCF analysis, we only report time lags with a significance p_τ < 0.1. The range of core positions in square brackets is due to the uncertainty on τ_peak, whereas the remaining uncertainties are the propagated errors on θ_obs and the jet proper motion μ.

Download table as:  ASCIITypeset image

Combining the core positions and ${d}_{\gamma ,r}$, we find that for 3C 454.3, the γ-ray emitting region is consistent with the position of the SMA millimeter core at $\sim {0.8}_{-0.5}^{+0.4}\,\mathrm{pc}$, also in agreement with the LCCF from ALMA. The OVRO result is only marginally in agreement with this result, suggesting instead that d_γ ≳ 3 pc. Given the larger uncertainty on the time lag and lower significance of the correlation, we deem the results at millimeter wavelengths as more robust. They also agree with the findings of Fuhrmann et al. (2014).

For CTA 102, we do not find a significant correlation with the OVRO light curve; therefore, we use the core position at 86 GHz, d_{core,86 GHz} ∼ 7 pc (Fromm et al. 2015). With the results on ${d}_{\gamma ,r}$, we also find that for CTA 102, the distance of the γ-ray emitting region is consistent with the millimeter core.

For 3C 273, we find that the derived core shift between 230 and 15 GHz of ∼4 pc is not consistent with the ν^−1/1.4 relation obtained by Vol'vach et al. (2013). Nevertheless, within the uncertainties, ${d}_{\gamma ,r}$ is also consistent with zero, and hence γ-ray and millimeter emission could be produced cospatially.

We conclude this section with a word of caution. One should note that a peak in the LCCF does not necessarily measure the lag in the intrinsic processes that generate the "injections" that produce flares. There will be an offset that depends on the relative shapes of the flare light curves at the different wavelengths. To estimate the size of this effect, we compute the autocorrelation functions (ACFs) for the γ-ray and radio/millimeter light curves for the sources for which we find a significant lag. For each source, the central peaks of the ACFs in radio/millimeter and γ-rays have comparable widths (within a factor of 2); the only exception is for 3C 273, for which the central peak in the ACF of the OVRO light curve is much broader than the ACF of the Fermi light curve. From this, we conclude that it is possible that the intrinsic flare shapes at the different wavelengths are similar. However, one has to keep the caveat in mind that different flare shapes can produce similar ACFs, as demonstrated in Figure 14 of Scargle (1981). A more quantitative exploration of this issue will be presented elsewhere.

From the weekly binned full 9.5 yr light curves, we are able to determine global temporal properties (Section 3) such as flux distributions and PSDs. The former are well described with BPLs or lognormal distributions, while the latter, derived following Max-Moerbeck et al. (2014b), indicate that the light curves show power-law-type PSDs with indices β ∼ 1 indicating flicker noise. With the exception of 3C 279, our slopes are also compatible within 1σ–2σ with the slopes found in the optical R band for the FSRQs considered by Chatterjee et al. (2012; 3C 273, 3C 279, PKS 1510–089, 3C 454.3). This could indicate that the emission in the γ-ray band and at optical wavelengths is produced by the same underlying electron population through external Compton scattering and synchrotron emission, respectively (Finke & Becker 2014, 2015). We also developed a novel objective algorithm to determine the flux level of the QB (see Section 3.1). The determination of the quiescent flux can have important implications for emission models for large-scale blazar jet components (Meyer & Georganopoulos 2014). A thorough analysis of the quiescent state of the studied sources will be provided elsewhere.

We use the light curves with one flux bin per orbit to derive local temporal flare properties by fitting the light curves with exponential profiles (Section 4). The obtained rise and decay times show that rapid variability at timescales of the order of the horizon-crossing timescale, which range from ∼0.5 to ∼2 hr for the considered black hole masses, are a common feature of all sources. In general, we find a large variety of flares showing both FRED behavior and the opposite, similar to the results found by Roy et al. (2019; see their Figure 10). No clear trend between flare asymmetry and other flare parameters, such as peak flux, integrated flux, or flare duration, is found. No apparent evolution of these quantities with time is observed either. This variety of flare profiles could be explained in the scheme of magnetic reconnection with different orientations of the reconnection layers leading to a variety of Doppler factors of the injected plasmoids (e.g., Petropoulou et al. 2016; Christie et al. 2019). With our novel approach of identifying flares, it will be possible in the future to build large statistical samples of flares (selected not only by their peak flux but also, e.g., by their integrated flux) whose properties could be compared in more detail to predictions of the reconnection scenario. Small and fast plasmoids injected close to the line of sight could also explain minute-scale variability (Petropoulou et al. 2016) for which we find evidence at the 2σ significance level (posttrial) in suborbital light curves of at least two sources, 3C 279 and CTA 102 (Table 4). For 3C 454.3 and PKS 1510–089, the BBs also indicate variability on such short timescales; however, a fit with a constant flux to the light curves of these orbits cannot be rejected beyond the 2σ (posttrial) level. Other possible explanations of such short variability include an energy-dependent kinetic beaming of particles during reconnection (Cerutti et al. 2012), radiative cooling of a plasma accelerated by recollimation shocks (Bodo & Tavecchio 2018), or synchrotron radiation by electrons accelerated to energies beyond γ' ≳ 10⁶ (Ackermann et al. 2016b), a scenario motivated by the γ-ray flares of the Crab Nebula (Abdo et al. 2011b) and discussed for the flare of PKS B1222+216 (Nalewajko et al. 2012). We have proposed two alternative scenarios in which the γ-rays are shielded from low-energy photons by a plasma sheath (Section 5.1.2) and γ-rays are produced by synchrotron emission of electron–positron pairs created by the interaction of protons with low-energy photons (Section 5.2.2). Our discussion is of a qualitative nature only, and a full treatment will be presented elsewhere.

We have also investigated the location of the γ-ray emitting region through three approaches: searches for a spectral cutoff, a comparison of decay times with predictions of radiative cooling times, and a correlation between γ-ray and radio light curves. Using the BLR model of Finke (2016), we find no significant spectral cutoff due to the interaction of γ-rays with BLR photons, which places the γ-ray emission region outside or on the edge of the BLR, $r\gtrsim {R}_{\mathrm{Ly}\alpha }$ or $\gtrsim {10}^{3}{r}_{g}$ (see Table 5). These lower limits are conservative in the sense that more sophisticated models of the BLR (e.g., Abolmasov & Poutanen 2017), which include continuum emission and a more realistic geometry, predict even larger optical depths than the model used here. We also do not account for a possible brightening of the BLR emission during γ-ray flares, as found by León-Tavares et al. (2013). The observed spectra over different time periods are provided for completeness in the Appendix.

The observed decay times of the brightest flares are compatible with the radiative cooling times predicted from IC scattering of electrons with BLR photons for distances of the γ-ray emitting regions up to r ∼ 1 pc. The IC scattering with photons of the dust torus yields cooling times on the order of hours for values of r up to the sublimation radius of ∼3 pc. This is only compatible with a subset of the analyzed flares; see Table 5. In order to reconcile cooling times with minute-scale variability, the emission region would have to be close to the obtained lower limits. At the same time, the value of r needs to be compatible with the amount of observed γ-ray emission. The IC emission scales as ${\delta }_{{\rm{D}}}^{3}/{x}^{2}$, where ${x}^{2}={R}_{\mathrm{li}}^{2}+{r}^{2}$ (see, e.g., Equation (87) in Finke 2016). A future comparison of the distance derived from IC emission predicted from multiwavelength modeling and cooling times could provide further insight into this issue.

In principle, the energy dependence of the observed decay times can be used to distinguish cooling with BLR and dust torus photons, as proposed by Dotson et al. (2012). At high electron energies, cooling with BLR photons occurs in the Klein–Nishina regime, whereas cooling with dust torus photons occurs in the Thompson regime (see Figure 13). However, the fit of the exponential flare profiles to light curves in different energy bands yields inconclusive results, since the photon statistics are not sufficient to distinguish between the two scenarios.

We find correlations significant beyond the 2.1σ level between the γ-ray and radio light curves for 3C 273 (correlations found with OVRO and SMA light curves), CTA 102 (with ALMA), and 3C 454.3 (with OVRA, ALMA, and SMA). The time lags between the γ-ray and millimeter light curves of ALMA and SMA are consistent with zero, which could indicate a cospatial production. This is consistent with the picture of superluminal knots passing through a standing shock associated with the radio core, as argued for the flares in PKS 1510–089 and 3C 454.3 (Marscher et al. 2010; Wehrle et al. 2012; Orienti et al. 2013). This would entail values of r on the order of parsecs. One has to keep in mind, however, that the inferred uncertainties on the time lags are large. Taken together with the uncertainties on the Doppler boost and bulk Lorentz factor, smaller values of r cannot be ruled out by our analysis.

Our three approaches to constrain the location of the γ-ray emitting region are all consistent with rather large distances from the central black hole, r ∼ 1 pc. If the distances are even larger, the radiative cooling with BLR photons becomes inefficient, and cooling through IC scattering with dust torus photons does not reproduce the observed flare decay times. Such large distances are at odds with the evidence with minute-scale variability, which we, however, only observe at ∼2σ posttrial significance in the rising parts of the flares. Densely sampled light curves at other wavelengths could provide further insight if this short variability is indeed present and possibly connected to the injection of relativistic particles. As noted by Costamante et al. (2018), the fact that we do not find significant absorption provides more promising prospects to detect these sources with future observations with the Cerenkov Telescope Array (CTA). The improved point-source sensitivity of the CTA, together with its energy range between 20 GeV and 300 TeV (Cerenkov Telescope Array Consortium et al. 2017), could lead to the detection of spectral absorption features during FSRQ flares due to the interaction of γ-rays with infrared photons from the dust torus that should become important at ∼TeV energies (see, e.g., Figure 14 in Finke 2016). The CTA observations could further effectively probe the shortest-variability timescales at very high γ-ray energies.