Multiwavelength observations of the blazar BL Lacertae: a new fast TeV gamma-ray flare
A. U. Abeysekara W. Benbow R. Bird T. Brantseg R. Brose M. Buchovecky J. H. Buckley V. Bugaev M. P. Connolly W. Cui M. K. Daniel A. Falcone Q. Feng J. P. Finley L. Fortson A. Furniss G. H. Gillanders I. Gunawardhana M. Hütten D. Hanna O. Hervet J. Holder G. Hughes T. B. Humensky C. A. Johnson P. Kaaret P. Kar M. Kertzman F. Krennrich M. J. Lang T. T. Y. Lin S. McArthur P. Moriarty R. Mukherjee S. O'Brien R. A. Ong A. N. Otte N. Park A. Petrashyk M. Pohl E. Pueschel J. Quinn K. Ragan P. T. Reynolds G. T. Richards E. Roache C. Rulten I. Sadeh M. Santander G. H. Sembroski K. Shahinyan S. P. Wakely A. Weinstein R. M. Wells P. Wilcox D. A. Williams B. Zitzer S. G. Jorstad A. P. Marscher M. L. Lister Y. Y. Kovalev A. B. Pushkarev T. Savolainen I. Agudo S. N. Molina J. L. Gómez V. M. Larionov G. A. Borman A. A. Mokrushina M. Tornikoski A. Lähteenmäki W. Chamani S. Enestam S. Kiehlmann T. Hovatta P. S. Smith P. Pontrelli




Qi Feng, Olivier Hervet, Svetlana Jorstad

qifeng@nevis.columbia.edu, ohervet@ucsc.edu, jorstad@bu.edu





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


Fred Lawrence Whipple Observatory, Harvard-Smithsonian Center for Astrophysics, Amado, AZ 85645, USA


Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA


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


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


Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA


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


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


School of Physics, National University of Ireland Galway, University Road, Galway, Ireland


Department of Physics and Astronomy, Purdue University, West Lafayette, IN 47907, USA
Department of Physics and Center for Astrophysics, Tsinghua University, Beijing 100084, China.


Fred Lawrence Whipple Observatory, Harvard-Smithsonian Center for Astrophysics, Amado, AZ 85645, USA


Department of Astronomy and Astrophysics, 525 Davey Lab, Pennsylvania State University, University Park, PA 16802, USA



Physics Department, McGill University, Montreal, QC H3A 2T8, Canada
Physics Department, Columbia University, New York, NY 10027, USA



Department of Physics and Astronomy, Purdue University, West Lafayette, IN 47907, USA


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


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


School of Physics, National University of Ireland Galway, University Road, Galway, Ireland


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


DESY, Platanenallee 6, 15738 Zeuthen, Germany


Physics Department, McGill University, Montreal, QC H3A 2T8, Canada


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


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


Fred Lawrence Whipple Observatory, Harvard-Smithsonian Center for Astrophysics, Amado, AZ 85645, USA


Physics Department, Columbia University, New York, NY 10027, USA


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


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


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


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


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


School of Physics, National University of Ireland Galway, University Road, Galway, Ireland


Physics Department, McGill University, Montreal, QC H3A 2T8, Canada


Department of Physics and Astronomy, Purdue University, West Lafayette, IN 47907, USA


School of Physics, National University of Ireland Galway, University Road, Galway, Ireland


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


School of Physics, University College Dublin, Belfield, Dublin 4, Ireland


Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA


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


Enrico Fermi Institute, University of Chicago, Chicago, IL 60637, USA


Physics Department, Columbia University, New York, NY 10027, USA


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


DESY, Platanenallee 6, 15738 Zeuthen, Germany


School of Physics, University College Dublin, Belfield, Dublin 4, Ireland


Physics Department, McGill University, Montreal, QC H3A 2T8, Canada


Department of Physical Sciences, Cork Institute of Technology, Bishopstown, Cork, Ireland


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


Fred Lawrence Whipple Observatory, Harvard-Smithsonian Center for Astrophysics, Amado, AZ 85645, USA


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


DESY, Platanenallee 6, 15738 Zeuthen, Germany


Department of Physics and Astronomy, Barnard College, Columbia University, NY 10027, USA
Department of Physics and Astronomy, University of Alabama, Tuscaloosa, AL 35487, USA


Department of Physics and Astronomy, Purdue University, West Lafayette, IN 47907, USA


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


Enrico Fermi Institute, University of Chicago, Chicago, IL 60637, USA


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


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


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


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


Physics Department, McGill University, Montreal, QC H3A 2T8, Canada


(The VERITAS Collaboration)


 
Institute for Astrophysical Research, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA 
Astronomical Institute, St.Petersburg State University, Universitetskij Pr. 28, Petrodvorets, 198504 St.Petersburg, Russia

 
Institute for Astrophysical Research, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA 

 
Purdue University, 525 Northwestern Avenue, West Lafayette, IN 47907, USA 

 
Astro Space Center of Lebedev Physical Institute, Profsoyuznaya 84/32, 117997 Moscow, Russia 
Moscow Institute of Physics and Technology, Dolgoprudny, Institutsky per., 9, Moscow region, 141700, Russia
Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany 

 
Crimean Astrophysical Observatory, 98409 Nauchny, Crimea, Russia 
Astro Space Center of Lebedev Physical Institute, Profsoyuznaya 84/32, 117997 Moscow, Russia 

 
Aalto University Metsähovi Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä, Finland
Aalto University Department of Electronics and Nanoengineering, PL 15500, FI-00076 Aalto, Finland 
Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany 

 
Instituto de Astrofísica de Andalucía (CSIC), Apartado 3004, E–18080 Granada, Spain 

  
Instituto de Astrofísica de Andalucía (CSIC), Apartado 3004, E–18080 Granada, Spain 

  
Instituto de Astrofísica de Andalucía (CSIC), Apartado 3004, E–18080 Granada, Spain 

  
Astronomical Institute, St.Petersburg State University, Universitetskij Pr. 28, Petrodvorets, 198504 St.Petersburg, Russia

  
Crimean Astrophysical Observatory, 98409 Nauchny, Crimea, Russia


Astronomical Institute, St.Petersburg State University, Universitetskij Pr. 28, Petrodvorets, 198504 St.Petersburg, Russia



Aalto University Metsähovi Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä, Finland




Aalto University Metsähovi Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä, Finland


Aalto University Department of Electronics and Nanoengineering, PL 15500, FI-00076 Aalto, Finland 


Tartu Observatory Observatooriumi 161602 Tõravere, Estonia




Aalto University Metsähovi Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä, Finland


Aalto University Department of Electronics and Nanoengineering, PL 15500, FI-00076 Aalto, Finland 





Aalto University Metsähovi Radio Observatory, Metsähovintie 114, FI-02540 Kylmälä, Finland


Aalto University Department of Electronics and Nanoengineering, PL 15500, FI-00076 Aalto, Finland 



Owens Valley Radio Observatory, California Institute of Technology, Pasadena, CA 91125, USA


Tuorla Observatory, Department of Physics and Astronomy, University of Turku, Väisäläntie 20, 21500 Kaarina, Finland


Steward Observatory, University of Arizona, Tucson, AZ 85716, USA


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




