Variability and Spectral Characteristics of Three Flaring Gamma-ray Quasars Observed by VERITAS and Fermi-LAT

Aryeh Brill

Olivier Hervet

Janeth Valverde

C. B. Adams
Department of Physics and Astronomy, Barnard College, Columbia University, NY 10027, USA

J. Batshoun
Department of Physics, California State University - East Bay, Hayward, CA 94542, USA

W. Benbow
Center for Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA 02138, USA

A. Brill
Physics Department, Columbia University, New York, NY 10027, USA
Now at NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

J. H. Buckley
Department of Physics, Washington University, St. Louis, MO 63130, USA

M. Capasso
Department of Physics and Astronomy, Barnard College, Columbia University, NY 10027, USA

B. Cavins
Santa Cruz Institute for Particle Physics and Department of Physics, University of California, Santa Cruz, CA 95064, USA

J. L. Christiansen
Physics Department, California Polytechnic State University, San Luis Obispo, CA 94307, USA

P. Coppi260 Whitney Avenue, Yale University, New Haven, CT 06511, USA

M. Errando
Department of Physics, Washington University, St. Louis, MO 63130, USA

K. A Farrell
School of Physics, University College Dublin, Belfield, Dublin 4, Ireland

Q. Feng
Department of Physics and Astronomy, Barnard College, Columbia University, NY 10027, USA

J. P. Finley
Department of Physics and Astronomy, Purdue University, West Lafayette, IN 47907, USA

G. M. Foote
Department of Physics and Astronomy and the Bartol Research Institute, University of Delaware, Newark, DE 19716, USA

L. Fortson
School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA

A. Furniss
Department of Physics, California State University - East Bay, Hayward, CA 94542, USA

A. Gent
School of Physics and Center for Relativistic Astrophysics, Georgia Institute of Technology, 837 State Street NW, Atlanta, GA 30332-0430

C. Giuri
DESY, Platanenallee 6, 15738 Zeuthen, Germany

D. Hanna
Physics Department, McGill University, Montreal, QC H3A 2T8, Canada

T. Hassan
DESY, Platanenallee 6, 15738 Zeuthen, Germany

O. Hervet
Santa Cruz Institute for Particle Physics and Department of Physics, University of California, Santa Cruz, CA 95064, USA

J. Holder
Department of Physics and Astronomy and the Bartol Research Institute, University of Delaware, Newark, DE 19716, USA

M. Houck
Davidson College, 405 N Main St, Davidson, NC 28035, USA
Department of Civil and Environmental Engineering, University of California, Irvine, CA 92617, USA

T. B. Humensky
Physics Department, Columbia University, New York, NY 10027, USA

W. Jin
Department of Physics and Astronomy, University of Alabama, Tuscaloosa, AL 35487, USA

P. Kaaret
Department of Physics and Astronomy, University of Iowa, Van Allen Hall, Iowa City, IA 52242, USA

M. Kertzman
Department of Physics and Astronomy, DePauw University, Greencastle, IN 46135-0037, USA

D. Kieda
Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA

F. Krennrich
Department of Physics and Astronomy, Iowa State University, Ames, IA 50011, USA

S. Kumar
Physics Department, McGill University, Montreal, QC H3A 2T8, Canada

M. Lundy
Physics Department, McGill University, Montreal, QC H3A 2T8, Canada

G. Maier
DESY, Platanenallee 6, 15738 Zeuthen, Germany

C. E McGrath
School of Physics, University College Dublin, Belfield, Dublin 4, Ireland

P. Moriarty
School of Physics, National University of Ireland Galway, University Road, Galway, Ireland

R. Mukherjee
Department of Physics and Astronomy, Barnard College, Columbia University, NY 10027, USA

D. Nieto
Institute of Particle and Cosmos Physics, Universidad Complutense de Madrid, 28040 Madrid, Spain

M. Nievas-Rosillo
DESY, Platanenallee 6, 15738 Zeuthen, Germany

S. O'Brien
Physics Department, McGill University, Montreal, QC H3A 2T8, Canada

R. A. Ong
Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA

A. Oppenheimer
Hastings High School, 1 Mt Hope Blvd, Hastings-On-Hudson, NY 10706, USA
University of Chicago, Chicago, IL 60637, USA

A. N. Otte
School of Physics and Center for Relativistic Astrophysics, Georgia Institute of Technology, 837 State Street NW, Atlanta, GA 30332-0430

S. Patel
Department of Physics and Astronomy, University of Iowa, Van Allen Hall, Iowa City, IA 52242, USA

K. Pfrang
DESY, Platanenallee 6, 15738 Zeuthen, Germany

M. Pohl
Institute of Physics and Astronomy, University of Potsdam, 14476 Potsdam-Golm, Germany and DESY, Platanenallee 6, 15738 Zeuthen, Germany

R. R. Prado
DESY, Platanenallee 6, 15738 Zeuthen, Germany

E. Pueschel
DESY, Platanenallee 6, 15738 Zeuthen, Germany

J. Quinn
School of Physics, University College Dublin, Belfield, Dublin 4, Ireland

K. Ragan
Physics Department, McGill University, Montreal, QC H3A 2T8, Canada

P. T. Reynolds
Department of Physical Sciences, Munster Technological University, Bishopstown, Cork, T12 P928, Ireland

A. Rhatigan
Physics Department, California Polytechnic State University, San Luis Obispo, CA 94307, USA

D. Ribeiro
Physics Department, Columbia University, New York, NY 10027, USA

E. Roache
Center for Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA 02138, USA

J. L. Ryan
Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA

M. Santander
Department of Physics and Astronomy, University of Alabama, Tuscaloosa, AL 35487, USA

G. H. Sembroski
Department of Physics and Astronomy, Purdue University, West Lafayette, IN 47907, USA

D. A. Williams
Santa Cruz Institute for Particle Physics and Department of Physics, University of California, Santa Cruz, CA 95064, USA

T. J Williamson
Department of Physics and Astronomy and the Bartol Research Institute, University of Delaware, Newark, DE 19716, USA

(VERITAS Collaboration)

J. Valverde
Center for Space Sciences and Technology, University of Maryland Baltimore County, Baltimore, MD 21250, USA
NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
Laboratoire Leprince-Ringuet, École Polytechnique, CNRS/IN2P3, 91128 Palaiseau, France
Department of Astronomy, Columbia University, New York, NY 10027, USA

D. Horan
Laboratoire Leprince-Ringuet, École Polytechnique, CNRS/IN2P3, 91128 Palaiseau, France

S. Buson
Julius-Maximilians-Universität, 97070, Würzburg, Germany

C. C. Cheung
Space Science Division, Naval Research Laboratory, Washington, DC 20375, USA

S. Ciprini
Istituto Nazionale di Fisica Nucleare (INFN) Sezione di Roma Tor Vergata, Via della Ricerca Scientica 1, 00133, Roma, Italy
ASI Space Science Data Center (SSDC), Via del Politecnico, 00133, Roma, Italy

D. Gasparrini
Istituto Nazionale di Fisica Nucleare, Sezione di Roma “Tor Vergata”, I-00133 Roma, Italy
Space Science Data Center - Agenzia Spaziale Italiana, Via del Politecnico, snc, I-00133, Roma, Italy

R. Ojha
NASA HQ, Washington, DC 20546, USA

P. van Zyl
Hartebeesthoek Radio Astronomy Observatory (HartRAO)
South African Radio Astronomy Observatory (SARAO)

(Fermi-LAT Collaboration)

L. Sironi
Department of Astronomy and Columbia Astrophysics Laboratory, Columbia University, New York, NY 10027, USA

Flat spectrum radio quasars (FSRQs) are the most luminous blazars at GeV energies, but only rarely emit detectable fluxes of TeV gamma rays, typically during bright GeV flares. We explore the gamma-ray variability and spectral characteristics of three FSRQs that have been observed at GeV and TeV energies by Fermi-LAT and VERITAS, making use of almost 100 hours of VERITAS observations spread over 10 years: 3C 279, PKS 1222+216, and Ton 599. We explain the GeV flux distributions of the sources in terms of a model derived from a stochastic differential equation describing fluctuations in the magnetic field in the accretion disk, and estimate the timescales of magnetic flux accumulation and stochastic instabilities in their accretion disks. We identify distinct flares using a procedure based on Bayesian blocks and analyze their daily and sub-daily variability and gamma-ray energy spectra. Using observations from VERITAS as well as Fermi, Swift, and the Steward Observatory, we model the broadband spectral energy distributions of PKS 1222+216 and Ton 599 during VHE-detected flares in 2014 and 2017, respectively, strongly constraining the jet Doppler factors and gamma-ray emission region locations during these events. Finally, we place theoretical constraints on the potential production of PeV-scale neutrinos during these VHE flares.

§ INTRODUCTION

Blazars are a class of active galactic nuclei (AGN) with jets oriented nearly along our line of sight. This alignment produces beamed emission, so that many blazars  show superluminal motion in their jets <cit.> and have a gamma-ray luminosity dominating their bolometric power. In jet models, high-energy electrons in a relativistically outflowing jet, ejected from an accreting supermassive black hole, are responsible for the synchrotron radiation seen as the radio to UV continuum from blazars [Blandford & Königl, 1979]. Blazars have a spectral energy distribution (SED) exhibiting a double-humped structure, with low-energy synchrotron and high-energy gamma-ray-peaked components.

Blazars are the most common gamma-ray-emitting objects in the extragalactic sky. Observationally, they can be divided into two classes: BL Lacertae (BL Lac) objects, the aligned counterparts to Fanaroff-Riley I radio galaxies, and flat spectrum radio quasars (FSRQs), the counterparts to Fanaroff-Riley II radio galaxies [Fanaroff & Riley, 1974].

