Multiwavelength Observations of the Blazar BL Lacertae: A New Fast TeV Gamma-Ray Flare

VERITAS is an array of four imaging atmospheric-Cherenkov telescopes located in southern Arizona (30°40' N, 110°57' W, 1.3 km above sea level; Holder 2011). 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 five 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 October 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 minutes after data-quality selection, with zenith angles ranging between 11° and 30°. The data were analyzed using two independent analysis packages (Cogan 2008; Daniel 2008) and a predetermined set of cuts optimized for lower-energy showers (see e.g., Archambault et al. 2014). 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⁻⁶ photon m⁻² s⁻¹ (or ∼0.95 C.U.).

2.1.1. VHE Gamma-Ray Flux Variability and the Modeling of the Flare Profile

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

Zoom In Zoom Out Reset image size

We first fitted the 4-minute-binned VERITAS light curve with a constant-flux model, obtaining a χ² value of 170.8 for 45 degrees of freedom (DOFs), corresponding to a p-value of 1.1 × 10⁻¹⁶ 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:

$$\begin{eqnarray}F(t)=\left\{\begin{array}{ll}{F}_{0}{e}^{(t-{t}_{\mathrm{peak}})/{t}_{\mathrm{rise}}}, & t\leqslant {t}_{\mathrm{peak}};\\ {F}_{0}{e}^{-(t-{t}_{\mathrm{peak}})/{t}_{\mathrm{decay}}}, & t\gt {t}_{\mathrm{peak}};\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where F₀ 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 emcee (Foreman-Mackey et al. 2013) 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⁵ 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 2. 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 1. The rise and decay times of the flare are determined to be ${140}_{-11}^{+25}\,\mathrm{minutes}$ and ${36}_{-7}^{+8}\,\mathrm{minutes}$, respectively. The best-fit peak time and flux are ${130}_{-3}^{+5}\,\mathrm{minutes}$ (after MJD 57666.165) and ${3.4}_{-0.2}^{+0.2}\times {10}^{-6}\,\mathrm{photon}\,{{\rm{m}}}^{-2}\,{{\rm{s}}}^{-1}$, respectively.

Zoom In Zoom Out Reset image size

Further VERITAS observations of BL Lacertae were made on 2016 October 6 with 37.6 minutes live exposure, and from October 22 to November 19 with 294.6 minutes 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 upper limits on integral flux (shown in the first panel of Figure 5) between 0.2 and 30 TeV at 99% confidence level from the observations on October 6 and between October 22 and November 19 were obtained as 2.0 × 10⁻⁷ photon m⁻² s⁻¹ and 2.8 × 10⁻⁸ photon m⁻² s⁻¹, respectively, assuming a power-law spectrum with a photon index of 3.3 (see Section 2.1.2).

Motivated by the existence of multiple radio emission zones identified in VLBA data (see Section 2.5) and several multi-zone models for BL Lacertae that are consistent with past observations (e.g., Raiteri et al. 2013; Hervet et al. 2016), 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}_{-0.8}^{+6.7}\,\mathrm{minutes}$, with a baseline flux of ${1.2}_{-0.2}^{+0.1}\times {10}^{-6}\,\mathrm{photon}\,{{\rm{m}}}^{-2}\,{{\rm{s}}}^{-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 2.2). However, we would like to highlight that it is not possible to unambiguously reject either model based on the statistics.

2.1.2. 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 χ² value χ²/DOF = 30.6, indicating that a simple power law does not adequately describe the spectrum.

A log-parabola model fits the VERITAS spectrum better:

$$\begin{eqnarray}\displaystyle \frac{{dN}}{{dE}} & = & (2.22\pm 0.07)\times {10}^{-5}\\ & & \times {\left(\displaystyle \frac{E}{0.2\mathrm{TeV}}\right)}^{\left[-(2.4\pm 0.1)-(1.8\pm 0.3){\mathrm{log}}_{10}(\tfrac{E}{0.2\mathrm{TeV}})\right]}\\ & & \times {{\rm{m}}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{TeV}}^{-1},\end{eqnarray} \tag{ 2 }$$

with χ²/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 model of extragalactic background light in Domínguez et al. (2011), the best-fit log-parabola model becomes

$$\begin{eqnarray}\displaystyle \frac{{dN}}{{dE}} & = & (2.36\pm 0.07)\times {10}^{-5}\\ & & \times {\left(\displaystyle \frac{E}{0.2\mathrm{TeV}}\right)}^{\left[-(2.2\pm 0.1)-(1.4\pm 0.3){\mathrm{log}}_{10}(\tfrac{E}{0.2\mathrm{TeV}})\right]}\\ & & \times {{\rm{m}}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{TeV}}^{-1},\end{eqnarray} \tag{ 3 }$$