Combined with very-long-baseline interferometry measurements, the observations of fast TeV gamma-ray flares probe the structure and emission mechanism of blazar jets. 
However, only a handful of such flares have been detected to date, and only within the last few years have these flares been observed from lower-frequency-peaked BL Lac objects and flat-spectrum radio quasars. 

We report on a fast TeV gamma-ray flare from the blazar BL Lacertae observed by VERITAS, with a rise time of ∼2.3 hr and a decay time of ∼36 min. 
The peak flux above 200 GeV is (4.2 ± 0.6) × 10^-6 photon m^-2 s^-1 measured with a 4-minute-binned light curve, corresponding to ∼180% of the flux which is observed from the Crab Nebula above the same energy threshold. 
Variability contemporaneous with the TeV gamma-ray flare was observed in GeV gamma-ray, X-ray, and optical flux, as well as in optical and radio polarization. 
Additionally, a possible moving emission feature with superluminal apparent velocity was identified in VLBA observations at 43 GHz, potentially passing the radio core of the jet around the time of the gamma-ray flare. 
We discuss the constraints on the size, Lorentz factor, and location of the emitting region of the flare, and the interpretations with several theoretical models which invoke relativistic plasma passing stationary shocks. 








§ INTRODUCTION
 

BL Lac objects belong to a subclass of radio-loud active galactic nuclei (AGNs), known as blazars. They are characterized by featureless optical spectra, non-thermal broadband spectra, and rapid variability, which jointly suggest that their emission originates in relativistic jets closely aligned to our line of sight <cit.>. 

Fast variability at very high energies (100 GeV ≲ E_γ≲ 100 TeV; VHE), with timescales as short as a few minutes, has been observed in several blazars <cit.>, including the prototypical BL Lacertae <cit.> located at redshift z=0.069 <cit.>. 

Long-term monitoring of BL Lacertae 
has led to no detection of the source in the TeV gamma-ray band by the current generation of instruments except during flaring episodes, when its flux has been observed to reach >100% of the Crab Nebula flux (C. U.) above 1 TeV in 1998 <cit.>, ∼0.03 C. U. above 200 GeV in 2005 <cit.>, and most recently ∼1.25 C. U. above 200 GeV with a short variability timescale of 13±4 min in 2011 <cit.>. 


The rapid gamma-ray variability observed in TeV blazars implies very compact emitting regions, as well as low gamma-ray attenuation by pair production on infrared/optical photons near the emission zone. 

While a one-zone synchrotron self-Compton (SSC) model, one of the simplest blazar models <cit.>, has 
been effective at explaining emission from high-frequency-peaked BL Lac (HBL) objects, the intrinsic pair-production opacity of a relativistic emission zone in such a model depends on its size and Doppler factor, and on the density of lower-energy photons. 

Therefore, if the synchrotron photons are the main source of the lower-energy radiation, the emitting region must have a small size and/or a large Doppler factor so that the gamma rays can escape pair production. 

Alternatively, if an external photon field (e.g., the broad-line region; BLR) dominates the lower-energy radiation, it can cause substantial gamma-ray absorption. As a result, the emitting region is generally expected to be far away from the central region of the AGN, especially for flat-spectrum radio quasars (FSRQs) 
whose broad-line emission is relatively strong. 

BL Lacertae was first classified as a low-frequency-peaked BL Lac (LBL) object as the synchrotron peak frequency of its spectral energy distribution (SED) was measured to be 2.2 × 10^14 Hz <cit.>, but was later reclassified as an intermediate-frequency-peaked BL Lac (IBL) object <cit.>. 
It has been reported that the SEDs of several IBLs/LBLs cannot be well described by a one-zone SSC model <cit.>, and more complex models such as multi-zone SSC models or external-radiation Compton (ERC) models are needed. 





The large Doppler factor and/or distant downstream emitting region required by the observed fast TeV variability of blazars, together with the knotty jet structures (both moving and stationary) identified with high-resolution radio observations <cit.>, can be explained consistently by theoretical models with multiple emitting zones that are either spatially or temporally separated, 
e.g., structured jets <cit.>, jet deceleration <cit.>, jets in a jet <cit.>, and plasma passing a standing shock <cit.>. 


However, the details regarding the location and the mechanism of blazar emission are still not well understood <cit.>. 
Simultaneous multiwavelength (MWL) observations can provide insights into the flaring mechanisms (e.g., leptonic or hadronic processes) 
of these objects, particularly at the wavelengths where SEDs often peak. 


In practice, such observations are limited in the case of fast flares at sub-hour timescales,
even with dedicated strategies <cit.>. 

Nevertheless, contemporaneous radio data are often relevant because the radio variability timescale is usually much longer <cit.>. 
In particular, the evolution of polarization (both radio and optical) before and after a gamma-ray flare provides information about magnetic field structures of the jet, and therefore the activity of possible gamma-ray emitting regions <cit.>. 

BL Lacertae exhibits both stationary radio cores/knots and superluminal radio knots <cit.>. 
Possible associations between the variability of superluminal radio knots and gamma-ray flares have been investigated for BL Lacertae <cit.> and other blazars <cit.>. 


On 2016 Oct 5, we observed BL Lacertae at an elevated flux level with sub-hour variability with the Very Energetic Radiation Imaging Telescope Array System (VERITAS; 
see Section <ref>). 
A series of observations with the Very Long Baseline Array (VLBA) at 43 GHz and 15.4 GHz was performed 
over a few months before and after the gamma-ray flare, revealing a possible knot structure emerging around the time of the TeV flare (see Section <ref>).
In this work, we report on the results of the VERITAS, VLBA, and other MWL observations and discuss their implications. 
The cosmological parameters assumed throughout this paper are Ω_m=0.27, Ω_Λ=0.73, and H_0=70 km s^-1 Mpc^-1 <cit.>. At the redshift of BL Lacertae, the luminosity distance and the angular size distance are 311 Mpc and 273 Mpc, respectively, and the angular scale is 1.3 pc/mas. 






§ OBSERVATIONS, DATA ANALYSIS, AND RESULTS






 §.§ VERITAS





VERITAS is an array of four imaging atmospheric-Cherenkov telescopes located in southern Arizona <cit.>. It is sensitive to gamma rays in the energy range from 85 GeV to >30 TeV with an energy resolution of ∼15% (at 1 TeV) and is capable of making a detection with a statistical significance of 5 standard deviations (5σ) of a point source of 0.01 C. U. in ∼25 hr. 




BL Lacertae was observed at an elevated TeV gamma-ray flux by VERITAS on 2016 Oct 5 as part of an ongoing monitoring program, 
and follow-up observations were immediately instigated 
based on a real-time analysis. The total exposure of these observations amounts to 153.5 min after data-quality selection, with zenith angles ranging between 11^∘ and 30^∘. 


The data were analyzed using two independent analysis packages <cit.> 
and a pre-determined set of cuts optimized for lower-energy showers <cit.>. A detection with a statistical significance of 71σ was made from the data on the night of the flare, with a time-averaged integral flux above 200 GeV of (2.24 ± 0.06) × 10^-6 photon m^-2 s^-1 (or ∼0.95 C. U.).





  §.§.§ VHE gamma-ray flux variability and the modelling of the flare profile