FSRQs are low-synchrotron-peaked (LSP) blazars, with synchrotron peak frequency less than $10^{14}$ Hz. The bolometric luminosity of FSRQs is typically greater than that of BL Lac objects. The anti-correlation of synchrotron luminosity with peak frequency is an empirical relationship known as the blazar sequence [Fossati et al., 1998, Nieppola et al., 2008], though its intrinsic validity has been disfavored by more recent work [Keenan et al., 2021]. Accordingly, while FSRQs make up only 8 of the 79 AGN that have been detected in the TeV band to date[<http://tevcat.in2p3.fr/>], they are more commonly detected at GeV energies, comprising 650 of 2863 AGN detected by the Large Area Telescope (LAT) on board the Fermi Gamma-Ray space telescope (Fermi-LAT) [Ajello et al., 2020] and dominating the blazar population detected by the Energetic Gamma Ray Experiment Telescope (EGRET) on the Compton Gamma-Ray Observatory [Mukherjee, 2001].

The SED of an FSRQ is generally dominated by the gamma-ray emission component, which peaks in the high-energy (HE; ${\sim}$GeV) band. FSRQs are believed to possess several structures producing radiation fields external to the jet, including a broad line region (BLR) and a dust torus. TeV detections of FSRQs are particularly interesting because the external radiation fields might be expected to produce increased Compton cooling of electrons and to absorb energetic gamma rays by pair production, leading to a cutoff in the gamma-ray spectrum above the GeV band <cit.>.

Blazars have been observed to be variable at all wavelengths and at timescales down to several minutes in both the GeV and TeV bands [Ackermann et al., 2016, Aharonian et al., 2007]. However, the physical mechanisms that drive this variability are unclear. Different processes, possibly originating at different locations in the AGN, may drive variable emission occurring at different timescales. By providing an upper bound on the light crossing time, the timescale of variability constrains the apparent size of the emission region, giving information on the location and mechanism of the gamma-ray emission. While short variability timescales observed in blazars suggest that the emission may be connected to processes in the central engine or accretion disk, the ability of very-high-energy (VHE; ${\gtrsim}100$ GeV) emission to escape the AGN implies an origin further out in the jet, where absorption is reduced [Abeysekara et al., 2015].

Over longer timescales, blazar variability can be studied through the flux distribution describing the relative frequencies of different flux levels. Blazar flux distributions exhibit long tails, and have been described using log-normal models <cit.>, which could indicate evidence of an underlying multiplicative physical process. Meyer et al., 2019 fit the flux distributions of six bright FSRQs with a broken power law, though a log-normal distribution was also compatible with their data, and recently, Tavecchio et al., 2020 have described the gamma-ray flux variability of those same objects using a model based on a stochastic differential equation (SDE) including both deterministic and stochastic components.

The physical structure and multiwavelength emission mechanisms of a blazar can be further understood by modeling its SED.

In leptonic models, the gamma-ray SED component is explained by relativistic electrons and positrons scattering via the inverse Compton process off of a population of lower-energy seed photons, which may be their own emitted synchrotron photons, as in the synchrotron self-Compton process <cit.>, or radiation from an external structure, as in the external inverse Compton process <cit.>. The EIC seed photons are commonly taken to be radiation fields in the BLR, although this picture has been challenged by the lack of characteristic BLR absorption features in the average gamma-ray spectra of Fermi-LAT FSRQs [Costamante et al., 2018].

In hadronic models, however, some or all of the gamma-ray emission is due to relativistic protons emitting via photohadronic processes, proton synchrotron radiation, or other mechanisms, so that relativistic neutrino emission may occur as well. For example, the blazar high-energy emission may be dominated by synchrotron radiation losses of high-energy protons <cit.>. Alternatively, neutrinos may be produced by the photohadronic interaction of a proton with a photon, producing pions that quickly decay to gamma rays and neutrinos, that is, $p\gamma \rightarrow \Delta^+ \rightarrow p\pi^0~\mathrm{or}~n\pi^+$ [Dermer & Menon, 2009]. In this case, production of PeV-scale neutrinos requires a target photon population in the X-ray band. The $p\gamma$ process may co-occur with leptonic gamma-ray emission. Under this scenario, FSRQs may be sources of relativistic neutrinos at PeV or even EeV energies <cit.>. High-energy neutrinos have been detected coming from a direction compatible with the blazar TXS 0506+056 [IceCube Collaboration et al., 2018], which may be an FSRQ masquerading as a BL Lac object [Padovani et al., 2019].

In this paper, we investigate strong gamma-ray flares from three FSRQs at intermediate redshifts. These three sources were continuously monitored by Fermi-LAT (Section <ref>) during the ten-year period from 2008 to 2018, and observed during periods of high gamma-ray activity by the Very Energetic Radiation Imaging Telescope Array System (VERITAS;  Section <ref>). Table <ref> provides an overview of the gamma-ray data analyzed in this work.

3C 279, at $z=0.536$ [Lynds et al., 1965], is one of the most well-studied blazars. It is among the brightest and most variable extragalactic objects in the gamma-ray sky, giving rise to one of the first large amplitude gamma-ray flares measured by EGRET in 1996 [Wehrle et al., 1998]. In recent times, it underwent multiple bright gamma-ray flares in 2014, 2015, and 2018. Notably, during a flare beginning on June 16, 2015, it was detected by High Energy Stereoscopic System (H.E.S.S.), and Fermi-LAT observed minute-scale variability [Romoli et al., 2017, Ackermann et al., 2016]. H.E.S.S. again detected 3C 279 during the flaring states in January and June 2018 [Emery et al., 2019].

PKS 1222+216, at $z=0.432$ [Osterbrock & Pogge, 1987] and also known as 4C +21.35, has exhibited periods of extreme variability in the VHE gamma-ray band, with VHE detections occurring during gamma-ray flares in June 2010 [Aleksić et al., 2011] and February and March 2014 [Holder, 2014].

Finally, Ton 599, at $z=0.725$ (Schneider et al., 2010; see also Burbidge, 1968) and also known as 4C +29.45 and B1156+295, entered a months-long GeV high state in October 2017 [Cheung et al., 2017], leading to VHE detections on the nights of December 15 and 16 2017 [Mirzoyan, 2017, Mukherjee, 2017].

We describe the observations and data analysis of these sources in Section <ref>. In Section <ref>, we examine the long-term variability of these FSRQs and connect it to processes in the accretion disk. Next, we select gamma-ray flares (Section <ref>) and analyze the short-timescale variability (Section <ref>) and spectra (Section <ref>) during these events, focusing primarily on 3C 279, the brightest of the three sources in Fermi-LAT. In Section <ref>, we model the SEDs of PKS 1222+216 and Ton 599 during their respective VHE detections by VERITAS. The observed VHE emission places constraints on the Doppler factor and gamma-ray emission region location during these flares, which we confirm using an independent method in Section <ref>. In Section <ref>, we place theoretical constraints on the potential production of PeV-scale neutrinos during these VHE flares. We summarize our conclusions in Section <ref>. Throughout this paper, a flat $\Lambda$CDM cosmology was used, with $H_0 = 69$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_M = 0.286$, and $\Omega_\Lambda = 0.714$ [Bennett et al., 2014].

§ OBSERVATIONS AND DATA ANALYSIS

Overview of the datasets presented in this paper.

2c| 0pt4ex     5c|Fermi-LAT (HE gamma-ray)     3cVERITAS (VHE gamma-ray)

Source     1c|z
Date Range     11cmEnergy Range     11cmTime Binning     10.85cmNo. Bins     1c|12.1cmFlare Threshold (No. Flares)
11.15cmEnergy Threshold*          Exposure            10.75cmNo. Obs.

2c|     5c|

2c|     [UT]      [GeV]     [day]          1c|[ph cm$^{-2}$ s$^{-1}$]     [GeV]     [hr]
3C 279                0.5362        2008-08-04 – 2018-12-07         0.1-500       1         3471        $4\times10^{-6}$ (10)      200       54.4       139

PKS 1222+216          0.432      2008-08-04 – 2018-12-07         0.1-500       3        1158     $5\times10^{-7}$  (11)     110       34.7       95

Ton 599               0.725      2008-08-04 – 2018-12-12         0.1-500       7        512       $5\times10^{-7}$ (5)     140       8.8        20

*The energy threshold varies for different observations. A typical value is quoted for 3C 279 and the values during the VHE-detected flares are quoted for PKS 1222+216 and Ton 599.

§.§ VERITAS

VERITAS is an array of four imaging atmospheric Cherenkov telescopes located in southern Arizona <cit.>.

VERITAS preferentially performs observations of FSRQs when they exhibit an elevated flux in other wavebands, as a flare at TeV energies might also be occurring. The VERITAS observations of 3C 279, PKS 1222+216 and Ton 599 that were simultaneous with the HE flares considered here were taken in response to the elevated fluxes reported by Fermi-LAT. VERITAS also carries out short monitoring observations of FSRQs. Because these sources are not believed to be strong emitters of TeV gamma rays except during flares, the primary aim of this monitoring is to self-trigger on serendipitous flares. For 3C 279 and PKS 1222+216, these observations provide VERITAS data corresponding to low states observed by Fermi-LAT.

The total exposure of the VERITAS observations on each of the sources is reported in Table <ref>. The data were analyzed using a standard VERITAS data analysis package [Maier & Holder, 2017] and cross-checked using an independent package [Cogan, 2007]. Boosted decision trees with soft selection cuts (appropriate for sources with a photon spectral index softer than $\Gamma \approx 3.5$) were used for separating gamma rays from background cosmic rays [Krause et al., 2017]. Preliminary analysis results of the VERITAS observations of 3C 279 and PKS 1222+216 in 2013 and 2014 were reported by [Errando, 2014]. These are superseded by the more updated analyses reported here.

§.§ Fermi-LAT

Fermi-LAT detects gamma rays from 20 MeV to above 500 GeV using a pair-conversion technique [Atwood et al., 2009]. Fermi-LAT primarily operates in survey mode, during which it scans the entire sky every three hours.

We analyzed the Pass 8 data [Atwood et al., 2013, Bruel et al., 2018] in the 10.3 year period starting on August 4, 2008 (MJD 54682.7), the start of the Fermi-LAT all-sky survey, as reported in Table <ref>. We performed an unbinned likelihood analysis of the data using the LAT  1.0.3 and instrument response functions . The energy range from 0.1 GeV to 500 GeV was analyzed, and photons with zenith angle $> 90\degree$ were excluded to reduce contributions from the Earth's limb. For each source, the region of interest (ROI) considered was the circle of radius $10\degree$ surrounding the catalog source position. The background model consisted of, along with galactic () and isotropic () diffuse emission models, all sources in the FL8Y catalog[<https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/>] within a $20\degree$ circle surrounding the source. This is to ensure that the model would include gamma-ray emission from sources outside the ROI that could extend into the ROI due to the size of the point spread function of the LAT, especially at low energies. We excluded time ranges corresponding to solar flares and gamma-ray bursts in the ROI from the analysis.

When performing the likelihood fit, we iteratively fixed the parameters of the least significant sources (in increasing square powers of natural numbers up to TS equal to 25) until convergence was reached. Sources with TS values outside the allowed range, usually associated with flux parameter values close to zero, were removed from the model.

When fitting individual light curve and SED points, the spectral parameters were kept fixed, either to their catalog values for global analyses or to values derived from an analysis of the full flare period for flare analyses, with the diffuse background model normalization parameters left free.
We checked that the background model we used is consistent with the 4FGL-DR2 catalog [Abdollahi et al., 2020, Ballet et al., 2020] by comparing these two catalogs, finding no new bright, variable sources in the ROI of each of the three FSRQs that could significantly impact the analysis of our sources.

Since 3C 279 lies close to the ecliptic, the Sun and Moon contribute diffuse foreground emission in the ROI of this source during certain periods. This is at the level of $\sim 0.5\times 10^{-6}$ cm$^{-2}\,$s$^{-1}$ for the Sun [Abdo et al., 2011] and $\sim 1\times 10^{-6}$ cm$^{-2}\,$s$^{-1}$ for the Moon [Abdo et al., 2012]. The Sun's quiescent gamma-ray emission extends over a $20\degree$ radius, so this emission is partially degenerate with the diffuse backgrounds modeled in the likelihood fit. The Moon moves about $13\degree$ per day, so it appears within a time bin of a day as a strip. No contamination is expected during any of the flare states identified in Section <ref>, since both the Sun and Moon were more than $20\degree$ from 3C 279 during these periods.

To our knowledge, the maximum level of contamination due to the proximity of 3C 279 to the ecliptic has not previously been quantified. To do so, we generated extended templates for the Sun and the Moon for a 1-day time bin containing the closest approach of the Sun to 3C 279 during the period considered for our analysis [Johannesson & Orlando, 2013], that is, the bin containing October 9, 2018 when its annual occultation occurred [Barbiellini et al., 2014].
During this time bin, the Sun reached a distance of $\approx 0.2^{\circ}$ from 3C 279. By coincidence, the Moon passed within $5\degree$ of the source during the same interval. The templates accounted for the expected extended emission from the Sun and Moon during a 1-day time bin. When we include these templates in the model file for the likelihood analysis for this bin we find that the flux of 3C 279 decreases by approximately $28\%$ with respect to that obtained when only the point sources in the ROI and the galactic and isotropic diffuse backgrounds are included.

The gamma-ray emission from the quiescent Sun and Moon is expected to vary with the solar cycle. In order to estimate the worst-case contamination, we also chose a selection of time bins during which only the Moon was present. When its diffuse template is included in the analysis, this results in a decrease of the 3C 279 flux by up to 49% with respect to when the template is omitted.

For time bins in which the Sun or Moon is more than 5$^\circ$ from 3C 279, the flux of 3C 279, as returned by the likelihood analysis, does not change significantly when the solar and lunar templates are included. We conclude, therefore, that these contributions show no evidence of being statistically significant when deriving the spectral properties of 3C 279 for the time periods studied in this work.

The Sun and Moon each come within $5\degree$ of 3C 279 for approximately 11-13 days per year. A more complete treatment, beyond the scope of this paper, could include the Sun and Moon as extended sources in the likelihood fits for these time bins to fully account for their emission.

§.§ Swift-XRT

The X-Ray Telescope (XRT) on the Neil Gehrels Swift observatory is a grazing-incidence focusing X-ray telescope, and is sensitive to photons with energies between 0.2 and 10 keV [Gehrels et al., 2004, Burrows et al., 2005]. Swift-XRT observed PKS 1222+216 and Ton 599 during the VHE flares of those sources.

The Swift-XRT data were extracted from the Swift data archive and analyzed using . The fluxes and flux errors were deabsorbed using the fixed total column density of Galactic hydrogen $N_H = 2.29 \times10^{20} \;\text{cm}^{-2}$ for PKS 1222+216 and $1.89 \times10^{20} \;\text{cm}^{-2}$ for Ton 599 [Kalberla et al., 2005, Willingale et al., 2013]
and the photoelectric cross section $\sigma(E)$ to account for the effects of neutral hydrogen absorption. The deabsorbed X-ray spectrum was fitted with a broken power law model for PKS 1222+216 and a power law model for Ton 599.

§.§ Swift-UVOT

The ultraviolet/optical telescope (UVOT) on the Neil Gehrels Swift observatory is a photon counting telescope sensitive to photons with energies ranging from about 1.9 to 7.3 eV or 170 to 550 nm [Roming et al., 2005].
Swift-UVOT observed PKS 1222+216 and Ton 599 approximately concurrently with Swift-XRT.

The UVOT data were extracted from the Swift data archive and analyzed using . The counts from the sources and the background were extracted from regions of a radius of $5.0''$ centered on the position of the sources and nearby positions without any bright sources, respectively. The magnitude values of the sources were computed using , and converted to fluxes using the zero-points given by Poole et al., 2008. Extinction corrections were applied following Roming et al., 2009, using the reddening values $E(B-V)=$ 0.0199 and 0.0171 [Schlafly & Finkbeiner, 2011] for PKS 1222+216 and Ton 599, respectively.

§.§ Steward Observatory

During the first decade of the Fermi mission, the Steward Observatory of the University of Arizona obtained optical polarimetry, photometry, and spectra of the LAT-monitored blazars and Fermi targets of opportunity (ToOs) using the SPOL CCD Imaging/Spectropolarimeter [Smith et al., 2009]. We downloaded the spectrophotometric Johnson V and R band magnitudes from the Steward Observatory public archive[<http://james.as.arizona.edu/ psmith/Fermi/>]. These magnitudes were obtained by convolving the flux spectra between 4000 and 7550 Å with a synthetic filter bandpass for the V or R band, summing the flux, and computing the magnitude difference with a comparison star. Smith et al., 2009 give the full details of the observations and data reduction. We then converted the magnitude for each bandpass to its equivalent energy flux. Six observations were taken of Ton 599 and two of PKS 1222+216 during their respective VHE flares. There was no significant variability during either event.

§ FERMI-LAT FLUX DISTRIBUTIONS

The Fermi-LAT light curves of the three sources and the periods of the VERITAS observations are shown in Figure <ref>. The LAT time binnings, reported in Table <ref>, were chosen for each source depending on its typical strength to avoid having an excessive number of bins with no detection.

Fermi-LAT light curves of 3C 279 (top), PKS 1222+216 (middle), and Ton 599 (bottom). The flux points (black circles) are shown for 1, 3, and 7 day time bins for the three sources, respectively. $5\,\sigma$ Bayesian blocks are shown with blue lines. The time intervals in which VERITAS observed the sources are marked in magenta. For 3C 279, time intervals in which the Sun is less than $20^{\circ}$ from the source are shown in orange and a Fermi-LAT low state from MJD 56230 – 56465 (see Section <ref>) is marked with a dashed line.

Flux distributions of the three FSRQs, scaled as probability densities. The distributions are fit with a log-normal PDF (dashed orange) and the stationary-state PDF corresponding to the SDE of Tavecchio et al., 2020 (solid blue). In all three cases, the SDE provides a better fit.

The distribution of the LAT fluxes observed from each of these FSRQs may provide a clue to the origin of the gamma-ray emission. The observed flux distributions of the three sources (scaled to form probability density histograms) are shown in Figure <ref>. Time bins that have a test statistic (TS) less than 9 or that occur when the Sun is less than $20^{\circ}$ from the source were excluded.

To account for uncertainties from both the flux binning and the finite observation length, the flux histogram bin errors were calculated using a bootstrapping approach. 2,500 bootstrap samples were used, each consisting of the same number of flux points as the actual light curve. Each bootstrap sample was obtained by sampling from the set of actual flux points with replacement, so that a given flux point might be sampled multiple times or not at all.

In order to include the uncertainties of the individual flux points, an error term was added to each sampled point in each bootstrap sample, determined by randomly sampling from a Gaussian distribution with standard deviation equal to the measurement uncertainty of the respective sampled point. The bin errors were then defined as the standard deviations of the bin fluxes over all of the bootstrap samples binned using the same bins as the original dataset.

One form of flux distribution often used to describe blazars is the log-normal distribution [Giebels & Degrange, 2009, Sinha et al., 2017, Shah et al., 2018]. Log-normal distributions are of interest because they indicate the presence of an underlying multiplicative rather than additive physical process [Aitchison & Brown, 1973]. Light curves with a log-normal flux distribution have an amplitude of variability linearly proportional to their mean flux.

On the other hand, Tavecchio et al., 2020 have proposed an alternative model based on an SDE with two terms modeling a deterministic tendency to return to equilibrium and stochastic fluctuations with amplitude proportional to the absolute flux level. The form of the SDE is motivated by an astrophysical scenario of stochastic disturbances perturbing a magnetically arrested accretion disk. In this model, the flux distribution is asymmetrical about a peak, falling off as a power law at high fluxes and exponentially at low fluxes, with the relative importance of the deterministic and stochastic components dictating the shape of the distribution. Figure <ref> shows a comparison between the best-fit probability density functions (PDFs) corresponding to a log-normal distribution and the stationary state of the SDE proposed by [Tavecchio et al., 2020].

Best-fit parameters and goodness of fit ($\chi^2_{\rm{red}}$) for the log-normal and SDE PDF fits to the LAT flux distributions. The mean $\mu$ and standard deviation $\sigma$ of the exponentiation of the log-normal distribution and the equilibrium flux $\mu$ of the SDE distribution are normalized to $1\times10^{-7}$ ph cm$^{-2}$ s$^{-2}$. The peak flux of the SDE PDF is given by $X_\mathrm{max} = \mu \lambda / (\lambda + 2)$.

Source     3cLog-normal

$\mu$     $\sigma$     $\chi^2_{\rm{red}}$     $\mu$     $\lambda$     $X_\mathrm{max}$     $\chi^2_{\rm{red}}$
3C 279       1.65$\pm$0.02        0.73$\pm$0.01        14.7      8.63$\pm$0.28      1.01$\pm$0.06     2.90$\pm$0.21     2.95

PKS 1222+216        1.07$\pm$0.03        0.92$\pm$0.03        4.44       6.27$\pm$0.60      0.54$\pm$0.08     1.33$\pm$0.23      1.75

Ton 599       0.21$\pm$0.04        0.75$\pm$0.03        1.76       2.11$\pm$0.19      0.94$\pm$0.15     0.68$\pm$0.13     0.73

The stationary-state PDF corresponding to the SDE <cit.> is

\begin{equation}
p(X) = \frac{(\lambda \mu)^{1 + \lambda}}{\Gamma(1 + \lambda)}\frac{e^{-\lambda \mu / X}}{X^{\lambda + 2}},
\end{equation}

where $X$ is a dimensionless random variable proportional to the flux, $\mu$ is a parameter representing the equilibrium value of $X$, $\lambda$ is a parameter representing the relative weight of the deterministic and stochastic terms, and $\Gamma$ is the gamma function. Here, $X$ was related to the flux by a proportionality constant of $1\times10^{-7}$ ph cm$^{-2}$ s$^{-1}$. The distribution peaks at $X_\mathrm{max} = \mu \lambda / (\lambda + 2)$. The stationary-state PDF is valid on timescales much longer than the timescale for the system to return to equilibrium, which is clearly the case for the ten-year periods considered here.

The PDFs were fit to the histogram bins using a nonlinear least-squares algorithm. The best-fit parameters and reduced $\chi^2$ values of the two models are reported in Table <ref>. In all three cases, the SDE PDF provides a better fit than the log-normal PDF. Both models have two free parameters. We verified that the preference for the SDE model is preserved if the histogram bins at the lowest fluxes, which might be affected by requiring light curve bins to have TS $>9$, are excluded from the fit.

The SDE model PDF is parameterized by the shape parameter $\lambda \equiv 2\theta/\sigma^2$, where $\theta$ and $\sigma$ are the coefficients of the deterministic and stochastic terms. These parameters can be interpreted by associating $1/\theta$ with the timescale of magnetic field accumulation in the accretion disk, while $\sigma$ is related to the dynamics of the perturbative processes. A large value of $\lambda$ therefore represents a high relative importance of the deterministic variability component compared to the stochastic one, while a small value indicates the opposite. To relate these timescales to the gravitational radii of the central supermassive black holes, $r_g = GM/c^2$, we adopt values of ${\sim}5\times10^8$, $6\times10^8$, and $3.5\times10^8 M_\odot$ for the black hole masses of 3C 279, PKS 1222+216, and Ton 599, respectively [Hayashida et al., 2015, Farina et al., 2012, Liu et al., 2006].

One can estimate $\sigma^2$ from the light curve using the expression [Tavecchio et al., 2020]:

\begin{equation}
\sigma^2 \simeq \frac{1}{n} \sum_{i=0}^{n} \frac{(X_i - X_{i - 1})^2}{X_{i - 1}^2 (t_i - t_{i - 1})},
\end{equation}

where $X_i$ is the scaled flux at time step $i$. Using this expression, we obtain $\sigma^2$ equal to 0.35, 0.16, and 0.062 day$^{-1}$, or 100, 200, and 800 $r_g / c$, for 3C 279, PKS 1222+216, and Ton 599, respectively. These values are consistent with the ${\gtrsim}100~r_g / c$ variability timescale injected into the jet by magneto-rotational instability in the accretion disk estimated in theoretical work [Giannios & Uzdensky, 2019]. Using the relation $1/\theta = 2/\lambda \sigma^2$, we can then constrain the physics of accretion flow in 3C 279, PKS 1222+216, and Ton 599 by estimating their magnetic flux accumulation timescales to be 200, 700, and 1800 $r_g / c$, respectively, within the magnetically-arrested disk scenario.

§ FLARE SELECTION

Flare states were identified in the Fermi-LAT data using the following procedure:

* Segment the data using Bayesian blocks. We set the false positive rate $p_0$ to the value equivalent to $5\,\sigma$ using Equation 13 of Scargle et al., 2013.

* Choose a flux threshold above which the blocks are designated as flaring.

* Designate each contiguous set of flare blocks as an individual flare state and all non-flare blocks as the quiescent state.

This empirical procedure reflects a picture of individual flares superimposed on a constant quiescent background, but identifies them purely as states of elevated flux, making no explicit assumptions about the flares' shape or spectra. Due to its basis on Bayesian blocks, it guarantees that states identified as flares have flux significantly greater than the states surrounding them.

The flux threshold to identify flares must be tuned on a source-by-source basis. Choosing the flux threshold to identify flares involves a trade-off between ensuring that dimmer flares are selected and avoiding misidentifying fluctuations in the quiescent background as flares. In addition, because the sources differ in average flux, the threshold must necessarily vary on an absolute level from source to source. Performing the flare selection procedure with the flare selection thresholds listed in Table <ref> results in 10 flares selected for 3C 279, 11 for PKS 1222+216, and 5 for Ton 599, listed in Table <ref>. We set the threshold low enough for each source to ensure that all flares that triggered VERITAS observations were selected.

Fermi-LAT flares selected using the algorithm given in Section <ref>. For each enumerated flare, the date range in MJD, approximate calendar date, number of Bayesian blocks, and amount of VHE gamma-ray exposure taken by VERITAS (if any) are provided. All of the times in the date ranges given in Table <ref> but not listed here are considered to be quiescent.

#     Date Range (MJD)     Approx. Date     Blocks     VHE Exp.
5c3C 279

1     56645.66 - 56647.66     Dec 2013     1     -

2     56717.66 - 56718.66     Mar 2014     1     -

3     56749.66 - 56754.66     Apr 2014     1     6.79 hr

4     57186.66 - 57190.66     Jun 2015     3     1.00 hr

5     58116.66 - 58119.66     Dec 2017     1     -

6     58130.66 - 58141.66     Jan 2018     4     1.38 hr

7     58168.66 - 58173.66     Feb 2018     1     -

8     58222.66 - 58230.66     Apr 2018     5     0.83 hr

9     58239.66 - 58247.66     May 2018     1     -

10     58268.66 - 58275.66     Jun 2018     2     3.95 hr

5cPKS 1222+216

1     55096.66 - 55114.66     Sep-Oct 2009     3     -

2     55144.66 - 55201.66     Nov-Dec 2009     5     -

3     55231.66 - 55594.66     2010     27     -

4     55603.66 - 55639.66     Feb-Mar 2011     4     5.38 hr

5     55777.66 - 55783.66     Aug 2011     1     -

6     56494.66 - 56500.66     Jul 2013     1     -

7     56536.66 - 56665.66     Sep 2013     5     -

8     56680.66 - 56752.66     Jan-Apr 2014     3     15.53 hr

9     56926.66 - 57004.66     Sep-Dec 2014     5     -

10     58243.66 - 58249.66     May 2018     1     -

11     58321.66 - 58327.66     Jul 2018     1     -

5cTon 599

1     55417.66 - 55445.66     Aug-Sep 2010     1     -

2     57342.66 - 57356.66     Nov 2015     1     -

3     57944.66 - 57958.66     Jul 2017     1     -

4     58042.66 - 58140.66     Oct 2017 - Jan 2018     5     8.30 hr

5     58217.66 - 58266.66     Apr-May 2018     1     2.00 hr

Because the flux distributions are best fit by the single-component SDE model PDF, it is not natural to calculate a duty cycle of flares based on a division into baseline and flaring components <cit.>. The amount of time spent in the highest-flux states can be estimated directly from the flux distribution by defining the “typical" flux as the peak of the PDF, given in Table <ref>. 3C 279, PKS 1222+216, and Ton 599 have flux greater than 5 (10) times the typical flux 12% (4%), 19% (8%), and 13% (4%) of the time, respectively.

Our flare selection flux thresholds for 3C 279 and Ton 599 are comparable at 13.8 and 11.8 times their typical fluxes, consistent with their similar values of the PDF shape parameter $\lambda \approx 1$. For PKS 1222+216, our threshold is 3.8 times the typical flux. This source has a lower value of $\lambda \approx 0.5$, with a correspondingly harder power law of the flux distribution at high fluxes. This is perhaps reflected in the long epochs of high flux seen in this source's light curve, such as its Flare 3 in 2010 which is approximately a year in duration (Table <ref>). A relatively low threshold was therefore needed to also capture the smaller flares of the approximately weekly timescales that typically trigger VERITAS observations, consistent with the other two sources.

§ DAILY AND SUB-DAILY VARIABILITY

In order to deduce the smallest variability time around the rising and decaying periods of each flare selected according to the algorithm described in Section <ref>, we extracted sub-daily light curves of the three sources in time bins ranging from 12 hours down to 1.5 hours for the brightest source, 3C 279. Starting from daily time bins, we refined the light curve iteratively by splitting the time bin duration until each bin had a TS of $\gtrsim50$ or until further refinement would not change the local trend of the light curve. For PKS 1222+216 and Ton 599, the minimum bin sizes were 12 and 6 hours, respectively.

To characterize the flares with multiple peaks we used a sum of exponential profiles [Valtaoja et al., 1999, Abdo et al., 2010], $F_i$, where each one has the form:

\begin{equation} \label{eq:flareprofile}
F_i(t) = \begin{cases}
F_{0_i}\, e^{(t-t_{\rm{peak}_i})/t_{\rm{rise}_i}}, & t\leq t_{\rm{peak}_i}\\
F_{0_i}\, e^{-(t-t_{\rm{peak}_i})/t_{\rm{decay}_i}}, & t> t_{\rm{peak}_i}.
\end{cases}
\end{equation}

For flares with a single peak we used:

\begin{equation} \label{eq:flareprofileconst}
F(t) = \begin{cases}
F_0\, e^{(t-t_{\rm{peak}})/t_{\rm{rise}}}+F_\mathrm{const}, & t\leq t_{\rm{peak}}\\
F_0\, e^{-(t-t_{\rm{peak}})/t_{\rm{decay}}}+F_\mathrm{const}, & t> t_{\rm{peak}},
\end{cases}
\end{equation}

including a constant term to avoid having a bias towards large rise and decay timescales, which is minimal when multiple peaks are included.

The fitting procedure started by considering a single peak characterized by Equation (<ref>). In order to evaluate the possibility of adding a second peak, a fit to the sum of two exponential profiles, as given by Equation (<ref>), was performed and compared against the one-peak scenario using the reduced $\chi^2$ method. The two-peak model was then taken when an improvement was observed over the one-peak function. More peaks were then added following a similar procedure until a reasonable reduced $\chi^2$ value was reached, or when the best fit values obtained no longer provided relevant information for constraining the variability timescales of the sources under study. At each step, human judgment was used to initialize the profile positions and determine by eye that the fits made sense. The peaks were not required to match the Bayesian blocks used for flare selection, which were defined using the coarsely binned light curves.

The flare profiles of the three sources are shown in Figure <ref>. Two selected flares of 3C 279 are shown, as are the two flares each of PKS 1222+216 and Ton 599 that were observed by VERITAS. Profiles of all ten flares of 3C 279 are provided in Appendix <ref>. In order to illustrate when VERITAS observed the source relative to the LAT flare peaks, the VERITAS daily-binned light curves for each of the flares are also shown in Figure <ref>.

LAT daily and sub-daily light curves (blue points) around selected flaring episodes (light shaded areas). The dotted blue lines show the fitted exponential profiles, with their sums shown in solid blue. The dark shaded areas indicate the periods considered for the SED modeling (Section <ref>). The VERITAS data points and 95% upper limits are shown as black squares and downwards arrows.

The fit results for the three sources are reported in Tables <ref>, <ref>, and <ref>. For 3C 279, each flare lasted between one and eleven days and consisted of between one and four separately resolved components, modeled using exponential profiles. Twenty-four distinct components are resolved within the ten flares. The rise and decay times range from timescales of days to less than one hour. The smallest resolved variability timescale was $36\pm13$ minutes, which occurred around MJD 58227.945, during the rising period of Flare 8 (MJD 58222.655 –  58230.655), indicated in boldface in Table <ref> and Table <ref>.

Results of the LAT flare profile fits for 3C 279, with flare timescales in minutes.

Amplitude ($F_0$)     $t_{\rm{peak}}$      $t_{\rm{rise}}$     $t_{\rm{decay}}$

($10^{-9} $ erg cm$^{-2}\,$s$^{-1}$)        (MJD)       (min)      (min)

5cFlare 4 (MJD 57186.655 – 57190.655): $\chi^2/$d.o.f.= 77.31/19 = 4.07
12.07 $\pm$ 0.67     57187.446 $\pm$ 0.031     378 $\pm$ 46     1784 $\pm$ 147

9.79 $\pm$ 2.29     57188.425 $\pm$ 0.028     216 $\pm$ 101     155 $\pm$ 64

21.72 $\pm$ 1.59     57189.069 $\pm$ 0.008     137 $\pm$ 18     512 $\pm$ 55

12.41 $\pm$ 1.30     57189.532 $\pm$ 0.010     220 $\pm$ 63     77 $\pm$ 25
5cFlare 8 (MJD 58222.655 – 58230.655): $\chi^2/$d.o.f.= 177.25/106 = 1.67
5.29 $\pm$ 1.29     58224.773 $\pm$ 0.105     1996 $\pm$ 716     5899 $\pm$ 4035

17.70 $\pm$ 2.01     58227.945 $\pm$ 0.004     36 $\mathbf{\pm}$ 13     329 $\pm$ 131

16.42 $\pm$ 1.87     58228.323 $\pm$ 0.012     140 $\pm$ 54     115 $\pm$ 48

5.59 $\pm$ 1.69     58227.139 $\pm$ 0.133     3816 $\pm$ 1450     4077 $\pm$ 2080
The smallest variability time found is indicated in boldface.

Results of the LAT flare profile fits for PKS 1222+216 and Ton 599, with flare timescales in days.

Amplitude ($F_0$)     $t_{\rm{peak}}$      $t_{\rm{rise}}$     $t_{\rm{decay}}$

($10^{-9} $ erg cm$^{-2}\,$s$^{-1}$)        (MJD)       (days)      (days)

5cPKS 1222+216
5cFlare 4 (MJD 55603.7 – 55639.7): $\chi^2/$d.o.f.= 102.25/69 = 1.48
0.56 $\pm$ 0.09     55607.1 $\pm$ 0.3     1.5 $\pm$ 0.5     2.7 $\pm$ 0.8

0.48 $\pm$ 0.03     55629.9 $\pm$ 1.2     22.3 $\pm$ 3.5     12.1 $\pm$ 2.2
5cFlare 8 (MJD 56680.7 – 56752.7): $\chi^2/$d.o.f.= 166.40/104 = 1.60
0.72 $\pm$ 0.09     56692.9 $\pm$ 0.1     13.6 $\pm$ 2.7     0.4 $\pm$ 0.3

1.75 $\pm$ 0.17     56702.8 $\pm$ 0.2     3.6 $\pm$ 0.6     1.5 $\pm$ 0.3
0.43 $\pm$ 0.05     56721.9 $\pm$ 1.6     19.9 $\pm$ 8.5     10.4 $\pm$ 6.2

0.41 $\pm$ 0.10     56732.0 $\pm$ 0.4     0.9 $\pm$ 0.8     8.6 $\pm$ 3.6

0.44 $\pm$ 0.06     56746.9 $\pm$ 0.1     0.3 $\pm$ 0.3     10.3 $\pm$ 2.5
5cTon 599
5cFlare 4 (MJD 58042.7 – 58140.7): $\chi^2/$d.o.f.= 456.78/296 = 1.54
1.89 $\pm$ 0.29     58057.1 $\pm$ 0.2     1.1 $\pm$ 0.3     1.9 $\pm$ 0.6

1.06 $\pm$ 0.11     58065.4 $\pm$ 1.0     11.9 $\pm$ 2.6     9.0 $\pm$ 2.1
1.37 $\pm$ 0.06     58103.5 $\pm$ 0.8     47.0 $\pm$ 4.7     11.8 $\pm$ 1.1
5cFlare 5 (MJD 58217.7 – 58266.7): $\chi^2/$d.o.f.= 153.01/96 = 1.59
0.57 $\pm$ 0.12     58219.2 $\pm$ 1.3     3.5 $\pm$ 2.8     7.3 $\pm$ 2.1

0.48 $\pm$ 0.03     58246.3 $\pm$ 0.7     6.6 $\pm$ 1.7     34.9 $\pm$ 5.6

The flare components coincident with VHE flares are marked with a *, with corresponding smallest variability times indicated in boldface.

For PKS 1222+216 and Ton 599, the variability timescales were of the order of days. Notably, for both sources, the fastest variability did not occur during the detected VHE flares. The shortest variability timescale observed by LAT during the VHE flare of PKS 1222+216 was $10.4 \pm 6.2$ days, which was the decay timescale of the coincident flare component. The shortest variability timescale of Ton 599 observed by LAT during its VHE flare was $11.8 \pm 1.1$ days, which also was the coincident flare component's decay timescale. In the case of Ton 599, the VERITAS detection occurred over a period of 2 days, after which the observed VHE flux became insignificant. No significant intra-flare variability was observed by Fermi-LAT or VERITAS during either event. Therefore, in the remainder of this work, we take the most constraining variability timescales during the VHE flares of PKS 1222+216 and Ton 599 to be 10 and 2 days, respectively.

Decay time vs. rise time for each of the resolved exponential components in the flares of 3C 279. The points corresponding to all of the components are shown in gray. The dashed reference line shows where the rise and decay times are equal.

The symmetry or asymmetry of flares can provide information on the timescales of the particle acceleration and cooling processes in the emission region <cit.>. If the cooling time is longer than the light travel time through the emission region, the decay time will be longer than the rise time, producing an asymmetric flare. If the cooling time is shorter than the light travel time, the flare will appear more symmetrical. Flares with a slow rise and fast decay may be produced by relativistic magnetic reconnection [Petropoulou et al., 2016].

Fluence distributions of the twenty-four resolved flare components and ten flares of 3C 279.

Figure <ref> shows the fitted rise and decay times for each of the exponential flare components of 3C 279. No clear trend in the flare asymmetry is observable, whether overall, among components within a single flare, or between the components belonging to different flares. Both longer decay times and longer rise times are observed, and many flares appear symmetric. A Wilcoxon signed-rank test [Wilcoxon, 1945] finds no significant preference ($p=0.178$) for flares to have a faster rise time than decay time rather than the reverse. These findings are consistent with previous studies of gamma-ray flares in bright Fermi blazars <cit.>.

Models of blazar flares powered by relativistic reconnection predict that flare components produced by large, non-relativistic plasmoids should have similar fluences to components produced by small, relativistic ones, so that flare components should have similar fluence regardless of their variability timescales [Petropoulou et al., 2016]. Figure <ref> shows the distributions of fluences of the components of the ten flares and the twenty-four individual flare components of 3C 279. The fluence $\mathcal{F}$ of a flare with exponential components $F_i$ is given by:

\begin{equation}
\mathcal{F} = \sum_i F_{0_i} (t_\mathrm{rise} + t_\mathrm{decay}).
\end{equation}

For 3C 279 Flares 1, 2, 5, and 7, the best fit is given by a single component plus a constant baseline flux. In these cases, the baseline flux is included in the fluence estimate for consistency with the other flares, approximating the flare duration as $t_\mathrm{rise} + t_\mathrm{decay}$, so that the fluence is given by:

\begin{equation}
\mathcal{F} = (F_0 + F_\mathrm{const}) (t_\mathrm{rise} + t_\mathrm{decay}).
\end{equation}

The median flare fluence is $2.1\times10^{-3}$ erg cm$^{-2}$ and the median component fluence is $0.85\times10^{-3}$ erg cm$^{-2}$. The observed component fluences range over about one order of magnitude, as do the flare amplitudes, while the rise and decay timescales span about two orders of magnitude. These dynamic ranges are generally compatible with the expectations for plasmoid-powered flares derived from particle-in-cell simulations of relativistic magnetic reconnection [Petropoulou et al., 2016].

Left: Global Fermi-LAT spectra for 3C 279, PKS 1222+216, and Ton 599. The LAT spectra are extrapolated to the VERITAS energy range, incorporating EBL absorption. Right: Two baseline states of 3C 279. The Fermi-LAT spectrum corresponding to the multiwavelength low state (MJD 56230–56465; see Figure <ref>), is shown by filled squares. The strictly simultaneous Fermi-LAT spectrum and VERITAS spectral upper limits during the quiescent state are shown by black filled and unfilled circles and contours. The symbol “$^{\dagger}$" indicates that the LAT spectrum corresponds to data strictly simultaneous with VERITAS observations. Downward arrows show 95% confidence level upper limits.

The long-term gamma-ray variability study of the three FSRQs presented here is compatible with the extensive flare characteristics study carried out recently by [Meyer et al., 2019] on the brightest flares detected by Fermi-LAT. A similar Bayesian blocks analysis was carried out to identify flares and look for variability on sub-hour timescales. Consistent with their findings, we find sub-hour-scale variability in 3C 279, where it was possible to resolve flares in finer time bins, suggesting that extremely compact emission regions may be present within the jet.

§ GAMMA-RAY SPECTRA

Figure <ref> shows the LAT energy spectra corresponding to the entire data sets of each of the three sources, along with the VERITAS spectral upper limits for 3C 279. The best-fit spectral parameters are reported in Appendix <ref>. Since all three sources were best fit by a log-parabola model in the 4FGL catalog [Abdollahi et al., 2020], we fit the LAT spectra with this model, parametrized as

\begin{equation}
\frac{dN}{dE} = N_0 \left( \frac{E}{E_b} \right)^{-(\alpha + \beta(\log(E/E_b))},
\end{equation}

where $E_b$ was fixed to the FL8Y catalog value of 457.698 MeV.

We checked that the log-parabola model provides a better fit than a power-law model using the likelihood ratio test. A power-law sub-exponential cutoff model was also preferred over a power law, but we could not establish whether this model is significantly preferred with respect to the log-parabola model using a likelihood ratio test. This is because the two curved models are non-nested, i.e. neither is a special case of the other, and therefore it is not possible to calculate the statistical significance of a preference for one over the other. We assumed a log-parabola spectrum for all subsequent LAT analyses. To facilitate comparison with the VERITAS points, the LAT model fits and butterfly contours were extended beyond the LAT maximum energy of 500 GeV, and extragalactic background light absorption was applied to them using the model of [Franceschini & Rodighiero, 2017].

The global spectral shapes of the three sources are similar, with an index $\alpha$ of ${\sim}$2.1–2.3 and a curvature parameter $\beta$ of ${\sim}$0.04–0.06, and they differ primarily by their normalization.

Fermi-LAT spectrum and VERITAS spectral upper limits of 3C 279 during four flares, strictly simultaneous with the VERITAS observations. The quiescent spectrum (gray circles and contour) is shown for comparison. The strictly simultaneous quiescent state LAT data and VERITAS upper limits are represented with filled and unfilled markers, respectively. The LAT spectra are extrapolated to the VERITAS energy range, incorporating EBL absorption. Downward arrows show 95% confidence level upper  limits.

Using the data from 3C 279, we compared several methods to determine a baseline non-flaring spectrum. First, we defined a low state lasting from MJD 56230 to 56465 (see Figure <ref>), during which the flux was quiescent and stable in HE gamma rays, $R$-band optical, and X-rays. We checked publicly available Tuorla[<https://users.utu.fi/kani/1m/3C_279_jy.html>] data for the $R$-band light curve. For the X-rays, we analyzed the Swift-XRT light curve using the online data products generator[<https://www.swift.ac.uk/user_objects/>]. To ensure low, stable gamma-ray emission, we selected the interval to span the Bayesian blocks with the lowest flux while excluding intervals with the Sun in the ROI. The low-state LAT SED is shown in Figure <ref>. Only one VERITAS observation occurred during this interval, on MJD 56417, so the corresponding VERITAS upper limits are not constraining and are not shown.

Fermi-LAT spectra of 3C 279 during ten flares, for the intervals selected using the algorithm proposed in this work and described in Section <ref>. For comparison, the LAT low state spectrum is shown in gray squares in all of the panels. The four flares shown in color have a corresponding spectrum in Figure <ref>.

Next, using the algorithm proposed in this work and described in Section <ref>, we designated all epochs of the LAT light curve other than the flaring episodes as quiescent. From those epochs, we extracted those LAT data strictly simultaneous with the VERITAS observations, integrating a total of 43.6 hours of observations. The resulting strictly simultaneous LAT spectrum and VERITAS spectral upper limits are shown in Figure <ref>. We then performed the same procedure for four flaring epochs during which a significant Fermi-LAT detection could be obtained strictly simultaneous with the VERITAS observations, which occurred on the nights of April 3, 2014; June 16, 2015; April 19, 2018; and June 3, 2018. These strictly simultaneous LAT and VERITAS SEDs are shown in Figure <ref>.

The spectral shapes of the 3C 279 low and quiescent states are similar to each other and to the global state, although the uncertainties on their fit parameters are high due to the low significance. The spectra differ primarily in their flux normalization. The normalization of the low state is lower than that of the global state by design, while the normalization of the strictly simultaneous quiescent state is higher. This could result from the timing of the VERITAS monitoring and triggered observations which often follow up on Fermi-LAT flares and may tend to catch mildly elevated activity in Fermi-LAT even if the source is not actually flaring.

Finally, we derived LAT SEDs for all of the ten identified flares of 3C 279, using the entire flare time periods, irrespective of strict simultaneity with VERITAS, shown in Figure <ref>. The average flare spectrum is more strongly curved than the global spectrum, with $\alpha = 2.02 \pm 0.01$ and $\beta = 0.093 \pm 0.008$, compared to $\alpha = 2.228 \pm 0.004$ and $\beta = 0.061 \pm 0.003$ for the global state.

§ SED MODELING

Multiwavelength SED modeling can shed light on the mechanisms of gamma-ray production during VHE flares. For 3C 279, we refer the reader to those works in the literature in which multiwavelength SED modeling of the epochs considered here has been performed, and we do not perform any additional modeling <cit.>.

PKS 1222+216 was first detected at TeV energies by MAGIC during a flaring event in June 2010 [Aleksić et al., 2011], and multiwavelength SED modeling of this event has been performed by e.g. [Tavecchio et al., 2011]. We therefore restricted our SED modeling of the source to the duration of the second VHE detection by VERITAS in February and March 2014. We considered data from all instruments taken from UT 2014-02-26 to 2014-03-10, inclusive.

Ton 599 has not been studied as extensively as the other two sources. [Prince, 2019] and [Patel & Chitnis, 2020] model its variability characteristics and multiwavelength SED, respectively, during the high state in December 2017, but do not have access to TeV data. We therefore modeled the multiwavelength SED of Ton 599 during the VERITAS detection in December 2017. We considered data from all instruments taken from UT 2017-12-15 to 2017-12-16, inclusive.

To assemble our multiwavelength SEDs, in addition to the gamma-ray data from VERITAS and Fermi-LAT, we incorporated X-ray and ultraviolet data from the XRT and UVOT instruments aboard the Swift satellite and optical observations from the Steward Observatory.

We described the multiwavelength SEDs of the two FSRQs using a multi-component synchrotron self-Compton (SSC) blob-in-jet model, implemented using the framework of the “" code, developed by [Hervet et al., 2015] and based on [Katarzyński et al., 2001]. We modeled the radiative interactions of a compact leptonic emission zone (a blob), including an EIC emission component resulting from the interactions of the blob particles with the thermal accretion disk emission reprocessed by the BLR. Figure <ref> shows a schematic illustration of the components producing the emission in this model.

(-235,35)Accretion disk
A schematic illustration of the emission model used in this work (not to scale).
The green dashed arc represents the nominal BLR radius $r_\mathrm{BLR}$ corresponding to the region of the maximal BLR density.
The observer measures the beamed emission from the blob interacting with the BLR as well as the accretion disk's thermal emission.

The accretion disk is assumed to be a point source.

We consider a simplified BLR model with a normalized density profile, based on [Nalewajko et al., 2014], where $\rho_\mathrm{BLR}(r)$ is at a maximum at the characteristic BLR radius $r = r_\mathrm{BLR}$ and decreasing as $r^{-2}$ with the distance to the core such that

\begin{equation}
\rho_\mathrm{BLR}(r) = \frac{(r/r_\mathrm{BLR})^2}{1+(r/r_\mathrm{BLR})^4}~,
\end{equation}

with $r_\mathrm{BLR}$ scaled to the bolometric disk luminosity $L_d$ as $r_\mathrm{BLR} = 0.1 \sqrt{L_d/1\times10^{46}~\mathrm{erg}~\mathrm{s}^{-1}}$ pc [Sikora et al., 2009, Ghisellini & Tavecchio, 2009]. From SED modeling of PKS 1222+216 and Ton 599 we deduce a BLR radius of 0.17 pc and 0.15 pc respectively.
We assume an isotropic diffusion of the disk light by the BLR, where the specific intensity of this field can be expressed as

\begin{equation}
I_\mathrm{BLR}(\nu, T_{d}, r) = \epsilon_\mathrm{BLR} \rho_\mathrm{BLR}(r) \frac{L_{d}}{4\pi r^2}\frac{I_{p}(\nu, T_{d})}{(\sigma_\mathrm{SB}/\pi)T_{d}^4}

where $\sigma_\mathrm{SB}$ is the Stefan-Boltzmann constant, $I_{p}$ is the Planck intensity, and $\epsilon_\mathrm{BLR}$ is the covering factor. This equation is similar to Eq. 12 in [Hervet et al., 2015] with the addition of the BLR density profile.
Only the extension of the BLR in front of the blob plays a significant role in our modeling since it drives the number of gamma rays produced by the blob that will be absorbed by pair production.
The BLR is by default defined between $r = 0$ and $r = 100~r_\mathrm{BLR}$. Given the fast convergence of the BLR opacity ($I_\mathrm{BLR} \propto r^{-4}$), the maximum extension of the BLR does not play a significant role in the model. Although we assume for simplicity that the BLR is isotropic, any anisotropy should have a small effect on the opacity <cit.>.

Figures <ref> and <ref> show the multiwavelength SED models of PKS 1222+216 and Ton 599. In these figures, the synchrotron and SSC emission are shown by solid blue lines. The subdominant second-order self-Compton emission caused by the interactions of the electrons with the self-Compton photons is shown by a dotted blue line. The thermal emission from the accretion disk is modeled as a point source radiating as a black body, and is shown by a heavy dashed green line. The inverse Compton emission due to the interaction of the electrons with the disk photons reprocessed in the BLR is shown by a dashed green line. Table <ref> gives the parameters characterizing the SED models.

Our model does not include any secondary radiation from pair cascades produced by the absorption of gamma rays in the BLR. While detailed modeling of this effect is beyond the scope of this paper, we estimate that the potential contribution of such cascades would be $\ll 1\%$ of the total bolometric luminosity for PKS 1222+216 and $\lesssim 1\%$ for Ton 599, given the respective levels of absorption in our models, which are described below. We evaluated these contributions by comparing the radiative output of each model with that from the same model with the BLR opacity set to zero. This effect may be noted as a source of systematic uncertainty when interpreting our results.

We note that our model requires that the dust-torus luminosity be negligible compared to the disk luminosity. As evidence of far-infrared dust-torus thermal emission is lacking in the SED, we consider this assumption to be reasonable in our study. Observing campaigns with good microwave to IR coverage would be needed to fully confirm this approach. The presence of strong dust-torus emission would require that the gamma-ray emission zone be farther downstream in the jet so as not to produce too large an opacity by pair production.

§.§ PKS 1222+216 modeling

Broadband SED of PKS 1222+216 during the VERITAS detection from UT 2014-02-26 to 2014-03-10. Left: Pure synchrotron self-Compton model. Right: Model considering an external inverse Compton component at high energy from the interaction of blob particles with the thermal accretion disk emission reprocessed by the BLR. The solid blue lines show synchrotron and SSC emission; the dotted blue line shows second-order self-Compton emission; the heavy dashed green line shows thermal emission from the accretion disk; and the dashed green line shows inverse Compton emission from the BLR. The EBL absorption is taken into account considering the model of [Franceschini & Rodighiero, 2017].

In order to investigate the necessity of including an EIC component, we represented the multiwavelength SED of PKS 1222+216 with a one-zone pure SSC model, shown in Figure <ref> (left). As can be seen by the similar amplitudes of the synchrotron and inverse Compton peaks in the figure, the SED is only weakly Compton dominated, with the inverse Compton luminosity about twice the synchrotron luminosity. The Swift-XRT spectrum contains a well-resolved break showing the transition between synchrotron and inverse Compton dominated emission, which sets a strong constraint on the model.
Our best attempt does not provide a satisfying representation of the observed SED. The main issue is that the optical-to-X-ray components of the SED have steep slopes which would require a narrow, sharp synchrotron bump to achieve a good representation, while the X-ray-to-VHE needs a wide, flat inverse Compton bump. This is not compatible with the usual simple SSC framework, especially when the SED is not heavily Compton dominated.

In our EIC model, the IR-to-UV SED is dominated by the blackbody big-blue-bump emission of an accretion disk (see Figure <ref>, right), which resolves the tension by eliminating the constraint on the synchrotron spectral shape. This allows for a broad SSC component matching the spectral break observed in the X-ray band. In this scenario, the VHE emission is produced by the EIC process between a relativistic blob and the disk thermal emission reprocessed by the BLR. The blob is set to a distance of 3.56 pc from the SMBH, corresponding to $21.3~r_\mathrm{BLR}$.
It should be noted that a thermal EIC process was also favored in previous models of PKS 1222+216 where clear disk emission and a strongly Compton-dominated SED were observed during a major outburst in 2010 [Tavecchio et al., 2011].

Because the peak frequency of the EIC emission is directly proportional to the blob Lorentz factor, this scenario imposes a strong constraint on the jet parameters. For PKS 1222+216, in order to match the VHE spectrum, the bulk Lorentz factor needs to be above approximately 23, which was achieved by assuming a Doppler factor $\delta = 40$ and an angle with the line of sight $\theta_\mathrm{obs} = 1\degree$.

This assumption is consistent with the jet constraints derived by Hervet et al., 2016 from the fastest motion observed in the radio jet of PKS 1222+216, which led to estimations of $\theta_\mathrm{obs} = 1.3\degree$, $\delta = 41.3$ and $\Gamma = 29.2$.

Because no significant variability was observed in any waveband during the time period selected for modeling for either source, we considered a stationary model giving a snapshot of the observed activity.
As a consistency check, we compared the expected radiative cooling time from the model with the observed flare decay timescale. The cooling time associated with the full radiative output (synchrotron, SSC and EIC emissions) can be expressed in the Thomson regime as

\begin{equation}
T_{\mathrm{cool}} (\gamma) = \frac{3 m_e c}{4 \sigma_T \gamma  (U'_B + U'_\mathrm{syn} + U'_\mathrm{blr})},
\end{equation}

with $m_e$ the electron mass, $\sigma_T$ the Thomson cross section, $\gamma$ the Lorentz factor of the emitting particle, and  $U'_B$, $U'_\mathrm{syn}$, $U'_\mathrm{blr}$ respectively the energy density in the blob frame of the magnetic field, synchrotron field, and external BLR field <cit.>.
One can associate the energy at the break of the spectral particle distribution $\gamma_{\mathrm{brk}}$ with the emission at the peaks of the SED. The Fermi-LAT energy range being mostly above this peak, we can deduce $T_\mathrm{cool}(\mathrm{\textit{Fermi}}) \lesssim 17$ days. This is consistent with the observed Fermi-LAT flare decay of $10.4\pm6.2$ days.

The minimum possible variability predicted by our model is 18 h, given by the blob's radius and Doppler factor such that $\tau_{\mathrm{min}} = R (1+z)/(c \delta)$.

The total power of the jet is approximately $3.4\times 10^{45}$ erg s$^{-1}$, in a particle-dominated regime with the equipartition parameter $U_B/U_e =1.7\times 10^{-3}$.

§.§ Ton 599 modeling

Contrary to PKS 1222+216, the SED of Ton 599 is heavily Compton dominated, with a ratio of inverse Compton to synchrotron luminosity of approximately one order of magnitude. This is a usual signature of an EIC component dominating the gamma-ray emission. We therefore consider the same scenario as for PKS 1222+216. As shown in Figure <ref>, the model describes the data well.

Broadband SED of Ton 599 during the VERITAS detection from UT 2017-12-15 to 2017-12-16. The solid blue lines show synchrotron and SSC emission; the dotted blue line shows second-order self-Compton emission; the heavy dashed green line shows thermal emission from the accretion disk; and the dashed green line shows inverse Compton emission from the BLR.

As in the case of PKS 1222+216, the thermal EIC emission imposes strong constraints on the properties of the emitting region. The largest constraint comes from the gamma-ray opacity by pair creation from the luminous thermal field surrounding the blob. We found that only for a Doppler factor of ${\gtrsim} 50$ is the EIC emission at VHE strong enough to produce the observed VHE gamma rays, given the BLR opacity.

Parameters of the SED models.

Parameter     PKS 1222+216     Ton 599     Unit
$\theta_\mathrm{obs}$     $1.0$     $1.0$      deg

$\delta$     $40$     $53$

$N_{e}^{(1)}$     $2.0\times 10^{4}$     $2.7\times 10^{5}$     cm$^{-3}$

$n_1$     $2.1$     $2.5$     $-$

$n_2$     $3.9$     $3.0$     $-$

$\gamma_{\mathrm{min}}$     $5.5\times 10^{2}$     $3.0\times 10^{2}$     $-$

$\gamma_{\mathrm{max}}$     $3.0\times 10^{5}$     $7.0\times 10^{4}$     $-$

$\gamma_{\mathrm{brk}}$     $5.0\times 10^{3}$     $1.5\times 10^{4}$     $-$

$B$     $3.0\times 10^{-2}$     $3.0\times 10^{-2}$     G

$R$     $5.5\times 10^{16}$     $6.0\times 10^{16}$     cm

$D_\mathrm{BH}$*     $3.56$     $2.33$     pc

$L_\mathrm{disk}$     $2.8\times 10^{46}$     $2.2\times 10^{46}$     erg s$^{-1}$

$T_\mathrm{disk}$     $2.8\times 10^{4}$      $1.1\times 10^{4}$     K

$\epsilon_\mathrm{BLR}$     $2.0\times 10^{-2}$     $2.0\times 10^{-2}$     $-$

$\theta_\mathrm{obs}$ is the angle of the blob direction of motion with respect to the line of sight. The electron energy distribution between Lorentz factors $\gamma_\mathrm{min}$ and $\gamma_\mathrm{max}$ is given by a broken power law with indices $n_1$ and $n_2$ below and above $\gamma_\mathrm{brk}$, with $N_e^{(1)}$ the normalization factor at $\gamma = 1$. The blob Doppler factor, magnetic field, radius, and distance to the black hole are given by $\delta$, $B$, $R$, and $D_\mathrm{BH}$, respectively. The disk luminosity and temperature are given by $L_\mathrm{disk}$ and $T_\mathrm{disk}$, while $\epsilon_\mathrm{BLR}$ is the covering factor of the broad line region.
* Host galaxy frame.

The solution presented in Figure <ref>, with $\delta = 53$, is consistent with a maximum VHE emission undergoing strong BLR absorption ($E_\mathrm{max} = 630$ GeV), with an opacity of $\tau_{\gamma \gamma,E_\mathrm{max}} = 2.8$. In this scenario we set the blob at a distance of 2.33 pc from the SMBH, corresponding to $15.7~r_\mathrm{BLR}$.

By applying the same consistency check for variability as PKS 1222+216, we found a minimal possible variability timescale predicted by the model of 18 h (coincidentally the same as PKS 1222+216), and a radiative cooling time
$T_\mathrm{cool}(\mathrm{\textit{Fermi}}) \simeq 8.7$ days, in good agreement with the observed Fermi-LAT flare decay of $11.8\pm1.1$ days.
The VERITAS observed variability can be associated with the cooling time $T_\mathrm{cool}(\gamma_{\mathrm{max}})$, which leads to $T_\mathrm{cool}(\mathrm{VERITAS}) = 45$ h, fully compatible with the observed variability of $\sim$ 2 days.

The blob is estimated to have a total power of approximately $1.2\times 10^{46}$ erg s$^{-1}$, and to be extremely particle-dominated with equipartition parameter $U_B/U_e =3.8\times 10^{-4}$.

Parameters used to calculate constraints on the parameter space.

Source     z      $D_L$     $t_\mathrm{var}$      $L_\mathrm{syn}$     $L_\mathrm{gamma}$     $L_d$     $M_\mathrm{BH}$a     $E_\mathrm{cool}$     $\epsilon_\mathrm{BLR}$b     $\epsilon_\mathrm{IR}$b

Gpc     day      erg s$^{-1}$     erg s$^{-1}$     erg s$^{-1}$     $M_\odot$     GeV

PKS 1222+216     0.434     2.44     10.0     $3.5 \times 10^{46}$     $7.8 \times 10^{46}$     $2.8 \times 10^{46}$     $3.47 \times 10^8$     7.07     0.02     0.2

Ton 599     0.725     4.54     2.0     $4.4 \times 10^{46}$     $1.2 \times 10^{48}$     $2.2 \times 10^{46}$     $6.8 \times 10^8$     326     0.02     0.2
$z$ and $D_L$ are the redshift and luminosity distance of the source. $t_\mathrm{var}$ is the variability timescale of cooling derived from each flare's fitted exponential decay. $L_\mathrm{syn}$, $L_\mathrm{gamma}$, and $L_d$ are the synchrotron luminosity, gamma-ray luminosity, and disk luminosity from the SED model. $M_\mathrm{BH}$ is the black hole mass. $E_\mathrm{cool}$ is the maximum photon energy due to the external Compton cooling of relativistic electrons. $\epsilon_\mathrm{BLR}$ and $\epsilon_\mathrm{IR}$ are the covering factors of the broad line region and IR-emitting torus region, respectively.
aFarina et al., 2012, Liu et al., 2006
bTavecchio et al., 2011

§ LORENTZ FACTORS AND LOCATIONS OF THE GAMMA-RAY-EMITTING REGIONS

(-155,190)(0,0)[lb]PKS 1222+216

(-140,190)(0,0)[lb]Ton 599
Constraints on the Lorentz factor $\Gamma$ and distance $r$ between the gamma-ray emission location and central black hole, adapted from [Nalewajko et al., 2014]. The allowed region is filled in purple. The black vertical line shows the opacity constraint on $r$ derived from the BLR modeling. The values of $\Gamma$ and $r$ derived from the SED modeling are indicated with dashed black lines.

We determined constraints on the Lorentz factor $\Gamma$ and distance $r$ of the gamma-ray emission region from the central black hole for PKS 1222+216 and Ton 599 following the method and assumptions of Nalewajko et al., 2014. The constraints are plotted in Figure <ref>. The parameters used to determine the constraints are given in Table <ref>. In order to obtain a conservative SSC constraint, we set the SSC luminosity equal to the observed gamma-ray luminosity $L_\mathrm{gamma}$. For PKS 1222+216 the fastest variability is observed with Fermi-LAT, while for Ton 599 it is observed with VERITAS. We therefore set the maximum energy $E_\mathrm{cool}$ for the EIC cooling constraint equal to the geometric mean of the energy ranges observed by Fermi-LAT and VERITAS for the two sources, respectively.

Three constraints on $\Gamma$ and $r$ are calculated. The collimation constraint requires that the size of the emission region be less than the size of the jet at the emission region location such that $\Gamma \theta \leq 1$, where $\theta$ as defined by <cit.> is the angle subtended by the blob expanding while propagating.

A caveat of the collimation constraint is the underlying assumption that the blob size is defined by the observed variability such that $R = c \delta t_\mathrm{var,obs} /(1+z)$. However, the observed variability gives only an upper limit on the blob radius, meaning that the actual size of the emission zone is likely smaller than that extrapolated from the observed variability. Indeed, our modeling of PKS 1222+216 and Ton 599 predicts observed variability much shorter than the one observed within the reconstructed SED periods. This discrepancy explains why the parameters predicted by our model lie only just below the line $\Gamma \theta = 1$.

The SSC constraint requires that the SSC luminosity should not exceed the total gamma-ray luminosity, which includes contributions from external radiation fields <cit.>. The cooling constraint requires that electrons radiatively emitting gamma rays at energies above $E_\mathrm{cool}$ cool through interactions with external radiation fields faster than the flare decay timescale <cit.>.

These parameter limits do not take into account the constraints given by the BLR and dust-torus opacity on the gamma-ray emission. We show with black vertical lines the minimum distance $r$ from the black hole in the SED models where the BLR would become fully opaque for the maximum observed energy $E_\mathrm{max}$ (370 GeV for PKS 1222+216 and 630 GeV for Ton 599). We consider the BLR opaque when $\tau_{\mathrm{BLR}, E_\mathrm{max}} > 5$, meaning that less than 1% of the gamma rays can escape.
We can clearly see that considering the BLR opacity significantly tightens the constraints on the gamma-ray emission location in Ton 599, as mentioned in the previous section. The opacity constraint on PKS 1222+216 is weaker, as in that case the blob does not have to be as deep inside the BLR to reproduce the observed EIC emission.

§ NEUTRINO EMISSION DURING VHE FLARES

Luminous gamma-ray flares of FSRQs are potential sources of PeV-scale (${\sim}100~\mathrm{TeV}$ – ${\sim}10~\mathrm{PeV}$) neutrino emission <cit.>. While the lack of point sources observed in IceCube data suggests that FSRQs are not the dominant population of neutrino sources, the possibility of neutrino emission during rare, bright flares has not been excluded [Murase & Waxman, 2016]. While Righi et al., 2020 have suggested that the bulk of the average neutrino emission from FSRQs occurs in the sub-EeV – EeV energy range, their results do not exclude PeV-scale neutrino emission during outlier states. In the SED modeling of the VHE flares of PKS 1222+216 and Ton 599 presented here, a purely leptonic model gives an adequate representation of the data, and performing full hadronic modeling is beyond the scope of this work.

However, we can place analytic constraints on the potential PeV-scale neutrino flux produced during these events by considering a lepto-hadronic scenario in which synchrotron emission from secondary electrons produced by pion decay contributes a subdominant component to the second peak of the SED, similar to models used to describe the flaring emission of TXS 0506+056 coincident with the detection of a neutrino by IceCube <cit.>. In this section, we make use of the assumptions and methods of Gao et al., 2017, particularly Appendix A of that work. All quantities in the following equations are in the comoving frame of the blob, unless explicitly noted with the superscript “ob”.

We consider neutrinos produced by the $p\gamma$ interaction via the $\Delta^+$(1232) resonance with threshold energy $\epsilon_{p\gamma,\mathrm{th}} \sim 0.3$ GeV. The characteristic proton energy is $E_{p,\mathrm{char}} \sim E_\nu/K_\nu \sim 2~\mathrm{PeV}$, where $E_\nu = E_\nu^\mathrm{ob}(1+z)/\Gamma$ and $K_\nu \sim 0.05$ [Murase et al., 2014]. Therefore, to check whether these sources could in principle support PeV neutrino emission, we first estimate the maximum energy to which protons can be accelerated in the source without escaping, following Hillas, 1984, as

\begin{equation}
E_{p,\mathrm{max}} = Z e \beta c B R,
\end{equation}

where the atomic number $Z = 1$ for protons, $e$ is the elementary charge, $\beta = v / c \sim 1$ for highly relativistic particles, $c$ is the speed of light, $B$ is the magnetic field in the source, and $R$ is the size of the source. Using the values in Table <ref>, the maximum energy to which protons could have been accelerated in the gamma-ray emission regions for the flares of PKS 1222+216 and Ton 599 is $E_p \sim 500$ PeV, equivalent to an upper limit on the neutrino energy of $E_\nu^\mathrm{ob} \sim 400$ PeV, so PeV-scale neutrino emission is certainly feasible.

A limit on the neutrino flux can be imposed by considering a steady state within the time period of the flare in which the synchrotron luminosity of the secondary electrons equals the power injected by pion decay. Because the lifetime of the ultra-high-energy protons cooling by photo-pion production may be longer than that of the electrons, the resulting model can be considered to provide an upper limit on the neutrino luminosity. The steady-state proton energy density at $E_{p,\mathrm{char}}$ is given by [Gao et al., 2017]:

\begin{equation}\label{eq:proton_energy_density}
u_p(E_{p,\mathrm{char}}) = \frac{\alpha_\mathrm{fs} \epsilon_{p\gamma,\mathrm{th}}}{c \sigma_{p\gamma} K_{p\gamma}K_{\pi \rightarrow e}} \frac{m_p c^2}{E_{p,\mathrm{char}}} \frac{\nu F_{\nu,\mathrm{2}}^\mathrm{ob}}{\nu F_{\nu,\mathrm{t}}^\mathrm{ob}},
\end{equation}

where $\alpha_\mathrm{fs} = 4c/(3R)$ is the free-streaming escape rate, $\sigma_{p\gamma} = 5.0\times10^{-28}~\mathrm{cm}^2$ is the $p\gamma$ cross-section, $K_{p\gamma} \sim 0.2$ is the average inelasticity for the proton in the $p\gamma$ interaction, $K_{\pi \rightarrow e} \sim 1/8$ is the fraction of energy transferred to $e^\pm$ pairs from pion decay, $\nu F_{\nu,\mathrm{2}}^\mathrm{ob}$ is the observed flux due to synchrotron emission from the secondary electrons, and $\nu F_{\nu,\mathrm{t}}^\mathrm{ob}$ is the observed flux of the target photons of the $p\gamma$ interaction at $E_\mathrm{t} \sim \epsilon_{p\gamma,\mathrm{th}} m_p c^2 / E_{p,\mathrm{char}}$, which is directly constrained by the Swift-XRT measurement at $E_\mathrm{t}\Gamma/(1+z) \sim 2~\mathrm{keV}$. We can estimate $\nu F_{\nu,\mathrm{2}}^\mathrm{ob}$ from the SED at the synchrotron peak frequency of the secondary electrons at

\begin{equation}
\begin{split}
\nu_2^\mathrm{ob} & = \frac{c e (K_{p\gamma}K_{\pi \rightarrow e})^2}{2\pi (m_e c^2)^3}\frac{\Gamma}{1 + z} B E_{p,\mathrm{char}}\\[1ex]
& \approx 10^{22}~\left(\frac{\Gamma}{23}\right) \left(\frac{B}{30 ~\mathrm{mG}}\right) \left(\frac{E_{p,\mathrm{char}}}{2~\mathrm{PeV}}\right)~\mathrm{Hz},
\end{split}
\end{equation}

for redshift $z \sim 0.5$. The corresponding power in protons can be estimated as

\begin{equation}
L_p \sim \Gamma^2 u_p(E_{p,\mathrm{char}}) \alpha_\mathrm{esc} V(R),
\end{equation}

where for simplicity we assume the proton escape time $\alpha_\mathrm{esc} = 0.1\alpha_\mathrm{fs}$ and let $V(R) = 4\pi R^3/3$ for a spherical blob. Parameterizing this power by the Eddington luminosity boosted into the jet frame yields

\begin{equation}
\begin{split}
L_p \sim 1 &\left(\frac{\Gamma}{23}\right)^2\left(\frac{R}{6\times10^{16}~\mathrm{cm}}\right) \left(\frac{M}{5\times10^{10}~\mathrm{M}_\odot}\right)^{-1}\\[1ex]
\times &\left(\frac{E_{p,\mathrm{char}}}{2~\mathrm{PeV}}\right)^{-1} \left(\frac{\nu F_{\nu,2}^\mathrm{ob} / \nu F_{\nu,\mathrm{t}}^\mathrm{ob}}{0.5}\right)~L_\mathrm{Edd},
\end{split}
\end{equation}

where $\nu F_{\nu,2}^\mathrm{ob} \lesssim 0.5~\nu F_{\nu,\mathrm{t}}^\mathrm{ob} \sim 1\times10^{-12}~\mathrm{erg}~\mathrm{cm}^{-2}~\mathrm{s}^{-1}$ is a conservative estimate of the largest energetically reasonable contribution[This assumption requires about $5\times10^{-3}$ protons for every electron, from Eq. <ref> and the electron energy distributions reported in Table <ref>.] to the SED at ${\sim}10^{22}$ Hz. The contribution is clearly subdominant. We can then estimate the observed neutrino energy flux, where $\dot{N}_{p\gamma}$ is the $p\gamma$ event rate per physical volume, using the relation

\begin{equation}
\dot{N}_{p\gamma} \sim c\sigma_{p\gamma}\frac{u_\mathrm{ph}(E_\mathrm{t})}{E_\mathrm{t}}\frac{u_p(E_{p,\mathrm{char}})}{E_{p,\mathrm{char}}} = \frac{\alpha_\mathrm{fs}u_\nu(E_\nu)}{E_\nu},
\end{equation}

where $u_\mathrm{ph}(E_\mathrm{t})$ is the energy density of the target photons. Since $u_\nu(E_\nu)/u_\mathrm{ph}(E_\mathrm{t}) = \nu F_{\nu,\nu}^\mathrm{ob} / \nu F_{\nu,\mathrm{t}}^\mathrm{ob}$, we obtain the simple relation

\begin{equation}
\nu F_{\nu,\nu}^\mathrm{ob} = \frac{K_\nu}{K_{p\gamma} K_{\pi\rightarrow e}} \nu F_{\nu,2}^\mathrm{ob} \sim 2~\nu F_{\nu,2}^\mathrm{ob}.
\end{equation}

The number of PeV-scale neutrinos of any flavor expected to be detected by IceCube during the VHE flare of PKS 1222+216 or Ton 599 is then

\begin{equation}
\begin{split}
N_\nu \lesssim 0.001 &\left( \frac{\nu F_{\nu,\nu}^\mathrm{ob}}{2\times10^{-12}~\mathrm{erg}~\mathrm{cm}^{-2}~\mathrm{s}^{-1}}\right)\\[1ex]
\times &\left(\frac{\Delta T}{5~\mathrm{day}}\right) \left(\frac{\mathcal{A}_\mathrm{eff}}{10^6~\mathrm{cm}^2}\right),
\end{split}
\end{equation}

where $\Delta T$ is the duration of the VHE flare, $\mathcal{A}_\mathrm{eff} \sim 10^6~\mathrm{cm}^2$ is the IceCube effective area for extremely high-energy real-time alerts in the PeV range [Aartsen et al., 2017], and $\Delta \nu \sim \ln(10)$ is assumed for the width of the neutrino spectrum. We conclude that it is plausible that PKS 1222+216 and Ton 599 could have produced PeV-scale neutrinos during their VHE flaring activity at a flux consistent with a null detection by current instruments.

To reduce the model-dependence of our constraints, $R$ could also be estimated using the timescale of gamma-ray flare variability,

\begin{equation}
R \sim \frac{\delta}{1 + z} c \Delta T,
\end{equation}

from which estimates of $R \sim 7\times10^{17}~\mathrm{cm}~\mathrm{and}~1\times10^{17}~\mathrm{cm}$ are obtained for PKS 1222+216 and Ton 599. For the two sources, the constraints on the maximum neutrino energy are loosened to $E_\nu^\mathrm{ob} \sim 5~\mathrm{EeV}~\mathrm{and}~700~\mathrm{PeV}$, respectively, and the required proton luminosities are increased by a factor of ${\sim}10$ and ${\sim}2$, or within a few times the Eddington luminosity for both sources.

§ CONCLUSIONS

In this paper, we present an analysis of the gamma-ray and multiwavelength emission of three bright, variable FSRQs observed by Fermi-LAT and VERITAS: 3C 279, PKS 1222+216, and Ton 599, making use of almost 100 hours of VERITAS observations of these sources. No VHE gamma-ray activity was observed during multiple flares of 3C 279, which is the brightest of the three sources as observed with Fermi-LAT, but VERITAS detected flares of both PKS 1222+216 and Ton 599.

The flux distributions of the Fermi-LAT light curves of all three sources are consistent with the PDF derived from the stochastic differential equation proposed by Tavecchio et al., 2020, in which the timescale of variability is controlled by processes in the accretion disk. The timescales associated with magneto-rotational instabilities and magnetic flux accumulation in the accretion disk obtained from this model are consistent with theoretical estimates.

We selected gamma-ray flaring states from the Fermi-LAT light curves of the three sources using a procedure based on Bayesian blocks. Daily to sub-daily variability was observed by Fermi-LAT during the flaring states. No pattern of asymmetry was found in the rise and decay times of the exponential components of the ten identified flares of 3C 279, while the fluence distribution of the flare components extended over one order of magnitude.

All three sources have similar gamma-ray spectra consistent with a log-parabola model. The average flaring spectrum of 3C 279 was found to exhibit stronger curvature than the baseline state, consistent with a null VHE detection even during extremely bright flares.

The SEDs of VHE flares of PKS 1222+216 and Ton 599 are described well by a purely leptonic emission model including an external inverse Compton emission component. For both sources, a strong constraint was placed on the Doppler factor, which must be ${\gtrsim}40$ for PKS 1222+216 and ${\gtrsim}50$ for Ton 599, to produce the observed gamma-ray emission at up to TeV energies despite internal absorption. We constrained the jet Lorentz factor and distance of the gamma-ray emission region from the central black hole using the independent method of Nalewajko et al., 2014, which we augmented using an opacity constraint derived from our SED modeling.
We found that both sources are operating in a similar regime, with both of them having strongly matter-dominated energetics and a gamma-ray emission zone located a few parsecs from the SMBH. It would be interesting in future work to extend this study to a large sample of TeV FSRQs in order to define the common physical properties of this blazar sub-class.

We calculated analytic constraints on a supposed subdominant hadronic component to the leptonic SED model, from which we estimated upper limits on potential PeV-scale neutrino emission during the TeV flares. We found that neutrino emission is energetically plausible at a flux consistent with a null detection by IceCube during the TeV flares.

This research is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation and the Smithsonian Institution, by NSERC in Canada, and by the Helmholtz Association in Germany. This research used resources provided by the Open Science Grid, which is supported by the National Science Foundation and the U.S. Department of Energy's Office of Science, and resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the instrument.

The Fermi LAT Collaboration acknowledges generous ongoing support
from a number of agencies and institutes that have supported both the
development and the operation of the LAT as well as scientific data analysis.
These include the National Aeronautics and Space Administration and the
Department of Energy in the United States, the Commissariat à l'Energie Atomique
and the Centre National de la Recherche Scientifique / Institut National de Physique
Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana
and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education,
Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research
Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and
the K. A. Wallenberg Foundation, the Swedish Research Council and the
Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully
acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre
National d'Études Spatiales in France. This work performed in part under DOE
Contract DE-AC02-76SF00515.

Data from the Steward Observatory spectropolarimetric monitoring project were used. This program is supported by Fermi Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX15AU81G. L.S. acknowledges support from the Sloan Fellowship, the Cottrell Scholars Award, DoE DE-SC0016542, NASA 80NSSC18K1104 and NSF PHY-1903412.
This material is based upon work supported in part by the National Science Foundation under Grant PHY-1659528.

A.B.'s work was supported by the NSF Grant PHY-1806554 at Barnard College, Columbia University. O.H. thanks NSF for support under grants PHY-1707432 and PHY-2011420.
J.V. was partially supported by the Alliance program, a partnership between Columbia University, NY, USA, and three major French institutions: École Polytechnique, Paris 1 Panthéon-Sorbonne University and Science Po. This material is based upon work supported by NASA under award number 80GSFC21M0002. S.B. acknowledges financial support by the European Research Council for the ERC Starting grant MessMapp, under contract no. 949555.

This work made use of the online cosmological calculator of Wright, 2006 and of NASA's Astrophysics Data System. We would like to thank the anonymous referee for the helpful and timely comments that helped improve the quality of the paper.

Astropy [Astropy Collaboration et al., 2013, Astropy Collaboration et al., 2018],
Fermitools [Fermi Science Support Development Team, 2019],
NumPy [Harris et al., 2020],
Matplotlib [Hunter, 2007],
SciPy [Virtanen et al., 2020].

[Aartsen et al., 2017]
Aartsen, M., Ackermann, M., Adams, J., et al. 2017, Astroparticle Physics,
92, 30

[Abdo et al., 2010]
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, , 722, 520

[Abdo et al., 2011]
Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, , 734, 116

[Abdo et al., 2012]
—. 2012, , 758, 140

[Abdollahi et al., 2020]
Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, , 247, 33

[Abeysekara et al., 2015]
Abeysekara, A. U., Archambault, S., Archer, A., et al. 2015, ,
815, L22

[Abolmasov & Poutanen, 2017]
Abolmasov, P., & Poutanen, J. 2017, , 464, 152

[Ackermann et al., 2016]
Ackermann, M., Anantua, R., Asano, K., et al. 2016, , 824, L20

[Aharonian et al., 2007]
Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2007, ,
664, L71

[Aharonian, 2000]
Aharonian, F. A. 2000, , 5, 377

[Aitchison & Brown, 1973]
Aitchison, J., & Brown, A. C. 1973, Cambridge University Press, 60, 104

[Ajello et al., 2020]
Ajello, M., Angioni, R., Axelsson, M., et al. 2020, , 892, 105

[Aleksić et al., 2011]
Aleksić, J., Antonelli, L. A., Antoranz, P., et al. 2011, ,
730, L8

[Astropy Collaboration et al., 2013]
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,
, 558, A33

[Astropy Collaboration et al., 2018]
Astropy Collaboration, Price-Whelan, A. M., SipHocz, B. M., et al.
2018, , 156, 123

[Atwood et al., 2013]
Atwood, W., Albert, A., Baldini, L., et al. 2013, arXiv e-prints,

[Atwood et al., 2009]
Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, The
Astrophysical Journal, 697, 1071

[Ballet et al., 2020]
Ballet, J., Burnett, T. H., Digel, S. W., & Lott, B. 2020,

[Barbiellini et al., 2014]
Barbiellini, G., Bastieri, D., Bechtol, K., et al. 2014, , 784, 118

[Bennett et al., 2014]
Bennett, C. L., Larson, D., Weiland, J. L., & Hinshaw, G. 2014, ,
794, 135

[Blandford & Königl, 1979]
Blandford, R. D., & Königl, A. 1979, , 232, 34

[Bruel et al., 2018]
Bruel, P., Burnett, T. H., Digel, S. W., et al. 2018, arXiv e-prints,

[Burbidge, 1968]
Burbidge, E. M. 1968, , 154, L109

[Burrows et al., 2005]
Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, , 120, 165

[Cerruti et al., 2019]
Cerruti, M., Zech, A., Boisson, C., et al. 2019, , 483, L12

[Cheung et al., 2017]
Cheung, C. C., Gasparrini, D., & Buson, S. 2017, The Astronomer's
Telegram, 10931, 1

[Cogan, 2007]
Cogan, P. 2007, in Proc. 30th ICRC (Mérida), Vol. 3, 1385–1388

[Costamante et al., 2018]
Costamante, L., Cutini, S., Tosti, G., Antolini, E., & Tramacere, A.
2018, , 477, 4749

[Dermer & Menon, 2009]
Dermer, C. D., & Menon, G. 2009, High Energy Radiation from Black Holes:
Gamma Rays, Cosmic Rays, and Neutrinos (Princeton University Press)

[Dermer et al., 2014]
Dermer, C. D., Murase, K., & Inoue, Y. 2014, Journal of High Energy
Astrophysics, 3, 29

[Emery et al., 2019]
Emery, G., Cerruti, M., Dmytriiev, A., et al. 2019, in Proc. 36th ICRC
(Madison), Vol. 36, 668

[Errando, 2014]
Errando, M. 2014, in AAS/High Energy Astrophysics Division, Vol. 14, AAS/High
Energy Astrophysics Division #14, 106.11

[Fanaroff & Riley, 1974]
Fanaroff, B. L., & Riley, J. M. 1974, , 167, 31P

[Farina et al., 2012]
Farina, E. P., Decarli, R., Falomo, R., Treves, A., & Raiteri, C. M. 2012,
, 424, 393

[Fermi Science Support Development Team, 2019]
Fermi Science Support Development Team. 2019, Fermitools: Fermi Science
Tools,  Astrophysics Source Code Library, ascl:1905.011

[Fossati et al., 1998]
Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini, G. 1998,
, 299, 433

[Franceschini & Rodighiero, 2017]
Franceschini, A., & Rodighiero, G. 2017, , 603, A34

[Gao et al., 2019]
Gao, S., Fedynitch, A., Winter, W., & Pohl, M. 2019, Nature Astronomy, 3, 88

[Gao et al., 2017]
Gao, S., Pohl, M., & Winter, W. 2017, , 843, 109

[Gehrels et al., 2004]
Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, , 611, 1005

[Ghisellini et al., 1998]
Ghisellini, G., Celotti, A., Fossati, G., Maraschi, L., & Comastri,
A. 1998, , 301, 451

[Ghisellini & Madau, 1996]
Ghisellini, G., & Madau, P. 1996, , 280, 67

[Ghisellini & Tavecchio, 2009]
Ghisellini, G., & Tavecchio, F. 2009, , 397, 985

[Giannios & Uzdensky, 2019]
Giannios, D., & Uzdensky, D. A. 2019, , 484, 1378

[Giebels & Degrange, 2009]
Giebels, B., & Degrange, B. 2009, , 503, 797

[Harris et al., 2020]
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585,

[Hayashida et al., 2015]
Hayashida, M., Nalewajko, K., Madejski, G. M., et al. 2015, , 807,

[Hervet et al., 2015]
Hervet, O., Boisson, C., & Sol, H. 2015, , 578, A69

[Hervet et al., 2016]
—. 2016, , 592, A22

[Hillas, 1984]
Hillas, A. M. 1984, , 22, 425

[Holder, 2011]
Holder, J. 2011, in Proc. 32nd ICRC (Beijing), Vol. 11, 137

[Holder, 2014]
Holder, J. 2014, The Astronomer's Telegram, 5981

[Hunter, 2007]
Hunter, J. D. 2007, Computing in Science Engineering, 9, 90

[IceCube Collaboration et al., 2018]
IceCube Collaboration, Aartsen, M. G., Ackermann, M., et al. 2018,
Science, 361, eaat1378

[Inoue & Takahara, 1996]
Inoue, S., & Takahara, F. 1996, , 463, 555

[Johannesson & Orlando, 2013]
Johannesson, G., & Orlando, E. 2013, in International Cosmic Ray
Conference, Vol. 33, International Cosmic Ray Conference, 3106

[Jorstad et al., 2001]
Jorstad, S. G., Marscher, A. P., Mattox, J. R., et al. 2001, , 134, 181

[Kadler et al., 2016]
Kadler, M., Krauß, F., Mannheim, K., et al. 2016, Nature Physics,
12, 807

[Kalberla et al., 2005]
Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, ,
440, 775

[Katarzyński et al., 2001]
Katarzyński, K., Sol, H., & Kus, A. 2001, , 367, 809

[Keenan et al., 2021]
Keenan, M., Meyer, E. T., Georganopoulos, M., Reddy, K., & French,
O. J. 2021, , 505, 4726

[Keivani et al., 2018]
Keivani, A., Murase, K., Petropoulou, M., et al. 2018, , 864, 84

[Krause et al., 2017]
Krause, M., Pueschel, E., & Maier, G. 2017, Astroparticle Physics, 89, 1

[Liu et al., 2006]
Liu, Y., Jiang, D. R., & Gu, M. F. 2006, , 637, 669

[Lynds et al., 1965]
Lynds, C. R., Stockton, A. N., & Livingston, W. C. 1965, , 142, 1667

[Maier & Holder, 2017]
Maier, G., & Holder, J. 2017, in Proc. 35th ICRC (Busan), Vol. 301, 747

[Mannheim, 1993]
Mannheim, K. 1993, , 269, 67

[Maraschi et al., 1992]
Maraschi, L., Ghisellini, G., & Celotti, A. 1992, , 397, L5

[Meyer et al., 2019]
Meyer, M., Scargle, J. D., & Blandford, R. D. 2019, , 877, 39

[Mirzoyan, 2017]
Mirzoyan, R. 2017, The Astronomer's Telegram, 11061

[Mücke & Protheroe, 2000]
Mücke, A., & Protheroe, R. J. 2000, in AIP Conference Series, ed.
B. L. Dingus, M. H. Salamon, & D. B. Kieda, Vol. 515, 149–153

[Mukherjee, 2001]
Mukherjee, R. 2001, AIP Conference Proceedings, 558, 324

[Mukherjee, 2017]
Mukherjee, R. 2017, The Astronomer's Telegram, 11075

[Murase et al., 2014]
Murase, K., Inoue, Y., & Dermer, C. D. 2014, , 90, 1

[Murase & Waxman, 2016]
Murase, K., & Waxman, E. 2016, , 94, 103006

[Nalewajko et al., 2014]
Nalewajko, K., Begelman, M. C., & Sikora, M. 2014, , 789, 161

[Nieppola et al., 2008]
Nieppola, E., Valtaoja, E., Tornikoski, M., Hovatta, T., &
Kotiranta, M. 2008, , 488, 867

[Osterbrock & Pogge, 1987]
Osterbrock, D. E., & Pogge, R. W. 1987, , 323, 108

[Padovani et al., 2019]
Padovani, P., Oikonomou, F., Petropoulou, M., Giommi, P., & Resconi, E. 2019,
, 484, L104

[Patel & Chitnis, 2020]
Patel, S. R., & Chitnis, V. R. 2020, , 492, 72

[Petropoulou et al., 2016]
Petropoulou, M., Giannios, D., & Sironi, L. 2016, , 462, 3325

[Poole et al., 2008]
Poole, T. S., Breeveld, A. A., Page, M. J., et al. 2008, , 383,

[Prince, 2019]
Prince, R. 2019, , 871, 101

[Prince, 2020]
—. 2020, , 890, 164

[Reimer et al., 2019]
Reimer, A., Böttcher, M., & Buson, S. 2019, , 881, 46

[Resconi et al., 2009]
Resconi, E., Franco, D., Gross, A., Costamante, L., & Flaccomio, E. 2009,
, 502, 499

[Righi et al., 2020]
Righi, C., Palladino, A., Tavecchio, F., & Vissani, F. 2020, ,
642, A92

[Roming et al., 2005]
Roming, P. W. A., Kennedy, T. E., Mason, K. O., et al. 2005, , 120,

[Roming et al., 2009]
Roming, P. W. A., Koch, T. S., Oates, S. R., et al. 2009, , 690,

[Romoli et al., 2017]
Romoli, C., Zacharias, M., Meyer, M., et al. 2017, in Proc. 35th ICRC
(Busan), Vol. 301, 649

[Roy et al., 2019]
Roy, N., Chatterjee, R., Joshi, M., & Ghosh, A. 2019, , 482, 743

[Scargle et al., 2013]
Scargle, J. D., Norris, J. P., Jackson, B., & Chiang, J. 2013, ,
764, 167

[Schlafly & Finkbeiner, 2011]
Schlafly, E. F., & Finkbeiner, D. P. 2011, , 737, 103

[Schneider et al., 2010]
Schneider, D. P., Richards, G. T., Hall, P. B., et al. 2010, , 139,

[Shah et al., 2018]
Shah, Z., Mankuzhiyil, N., Sinha, A., et al. 2018, Research in
Astronomy and Astrophysics, 18, 141

[Sikora et al., 2009]
Sikora, M., Stawarz, Ł., Moderski, R., Nalewajko, K., & Madejski,
G. M. 2009, , 704, 38

[Sinha et al., 2017]
Sinha, A., Sahayanathan, S., Acharya, B. S., et al. 2017, , 836, 83

[Smith et al., 2009]
Smith, P. S., Montiel, E., Rightley, S., et al. 2009, 2009 Fermi
Symposium, eConf Proceedings, C091122, arXiv:0912.3621

[Tavecchio et al., 2011]
Tavecchio, F., Becerra-Gonzalez, J., Ghisellini, G., et al. 2011, , 534,

[Tavecchio et al., 2020]
Tavecchio, F., Bonnoli, G., & Galanti, G. 2020, , 497, 1294

[Valtaoja et al., 1999]
Valtaoja, E., Lähteenmäki, A., Teräsranta, H., & Lainela,
M. 1999, , 120, 95

[Virtanen et al., 2020]
Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17,

[Wehrle et al., 1998]
Wehrle, A. E., Pian, E., Urry, C. M., et al. 1998, , 497, 178

[Wilcoxon, 1945]
Wilcoxon, F. 1945, Biometrics Bulletin, 1, 80

[Willingale et al., 2013]
Willingale, R., Starling, R. L. C., Beardmore, A. P., Tanvir, N. R., &
O'Brien, P. T. 2013, , 431, 394

[Wright, 2006]
Wright, E. L. 2006, , 118, 1711

[Yoo & An, 2020]
Yoo, S., & An, H. 2020, , 902, 2

§ COMPLETE SET OF FERMI-LAT FLARE PROFILES FOR 3C 279

3C 279 LAT sub-daily light curves (blue points) around the flaring episodes selected as described in Section <ref> (shaded areas). The dotted blue lines show the fitted exponential profiles, with their sums shown in solid blue. The VERITAS 95% upper limits are shown as black downwards arrows. For Flares 1, 5, 7, and 9, no VERITAS observations were taken around the time of the flare.

Results of the LAT flare profile fits for 3C 279.

Amplitude ($F_0$)     $t_{\rm{peak}}$      $t_{\rm{rise}}$     $t_{\rm{decay}}$      Constant ($F_{\rm{const.}}$)

($10^{-9} $ erg cm$^{-2}\,$s$^{-1}$)        (MJD)       (min.)      (min.)     ($10^{-9} $ erg cm$^{-2}\,$s$^{-1}$)

5cFlare 1 (MJD 56645.655 – 56647.655): $\chi^2/$d.o.f.= 12.05/8 = 1.51
9.56 $\pm$ 1.07     56646.330 $\pm$ 0.033     130 $\pm$ 45     674 $\pm$ 73     0.28 $\pm$ 0.06
5cFlare 2 (MJD 56717.655 – 56718.655): $\chi^2/$d.o.f.= 16.63/10 = 1.66
4.40 $\pm$ 0.70     56718.142 $\pm$ 0.043     445 $\pm$ 95     307 $\pm$ 86     0.55 $\pm$ 0.07
5cFlare 3 (MJD 56749.655 – 56754.655): $\chi^2/$d.o.f.= 69.98/34 = 2.06
7.27 $\pm$ 0.64     56750.382 $\pm$ 0.015     229 $\pm$ 25     267 $\pm$ 42     N.A.

2.78 $\pm$ 1.15     56751.238 $\pm$ 0.024     140 $\pm$ 82     69 $\pm$ 47     N.A.

4.80 $\pm$ 0.37     56752.532 $\pm$ 0.067     2001 $\pm$ 116     631 $\pm$ 136     N.A.
5cFlare 4 (MJD 57186.655 – 57190.655): $\chi^2/$d.o.f.= 77.31/19 = 4.07
12.07 $\pm$ 0.67     57187.446 $\pm$ 0.031     378 $\pm$ 46     1784 $\pm$ 147     N.A.

9.79 $\pm$ 2.29     57188.425 $\pm$ 0.028     216 $\pm$ 101     155 $\pm$ 64     N.A.

21.72 $\pm$ 1.59     57189.069 $\pm$ 0.008     137 $\pm$ 18     512 $\pm$ 55     N.A.

12.41 $\pm$ 1.30     57189.532 $\pm$ 0.010     220 $\pm$ 63     77 $\pm$ 25     N.A.
5cFlare 5 (MJD 58116.655 – 58119.655): $\chi^2/$d.o.f.= 54.70/29 = 1.89
3.72 $\pm$ 0.20     58118.171 $\pm$ 0.069     1278 $\pm$ 220     2521 $\pm$ 309     0.06 $\pm$ 0.11
5cFlare 6 (MJD 58130.655 – 58141.655): $\chi^2/$d.o.f.= 141.28/72 = 1.96
7.08 $\pm$ 1.01     58134.520 $\pm$ 0.055     3719 $\pm$ 390     421 $\pm$ 259     N.A.

10.95 $\pm$ 3.74     58135.229 $\pm$ 0.053     718 $\pm$ 232     3535 $\pm$ 1394     N.A.

11.78 $\pm$ 7.00     58136.266 $\pm$ 0.048     349 $\pm$ 160     1839 $\pm$ 1055     N.A.

3.44 $\pm$ 0.66     58139.546 $\pm$ 0.033     233 $\pm$ 175     6119 $\pm$ 1824     N.A.
5cFlare 7 (MJD 58168.655 – 58173.655): $\chi^2/$d.o.f.= 78.81/58 = 1.36
2.36 $\pm$ 0.50     58172.345 $\pm$ 0.242     8540 $\pm$ 4159     4458 $\pm$ 2319     0.45 $\pm$ 0.59
5cFlare 8 (MJD 58222.655 – 58230.655): $\chi^2/$d.o.f.= 177.25/106 = 1.67
5.29 $\pm$ 1.29     58224.773 $\pm$ 0.105     1996 $\pm$ 716     5899 $\pm$ 4035     N.A.

17.70 $\pm$ 2.01     58227.945 $\pm$ 0.004     36 $\mathbf{\pm}$ 13     329 $\pm$ 131     N.A.

16.42 $\pm$ 1.87     58228.323 $\pm$ 0.012     140 $\pm$ 54     115 $\pm$ 48     N.A.

5.59 $\pm$ 1.69     58227.139 $\pm$ 0.133     3816 $\pm$ 1450     4077 $\pm$ 2080     N.A.
5cFlare 9 (MJD 58239.655 – 58247.655): $\chi^2/$d.o.f.= 46.25/34 = 1.36
2.96 $\pm$ 0.40     58241.258 $\pm$ 0.149     2546 $\pm$ 595     2226 $\pm$ 1088     N.A.

3.27 $\pm$ 0.25     58245.648 $\pm$ 0.133     3080 $\pm$ 1384     3028 $\pm$ 303     N.A.
5cFlare 10 (MJD 58268.655 – 58275.655): $\chi^2/$d.o.f.= 75.80/55 = 1.37
6.23 $\pm$ 9.46     58269.171 $\pm$ 0.182     73 $\pm$ 236     177 $\pm$ 102     N.A.

8.81 $\pm$ 0.84     58270.137 $\pm$ 0.107     2392 $\pm$ 243     2449 $\pm$ 956     N.A.

4.46 $\pm$ 1.91     58271.223 $\pm$ 0.088     477 $\pm$ 431     5824 $\pm$ 862     N.A.

The smallest variability time found is indicated in boldface.

§ FERMI-LAT SPECTRAL FIT PARAMETERS

Fermi-LAT spectral fit parameters.

State     Epoch     TS     $N_0$     $\alpha$     $\beta$     Flux

(MJD)          ($10^{-10}\,$MeV$^{-1}\,$cm$^{-2}\,$s$^{-1}$)          ($\times\,10^{-2}$)     ($10^{-6}\,$ph$\,$cm$^{-2}\,$s$^{-1}$)
7c3C 279
Global     54682.66 – 58459.35     271945     3.33 $\pm$ 0.02     2.228 $\pm$ 0.004     6.1 $\pm$ 0.3     0.751 $\pm$ 0.004

Low state     56230.66 – 56465.66     1130     0.54 $\pm$ 0.03     2.38 $\pm$ 0.06     2.9 $\pm$ 3.1     0.14 $\pm$ 0.01

VER-LAT quiescent     Various     322     5.7 $\pm$ 0.8     2.2 $\pm$ 0.1     1.8 $\pm$ 7.0     1.4 $\pm$ 0.2
41.7cmLAT Flares simultaneous with VER obs.     56750.27 – 56750.34     578     69 $\pm$ 10     2.1 $\pm$ 0.1     31 $\pm$ 13     11 $\pm$ 2

57189.17 – 57189.23     1141     173 $\pm$ 15    2.07 $\pm$ 0.07     14 $\pm$ 6     32 $\pm$ 3

58227.22 – 58227.27     355     79 $\pm$ 14     1.8 $\pm$ 0.2     21 $\pm$ 13     12 $\pm$ 4

58272.18 – 58272.22     235     62 $\pm$ 15     1.5 $\pm$ 0.3     42 $\pm$ 2     7 $\pm$ 3
101.7cmLAT Flares     56645.66 – 56647.66     1633     22.0 $\pm$ 1.5     1.73 $\pm$ 0.07     9.6 $\pm$ 3.2     3.7 $\pm$ 0.3

56717.66 – 56718.66     900     23.2 $\pm$ 2.0     2.08 $\pm$ 0.08     11 $\pm$ 6     4.5 $\pm$ 0.5

56749.66 – 56754.66     7680     24.0 $\pm$ 0.8     2.20 $\pm$ 0.03     13 $\pm$ 2     5.0 $\pm$ 0.2

57186.66 – 57190.66     23623     77.9 $\pm$ 1.5     2.04 $\pm$ 0.02     11 $\pm$ 1     14.7 $\pm$ 0.5

58116.66 – 58119.66     3543     19.5 $\pm$ 0.9     2.06 $\pm$ 0.04     4.1 $\pm$ 2.2     4.0 $\pm$ 0.2

58130.66 – 58141.66     27256     53.3 $\pm$ 0.9     2.14 $\pm$ 0.02     8.5 $\pm$ 1.1     11.0 $\pm$ 0.2

58168.66 – 58173.66     4932     19.5 $\pm$ 0.7     2.10 $\pm$ 0.03     5.8 $\pm$ 2.0     4.1 $\pm$ 0.2

58222.66 – 58230.66     53745     59.4 $\pm$ 0.7     2.00 $\pm$ 0.01     9.6 $\pm$ 0.7     11.2 $\pm$ 0.1

58239.66 – 58247.66     7989     19.2 $\pm$ 0.6     1.90 $\pm$ 0.03     6.2 $\pm$ 1.4     3.6 $\pm$ 0.1

58268.66 – 58275.66     107456     47.8 $\pm$ 1.3     1.91 $\pm$ 0.02     14 $\pm$ 2     8.2 $\pm$ 0.2
7cPKS 1222+216
Global     54682.66 – 58459.64     94556     1.66 $\pm$ 0.01     2.305 $\pm$ 0.007     3.8 $\pm$ 0.4     0.337 $\pm$ 0.002
7cTon 599
Global     54682.66 – 58464.49     48176     6.55 $\pm$ 0.06     2.11 $\pm$ 0.01     5.5 $\pm$ 0.5     0.161 $\pm$ 0.002