with χ²/DOF = 1.7. The observed and de-absorbed TeV gamma-ray spectra are shown together with the GeV gamma-ray spectra (Section 2.2) in Figure 3 in the νF_ν representation.

Zoom In Zoom Out Reset image size

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 (Atwood et al. 2009).

An unbinned likelihood analysis was performed with the LAT ScienceTools v10r0p5 and Pass-8 P8R2_SOURCE_V6_v06 instrument response functions (Atwood et al. 2013). SOURCE class events with energy between 100 MeV and 300 GeV within 10° of 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 (gll_iem_v06) and isotropic (iso_P8R2_SOURCE_V6_v06) 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 3FGL 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 5, 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 3. In order to obtain the GeV gamma-ray SEDs, we used the user-contributed tool likeSED.py⁴⁰ 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 3FGL catalog (with an index of 2.25; Acero et al. 2015).

Both the GeV and TeV gamma-ray spectral indices of this flare in 2016 are comparable to those of the flare in 2011 (Arlen et al. 2013), and they suggest that the peak energy of the gamma-ray SED during the flare is between 5 GeV and 100 GeV.

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 (Gehrels et al. 2004; Burrows et al. 2005).

Follow-up observations of BL Lacertae were carried out with the Swift-XRT on 2016 October 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 October 27 and November 2. The XRT data, taken in the photon-counting (PC) mode, were analyzed using the HEAsoft package (v6.19). The data were first processed using xrtpipeline (v0.13.2) with calibration database (CALDB v20160706). The count rates in the PC mode were >0.5 count s⁻¹, and the effect of potential pile-up was checked for all observations by fitting a King function to the point-spread functions at >15 arcsec (Moretti et al. 2005). Those central pixels where the data fall below the model curve, indicating pile-up, were excluded.

For the observations on 2016 October 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. Annular source regions were used for the data taken on 2016 October 7 and 8, 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 October 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 October 7 consisted of two intervals of duration 486 s and 1422 s, separated by roughly one satellite orbital period (∼90 minutes). 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 minutes of this relatively short exposure were checked for data quality, and no anomaly was found.

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

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}={e}^{-{N}_{H}\sigma (E)}K{\left(\displaystyle \frac{E}{1\mathrm{keV}}\right)}^{-\alpha },\end{eqnarray} \tag{ 4 }$$

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 1. Note that the best-fit values of N_H are in agreement with the archival results from X-ray spectral fit (Madejski et al. 1999; Ravasio et al. 2003; Arlen et al. 2013), but are larger than the value N_H = 1.8 × 10²¹ cm⁻² from the Leiden/Argentine/Bonn (LAB) survey of Galactic H i (Kalberla et al. 2005). This difference is likely due to the additional contribution of the Galactic molecular gas (e.g., CO emission), because BL Lacertae is relatively close to the Galactic plane (with a Galactic latitude b = −1044) (Madejski et al. 1999). 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²¹ cm⁻² 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²¹ cm⁻² 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 October 6, 7, and 8 are shown in Figure 4. The X-ray emission from the source was stronger and harder on 2016 October 7 (two days after the TeV gamma-ray flare) compared to the day before and the day after (see Table 1). The energy flux values based on the best-fit absorbed-power-law model between 0.3 keV and 10 keV on 2016 October 6, 7, and 8 were (1.4 ± 0.1), (14.2 ± 0.9), and (1.1 ± 0.1) × 10⁻¹¹ erg cm⁻² s⁻¹, respectively, as shown in the third panel of Figure 5 along with the results from two more observations taken on 2016 October 27 and November 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Swift-XRT Spectral-fit Results Using the Absorbed-power-law Model Described in Equation (4), with N_H Free and Fixed

Date	α	K	NH	χ2/DOF
		(10−2 keV−1 cm−2 s−1)	(1021 cm−2)
Oct 6	2.5 ± 0.1	${0.62}_{-0.06}^{+0.07}$	${2.7}_{-0.3}^{+0.3}$	0.83
Oct 7	2.1 ± 0.1	${4.6}_{-0.5}^{+0.6}$	${3.1}_{-0.4}^{+0.5}$	1.07
Oct 8	2.3 ± 0.1	${0.43}_{-0.05}^{+0.06}$	${3.0}_{-0.4}^{+0.5}$	0.54
Oct 27	${1.4}_{-0.3}^{+0.4}$	${0.14}_{-0.04}^{+0.08}$	${1.2}_{-1.2}^{+2.2}$	0.48
Nov 2	1.3 ± 0.3	${0.11}_{-0.03}^{+0.04}$	${1.4}_{-0.8}^{+1.1}$	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

Note. The errors quoted denote 68% confidence intervals.

Download table as:  ASCIITypeset image