Figure <ref> shows the VERITAS TeV gamma-ray light curve of BL Lacertae above 200 GeV on 2016 Oct 5 with 4-minute and 30-minute bins. 
A gradual rise of the TeV flux by a factor of ∼2 followed by a faster decay was observed. 
The measured peak flux for the 30-minute-binned light curve is (3.0 ± 0.2) × 10^-6 photon m^-2 s^-1, corresponding to ∼1.25 C. U., and that for the 4-minute-binned light curve is (4.2 ± 0.6) × 10^-6 photon m^-2 s^-1, or ∼1.8 C. U.

We first fitted the 4-minute-binned VERITAS light curve with a constant-flux model, obtaining a χ^2 value of 170.8 for 45 degrees of freedom (DOF), corresponding to a p-value of 1.1×10^-16 and rejecting the constant-flux hypothesis. 

To quantify the rise and decay times of the TeV flare, we then fitted the VHE gamma-ray light curve with a piecewise exponential function as follows: 

    F(t) = 
    
    
    F_0 e^(t - t_peak)/t_rise,     t ⩽ t_peak;  
    
    F_0 e^-(t- t_peak)/t_decay,     t > t_peak;

where 
F_0 is the peak flux, t_peak is the time of the peak flux, and t_rise and t_decay are the rise and decay times, respectively, on which the flux varies by a factor of e. 


The optimal values of the parameters and their uncertainties were determined from the posterior distributions obtained from Markov chain Monte Carlo (MCMC) simulations, for which the Python package  <cit.> was used. 
The MCMC chain contains 100 random walkers in the parameter space initialized with a uniform random prior. Each random walker walks 4000 steps, the first 2000 steps of which are discarded as the “burn-in” samples. This amounts to 2×10^5 effective MCMC simulations. 

A proposal scale parameter was chosen so that the mean proposal acceptance fraction is 37%, ensuring an adequate yet efficient sampling of the posterior distributions. 
Note that the parameters are bounded to be positive, so that they are physically meaningful, and sufficiently large upper bounds were also provided for computational efficiency. 
After the posterior distributions were obtained, kernel density estimation with Gaussian kernels of bandwidths equal to 1% of the range of the corresponding parameter was used to estimate the most likely value (maximum a posteriori) and the 68% confidence interval of each parameter. 


The joint posterior distributions of the parameters from the MCMC sampling are shown in Figure <ref>. 
The diagonal plots show the posterior probability distributions of each parameter, some of which (e.g., t_peak) appear non-Gaussian. 
Correlations between t_peak and t_rise, as well as between t_peak and t_decay, are also apparent in the off-diagonal joint distributions. 
The best-fit model and the 99% confidence intervals from the MCMC sampling are shown in Figure <ref>. 
The rise and decay times of the flare are determined to be 140^+25_-11 min and 36^+8_-7 min, respectively. The best-fit peak time and flux are 130^+5_-3 min (after MJD 57666.165) and 3.4^+0.2_-0.2× 10^-6 photon m^-2 s^-1, respectively. 


Further VERITAS observations of BL Lacertae were made on 2016 Oct 6 with 37.6-minute live exposure, and from Oct 22 to Nov 19 with 294.6-minute live exposure, after data-quality selection; neither of these sets of observations led to a detection of the source (signal significances of only 2.6σ and 0.9σ, respectively). 
The integral flux upper limits (shown in the first panel of Figure <ref>) between 0.2 and 30 TeV at 99% confidence level from the observations on Oct 6 and between Oct 22 and Nov 19 were obtained as 2.0×10^-7 photon m^-2 s^-1 and 2.8×10^-8 photon m^-2 s^-1, respectively, assuming a power-law spectrum with a photon index of 3.3 (see Section <ref>). 

Motivated by the existence of multiple radio emission zones identified in VLBA data (see Section <ref>) and several multi-zone models for BL Lacertae that are consistent with past observations <cit.>, we also fitted the light curve with a model including an additional constant-flux baseline. 
In a multi-zone model, different zones can be of different sizes and vary independently on different timescales. Therefore, it is possible to have a larger emitting zone that varies slowly which can be adequately described by a constant-baseline component on the timescale considered, and a smaller, more energetic zone that is responsible for the rapid flare described by the exponential components. 
With the more complex model, the best-fit decay time is only 2.6^+6.7_-0.8 min, with a baseline flux of 1.2^+0.1_-0.2×10^-6 photon m^-2s^-1. This best-fit baseline flux is higher than the upper limit obtained from the observations on the next day, 
indicating that the slower component would be required to vary on timescales of ∼1 day, 

consistent with the GeV gamma-ray observations (Section <ref>). However, we would like to highlight that it is not possible to unambiguously reject either model based on the statistics. 







  §.§.§ The VHE spectrum



A power-law fit to the VERITAS spectrum of BL Lacertae 
yields a best-fit photon index of 3.28 and a reduced χ^2 value χ^2/DOF=30.6, indicating that a simple power law does not adequately describe the spectrum. 










A log-parabola model 
fits the VERITAS spectrum better: 

    dN/dE   = (2.22± 0.07)× 10^-5
       ×( E/0.2 TeV)^[ -(2.4 ± 0.1) -  ( 1.8 ± 0.3 ) log_10 (E/0.2 TeV)  ]
       m^-2 s^-1 TeV^-1,

with χ^2/DOF=1.6. 

After de-absorbing the VHE spectrum using the optical depths for a source at a redshift of 0.069 according to the extragalactic background light model in <cit.>, the best-fit log-parabola model becomes: 

    dN/dE   = (2.36± 0.07)× 10^-5
       ×( E/0.2 TeV)^[ -(2.2 ± 0.1) -  ( 1.4 ± 0.3 ) log_10 (E/0.2 TeV)  ]
       m^-2 s^-1 TeV^-1,

with χ^2/DOF=1.7. 
The observed and de-absorbed TeV gamma-ray spectra are shown together with the GeV gamma-ray spectra (Section <ref>) in Figure <ref> in the ν F_ν representation. 





 §.§ Fermi-LAT




The Large Area Telescope (LAT) on board the Fermi satellite is a pair-conversion gamma-ray telescope sensitive to energies from ∼20 MeV to >300 GeV <cit.>. 


An unbinned likelihood analysis was performed with the LAT ScienceTools  and Pass-8  instrument response functions <cit.>.  class events with energy between 100 MeV and 300 GeV within 10^∘ from the position of BL Lacertae were selected. 
For the short durations of interest to the TeV flare, a simple model containing BL Lacertae, another point source 3FGL J2151.6+4154 ∼2^∘ away from BL Lacertae, and the contributions from the Galactic () and isotropic () diffuse emission were included. A maximum zenith angle cut of 90^∘ was applied. We checked in the residual test-statistics map that no significant excess was left unaccounted for within the model. For the short durations, a power law was used to model BL Lacertae instead of the log-parabola model used in the  catalog. 
We verified with an analysis using a log-parabola spectral model and obtained consistent flux values.

For the light curve shown in the second panel of Figure <ref>, an unbinned likelihood analysis was performed on each one-day interval, leaving the normalizations and power-law indices of BL Lacertae and 3FGL J2151.6+4154 free, as well as the normalization of the diffuse components. The source was in an elevated GeV gamma-ray state when the TeV flare was observed, although the GeV flux varied on a much longer timescale. 
An exponential fit to a 15-day interval around the TeV gamma-ray flare yields a rise time of 2.1±0.2 days and a decay time of 7±2 days. 