BL Lacertae was monitored in the R-band at high cadence by a number of optical facilities, including the Steward Observatory⁴¹ (Smith et al. 2009), the AZT-8 reflector of the Crimean Astrophysical Observatory, the Perkins telescope,⁴² the LX-200 telescope in St. Petersburg, Russia, and the Calar Alto 2.2 m Telescope (with observations obtained through the MAPCAT⁴³ 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 (Fan et al. 1998) and the transformation r'− R = 0.19(V − R) + 0.13 (Smith et al. 2002). 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 5. 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 5. 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 (see, e.g., Abdo et al. 2010). 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 September to roughly parallel in late 2016 October, 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 October when the EVPA was again aligned with the jet.

BL Lacertae was observed throughout the period of interest at 43 GHz with the VLBA under the VLBA-BU-Blazar monitoring program (Jorstad & Marscher 2016) and at 15.4 GHz with the Monitoring Of Jets in AGNs with VLBA Experiments (MOJAVE) program (Lister et al. 2009). The 43 GHz and 15.4 GHz VLBA data calibration and imaging procedures were identical to those described by Jorstad et al. (2005) and Lister et al. (2009), respectively.

Figure 6 presents 43 GHz images of the parsec-scale jet of BL Lacertae at five epochs from 2016 September 5 to December 23. The second epoch, 2016 October 6, took place only one day after the VHE flare. The images are convolved with a circular Gaussian restoring beam with a 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 October 6 observation was affected by equipment failure at the Maunakea 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.4}^{+0.8}\,\mathrm{pc}$ if we adopt a viewing angle of 42 ± 13 between the jet axis and line of sight (Jorstad et al. 2017).

Zoom In Zoom Out Reset image size

As was the case in previous observations (Jorstad et al. 2005; Arlen et al. 2013; Gómez et al. 2016; Wehrle et al. 2016), 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 a 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 Jorstad et al. (2007) estimated to be low (−16°) between 43 GHz and 300 GHz. It is worth mentioning that Hovatta et al. (2012) measured a much lower (by an order of magnitude) Faraday rotation using 8–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 (Jorstad et al. 2007), because the core at 15 GHz is located further away from the black hole than 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 October 23 to ∼0.28 mas from A0 on December 23. This corresponds to an apparent speed of $6c$, within the range typically observed in BL Lacertae (Jorstad et al. 2005, 2017; Marscher et al. 2008; Arlen et al. 2013; Lister et al. 2013; Wehrle et al. 2016). Extrapolation back to October 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 (Jorstad et al. 2017). This implies that the VHE flare occurred as the moving knot crossed the stationary "core," which Marscher et al. (2008) 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 7, 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 December 26 and gradually increased, with a potentially bright feature with distinct polarization angle (consistent with the EVPAs measured at 43 GHz on 2016 December 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 Bach et al. (2006), 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.

Zoom In Zoom Out Reset image size

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 September. This implies that the magnetic field is toroidal/helical at the core, as we discuss in Section 3. Further downstream in the jet, the EVPAs become more perpendicular to the PA of the jet, as shown in Figure 7. Such location-dependent radio EVPAs help us to interpret the dominant optical component based on the optical EVPA data.

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

Zoom In Zoom Out Reset image size

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 dual-beam receiver of 1 GHz bandwidth centered at 36.8 GHz. The front end, a pseudomorphic transistor of high electron mobility, 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. The typical integration time to obtain one flux-density data point is between 1200 and 1600 s. The detection limit of the telescope at 37 GHz is of 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 H ii 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 Teräsranta et al. (1998). 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 drifts in receiver gain, 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 Baars et al. (1977). The systematic uncertainty in the flux-density scale is ∼5%, which is not included in the error bars in Figure 8. Complete details of the reduction and calibration procedures are given in Richards et al. (2011).

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.

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

$$\begin{eqnarray}&&R\leqslant \displaystyle \frac{{{ct}}_{\mathrm{decay}}\delta }{1+z},\end{eqnarray} \tag{ 5 }$$

where c is the speed of light, δ is the Doppler factor of the jet,⁴⁴ 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⁸M_⊙ by Wu et al. (2009) using the R-band absolute magnitude and the empirical correlation between black hole mass and bulge luminosity of the host galaxy (McLure & Dunlop 2002). The corresponding Schwarzschild radius R_s of the central black hole of BL Lacertae is ∼1.1 × 10¹² m (∼3.6 × 10⁻⁵ pc). It is worth noting that measuring the mass of a black hole is a challenging task, and M_BH values in the range (0.16–5.01) × 10⁸M_⊙ have been reported for BL Lacertae (see Gupta et al. 2012, and references therein).

The Doppler factor δ was estimated to be ∼24 according to the method described in Hervet et al. (2016) from the propagation of a possible perturbation in the radio jet observed with VLBA at 15 GHz (Lister et al. 2013), assuming a viewing angle of 22 based on radio measurements of apparent velocity. Taking the best-fit value of ${t}_{\mathrm{decay}}={36}_{-7}^{+8}\,\mathrm{minutes}$ (see Section 2.1.1) and using Equation (5), we estimate the upper limit on the size of the emitting region to be R ≲ 11.9R_s.

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 (see, e.g., Katarzyński et al. 2003; Petropoulou et al. 2016). 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 8) are consistent with the model proposed by Petropoulou et al. (2016). 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 (analogous to solar flares; see, e.g., Harra et al. 2016).

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 October 7 (two days after the TeV flare), about a factor of 10 stronger than the flux on October 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 (e.g., a radio core; see Kovalev et al. 2009), while the TeV gamma-rays were produced through magnetic reconnection in a localized region (e.g., a small plasmoid in a magnetic reconnection layer, possibly at the interface between a radio core and a moving knot; see Petropoulou et al. 2016).

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 (e.g., Madejski et al. 1999; Raiteri et al. 2013). 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 (e.g., Abdo et al. 2011; Rani et al. 2013). Such behavior is most frequently observed in FSRQs and can be interpreted as an ERC process acting on IR photons in the torus region (e.g., Ghisellini & Tavecchio 2009; Tagliaferri et al. 2015). Such an ERC process was also used to interpret the emission of BL Lacertae in a flaring state (Madejski et al. 1999; Ravasio et al. 2003).

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 Dondi & Ghisellini (1995) (see also Arlen et al. 2013). We found that δ ≳ 13 using the following observables: the best-fit decay time of the TeV gamma-ray flare (36 minutes), 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 (1.5; Allen et al. 1982). Assuming a viewing angle of 22, 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 Nalewajko et al. (2014). 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 an 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 2. Some of the parameters are constrained by observations, and the others 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⁴⁴ erg s⁻¹ (Abdo et al. 2011) and the relations ${r}_{\mathrm{BLR}}=1\times {10}^{15}\sqrt{({L}_{{\rm{d}}}/{10}^{45}\,\mathrm{erg}\ {{\rm{s}}}^{-1})}\,{\rm{m}}\approx 0.025\,\mathrm{pc}$ and ${r}_{\mathrm{IR}}=2\times {10}^{16}\sqrt{({L}_{{\rm{d}}}/{10}^{45}\,\mathrm{erg}\ {{\rm{s}}}^{-1})}\,{\rm{m}}\approx 0.5\,\mathrm{pc}$ (Ghisellini & Tavecchio 2015).