The gamma-ray SEDs measured by the Fermi-LAT and VERITAS on the night of the TeV flare are shown in Figure <ref>. 
In order to obtain the GeV gamma-ray SEDs, we used the user-contributed tool [https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/] to perform the unbinned likelihood analysis in several energy bands. 

The power-law index that gives the best fit to the LAT data completely simultaneous with VERITAS is 1.83±0.21, which is similar to that of the three-day binned LAT data, 1.85±0.07. These values indicate a harder GeV gamma-ray spectrum during the flare than that reported in the  catalog <cit.>. 

Both the GeV and TeV gamma-ray spectral indices of this flare in 2016 are comparable to those of the flare in 2011 <cit.>, 
and they suggest that the peak energy of the gamma-ray SED during the flare is between 5 GeV and 100 GeV. 




 §.§ Swift XRT



The X-Ray Telescope (XRT) on board the Swift satellite is a grazing-incidence focusing X-ray telescope, and is sensitive to photons in the energy range 0.2–10 keV <cit.>. 

Follow-up observations of BL Lacertae were carried out with the Swift-XRT on 

2016 Oct 6, 7, and 8; the only other XRT observations within a 45-day window around the time of the VHE flare were made on 2016 Oct 27 and Nov 2. 
The XRT data, taken in the photon-counting (PC) mode, were analyzed using the  package (v6.19). 

The data were first processed using  (v0.13.2) with calibration database (). 

The count rates in the PC mode were >0.5 count s^-1, and the effect of potential pile-up was checked for all observations by fitting a King function to the point-spread functions (PSFs) at >15 arcsecond <cit.>. Those central pixels where the data fall below the model curve, indicating pile-up, were excluded. 

For the observations on 2016 Oct 6, the King function agrees with the data even on the brightest pixels. Therefore, a circular source region of a radius of 20 pixels centered on BL Lacertae was used. 
For the data taken on 2016 Oct 7 and 8, annular source regions were used, with inner radii of four and two pixels, and an outer radius of 20 pixels. 
For all three observations, an annular background region with inner and outer radii of 70 and 120 pixels, respectively, was used. 
Note that source regions excluding the central two and four pixels were also tested for the observations on 2016 Oct 6, and consistent results were obtained. Therefore, we are confident that no bias was introduced by the different exclusion regions used for pile-up correction. 

The observations on Oct 7 consisted of two intervals of duration 486 s and 1422 s, separated by roughly one satellite orbital period (∼90 min). A sustained dark stripe (likely due to bad CCD columns) appears in the XRT image near the position of BL Lacertae, contaminating the second interval. Therefore, we conservatively chose to use only the data recorded during the first interval. The image and spectrum of each ∼3 min of this relatively short exposure were checked for data quality, and no anomaly was found. 

ccccc[h]
5

Swift-XRT spectral-fit results using the absorbed-power-law model described in Equation <ref>, with N_H free and fixed. The errors quoted denote 68% confidence intervals. 




Date     α     
K     N_H     χ^2/DOF 
 
           
10^-2 keV^-1cm^-2s^-1     10^21 cm^-2      
 
Oct 6		    2.5 ± 0.1 	     0.62^+0.07_-0.06 	    2.7^+0.3_-0.3	    0.83 	

       Oct 7 		    2.1 ± 0.1	     4.6^+0.6_-0.5 	    3.1^+0.5_-0.4	    1.07 	

       Oct 8 		    2.3 ± 0.1	     0.43^+0.06_-0.05	    3.0^+0.5_-0.4	    0.54 	
        
       Oct 27 	    1.4^+0.4_-0.3	     0.14^+0.08_-0.04	    1.2^+2.2_-1.2	    0.48 	
 
       Nov 2		    1.3 ± 0.3	     0.11^+0.04_-0.03	    1.4^+1.1_-0.8	    0.33 	
 Oct 6		    2.54 ± 0.07 	     0.65 ± 0.03 		    2.9 (fixed)    0.82 	

       Oct 7 		    2.08 ± 0.07	     4.34  ± 0.24 	    2.9 (fixed)    1.04	

       Oct 8 		    2.26 ± 0.08	     0.41 ± 0.02		    2.9 (fixed)    0.52 	
   
       Oct 27 	    1.64 ± 0.20	     0.19  ± 0.03 	    2.9 (fixed)    0.50	

       Nov 2 	    1.76 ± 0.19	     0.17  ± 0.02 	    2.9 (fixed)    0.63	

       
       
       




Ancillary response files were generated using the  task with the response matrix file . 
The spectrum was fitted with an absorbed-power-law model (): 

    dN/dE = e^-N_Hσ(E) K ( E/1 keV)^-α,

where N_H is the column density of neutral hydrogen, σ(E) is the photoelectric cross-section, and K and α are the normalization and index of the power-law component, respectively. 
The best-fit values of the parameters are shown in Table <ref>. Note that the best-fit values of N_H are in agreement with the archival results from X-ray spectral fit <cit.>, 
but are larger than the value N_H = 1.8  ×10^21 cm^-2 from the Leiden/Argentine/Bonn (LAB) survey of Galactic HI <cit.>. This difference is likely due to the additional contribution of the Galactic molecular gas (e.g., CO emission), as BL Lacertae is relatively close to the Galactic plane (with a Galactic latitude b=-10.44^∘) <cit.>. 


As the N_H value is expected to stay constant over the period of interest, we also fit the same model with N_H fixed at 
the average best-fit value N_H = 2.9  ×10^21 cm^-2 over the three nights of XRT observations, to better constrain the spectral index and normalization. 
We also investigated an absorbed-log-parabola model with N_H fixed at 2.9 ×10^21 cm^-2 to fit the X-ray spectra. We see no evidence for spectral curvature, as the best-fit log-parabola model reduces to a power law.
The X-ray SEDs of BL Lacertae measured on 2016 Oct 6, 7, and 8 are shown in Figure <ref>. 
The X-ray emission from the source was strongest and hardest on 2016 Oct 7 (two days after the TeV gamma-ray flare) compared to the day before and the day after (see Table <ref>). 
The energy flux values based on the best-fit absorbed-power-law model between 0.3 keV and 10 keV on 2016 Oct 6, 7, and 8 were (1.4± 0.1), (14.2± 0.9), and (1.1± 0.1) ×10^-11 erg cm^-2 s^-1, respectively, as shown in the third panel of Figure <ref> along with the results from two more observations taken on 2016 Oct 27 and Nov 2. Note that the






 §.§ Optical facilities