Table 2.  Parameters Used to Calculate the Constraints Shown in Figure 9

Lgamma	Lsyn	Lda	tvar	MBHb	Ecool	δ/Γ	BLR	IR	rBLR	rIR	EBLR	EIR	gSSC	gERC
(erg s−1)	(erg s−1)	(erg s−1)	(min)	(M⊙)	(TeV)				(pc)	(pc)	(eV)	(eV)
7.8 × 1045	7.8 × 1044	6.0 × 1044	36	3.8 × 108	1	1	0.1	0.1	0.025	0.5	10	0.3	0.75	0.5

Notes. L_gamma, L_syn, and L_d are the observed gamma-ray luminosity, the synchrotron luminosity, and the 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.

^aAbdo et al. (2011). ^bWu et al. (2009).

Download table as:  ASCIITypeset image

Zoom In Zoom Out Reset image size

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 (Jorstad et al. 2005, 2017) 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 than the bulk plasma flow (e.g., Lister et al. 2013; Hervet et al. 2016).

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 (e.g., Jorstad et al. 2005; Cohen et al. 2007; Lister et al. 2016). 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.

The 43 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 (e.g., Gómez et al. 2016). 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 (e.g., Kharb et al. 2008; Hervet et al. 2016).

Based on these radio observations, we can use the observed changes in the optical polarizations of BL Lacertae (as shown in Figure 5) to gain insights into the magnetic field structure and the location of the region that dominates the optical emission (e.g., D'arcangelo et al. 2009; Algaba et al. 2011). The optical EVPAs were observed roughly perpendicular to the PA of the jet in late 2016 September, 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 October, 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 (Zhang et al. 2014; Deng et al. 2016), 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 (e.g., Marscher et al. 2008). However, since the fractional polarization was relatively low during this period, it is also possible that the observed change in EVPA 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 Arlen et al. (2013).

In the model proposed by Marscher (2014), 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 (e.g., Marscher 2014; Zhang et al. 2014). This is consistent with the variation observed in the R-band polarization shortly before the VHE gamma-ray flare (see Figure 5), as well as the VLBA images at 15.4 GHz (see Figure 7). 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 6).

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 (Hervet et al. 2016). In this model, one or more recollimation shocks, of similar nature to those in Marscher (2014) (see also Mizuno et al. 2015; Fromm et al. 2016), 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 in 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.