BL Lacertae was monitored in the R-band at high cadence by a number of optical facilities, including the Steward Observatory[<http://james.as.arizona.edu/ psmith/Fermi>] <cit.>, the AZT-8 reflector of the Crimean Astrophysical Observatory, the Perkins telescope[<http://www.lowell.edu/research/research-facilities/1-8-meter-perkins/>], the LX-200 telescope in St. Petersburg, Russia, 

and the Calar Alto 2.2-m Telescope (with observations obtained through the MAPCAT[MAPCAT stands for the Monitoring of AGN with Polarimetry at the Calar Alto Telescopes, see: <http://www.iaa.es/ iagudo/_iagudo/MAPCAT.html>] program) in Almería, Spain. 

We also included the r'-band observations made with the 48-inch telescope at the Fred Lawrence Whipple Observatory (FLWO). We transformed the r'-band flux to R-band  using the color index V-R=0.73 ± 0.19 for BL Lacertae <cit.> and the 
transformation r'-R = 0.19 (V-R) + 0.13 <cit.>. 
The dominant uncertainty from the conversion to R-band comes from the variability of the V-R color index, resulting in an additional systematic uncertainty of ∼0.04 magnitude, which is included in the converted FLWO R-band magnitude shown in Figure <ref>. 
Any variability in the color index Δ(V-R)<0.19 during the epoch of observations shown would lead to a shift of the FLWO R-band magnitude within the error bars shown. 

The R-band flux and polarization measurements contemporaneous with the gamma-ray flare are shown in Figure <ref>. 
The lower three panels (from top to bottom) show the R-band magnitude, polarization fraction, and electric vector position angle (EVPA) of the source, respectively. 

A -180^∘ shift is applied to all the EVPA measurements before MJD 57662, so that the EVPA difference between MJD 57662 and 57658 is reduced to ∼80^∘ from ∼100^∘ before the shift was applied <cit.>. 

The measurements are reasonably consistent between the various instruments. 


The R-band flux from the source varied in a similar manner to the GeV flux, with an increase  
observed a few days before the VHE flare. 
The optical EVPA appeared to have rotated smoothly from roughly perpendicular to the position angle (PA) of the jet in late 2016 Sept to roughly parallel in late 2016 Oct, except for three days before the TeV gamma-ray flare when the optical EVPA was nearly aligned with the PA of the jet. This was followed by a sudden decrease in EVPA on the day before the TeV gamma-ray flare. 
The fractional polarization was relatively low around the time of the TeV flare, and increased to the highest value of the 45-day period in late Oct when the EVPA was again aligned with the jet. 











 §.§ Radio facilities













BL Lacertae was observed throughout the period of interest at 43 GHz with the Very Long Baseline Array under the VLBA-BU-BLAZAR monitoring program <cit.> and at 15.4 GHz with the Monitoring Of Jets in Active galactic nuclei with VLBA Experiments (MOJAVE) program <cit.>. The 43-GHz and 15.4-GHz VLBA data calibration and imaging procedures were identical to those described by <cit.> and <cit.>, respectively. 




Figure <ref> presents 43-GHz images of the parsec-scale jet of BL Lacertae at five epochs from 2016 Sept 5 to Dec 23. The second epoch, 2016 Oct 6, took place only one day after the VHE flare. The images are convolved with a circular Gaussian restoring beam with a full width at half maximum (FWHM) of 0.1 mas, 
which is similar to the angular resolution of the longest baselines along 
the (southern) direction of the jet. 
We note that the Oct 6 observation was affected by equipment failure at the Mauna Kea and Hancock antennas, at the extremities of the array, although this degraded the north-south angular resolution by only 14%. The corresponding linear resolution at the redshift of BL Lacertae is 0.13 pc in projection on the sky and 1.8^+0.8_-0.4 pc if we adopt a viewing angle of 4.2^∘±1.3^∘ between the jet axis and line of sight <cit.>. 

As was the case in previous observations <cit.>, the main structure of the jet consists of three quasi-stationary brightness peaks, designated as A0, A1 0.12 mas to the south of A0, and A2 0.30 mas to the south of A0. The locations of A1 and A2 appear to fluctuate as moving emission features (frequently referred to as “knots”) with superluminal apparent velocities pass through the region. Such combination of moving and stationary emission components complicates the interpretation of the changing structures of the total and polarized intensities. Because of this, the interpretation that we offer to explain the variations within the images is not unique. 

We ignore the effects of Faraday rotation on the polarization EVPA, which <cit.> estimated to be low (-16^∘) between 43 GHz and 300 GHz. 
It is worth mentioning that <cit.> measured a much lower (by an order of magnitude) Faraday rotation using 8-GHz to 15-GHz observations. 
This could be due to a combination of a possible variability in the rotation measure and a decrease of the rotation measure with distance from the central black hole <cit.>, as the core at 15 GHz is located further away from the black hole compared to the core at 43 GHz due to the effect of opacity. 

A knot of emission with enhanced polarization at 43 GHz, which we designate as K16, appears to propagate down the jet. 

Its centroid moves from ∼0.05 mas south of A0 on Oct 23 to ∼0.28 mas from A0 on Dec 23. 
This corresponds to an apparent speed of 6c, within the range typically observed in BL Lacertae <cit.>. Extrapolation back to Oct 6 places the knot K16 0.01 mas north of the centroid of A0, within the A0 emission region characterized by its angular size of 0.03±0.02 mas <cit.>. This implies that the VHE flare occurred as the moving knot crossed the stationary “core”, which <cit.> have interpreted as a standing shock located ∼1 pc from the central black hole.






The VLBA images at 15.4 GHz, as shown in Figure <ref>, reveal the evolution of the jet structures 
further away from the central source and on a larger spatial scale, as a result of optical depth and angular resolution, respectively, compared with the observations at 43 GHz. 
Therefore, a delay is expected between the measurements at these two frequencies. 
The polarized intensity of the stationary core of BL Lacertae at 15.4 GHz reached a minimum on 2016 Dec 26 and gradually increased, 
with a potentially bright feature with distinct polarization angle (consistent with the EVPAs measured at 43 GHz on 2016 Dec 23), 
which may correspond to the knot K16 observed at 43 GHz earlier, appearing at ∼1 mas southwest of the core. 
This is consistent with past observations of the same source with the VLBA at different frequencies reported by <cit.>, where new components of the jet were seen to fade as they separated from the core, disappearing at  ∼0.7 mas 
and reappearing at ∼1 mas. 


We show the fractional polarizations and EVPAs of the 43-GHz and 15-GHz core in the bottom two panels of Figure 5, along with the R-band results. 
The EVPAs of the core at 43 GHz and 15 GHz are roughly consistent with the PA of 10^∘ of the jet over the course of a few months since 2016 Sept. This implies that the magnetic field is toroidal/helical at the core, as we discuss in Section <ref>. 
Further downstream in the jet, the EVPAs become more perpendicular to the PA of the jet, as shown in Figure <ref>. Such location-dependent radio EVPAs help us to interpret the dominant optical component based on the optical EVPA data. 



We show in Figure <ref> the evolution of the total flux density of BL Lacertae, measured by the Metsähovi Radio Observatory (MRO) at 37 GHz and by Owens Valley Radio Observatory (OVRO) at 15 GHz, respectively, over about a year, as well as their z-transformed discrete cross-correlation <cit.>. 

The 37-GHz observations were made with the 13.7-m diameter Aalto
University Metsähovi radio telescope, which is a radome-enclosed
Cassegrain-type antenna situated in Finland. The measurements were made with a 1-GHz-bandwidth dual-beam
receiver centered at 36.8 GHz. The 
high electron mobility pseudomorphic transistor front end operates at room temperature. 


The 37-GHz observations are Dicke-switched ON–ON observations, alternating the
source and the sky in each feed horn to remove atmospheric and ground contamination. 
Typical integration time to
obtain one flux-density data point is between 1200 s and 1600 s. The
detection limit of the telescope at 37 GHz is on the order of 0.2 Jy
under optimal conditions. Data points with a signal-to-noise ratio < 4
are handled as non-detections.

The flux-density scale is set by observations of the HII region of DR 21, while NGC 7027, 3C 274 and 3C 84 are used as secondary
calibrators. A detailed description of the data reduction and analysis
can be found in <cit.>. The error estimate in
the flux density includes the contributions from the measurement rms
and the uncertainty of the absolute calibration.



The OVRO 40-m telescope uses off-axis dual-beam optics and a cryogenic pseudo-correlation receiver with a 
15.0-GHz center frequency and 3-GHz bandwidth. 
The source is alternated between the two beams in an ON–ON fashion to remove atmospheric and ground 
contamination. The fast gain variations are corrected using a 180^∘ phase switch. 

Calibration is achieved using a temperature-stable diode noise source to compensate for receiver gain drifts, and the flux-density scale is derived from observations of 3C 286 assuming the value of 3.44 Jy at 15.0 GHz reported by <cit.>. The systematic uncertainty in the flux density scale is ∼5%, which is not included in the error bars in Figure <ref>. Complete details of the reduction and calibration procedures are given in <cit.>.


The 37-GHz and 15-GHz light curves show that at the time of the TeV gamma-ray flare, BL Lacertae was transitioning from a steady radio flux state to a flaring state that lasted for about five months. 
The ZDCF shows no significant detection of any time lag between the fluxes at the two frequencies, suggesting that both observations are dominated by the flux from a region that is optically thin at 15 GHz. 









§ DISCUSSION




For the second time, VERITAS has detected a fast gamma-ray flare from BL Lacertae. 



While no information was obtained from the rising phase of the first VHE flare in 2011 <cit.>, the VERITAS measurements during the 2016 flare described in this work cover both the rise and decay phases of the flare. 





 §.§ On the size of the gamma-ray-emitting region

The fastest timescale of a flare (in this case the decay time) provides a constraint on the size R of the emitting region, as 

    R≤ct_decayδ/1+z,

where c is the speed of light, δ is the Doppler factor of the jet[δ = [Γ (1-βcosθ)]^-1, where
Γ is the bulk Lorentz factor of the jet, and θ is the angle between the axis of the jet and the line of sight.], and z is the redshift of the source. 

The mass of the central black hole (M_BH) of BL Lacertae was estimated to be ∼3.8× 10^8 M_⊙ by <cit.> using the R-band absolute magnitude and the empirical correlation between black hole mass and bulge luminosity of the host galaxy <cit.>. 

The corresponding Schwarzschild radius R_s of the central black hole of BL Lacertae is ∼1.1× 10^12 m (∼ 3.6× 10^-5 pc). 

It is worth noting that the mass measurement of a black hole is a challenging task, and M_BH values in the range (0.16 - 5.01)× 10^8 M_⊙ have been reported for BL Lacertae <cit.>. 

The Doppler factor δ was estimated to be ∼24 according to the method described in <cit.> from the propagation of a possible perturbation in the radio jet observed with VLBA at 15 GHz <cit.>, assuming a viewing angle of 2.2^∘ based on radio apparent-velocity measurements. Taking the best-fit value of t_decay=36^+8_-7 min (see Section <ref>) and using Equation <ref>, we estimate the upper limit on the size of the emitting region to be R ≲11.9 R_s. 





















 §.§ On the gamma-ray flare profile



An asymmetric profile with a faster decay of the VHE gamma-ray flux was observed in the flare, which would be caused by an abrupt cessation of the high-energy particle injection <cit.>. In this scenario, the flaring activity is attributed to fresh injection of high-energy particles into the emitting region instead of in situ acceleration of the particles. However, minimal variability in the radio band would be observable for this interpretation. 
Since strictly simultaneous radio observations were not performed, we cannot draw any conclusions regarding the radio variability at the time of the TeV gamma-ray flare. However, we note that the observed gamma-ray flare profile and the longer-term radio light curves (Figure <ref>) are consistent with the model proposed by <cit.>. In this model, a fast gamma-ray flare can be produced by a small plasmoid in the magnetic reconnection layer, with no concurrent radio flares from the single plasmoid but a delayed radio flare powered by the entire reconnection event. The delay timescale is expected to correspond to the duration of the reconnection event, typically a few weeks. 

The asymmetric flare is in contrast to the more frequently observed flaring profile, a fast rise followed by a slow decay, which can be the manifestation of in situ acceleration and/or a longer cooling time (i.e., longer than the acceleration time) associated with a steep particle-energy distribution <cit.>. 

BL Lacertae showed an enhancement in its GeV gamma-ray flux at the time of the TeV flare, but on a longer timescale of a few days. It also exhibited high X-ray flux on 2016 Oct 7 (two days after the TeV flare), 
about a factor of 10 stronger than the flux on Oct 6 and 8. 

These observations indicate efficient acceleration of relativistic particles in the jet to at least a few hundred GeV. We note, however, that the delayed X-ray flare may or may not be related to the TeV gamma-ray flare, 
since the lack of strictly simultaneous X-ray data precludes us from ruling out the possibility of an X-ray flare simultaneous with the TeV gamma-ray one. 
The different variability timescales of the observed TeV and GeV gamma rays give a hint that they may originate from different emitting zones. 
One possibility is that the GeV gamma rays were produced by particles injected into and accelerated in a large shock region <cit.>, while the TeV gamma rays were produced through magnetic reconnection in a localized region <cit.>. 


ccccccccccccccc[ht]
15
 The parameters used to calculate the constraints shown in Figure <ref>. 




L_gamma     L_syn     L_da     
t_var     M_BHb      E_cool      
δ/Γ      ϵ_BLR    ϵ_IR     
r_BLR     r_IR      
E_BLR     E_IR     
g_SSC     g_ERC  

(erg s^-1)     (erg s^-1)      (erg s^-1)      (min)      (M_⊙)     
(TeV)      
    (pc)   (pc)   (eV)   (eV)   
                7.8×10^45   7.8×10^44   6.0×10^44   36	   3.8×10^8   
                        1	   1	   0.1	   0.1 	    0.025 	   0.5	   10 	   0.3  	   0.75	   0.5	
                    L_gamma, L_syn, and L_d are the observed gamma-ray luminosity, the synchrotron luminosity, and disk luminosity, respectively; t_var is the observed variability time; M_BH is the mass of the central black hole; E_cool is the energy of the observed photons due to the external Compton cooling of relativistic electrons; δ/Γ is the ratio between the Doppler factor and Lorentz factor of the electrons; ϵ_BLR, r_BLR, and E_BLR are the covering factor, characteristic radius of the BLR, and the energy of BLR photons, respectively; ϵ_IR, r_IR, and E_IR are similar parameters for the IR-emitting torus region; g_SSC and g_ERC are the bolometric correction factors for SSC and ERC mechanisms. a<cit.>b<cit.>
                    
                    
                     §.§ On the Lorentz factor and the location of the gamma-ray-emitting region
                    
                    
                    Without simultaneous MWL observations with temporal resolution comparable to that of the TeV gamma-ray observations, we cannot construct a reliable broadband SED of the source during the TeV flaring state. Instead, we constrain the Lorentz factor (Γ) of the gamma-ray-emitting region 
                    based on the gamma-ray variability, assuming two different emission mechanisms, SSC and ERC. 
                    Both models have been used to describe the broadband SED of BL Lacertae in the past <cit.>. 
                    
                    However, we note that during the flare, the peak of the gamma-ray SED is located between ∼5 GeV and ∼100 GeV, higher than that in the lower flux state <cit.>. 
                    Such behaviour is most frequently observed in FSRQs and can be interpreted with ERC process on IR photons in the torus region <cit.>. Such ERC process was also used to interpret the emission of BL Lacertae in a flaring state <cit.>. 
                    
                    
                    
                    Assuming a one-zone SSC model, we can calculate an opacity constraint on the Doppler factor δ of the TeV gamma-ray emitting region by requiring the pair-production optical depth to be ⩽ 1, following Equation 3.7 and 3.8 in <cit.><cit.>. We found that δ≳ 13 using the following observables: the best-fit decay time of the TeV gamma-ray flare (36 min), the center of the highest-energy bin with significant excess of the TeV gamma-ray spectrum of the source during the flare (∼1.5 TeV), the R-band magnitude inferred from the FLWO observations on the same night (13.17), and the near-infrared spectral index <cit.>. Assuming a viewing angle of 2.2^∘, the constraint δ≳ 13 is equivalent to a constraint on the Lorentz factor of Γ≳ 7. 
                    
                    Assuming the gamma rays are emitted via an ERC process, we can constrain the Lorentz factor Γ and the distance r from the central black hole of the gamma-ray-emitting region
                    following the method described by <cit.>. 
                    
                    The collimation constraint was derived from the requirement Γθ≲1. 
                    Both SSC and ERC processes are considered in the calculation of the SSC constraint, while the majority of the gamma rays are assumed to be produced via ERC process. 
                    For the cooling constraint, only the ERC process on the thermal radiation fields close to the black hole is considered. 
                    
                    This does not take into account any possible inverse-Compton scattering of an external synchrotron field, 
                    which, as we consider below, would loosen the cooling constraint on Γ at large distances from the central black hole. 
                    We also assume that the emitting region is spherically symmetric. It is possible that the emitting region is not spherical (e.g., if it is passing a standing shock), and the constraints on Γ and r may change. 
                    
                    
                    The values of the parameters used for the calculation of the above three constraints on Γ and r are shown in Table <ref>. 
                    Some of the parameters are constrained by observations, and the other parameters are chosen so that a conservative constraint is derived. For example, we set the Compton dominance parameter q=L_gamma/L_syn=10 based on the observed R-band magnitude and the peak flux of the gamma-ray SED, the former of which should provide a good estimation of the peak of the synchrotron flux, considering that the source is a lower-frequency-peaked BL Lac object. The SSC luminosity was set equal to the observed gamma-ray luminosity L_SSC=L_gamma in order to obtain a conservative SSC constraint. We also used a relatively-high observed gamma-ray energy (1 TeV) for a conservative ERC cooling limit. 
                    
                    
                    We note that 
                    changes in the values of the parameters describing the geometry of the external radiation fields, namely the covering factor ϵ_BLR and characteristic radius r_BLR for the BLR, and similarly ϵ_IR and r_IR for the IR torus, which are poorly constrained by observations, could change the cooling constraint. The values of the radii used in this work are derived based on the disk luminosity L_d=6.0×10^44 erg s^-1<cit.> and the relations r_BLR = 1× 10^15√((L_d/10^45erg s^-1)) m≈ 0.025 pc and r_IR = 2× 10^16√((L_d/10^45erg s^-1)) m≈ 0.5 pc<cit.>. 
                    
                    
                    
                    In this analysis, the distance r between the central black hole and the gamma-ray-emitting region is constrained to be ≲ 12.4 pc. 
                    If we fix the Lorentz factor at Γ=24, then we constrain the distance to be 0.01≲ r/pc≲ 0.7. If we fix the distance at r=1 pc, the estimated distance between the core A0 and the central black hole, and assume that the gamma rays are produced as the knot K16 passes the core A0, then the Lorentz factor is only loosely constrained at 35 ≲Γ≲ 226. 
                    
                    At small r values (r≲0.68 pc), the SSC constraint on the lower limit of Γ is stricter than the cooling constraint. At r=r_BLR=0.025 pc (the smallest distance for the VHE-emitting region without heavy absorption from the radiation field in the BLR), we put a strong lower limit on the Lorentz factor Γ≳10.1, which is larger than the archival values (∼5–7) derived from radio observations <cit.>, and consistent with the value of 24 adopted in this work. A possible explanation for the lower values of Γ obtained from the radio observations is that they are calculated based on the apparent velocity of the superluminal features in the jet, which may travel at a lower speed compared to the bulk plasma flow <cit.>. 
                    
                    At large r values (r≳2 pc), the lower limit on the Lorentz factor Γ increases to >100, exceeding the typical range of Γ∼4-50 obtained from observations of blazars <cit.>. This indicates that another seed-photon population, such as an external synchrotron radiation field, is needed if the gamma-ray-emitting region lies beyond ∼2 pc. 
                    
                    
                    
                     §.§ On the radio and optical polarizations
                    
                    
                    The 43-GHz and 15-GHz observations reveal that the EVPAs at the core are mostly parallel to the PA of the jet. This implies that the magnetic field is likely toroidal or strongly helical near the core, consistent with earlier observations of BL Lacertae <cit.>. The 15-GHz EVPAs at larger distances away from the core become more perpendicular to the PA of the jet, indicating that the magnetic field may be more poloidal in the outer jet. Such a magnetic field configuration has been proposed for low-frequency-peaked BL Lac objects <cit.>. 
                    
                    Based on these radio observations, we can use the observed changes in the optical polarizations of BL Lacertae (as shown in Figure <ref>) to gain insights into the magnetic field structure and the location of the region that dominates the optical emission <cit.>. 
                    The optical EVPAs were observed roughly perpendicular to the PA of the jet in late 2016 Sept, indicating that the magnetic field is close to being aligned with the jet and likely dominated by the region downstream in the jet at that time. 
                    Similarly, the optical EVPAs became mostly parallel to the PA of the jet after late Oct, suggesting that the optical emission was then dominated by the core or the inner jet. 
                    We also observed the highest optical fractional polarization during this period, suggesting that 
                    the magnetic field of the core/inner jet is more ordered. 
                    During the three days preceding the TeV gamma-ray flare, the optical EVPA became (temporarily) nearly aligned with the PA of the jet, but on the day before the TeV gamma-ray flare it suddenly rotated back to a direction consistent with its direction prior to this quasi-alignment. 
                    Such abrupt changes in optical polarization associated with flares are found in numerical simulations for blazars and gamma-ray bursts <cit.>, 
                    
                    and can potentially be interpreted as resulting from the helical motion of an emitting component in a toroidal/helical magnetic field before that component reaches the shocked region <cit.>. 
                    However, since the fractional polarization was relatively low during this period, it is also possible that the observed EVPA change was a random fluctuation due to a turbulent magnetic field. 
                    
                    
                    A superluminal radio knot K16 was observed through a series of VLBA exposures on BL Lacertae at 43 GHz. Extrapolation of the knot position implies that the VHE gamma-ray flare happened as the knot K16 crossed the quasi-stationary radio core. 
                    This suggests a possible association between the fast VHE gamma-ray flare and the emergence of the superluminal radio knot for the source, similar to that reported by <cit.>. 
                    
                    
                     §.§ Interpretations of the TeV gamma-ray and the radio results
                    
                    
                    In the model proposed by <cit.>, the radio core is a Mach disk at the apex of a conical shock downstream in the jet, with a transverse orientation with respect to the jet axis. 
                    
                    
                    
                    
                    
                    When turbulent cells of plasma pass through the conical shock, relativistic 
                    electrons can be accelerated to higher energies in those cells where the magnetic 
                    field orientation relative to the shock normal is favorable. A fast gamma-ray flare 
                    can happen via inverse-Compton scattering as the relativistic plasma approaches the
                    Mach disk at the end of the conical shock, which provides a dense source of 
                    synchrotron and SSC seed photons. After the energized plasma passes the Mach disk, a 
                    conical rarefaction causes the flow to expand and accelerate, with the bright plasma 
                    appearing as a superluminal radio knot. 
                    
                    
                    
                    In some numerical simulations, the polarization fraction drops as the magnetic field direction changes, while the EVPA can rotate owing to random fluctuation of the field or the emergence of a new field component <cit.>. 
                    This is consistent with the variation observed in the R-band polarization shortly before the VHE gamma-ray flare (see Figure <ref>), as well as the VLBA images at 15.4 GHz (see Figure <ref>). The changing superposition of the magnetic fields as the moving knot (K16) passes the quasi-stationary knots (A0, A1, and A2) may also explain the change in the positions of A0, A1, and A2 between epochs (see Figure <ref>). 
                    
                    
                    An alternative hypothesis that can explain both the VHE gamma-ray flare and the superluminal radio knot of BL Lacertae 
                    is the breakout of a recollimation-shock zone 
                    <cit.>. 
                    In this model, one or more recollimation shocks, of similar nature to those in <cit.><cit.>, can form upstream in the jet where the magnetic energy density is high and appear as stationary radio knots; further downstream in the jet, particle kinetic energy becomes dominant, the magnetic field becomes unstable, and a stationary knot can be carried away by the underlying relativistic flow and become a superluminal knot. 
                    In the case of a compact region with large kinetic energy passing the recollimation-shock zone, 
                    a multi-component flare could be observed, with one component that varies slowly (i.e., on timescales of hours), thereby giving the appearance of a quasi-constant baseline in an intra-night light curve, as a result of the following sequence of events. 
                    First, in this scenario, an increase in the non-thermal emission of the shock region is expected, which leads to a flux increase on the timescale corresponding to the size of the entire shock region (as the baseline component). 
                    As the kinetic power of the jet increases at the shock zone, 
                    the magnetic field structure is subject to strong tearing instabilities, at which point a magnetic-reconnection event occurs, leading to the observed fast flare. 
                    Finally, the shock zone is dragged away by the flow and enters an adiabatic expansion and cooling phase, leading to a decrease in flux and a return to the low state of the source. 
                    In the case of the 2016 flare of BL Lacertae, there is no evidence for any disruption or breakout of a stationary knot, although it is possible that the recollimation zone reformed quickly between VLBA epochs and was therefore not sampled by the observations. 
                    Therefore, future observations of flares from gamma-ray blazars, with adequate coverage after the flux decreases, can potentially reduce the ambiguity in the interpretation. 
                    
                    
                    
                    
                    
                    VERITAS is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation and the Smithsonian Institution, and by NSERC in Canada. 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 VERITAS Collaboration is grateful to Trevor Weekes for his seminal contributions and leadership in the field of VHE gamma-ray astrophysics, which made this study possible.
                    
                    
                    The research at Boston University was supported in part by NASA Fermi Guest Investigator Program grant 80NSSC17K0694.
                    The VLBA is an instrument of the Long Baseline Observatory (LBO). The LBO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. 
                    
                    This research has made use of data from the MOJAVE database that is maintained by the MOJAVE team <cit.> and supported by NASA- Fermi grant NNX15AU76G. This work made use of the Swinburne University of Technology software correlator <cit.>, developed as part of the Australian Major National Research Facilities Programme and operated under licence. YYK and ABP are partly supported by the Russian Foundation for Basic Research (project 17-02-00197), the government of the Russian Federation (agreement 05.Y09.21.0018), and the Alexander von Humboldt Foundation. T.S. was funded by the Academy of Finland projects 274477 and 284495. 
                    
                    
                    This research has made use of data from the OVRO 40-m monitoring program <cit.> which is supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911.
                    
                    The monitoring of BL Lacertae and other blazars at the Steward Observatory is supported through NASA Fermi Guest Investigator grant NNX15AU81G.
                    
                    IA acknowledges support by a Ramón y Cajal grant of the Ministerio de Economía y Competitividad (MINECO) of Spain. Acquisition and reduction of the MAPCAT data was supported in part by MINECO through grants AYA2010-14844, AYA2013-40825-P, and AYA2016-80889-P, and by the Regional Government of Andalucía through grant P09-FQM-4784. The MAPCAT observations were carried out at the German-Spanish Calar Alto Observatory, which is jointly operated by the Max-Plank-Institut für Astronomie and the Instituto de Astrofísica de Andalucía-CSIC. 
                    
                    The St. Petersburg University team acknowledges support from Russian Science Foundation grant 17-12-01029.
                    
                    
                    VERITAS, Fermi(LAT), Swift(XRT), SO:Kuiper, Bok, CrAO:1.25m, CAO:2.2m, Perkins, LX-200, FLWO:1.2m, VLBA,  Metsähovi Radio Observatory, OVRO:40m
                    
                    
                    
                              Emcee <cit.>, 
                              NumPy <cit.>,
                              Matplotlib <cit.>,
                              SciPy <cit.>,
                              Seaborn <cit.>
                    
                    
                    
                    
                    § THE MODELLING OF THE VHE GAMMA-RAY FLARING LIGHT CURVE
                    
                    
                    
                    
                    ccccc ccccc cc[b!]
                    Comparisons between the fit to the VHE gamma-ray light curve with models I (Equation <ref>) and II (Equation <ref>) 12
                    0ptModel   Fitting method   F_0   t_peak   t_rise   t_decay   F_base   χ^2_red   AIC   BIC   F   2Δ L
                          (10^-6 photon m^-2 s^-1)   (min)   (min)   (min)   (10^-6 photon m^-2 s^-1)               
                    3*I    LM 	   3.4±0.2   133±4   145±18   31±6     -			    1.46     70.3	    76.6	   11.0	   13.0	
                        Soft-L1 	   3.4±0.2   128±6   146±25   40±9     -			    1.47     70.8	    77.2  	   1.7	   2.4		
                        MCMC	   3.4^+0.2_-0.2   130^+5_-3   140^+25_-11   36^+8_-7     -			    1.46     70.2 	    76.6 	   9.3	   11.4	
                    
                    				3*II    LM 	   2.7±0.3   140±1   63±12   3.7±1.5   1.2±0.1    1.18      59.9	    67.5	   11.0 	   13.0	
                        Soft-L1 	   2.6±0.4   131±5   79±29   16±10   1.0±0.4    1.45      70.9	    78.6	   1.7  	   2.4		
                        MCMC	   2.7^+0.3_-0.2   140^+3_-3   65^+19_-8   2.6^+6.7_-0.8   1.2^+0.1_-0.2    1.21     61.3	    68.9	   9.3	   11.4
