Multiwavelength Variability of BL Lacertae Measured with High Time Resolution

Blazars are a class of active galactic nuclei (AGNs) that possess extreme characteristics across the electromagnetic spectrum (Angel & Stockman 1980). They are the most common extragalactic sources of γ-ray photons with energies ≥0.1 GeV (VHE: E > 100 GeV; e.g., MAGIC Collaboration et al. 2019; Abdollahi et al. 2020). They are thought to be powered by relativistic jets of high-energy plasma flowing away from the central engine at nearly the speed of light (e.g., Lister et al. 2016; Jorstad et al. 2017), with trajectories closely aligned to the line of sight. The observed phenomena include ultraluminous emission (apparent luminosity as high as ∼10⁵⁰ erg s⁻¹; e.g., Abdo et al. 2010b, 2011a; Giommi et al. 2012; Şentürk et al. 2013), high amplitudes of variability on timescales as short as several minutes at various wave bands (e.g., H. E. S. S. Collaboration et al. 2010; Jorstad et al. 2013; Weaver et al. 2019), and high degrees of optical linear polarization (which can exceed 40%; e.g., Smith 2016). Both theoretical work (e.g., Königl 1981; Marscher 1987) and observations (e.g., Hartman et al. 1999; Jorstad et al. 2001; Lister et al. 2011) have found a tight connection between the high-energy and radio emission from the jets.

Blazars are split into two classes: flat-spectrum radio quasars (FSRQs) and BL Lacertae objects (BLs), based on their optical emission-line properties and compact radio morphologies (Wardle et al. 1984; Stickel et al. 1991; Weymann et al. 1991; Urry & Padovani 1995). The spectral energy distributions (SEDs) of both blazar types generally consist of two humps. The first, located between 10¹³ and 10¹⁷ Hz, is attributed to synchrotron radiation by relativistic electrons, and the second, peaking between 1 MeV and 100 GeV, is commonly interpreted as inverse Compton scattering of infrared/optical/UV photons by the same population of electrons (Sikora et al. 2009). BLs are further divided by their synchrotron peaks into high (HBLs), intermediate (IBLs), and low (LBLs) frequency peaking varieties, with ν_peak < 10¹⁴ Hz for LBLs, 10¹⁴–10¹⁵ Hz for IBLs, and >10¹⁵ Hz for HBLs (Padovani & Giommi 1995; Abdo et al. 2010a).

BL Lacertae (hereafter BL Lac; redshift z = 0.069; Miller et al. 1978) is the prototype of BLs. It is usually classified as an LBL (Nilsson et al. 2018), but is sometimes listed as an IBL (Ackermann et al. 2011; Hervet et al. 2016). The blazar has been a target of numerous multiwavelength observing campaigns (e.g., Hagen-Thorn et al. 2002; Agarwal & Gupta 2015; Gaur et al. 2015; Wierzcholska et al. 2015; Abeysekara et al. 2018; Bhatta & Webb 2018; MAGIC Collaboration et al. 2019), including several carried out under the Whole Earth Blazar Telescope (WEBT) GLAST-AGILE Support Program (GASP; e.g., Villata et al. 2002, 2004a, 2004b, 2009b; Böttcher et al. 2003; Bach et al. 2006; Raiteri et al. 2009, 2010).

Abdo et al. (2011b) described a campaign in which BL Lac was in a low, relatively quiescent γ-ray state. The low-level γ-ray emission was explained as inverse Compton (IC) scattering of photons originating outside the jet (external IC radiation, EIC) in addition to IC scattering of in-jet photons (synchrotron self-Compton emission, SSC; e.g., Madejski et al. 1999; Böttcher & Bloom 2000). Based on Fermi, Swift, Submillimeter Array, and WEBT observations from 2008 to 2012, as well as data from other studies, Raiteri et al. (2013) interpreted the variability of emission from BL Lac in terms of changes in the orientation of the emitting regions, possibly caused by a shock oriented perpendicular to the jet axis.

More recently, Wehrle et al. (2016) extended the period analyzed by Raiteri et al. (2013) by one year to include an extended interval of erratic changes in γ-ray flux. Their study filled gaps in the SED with Herschel far-infrared and Nuclear Spectroscopic Telescope Array (NuSTAR) hard X-ray observations. They described the flaring nature of the source in terms of turbulent plasma flowing across quasi-stationary shocks within 5 pc of the supermassive black hole, with high-energy electrons accelerated at the shock fronts.

A number of studies have been performed to analyze the optical polarization behavior of BL Lac (e.g., Hagen-Thorn et al. 2002; Sakimoto et al. 2013). A two-component system was proposed to describe changes of polarization parameters with flux and time: a long-lived underlying source of polarized radiation (perhaps variable on timescales of years) plus several short-lived components (associated with flares) with randomly oriented polarization directions and high degrees of polarization. Blinov & Hagen-Thorn (2009) have employed a Monte Carlo method to simulate simultaneous photometric and polarimetric data of BL Lac over a period of 22 yr within such a system. These authors have found that the observed photometric and polarimetric variability of BL Lac can be explained within a model containing a steady component with a high degree of polarization, ∼40%, and a position angle of polarization along the parsec-scale jet direction, plus 10 ± 5 components with variable polarization.

Exploration of the complex emission mechanisms and physical processes that operate in blazar jets requires observations of variability on both long and short timescales. For example, time-series studies of long-term light curves show that for many blazars there is a correlation between γ-ray and optical variations, with time delays ranging from zero to several days (e.g., Chatterjee et al. 2012; Jorstad et al. 2013; Raiteri et al. 2013; D'Ammando et al. 2019). However, the uncertainties are often comparable with the delays themselves. This limitation can be overcome with short-cadence observations in order to increase the precision of correlation analyses and search for patterns and characteristic timescales of variations at different wavelengths (e.g., Uttley et al. 2002).

Until recently, intensive monitoring campaigns to identify short-timescale variability have been limited to relatively brief time spans, usually with gaps in temporal coverage. The Kepler mission has provided optical short-cadence light curves over spans of weeks and without gaps for several AGNs (e.g., Mushotzky et al. 2011; Edelson et al. 2014; Smith et al. 2015, 2018; Aranzana et al. 2018). The resulting data allow for precise time-series analyses whose value can be amplified by the addition of simultaneous, well-sampled light curves at other wavelengths. The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) is capable of producing short-cadence, unbiased light curves of many more AGNs as it performs a nearly all-sky survey. It samples fluxes over a wide optical to near-IR band with a default cadence of 30 minutes, shortened to 2 minutes for selected objects.

In this paper, we report the results of an observing campaign that combines continuous monitoring of BL Lac with TESS along with broad multiwavelength coverage from other space- and ground-based facilities. These include TESS, the Fermi Large Area Telescope (LAT), NuSTAR, the Neil Gehrels Swift satellite, and ground-based telescopes within the WEBT collaboration. The paper is structured as follows. In Section 2 we describe the reduction of observations at all wavelengths. The light curves are presented in Section 3, first for the entire WEBT 3 months campaign and then for overlapping observations with the other telescopes. We investigate short-term variability observed with TESS at optical and NuSTAR at X-ray frequencies, as well as R-band optical polarization, in Section 4. We analyze the optical behavior in different bands in Section 5. The polarization behavior over the entire time span, observed with a subset of WEBT telescopes throughout the campaign, is presented in Section 6. In Section 7 we perform a correlation analysis between the TESS and other light curves to identify potential time lags. We review and discuss our results and offer theoretical interpretations in Section 8, and summarize our conclusions in Section 9.

In order to investigate the short-timescale variability of BL Lac, we organized a multiwavelength observing campaign around the TESS observations, which took place as part of observing sector 16 at a 2 minute cadence, from 2019 September 12 to October 6 (MJD: 58738–58762). We have retrieved ∼3 months of γ-ray data measured by Fermi-LAT, and obtained four-band, BVRI optical flux and R-band polarization data with numerous WEBT-affiliated ground-based telescopes, over the time period from 2019 August 5 to November 2 (MJD: 58700–58789). We also obtained observations at X-ray frequencies with the NuSTAR satellite and at X-ray, UV, and optical frequencies with the Swift satellite for 5 days during the TESS monitoring, 2019 September 14–19 (MJD: 58740–58745). In this section we discuss the processing of the various data.

2.1. Gamma-Ray Data

Fermi-LAT (Atwood et al. 2009) surveys the entire sky every ∼3 hr in the energy range 0.1–300 GeV, with data archived for public access. We retrieved P8R3 photon and spacecraft data centered on BL Lac (4FGL J2202.7+4216). The data were reduced using version v1.0.10 of the Fermi Science Tools, background models from the iso_P8R3_SOURCE_V2_v1.txt isotropic template, and the gll_iem_v07 Galactic diffuse emission model.²⁹ We utilized analysis cuts of evclass = 128, evtype = 3, and zmax = 90 for an unbinned likelihood analysis of the photon data, which we restricted to an energy range of 0.1–200 GeV.

The γ-ray emission from BL Lac and other point sources within a 25° radius of BL Lac were represented by spectral models listed in the 4FGL catalog of sources detected by LAT (The Fermi-LAT collaboration 2019). Specifically, the number of photons N per unit energy E of BL Lac was modeled as a log-parabola of the form

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}={N}_{0}{\left(\displaystyle \frac{E}{{E}_{{\rm{b}}}}\right)}^{-(\alpha +\beta \mathrm{log}E/{E}_{{\rm{b}}})}.\end{eqnarray} \tag{ 1 }$$

During the analysis, the spectral parameters of BL Lac were kept fixed at their 4FGL catalog values: α = 2.1755, β = 6.0062 × 10⁻², and break energy E_b = 7.47961 × 10² MeV. This was necessary because of the relatively low flux, which precludes an accurate spectral analysis. The prefactor N₀ was allowed to vary for BL Lac, as well as for all cataloged sources within 5° and bright (F_γ > 10⁻¹¹ erg cm⁻² s⁻¹) sources within 10°.

In the initial analysis, we integrated the observations over 6 hr. Contiguous periods of upper limits for the 6 hr binning were recalculated with 12 hr binning. This procedure yielded a γ-ray light curve with 322 measurements of BL Lac over the time period from 2019 August 5 to November 2. The source was considered detected if the test statistic (TS) provided by the maximum-likelihood analysis exceeded 10, which corresponds to approximately a 3σ detection level (Nolan et al. 2012). If TS < 10, we calculated 2σ upper limits with the Fermi Python script. Of the 322 measurements, 126 were detections and 196 were upper limits.

2.2. X-Ray Data

2.2.1. Swift X-Ray Data

The Neil Gehrels Swift Observatory (Swift) X-ray Telescope (XRT; Burrows et al. 2005) observes over the 0.3–10 keV band. We obtained 40 observations over 5 days from 2019 September 14 to 19, averaging one observation every three hours, for a total exposure time of ∼46 ks. All observations were made in photon counting mode.

We used v6.26.1 of the HEAsoft package and CALDB v20190412 to process the data. We defined a circular source region with a 70'' radius and an annular background region with inner radius 88'' and outer radius 118'' (selected to avoid contaminating sources), both centered on BL Lac. Using the standard reduction protocol, we first cleaned the data and created an exposure map with the xrtpipeline tool. The image and spectra were extracted using XSELECT, and the ancillary response file was generated with xrtmkarf. Because we used Cash statistics (Cash 1979; Humphrey et al. 2009) in the form of the modified C-statistic cstat to fit the data in XSPEC, we grouped our data by single photons in grppha. The spectrum was then fit in XSPEC and further evaluated using a χ² test.

The total hydrogen column density toward BL Lac consists of the atomic hydrogen column density, N_{H i} = 1.74 × 10²¹ cm⁻² (Kalberla et al. 2005) and molecular column density, N_H,mol (in atoms cm⁻²). The value of the latter ranges among various studies of the X-ray spectrum from 0.5 × 10²¹ cm⁻² (considered as a possibility by Madejski et al. 1999) to 1.7 × 10²¹ cm⁻² (Raiteri et al. 2009). An even higher value of N_H,mol = 2.8 × 10²¹ cm⁻², based on a CO emission line from a molecular cloud in the direction of BL Lac (Bania et al. 1991), has been proposed as well (see Raiteri et al. 2009, for more examples). In addition, Moore & Marscher (1995) have found changes by 14% in the equivalent width of H₂CO absorption lines along the line of sight to BL Lac on a timescale of ∼2 yr, which suggests that N_H,mol could be variable.

In order to estimate the value of N_H that best fits our data, we have combined our XRT observations over five sets of 24 hr each. These sets were obtained by summing the individual exposure maps using XIMAGE, followed by application of XSELECT to sum the individual event files. The ancillary response file for each set of observations was generated using xrtmkarf with the corresponding summed exposure map.

We modeled the daily combined data at 0.3–10 keV in XSPEC using an absorbed simple power law with all parameters free. The results are presented in Table 1, which yields $\langle {N}_{{\rm{H}}}\rangle =(2.21\pm 0.29)\times {10}^{21}$ cm⁻² over five days. We subsequently modeled each set with a power law for two different values of N_H (fixed during each model fit), 2.7 × 10²¹ cm⁻² and 3.4 × 10²¹ cm⁻². These are the closest values to those estimated from atomic and molecular line observations (see above), and correspond to those used for BL Lac by Madejski et al. (1999), Raiteri et al. (2009), and Wehrle et al. (2016). The results of the modeling are given in sections 2 and 3 of Table 1, from which we find that there is no statistically significant difference between the models as judged by the reduced χ². We then averaged the photon indices obtained over the five days and two fixed values of N_H, which resulted in Γ = 2.419. Fixing Γ at this value, we performed a search for the best-fit value of N_H. The results of this search are listed in section 4 of Table 1. The reduced χ² values are similar to those of the previous three model fits. The last procedure results in an average value of N_H = (2.80 ± 0.32) × 10²¹ cm⁻² over five days, which is in good agreement with the value adopted by Madejski et al. (1999). Based on these considerations, we have modeled the X-ray data presented below using the same fixed hydrogen column density as adopted by Madejski et al. (1999), N_H = 2.7 × 10²¹ cm⁻².

Table 1.  Summary of Swift 0.3–10 keV Modeling to Calculate N_H

MJD start	Exposure time	NH	Γ	Flux	D.o.F.	${\chi }_{\nu }^{2}$
	(s)	×1021 cm−2		×10−12 erg cm−2 s−1
58740.36	8910	${2.60}_{-0.27}^{+0.28}$	${2.474}_{-0.105}^{+0.108}$	${11.90}_{-0.62}^{+0.58}$	481	1.372
58741.28	9752	${2.31}_{-0.32}^{+0.34}$	${2.334}_{-0.127}^{+0.133}$	${7.54}_{-0.43}^{+0.37}$	454	1.063
58742.35	9382	${2.19}_{-0.48}^{+0.51}$	${1.943}_{-0.162}^{+0.171}$	${4.88}_{-0.40}^{+0.50}$	401	1.253
58743.27	9652	${1.79}_{-0.41}^{+0.44}$	${1.930}_{-0.145}^{+0.151}$	${5.64}_{-0.55}^{+0.45}$	406	1.002
58744.27	8296	${2.16}_{-0.33}^{+0.35}$	${2.194}_{-0.126}^{+0.132}$	${8.73}_{-0.69}^{+0.59}$	450	1.045
58740.36	8910	2.7	${2.504}_{-0.057}^{+0.056}$	${11.80}_{-0.50}^{+0.53}$	482	1.355
58741.28	9752	2.7	${2.462}_{-0.071}^{+0.072}$	${7.27}_{-0.37}^{+0.40}$	455	1.022
58742.35	9382	2.7	${2.081}_{-0.100}^{+0.101}$	${4.67}_{-0.36}^{+0.45}$	402	1.253
58743.27	9652	2.7	${2.181}_{-0.092}^{+0.093}$	${5.22}_{-0.35}^{+0.41}$	407	0.994
58744.27	8296	2.7	${2.362}_{-0.074}^{+0.075}$	${8.31}_{-0.38}^{+0.63}$	451	1.010
58740.36	8910	3.4	${2.730}_{-0.061}^{+0.062}$	${11.20}_{-0.43}^{+0.35}$	482	1.302
58741.28	9752	3.4	${2.687}_{-0.077}^{+0.078}$	${6.85}_{-0.32}^{+0.30}$	455	1.009
58742.35	9382	3.4	${2.254}_{-0.110}^{+0.111}$	${4.43}_{-0.31}^{+0.26}$	402	1.274
58743.27	9652	3.4	${2.360}_{-0.099}^{+0.101}$	${4.96}_{-0.35}^{+0.33}$	407	1.026
58744.27	8296	3.4	${2.567}_{-0.082}^{+0.083}$	${7.85}_{-0.39}^{+0.35}$	451	1.017
58740.36	8910	${2.49}_{-0.14}^{+0.15}$	2.419	${12.10}_{-0.36}^{+0.40}$	482	1.411
58741.28	9752	${2.49}_{-0.18}^{+0.18}$	2.419	${7.32}_{-0.24}^{+0.28}$	455	1.033
58742.35	9382	${3.37}_{-0.32}^{+0.36}$	2.419	${4.04}_{-0.21}^{+0.27}$	402	1.265
58743.27	9652	${2.97}_{-0.28}^{+0.30}$	2.419	${4.66}_{-0.22}^{+0.24}$	407	0.977
58744.27	8296	${2.67}_{-0.19}^{+0.20}$	2.419	${8.05}_{-0.41}^{+0.26}$	451	0.985

Note. In section 1 of this table, all parameters were allowed to vary. In sections 2 and 3, N_H was fixed to the listed values. In section 4, the photon index Γ was fixed at the average value of sections 2 and 3, while N_H was allowed to vary.

Download table as:  ASCIITypeset image

We have modeled the 40 Swift XRT observations at 0.3–10 keV with an absorbed simple power law. A summary of the results is provided in Table 2. Figure 1 reveals that the photon index becomes steeper at higher flux levels. The dependence is more apparent for the daily averaged X-ray data.

Zoom In Zoom Out Reset image size

Table 2.  Summary of Swift XRT 0.3–10 keV Observations

Statistic	Value
Number of observations	40
Total exposure time (s)	45,993
Averages per observation:
Count rate (cts s−1)	0.20
Counts per Obs.	233
Photon index	2.33
Minimum photon index	1.79
Maximum photon index	2.72
Flux (erg cm−2 s−1)	7.12 × 10−12
Min. flux (erg cm−2 s−1)	3.06 × 10−12
Max. flux (erg cm−2 s−1)	1.83 × 10−11

Download table as:  ASCIITypeset image

2.2.2. NuSTAR Data

NuSTAR observes in the 3–79 keV energy band (Harrison et al. 2013). Two independent, co-aligned telescopes (FPMA and FPMB) observe as photon counting modules, with each module consisting of a 2 × 2 array of four detectors. Observations of a source can be obtained throughout the satellite's ∼95 minute orbit, excluding dead time while the observatory is slewing, performing calibration or house-keeping activities, passing through the South Atlantic Anomaly (SAA), or occulted by the Earth. We obtained five days of continuous observations (ID 60501024002) from 2019 September 14 05:36:09 to September 19 06:01:09 UT, for a total dead-time-corrected exposure time of ∼197 ks spanning 75 orbits.

We processed the data using the NuSTAR Data Analysis Software (NuSTARDAS), downloaded as part of v6.25 of the HEAsoft package, with CALDB version 20190627. Upon examination of the provided SAA filtering reports, we chose to use the "strict" SAA calculation mode and SAA passage algorithm 1, with no tentacle correction. The nupipeline was run with the SCIENCE observing mode and required the creation of an exposure map. The results were processed through nuproducts. Using SAOImageDS9,³⁰ we defined a 70'' region for each telescope centered on the source and a 70'' background region on the same detector as the source (but sufficiently distant to avoid contamination). As for the Swift data, we used Cash statistics to fit the data in XSPEC, and so grouped our data by single photons.

We divided the observations first by orbit, defined to begin with the satellite's emergence from the SAA as indicated on the Good Time Intervals (GTI) file. We used XSELECT to generate the GTI that subsequently was fed into the nuproducts process. The data were then examined in XSPEC and modeled as an absorbed simple power law. We simultaneously fit the FPMA and FPMB files with a cross-normalization factor frozen to unity for FPMA and allowed to vary for FPMB. The fit was evaluated with the χ² test statistic.

The NuSTAR observations are summarized in Table 3. The average dead-time-corrected exposure time per orbit was 2626 s. Figure 2(a) shows the light curve for the full energy band, as well as two narrower, soft and hard, bands (3–10 and 10–79 keV, respectively). The full spectrum is dominated by flux from the lower energies, while the counts for the higher energies are too low to be analyzed by single orbits. Figure 2(b) shows the photon index for the full spectrum. There are enough counts in the soft band to allow the photon index to vary; this is not the case for the hard band.

Zoom In Zoom Out Reset image size

Table 3.  Summary of NuSTAR 3–79 keV Observations

Statistic	Value
Number of orbits	75
Total dead-time-corrected exoposure (s)	196,938
Averages per orbit:
FPMA count rate (cts s−1)	0.14
FPMA counts per orbit	366
FPMB count rate (cts s−1)	0.13
FPMB counts per orbit	332
Photon index	1.872
Minimum photon index	1.563
Maximum photon index	2.147
Flux (erg cm−2 s−1)	1.371 × 10−11
Min. flux (erg cm−2 s−1)	0.862 × 10−11
Max. flux (erg cm−2 s−1)	1.972 × 10−11

Download table as:  ASCIITypeset image

To analyze the hard energy band (10–79 keV), we divided the observations into groups of both five and 15 orbits. Figure 3 shows the counts per unit energy of an average group of five orbits (orbits 16–20). For clarity, only data from FPMA is shown; data from FPMB are similar. The source counts are significantly above the background level for the lower-energy portion of the spectrum, but steadily deteriorate toward harder energies. Based on the results of all groups of orbits, we have separated the hard spectrum into two bands, 10–22 and 22–79 keV. Figure 4(a) presents the 10–22 keV light curves, with data grouped over five and 15 orbits, while Figure 4(b) plots the photon index at 10–22 keV with data grouped over five orbits. For the 22–79 keV band grouped over five orbits, the counts are so few that it is necessary to fix the photon index at the value determined by the full 10–79 keV range, rather than allow the photon index to vary. Using this method permits us to use the five-orbit binned light curve to enhance the time resolution of the hard X-ray data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.2.3. Simultaneous Swift and NuSTAR Data Reduction

The photon indices for the Swift and NuSTAR observations provided in Tables 2 and 3, respectively, suggest that, in general, the photon index of the 0.3–10 keV emission is steeper than that of the 3–79 keV emission. This implies a break in the X-ray spectrum. We simultaneously fit the NuSTAR and Swift XRT data that are contemporaneous within a given day. This resulted in 15 FPMA/FPMB and eight XRT observations per 24 hr period, with five such daily sets over our observations.

We used XSPEC to simultaneously fit the FPMA, FPMB, and XRT data sets for each day. We employed two models: a single power law and a broken power law, each with photoelectric absorption, and fit the energy range from 0.3 to 79 keV. The results are given in Tables 4 and 5, respectively. There are no statistically significant differences in the reduced χ² values between the single and broken power-law models. Figure 5 gives an example of modeling the data by a broken power-law model. The results presented in Table 5 show that the spectral indices both before and after the break do not exhibit a dependence on flux and that, overall, Γ₁ = 2.40 ± 0.14 and Γ₂ = 1.72 ± 0.05 over the 5 days of observation. However, the break energy tends to increase with flux, with the best-fit models suggesting a break at the highest flux level, ∼6 keV, while E_b ∼ 2 keV at lower flux levels.

Zoom In Zoom Out Reset image size

Table 4.  Simultaneous Single Power-law Fits to 24 hr NuSTAR and Swift XRT Spectra

MJD starta	Exposurea	Exposureb	Γ	Fluxc	Count rate	d.o.f.	${\chi }_{\nu }^{2}$
	(ks)	(ks)			(counts s−1)
58740.24	39.6	8.9	${2.151}_{-0.021}^{+0.021}$	${18.02}_{-0.42}^{+0.52}$	0.678 ± 0.007	2387	1.110
58741.25	39.4	9.8	${1.973}_{-0.025}^{+0.024}$	${15.42}_{-0.40}^{+0.50}$	0.477 ± 0.006	2246	1.044
58742.26	39.3	9.4	${1.728}_{-0.029}^{+0.029}$	${14.24}_{-0.41}^{+0.63}$	0.327 ± 0.004	2191	1.056
58743.26	39.2	9.7	${1.804}_{-0.028}^{+0.028}$	${14.06}_{-0.45}^{+0.54}$	0.349 ± 0.005	2188	0.977
58744.27	39.4	8.3	${1.949}_{-0.024}^{+0.024}$	${17.33}_{-0.49}^{+0.55}$	0.526 ± 0.006	2305	0.978

Notes. Single power law: S(E) = KE^−Γ, where S is photon flux density in photons cm⁻² s⁻¹ keV⁻¹ and E is in keV.

^aOf the NuSTAR data. ^bOf the Swift data. ^cIn units of 10⁻¹² erg cm⁻² s⁻¹.

Download table as:  ASCIITypeset image

Table 5.  Simultaneous Broken Power-law Fits to 24 hr NuSTAR and Swift XRT Spectra

MJD starta	Γ1	Eb	Γ2	Fluxb	Count rate	d.o.f.	${\chi }_{\nu }^{2}$
		(keV)			(counts s−1)
58740.24	${2.363}_{-0.033}^{+0.038}$	${5.864}_{-0.678}^{+0.551}$	${1.716}_{-0.062}^{+0.067}$	${23.04}_{-1.10}^{+1.17}$	0.175 ± 0.002	2385	1.089
58741.25	${2.558}_{-0.126}^{+0.106}$	${2.430}_{-0.264}^{+0.475}$	${1.763}_{-0.036}^{+0.034}$	${17.43}_{-0.51}^{+0.76}$	0.138 ± 0.002	2244	1.011
58742.26	${2.146}_{-0.202}^{+0.771}$	${2.090}_{-1.053}^{+2.200}$	${1.651}_{-0.039}^{+0.061}$	${15.01}_{-0.45}^{+1.49}$	0.110 ± 0.002	2189	1.070
58743.26	${2.462}_{-0.188}^{+0.202}$	${1.770}_{-0.297}^{+0.392}$	${1.706}_{-0.036}^{+0.035}$	${14.96}_{-0.47}^{+0.96}$	0.117 ± 0.002	2186	0.976
58744.27	${2.449}_{-0.131}^{+0.167}$	${2.516}_{-0.483}^{+0.641}$	${1.787}_{-0.034}^{+0.038}$	${18.98}_{-0.65}^{+0.71}$	0.154 ± 0.002	2303	0.965

Notes. Broken power law: $S(E)={{KE}}^{-{{\rm{\Gamma }}}_{1}}$ if E ≤ E_b and $S(E)={{KE}}_{{\rm{b}}}^{({{\rm{\Gamma }}}_{2}-{{\rm{\Gamma }}}_{1})}{E}^{-{{\rm{\Gamma }}}_{2}}$ if E > E_b, with E in keV. Exposure times are the same as in Table 4.

^aOf the NuSTAR data. ^bIn units of 10⁻¹² erg cm⁻² s⁻¹.

Download table as:  ASCIITypeset image

2.3. UV and Optical Data

2.3.1. Swift Ultraviolet and Optical Data

The Swift satellite also provides UV and optical data via the UV/Optical Telescope (UVOT; Roming et al. 2005). We retrieved the data from the HEASARC Archive and reduced them with v6.26.1 of the HEAsoft software and CALDB v20170922. We defined a 5'' circular region centered on the source and a 20'' circular aperture on a source-free region of the image to represent the background. We ran the tool uvotunicorr if an aspect correction was not applied. Multiple extensions within a particular file were summed using uvotimsum, then processed with uvotsource, setting the detection significance to σ = 5. Images were retained if they had an exposure time ≥40 s and a magnitude error σ_mag < 0.2. None of the observations suffered a high coincidence loss. We used the count-rate-to-flux conversion factors reported by Breeveld et al. (2011) for γ-ray burst models, which correspond to continuum spectra similar to those of blazars. We have corrected for Galactic extinction using the values found by Raiteri et al. (2009). Our aperture size introduced a flux contamination from the host galaxy of ∼0.5 times the host-galaxy flux density (Raiteri et al. 2010), which we have subtracted from the source flux density. All Galactic extinction, magnitude-to-flux density conversion factors, and host-galaxy flux densities are given in Table 6.

Table 6.  UV and Optical Correction Factors Used in This Work

Filter	Extinction	Absolute flux	Host-galaxy
	(mag)	densitya	flux density
		$\left({10}^{-20}\right.$ erg	(mJy)
		cm−2 s−1 Hz−1)
UVW2	2.92	0.738	0.017
UVM2	3.04	0.689	0.020
UVW1	2.40	0.942	0.026
u	1.79	1.307	0.036
b	1.44	3.476	1.30
v	1.10	3.420	2.89
B	1.42	4.063	1.30
V	1.08	3.636	2.89
R	0.90	3.064	4.23
I	0.64	2.416	5.90

Note.

^aFor a zero-magnitude star.

References. Bessell et al. (1998), Raiteri et al. (2010), Wehrle et al. (2016).

Download table as:  ASCIITypeset image

2.3.2. WEBT Optical Data

WEBT was formed in 1997 as a network of optical, near-infrared, and radio observatories working together to obtain continuous well-sampled monitoring of the flux and polarization of blazars. In 2007 the WEBT started the GLAST-AGILE Support Program (GASP; e.g., Villata et al. 2008, 2009a). The GASP-WEBT data reported here correspond to four-band optical photometry (BVRI) and R-band polarimetry measured from 2019 August 5 to November 2. The data were checked for consistency between different observers and telescopes (following the standard WEBT prescription; e.g., Villata et al. 2002). Table 7 lists the observatories that participated in the campaign, while Table 8 gives the magnitudes of comparison stars used in the photometric analysis. The data were corrected for Galactic extinction and contamination from the host galaxy, assuming contamination of ∼60% of the total host flux density, as suggested by Raiteri et al. (2010) for a circular aperture with a radius of 8'' employed for the photometry of BL Lac. Galactic extinction along the line of sight to BL Lac was calculated according to Cardelli et al. (1989), with R_V = 3.1 and A_B = 1.42 from Schlegel et al. (1998).³¹ Table 6 gives the extinction, absolute flux density conversion coefficient, and host-galaxy total flux density for each filter.

Table 7.  WEBT-affiliated Ground-based Telescopes Used in This Work

Observatory	Bands	Number of	Marker
		observationsa	styleb
ARIES	BVRI	2, 2, 2, 2	blue $\bullet$
Abastumani	R	144	green $\bullet$
Belogradchik	VRI	14, 16, 15	red $\bullet$
Burke-Gaffney	R	1	cyan $\bullet$
Crimean (AZT-8; AP7p)	BVRI	61, 60, 30, 63	magenta $\bullet$
Crimeanc (AZT-8; ST-7)	BVRI	32, 31, 115 (55), 34	orange $\bullet$
Foggy Bottom	R	249	blue $\unicode{x02666}$
Las Cumbres	R	34	red $\blacktriangle$
Lulin	R	45	black $\bullet$
Mt. Maidanak	BVRI	133, 135, 136, 136	blue $\unicode{x025A0}$
Osaka Kyoiku	R	19	green $\unicode{x025A0}$
Perkinsc	BVRI	112, 116, 193 (193), 110	red $\unicode{x025A0}$
Rozhen (200 cm; 50/70 cm)	BVRI	7, 8, 23, 8	cyan $\unicode{x025A0}$
San Pedro Martirc	R	14 (14)	magenta $\unicode{x025A0}$
Sirio	R	2	orange $\unicode{x025A0}$
Skinakas	BVRI	124, 123, 124, 124	black $\unicode{x025A0}$
Skinakasc (Robopol)	R	5 (5)	blue $\blacktriangle$
St. Petersburgc (LX-200)	BVRI	15, 17, 48 (37), 47	green $\blacktriangle$
Tijarafe	R	219	cyan $\blacktriangle$
Vidojevica (140 cm; 60 cm)	BVRI	3, 3, 3, 3	magenta $\blacktriangle$
West Mountain	V	13	orange $\blacktriangle$

Notes.

^aListed for each filter. Number in parentheses refers to polarimetry measurements for that filter. ^bFor use in Figures 8–11. ^cPhotometry and polarimetry.

Download table as:  ASCIITypeset image

Table 8.  Magnitudes and Distances of Primary Comparison Stars in the Field of BL Lac

Star	Ba	Vb	Rb	Ib	ϖ	Distance
					(mas)	(pc)
B	14.68 ± 0.04	12.90 ± 0.04	11.99 ± 0.04	11.12 ± 0.05	0.3393 ± 0.0266	${2980}_{-343}^{+446}$
C	15.20 ± 0.05	14.26 ± 0.06	13.79 ± 0.05	13.32 ± 0.05	1.8078 ± 0.0213	${549}_{-8}^{+3}$
H	15.81 ± 0.06	14.40 ± 0.06	13.73 ± 0.06	13.07 ± 0.06	0.6922 ± 0.0928	${1453}_{-97}^{+111}$
K	16.36 ± 0.07	15.47 ± 0.07	15.00 ± 0.07	14.54 ± 0.07	0.8994 ± 0.0188	${1113}_{-37}^{+40}$

Notes.

^aFrom Bertaud et al. (1969). ^bFrom Fiorucci & Tosti (1996).

Download table as:  ASCIITypeset image

As in Raiteri et al. (2010), a comparison between the Swift UVOT b and v data and WEBT B and V data revealed an offset between the space-based and ground-based magnitudes. We used the offset determined by Raiteri et al. (2010), with B−b = 0.10, and V − v = −0.05.

During the campaign, the WEBT collaboration observed BL Lac 459 times in the B band, 492 times in the V band, 1417 times in the R band, and 507 times in the I band. During the same time period, a subset of the WEBT observatories measured the R-band polarization a total of 303 times.

2.3.3. Optical Polarization Observations

The R-band polarization observations were obtained at the five telescopes noted in Table 7. The Perkins telescope is equipped with the PRISM camera, which includes a polarimeter with a rotating half-wave plate. Each polarization observation consisted of four consecutive measurements at instrumental position angles 0°, 90°, 45°, and 135° of the wave plate to calculate the normalized Stokes parameters q and u. (For more detail see Jorstad et al. 2010.) Polarization observations at the LX-200 and AZT-8 telescopes were performed in the same manner, each using an identical photometer-polarimeter, with two Savart plates rotated by 45° relative to each other. Swapping the plates allows one to obtain a normalized Stokes parameter, either q or u (for more detail see Larionov et al. 2008). Several polarization observations were also performed at the San Pedro Martir Observatory and Skinakas Observatory (Robopol program). Details of these observations can be found in López & Hiriart (2011) and Ramaprakash et al. (2019), respectively. The degree, P_R, and position angle, χ_R, of the polarization in all cases are calculated from normalized Stokes q and u parameters. Throughout the paper we indicate the degree of polarization in percent. All polarization data have been corrected for the Rice statistical bias (Vinokur 1965), according to Wardle & Kronberg (1974), using the Modified Asymptotic estimator (MAS; Plaszczynski et al. 2014) in the case of the San Pedro Martir data. The instrumental polarization of each instrument has been estimated to be ≲0.5%, based on measurements of unpolarized calibration stars (e.g., Schmidt et al. 1992). We have calculated that the average uncertainty of a measurement of P_R is $\langle {\sigma }_{{P}_{{\rm{R}}}}\rangle \,=0.23 \%$.

As described above, a molecular cloud lies along the line of sight to BL Lac, accounting for a substantial fraction of the total hydrogen column density (e.g., Bania et al. 1991; Madejski et al. 1999). This molecular cloud could possibly contaminate the measured polarization due to dichroic absorption by aligned dust particles along the line of sight.

The data reduction methods applied to all of the polarization data reported here use field stars near the position of BL Lac to perform both interstellar and instrumental polarization corrections. If the field stars lie beyond the molecular cloud, then their polarization would include the dichroic absorption effects from the cloud, and the effects of the molecular cloud will be subtracted out from the polarization of the source. The distance to the cloud has been estimated to be ∼330 pc based on the Galactic latitude and average distance of molecular clouds in the solar neighborhood above and below the Galactic plane (Lucas & Liszt 1993; Moore & Marscher 1995).

We have obtained the Gaia DR2 parallaxes for the four main comparison stars listed in Table 8 (Gaia Collaboration et al. 2016, 2018). However, simple inversion of the parallax introduces known biases to the calculated distance, especially when the relative uncertainties are large. A proper calculation of distance requires a proper statistical treatment of the data. With the pyrallaxes program³² (Luri et al. 2018), we have used Bayesian inference to calculate the distances to each of the comparison stars. We have implemented the two recommended priors from Bailer-Jones (2015): a uniform distance prior (out to 100 kpc) and an exponentially decreasing space density prior (with a characteristic length scale L = 1.35 kpc). Both priors generate extremely similar distance measurements to the stars. We list the parallaxes and distances to the stars, with 90% uncertainty intervals, in Table 8. These distances are calculated with the exponentially decreasing space density prior.

All comparison stars lie beyond the molecular cloud, with the closest one still ∼200 pc beyond the cloud. If we assume that the emission from the stars is intrinsically unpolarized, then the measured polarization to the stars includes the dichroic absorption effects from the cloud. Using 98 of the polarization measurements from the Perkins telescope taken during the monitoring period (during good weather), we estimate the average ISM polarization parameters toward BL Lac to be P_R,ISM = 0.43% ± 0.08% and χ_R,ISM = 69° ± 5°. This level of polarization is well below the upper limit of the contribution of the ambient interstellar dust to polarization along the line of sight through the Milky Way, P_isp ≤ 9% ∗ E(B − V) (Serkowski et al. 1975).

In addition, contamination from the host galaxy can modify the observed degree of optical linear polarization. We assume that the emission from the host galaxy, dominated by starlight, is unpolarized. If the observed degree of linear polarization is P_obs, then the intrinsic degree of polarization, P_AGN, is found through a modulating factor

$$\begin{eqnarray}&&{P}_{\mathrm{AGN}}=\displaystyle \frac{{F}_{\mathrm{obs}}}{{F}_{\mathrm{obs}}-0.6{F}_{\mathrm{host}}}{P}_{\mathrm{obs}},\end{eqnarray} \tag{ 2 }$$

where F_obs is the observed flux density, F_host is the total flux density of the host galaxy, and the factor of 0.6 is the fraction of host-galaxy contamination. We have applied the correction for dilution of the polarization from the host-galaxy starlight to the values of the degree of polarization of BL Lac reported here. The position angle of polarization is unaffected by the light of the host galaxy.

2.3.4. TESS Data

TESS (Ricker et al. 2015) continuously monitors sectors of the sky at a wavelength band of 6000–10000 Å, shifting to different sectors after several weeks. This wavelength coverage nearly encompasses the R- and I-band WEBT coverage. A summary of the telescope specifications is provided in Table 9. TESS collects full-frame images (FFIs) of the entire field of view (FOV) every 30 minutes, with "postcard cutouts" of select targets obtained at a higher cadence of 2 minutes. A list of TESS identification numbers for sources is given by the TESS Input Catalogue (TIC; Stassun et al. 2018). Data products from the TESS mission are publicly available on the Mikulski Archive for Space Telescopes (MAST).³³

Table 9.  Summary of TESS Telescope Specifications

Attribute	Value
Single camera FOV	24° × 24°
Combined FOV	24° × 96° (3200 deg2)
Single camera aperture	10.5 cm
Focal ratio (f/#)	f/1.4
Wavelength range	6000–10000 Å
Pixel size on sky	21''
FFI exposure time	30 minutes
Orbital period	12–15 days

References. Ricker et al. (2015).

Download table as:  ASCIITypeset image

BL Lac (TIC 353622691) was observed by TESS from 2019 September 12 03:29:27 to October 6 19:39:27 UT, with its image located on camera 1 CCD chip 2. Figure 6 shows a cutout of a TESS FFI taken on 2019 October 4 07:15:36 UT, superposed on a Digitized Sky Survey image of the same field. The large pixel sizes and crowded field render photometry of faint targets difficult. However, BL Lac is usually one of the brightest objects in the field, and thus photometry is easier to perform.

Zoom In Zoom Out Reset image size

In order to carry out a preliminary exploration of the TESS data, we have used the eleanor software package (Feinstein et al. 2019) to reduce the FFI images for BL Lac and the four main WEBT-recommended comparison stars (see Table 8). The resulting light curves are displayed in Figure 7. We did not make any quality cuts while reducing the data, and for each source we defined a unique 3 or 4 pixel aperture that did not overlap with the aperture of another source of significant flux. The eleanor pipeline is able to remove the majority of systematic effects present in TESS data, as evidenced by the relatively constant flux of the four comparison stars. (Some small variations are seen in comparison stars close to BL Lac, but these variations are due to light from BL Lac leaking into the aperture of the comparison stars.)

Zoom In Zoom Out Reset image size

The removal of such artificial systematic trends from the data with 30 minute cadence is generally successful in stellar light curves made from TESS data. However, the long-term variability—months or longer—common in AGNs may be mistaken as instrumental effects and removed from the light curve by commonly used TESS data reduction software designed to search for evidence of exoplanets. Future studies, especially for sources in the CVZ, will need to be more carefully calibrated. We note that such a procedure is unnecessary for the present study, which focuses on short-term variability with 2 minute cadence data, as discussed below.

BL Lac was selected as a target for monitoring by TESS with a 2 minute cadence. These data were processed by the Science Processing Operations Center (SPOC; Jenkins et al. 2016). The pipeline performs both general CCD and pixel-level corrections, computes optimal apertures, completes a photometric analysis of the sources, and performs a "pre-search data conditioning" (PDC) procedure designed to take into account systematic effects. It also removes isolated outliers, corrects the flux of a source for crowding effects, and corrects for the aperture not containing all of the flux from a target source. The light curve obtained from this method is called the "PDCSAP" light curve.

There are several aspects of the PDCSAP light curve to check before using it in subsequent analyses. The pipeline calculates the optimal aperture, which contains only 66% of the total flux of the blazar (calculated using the pixel response function of the TESS detectors). Also, the flux from BL Lac represents only ∼30% of the total flux in the aperture from all sources.³⁴ We do not consider these issues as reasons to avoid using the PDCSAP light curve. The eleanor analysis in Figure 7 shows that the major bright, nearby stars are not variable, and that it is possible to separate the variable emission of BL Lac from the flux of stars within one TESS pixel. Also, we have normalized the light curve to its median value rather than convert the TESS electron counts to an energy flux density for the analysis presented in this paper.

We find no evidence of large-amplitude exponential decreases in the flux at the beginning of an orbit attributed to uneven heating of the telescope during data transmission modes (dubbed "thermal ramps") that are often present in TESS light curves. We thus make no cuts of the data near the beginning or end of an orbit.

We have checked the data quality raised by the pipeline for the PDCSAP light-curve data, but find an insignificant number of data points with quality issues. We have, however, eliminated 14 outlier points out of 16,006 observations of BL Lac.

The SPOC pipeline is subject to over-fitting similar to the eleanor program (see above). The PDC noise goodness metric (between 0 and 1), present in the header of the light-curve data product, is used as an indicator of over-fitting. The PDC noise goodness metric for BL Lac is 0.68, which implies a modest level of removal of intrinsic long-term trends. Since our study focuses on short-term variability of BL Lac, this over-fitting has an insignificant effect on our analysis.

The multiwavelength behavior of BL Lac over the entire WEBT monitoring campaign is shown in Figure 8. The optical data coverage is dense, especially in the R band.

Zoom In Zoom Out Reset image size

The γ-ray flux increased from ∼3 × 10⁻⁷ to 1.5 × 10⁻⁶ ph cm⁻² s⁻¹ over 10 days, peaking on September 29 (MJD: 58755), and then decayed quickly over the next two days. The optical light curves also rose to a peak in late September, although the details differ from the γ-ray behavior. In the R band, the underlying trend (defined by the minima of shorter-timescale variations) corresponded to an increase from ∼20 to ∼35 mJy before the flux density faded back down to 20 mJy. Shorter-timescale fluctuations, with durations of hours to days, occurred frequently during the monitoring period at optical wavelengths. Any similar events would be difficult to identify in the γ-ray light curve owing to the low-flux level and large number of non-detections. At least one event is identified at all wavelengths: the rapid brightening and decay on September 19 (MJD: 58745; see below).

The degree of linear polarization, P_R, fluctuated between 1% and 12% over the monitoring period. The average value, derived from the normalized q and u Stokes parameters, was $\langle {P}_{{\rm{R}}}\rangle =6.7 \%$ with a standard deviation of 2.1%. The mean electric-vector position angle was $\langle {\chi }_{{\rm{R}}}\rangle =-183^\circ \pm 15^\circ$, which is within 1σ uncertainty from the average 43 GHz radio jet direction of −173° (Jorstad et al. 2017; marked with a red dashed line in Figure 8(f)). Significant swings by up to ∼50° about this position angle were observed throughout the monitoring period.

Figure 9 shows the entire TESS light curve of BL Lac, along with the Fermi-LAT γ-ray and WEBT R-band light curves. The TESS count rates have been normalized to the median value of 1150.39 e⁻ s⁻¹. The WEBT and TESS light curves are very similar, despite the potential over-crowding, aperture, and over-fitting issues present in the PDCSAP light curve. The WEBT light curve, with its sparser sampling, is a reliable tracer of the major events and trends visible in the TESS light curve.

Zoom In Zoom Out Reset image size

In Figure 9 the peak of the γ-ray light curve on September 29 can be seen in greater detail. While the WEBT monitoring is sparse around this date, the TESS light curve shows a complicated structure during the γ-ray brightening. Of particular note is the large increase in optical flux density on September 19, clearly seen in both the WEBT and TESS light curves. This peak is also apparent in the γ-ray light curve as an increase in flux by a factor of ∼2.

Figures 10 and 11 display time variations of the flux or flux density and polarization of BL Lac over the five days of concurrent monitoring at all frequencies. All light curves show the same general trend of two periods of higher flux near September 15 and 19, labeled P1 and P2, respectively, with a period of lower flux in between. The high-flux state near P1 is most easily seen in the X-ray light curves (panels (b)–(d)). All of the light curves exhibit similar amplitudes of variability by a factor up to ∼2. The UV and optical light curves only show a moderate increase of flux density during P1 compared with the higher-amplitude increase of P2. The peak of P2 is very well sampled by both TESS and the WEBT observations, with a smooth rise and fall. The rise of P2 is also well sampled at higher energies; however, the observations ended before the decline could be detected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The optical linear polarization varied significantly near the periods P1 and P2, but was generally stable during the intervening low-flux plateau. During the plateau, P_R was high, near 9%, and χ_R was essentially parallel to the 43 GHz radio jet (Jorstad et al. 2017). The position angle χ_R was quite variable during P1, undergoing a ∼20° swing, but was relatively stable for several hours during P2.

Inspection of the light curves, especially those from TESS in Figures 9 and 11, reveals several periods of rapid changes in flux. In this section we calculate the shortest timescales of variability in the TESS and X-ray data. We discuss the variability of the optical polarization in Section 6.

4.1. Intraday Variability of TESS Data

Several statistical methods have been developed and applied to quasar variability in order to quantify the degree of short-timescale, including intraday, flux variations. de Diego (2010) compared several statistical tests using simulated light curves, and determined that a one-way analysis of variance (ANOVA) test is one of the most robust ways to identify statistically significant variability. An ANOVA test is a metric to judge the equivalence of measurements in a sample by breaking the sample into several groups and evaluating the means and variances of those groups. The null hypothesis for an ANOVA test is equality of all group means. Applied to blazar variability, this null hypothesis can be translated as non-variability of the source over the time period being tested. Two statistical measures are returned from the test, an F-statistic and a p-value. The null hypothesis can be rejected if (1) the returned F-statistic is greater than a critical value F_crit, calculated using the two degrees of freedom (dk1 = k − 1 and dk2 = N − k, where k is the number of groups and N is the total number of measurements in the sample) and a user-selected significance value, and (2) the returned p-value is smaller than the significance value. de Diego (2010) recommends using a number of groups k ≥ 5 in order for the test to have the most power to detect variability. In the following analyses, we label the critical value with the corresponding degrees of freedom as F_dk1,dk2.

An ANOVA test has been used to identify and characterize the optical flux density and polarization variability of the BL Mrk 421 (Fraija et al. 2017) and the FSRQ 3C454.3 (Weaver et al. 2019), with observed timescales of variability of ∼2 hr in both cases. The sampling rate of observations in these studies was on the order of several minutes between observations, for at most a few hours each night. In contrast, the TESS light curve of BL Lac, obtained at a 2 minute cadence over about 25 days, allows for a much more systematic survey of variability to be performed. This produces robust metrics to be used to test for variability (de Diego et al. 2015). Since the TESS light curves are evenly and densely sampled compared to the timescale of variations being investigated, a simple ANOVA test is sufficient (de Diego 2014).

We have broken the TESS PDCSAP light curve into hour-long sets, starting from the first TESS observation, each containing ∼30 data points.³⁵ In total, we perform tests on 516 hour long light-curve segments. All segments were normalized to the PDCSAP median value of 1150 e⁻ s⁻¹ in order to avoid the issues with the flux scaling present in TESS light curves (see Section 2.3.4). Following the recommendation of de Diego (2010), we passed each light curve through an ANOVA test, breaking these hourly segments into five groups (∼6 data points per group; df1 = 4 and average $\langle {df}2\rangle =26$). We have chosen a significance value of p < 0.001 (3σ) as the threshold for variability over the hour-long periods. Through this method, we have identified 107 hour long periods during which BL Lac was significantly variable (∼20% of the total observation period, after one accounts for the downlink break between orbits).

Figure 12 shows two examples of the analyzed hour-long light-curve segments. Panel (a) corresponds to the most statistically significant detected variability (starting on September 19 00:59:37 UT; MJD: 58745.0414), with F_4,26 =53.723 and p = 3.398 × 10⁻¹². During this hour, the flux increased by ∼10%. In contrast, panel (b) represents the hour-long period with the smallest F-value (starting on September 15 02:49:12 UT; MJD: 58741.1233), with F_4,26 = 0.046 and p = 0.996.

Zoom In Zoom Out Reset image size

We cannot conclusively determine whether the variability is slightly greater when the source is in a higher flux state. The average flux (normalized to the median value) of a variable hour was 1.036, with a standard deviation of 0.104. For the non-variable hours, the average flux was 0.998, with a standard deviation of 0.094. The difference is therefore not statistically significant.

With the ANOVA test showing that BL Lac is variable on sub-hour timescales during a significant fraction of the monitoring period, we now calculate the timescale of optical flux variability, τ_opt. We consider all pairs of flux measurements within the hour-long sets of data if, for a given pair, ${S}_{2}-{S}_{1}\gt \tfrac{3}{2}({\sigma }_{{S}_{1}}+{\sigma }_{{S}_{2}})$, where S_i and ${\sigma }_{{S}_{i}}$ refer to the flux and associated uncertainty of each measurement. Of the ∼50,000 possible pairs of observations, ∼11,000 met this uncertainty criterion. For these pairs, we calculate τ_opt using the formalism suggested by Burbidge et al. (1974): τ = Δt/ln(S₂/S₁) with S₂ > S₁, where ${\rm{\Delta }}t=| {t}_{2}-{t}_{1}|$ is the difference in the time of observation of each measurement. The average timescale is 15 hr, with a standard deviation of 7 hr. This is very similar to the derived minimum timescale of variability in the softer 3–10 keV band (see Section 4.2 below). The minimum calculated timescale of variability of the TESS data is 31 minutes. We plot the calculated timescales of variability in Figure 13, with histogram bins of 2 hr.

Zoom In Zoom Out Reset image size

4.2. X-Ray Variability

Analysis of multiwavelength IR/optical/UV data can identify separate components contributing to emission, each with its own continuum spectrum and variability properties. In order to isolate the contribution of the component of (likely synchrotron) radiation that is variable on the shortest timescales, we follow a method first suggested by Choloniewski (1981) and later developed by Hagen-Thorn (1997). A relative continuum spectrum can be constructed from essentially simultaneous flux density measurements in different bands by considering the slopes of the sets of cross-frequency flux density versus flux density (here shortened to "flux–flux") relations. This method has been successfully applied to a number of blazars (Hagen-Thorn et al. 2008; Larionov et al. 2008, 2010, 2016, 2020; Jorstad et al. 2010; Gaur et al. 2019). In the case of BL Lac (Larionov et al. 2010; Gaur et al. 2019), the relation between the optical and near-infrared flux densities over long timescales and major changes in flux density cannot be properly fit by a simple linear dependence. These authors obtained a second-order polynomial fit to the flux density of a given band i: ${F}_{\nu ,i}={a}_{i}{F}_{{\rm{R}}}^{2}+{b}_{i}{F}_{{\rm{R}}}+{c}_{i}$. They also found that the polynomial fits flatten toward higher frequencies, indicating that BL Lac exhibits a bluer-when-brighter trend, in agreement with other methods of determination of the spectral slope of the variable component (Villata et al. 2002, 2004b; Papadakis et al. 2007). The flux density range available to Hagen-Thorn et al. (2004) was not wide; hence, they did not detect any deviations from a linear dependence in the flux–flux plots.

Figure 14 shows the optical and UV flux–flux relations relative to the WEBT R band. To obtain this, we associated the UV data with the R-band observations that were nearest in time. For the BVI dependencies, only WEBT data from telescopes with quasi-simultaneous multiband observations of BL Lac were used from the entire time period of observations. The dependencies do not show any changes over time during the 3 months of observations. The optical (UBVI) behavior is shown in panel (a) of Figure 14, while the UV behavior is shown in panel (b).

Zoom In Zoom Out Reset image size

We have fit a straight line (F_ν = bF_R + c) and second-order polynomial (${F}_{\nu }={{aF}}_{{\rm{R}}}^{2}+{{bF}}_{{\rm{R}}}+c$) to the data. Table 11 gives the results of a χ² goodness of fit test for both fits for each band. In general, the χ² test indicates a slight preference for a second-order polynomial fit for almost every wave band. However, the difference between the χ² values for the two fits is small, as are the quadratic coefficients. We attribute this to the relatively modest amplitudes of variability during our observations, as was the case for Hagen-Thorn et al. (2004). Therefore, we use linear fits for subsequent analyses.

In Figure 15 we show the synthetic spectra of BL Lac for R-band flux densities F_R = 25, 40, and 50 mJy (in black, red, and blue, respectively). Table 12 lists the coefficients of power-law fits to the synthetic spectra of the form $\mathrm{log}({F}_{\nu })=\alpha \mathrm{log}(\nu )+\beta$, where α is the optical spectral index. Figure 15 shows that the synthetic spectra corresponding to different brightness levels have slightly different slopes (see also Table 12), with the slope flattening toward higher flux states. This bluer-when-brighter trend is also apparent in Figure 16, which displays the B−R evolution as a function of R-band brightness. Color indices made with combinations of the other available filters show similar trends.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 17 shows the SED of BL Lac at two different flux states, the low-flux plateau and peak P2. To describe the X-ray portion of the SEDs, we have calculated fluxes within seven energy intervals: 0.3–0.6, 0.6–1.2, 2.4–4.8, and 4.8–9.6 (Swift), plus 3–6, 6–12, and 12–24 keV (NuSTAR), for simultaneously measured data, as described in Section 2.2.3. BL Lac is brighter across all wave bands during event P2. The soft X-ray portion of the SEDs appears to include contributions from both steep-spectrum synchrotron and flatter-spectrum IC emission components. The synchrotron component becomes more prominent at higher flux states, while the IC component dominates at hard X-ray energies. As mentioned in Section 2.2.3, the break energy moves to higher energies as the X-ray flux increases.

Zoom In Zoom Out Reset image size

In the optical region, the spectral indices flatten from the B band toward the UV bands. Raiteri et al. (2009) have inferred a contribution from thermal emission in the optical/UV spectrum of BL Lac, which they attributed to an accretion disk component. They have modeled such a disk as blackbody emission with a temperature ≥20,000 K and a luminosity ≥6 × 10⁴⁴ erg s⁻¹. To investigate whether such a thermal component contributes to the SED in Figure 17, we have restricted the high- and low-flux SEDs to optical/UV frequencies in Figure 18(a). We apply a power-law model (drawn with dotted lines) for both states, adopting the spectral indices from Table 12. The small residuals of the model (Figure 18(b)) indicate that a single power law can adequately describe the optical/UV emission. (In all of our SEDs, the flux density in the B filter is low compared to the expected power-law model, regardless of whether the data is taken from the WEBT or Swift observations. We suggest that this low-flux density might arise from the wide filter bandpass and spectral shape of the source.) While a single power-law component provides a general description of the observed flux density, we do see an excess in the UV portion of the spectrum of the low-flux state compared to the high-flux state. This difference is clearly seen in Figure 18(c), with the difference between the high- and low-flux residuals increasing toward shorter wavelengths. This suggests that a UV excess occurs at the low-flux state. However, an accretion disk component with the same characteristics as reported by Raiteri et al. (2009) would significantly exceed the optical-UV SED when added to the synchrotron component. Taking into account that accretion disk emission can significantly change on timescales of months to years (Kaspi et al. 2000), this may indicate a weakening of the disk contribution in recent years.

Zoom In Zoom Out Reset image size

Over the entire WEBT campaign, we measured P_R at 303 times, with the densest sampling during the five days of intensive X-ray monitoring. The data are displayed in Figures 8(e) and (f). The degree of polarization P_R fluctuated rapidly between 1% and 12% over most of the observing period, although it was stable at ∼9% from mid-October to the beginning of November. The average uncertainty on a measurement of P_R is $\langle {\sigma }_{{P}_{{\rm{R}}}}\rangle =0.23 \%$.

Results of an ANOVA test can be considered reliable if the variable being tested is approximately normally distributed. However, as mentioned earlier, the degree of polarization follows a Rice distribution. Several Monte Carlo simulations have been performed to show that, as the sample size increases, the ANOVA test is robust against violations from normality (Donaldson 1966; Tiku 1971). As the number of measurements of P_R ∼ 300 and the polarization data have been corrected for the Rice bias, we have thus used an ANOVA test, as described above, on the values of P_R over the entire period. We have calculated the F-statistic for five groups to be F_4,298 = 13.47, with a p-value of p = 4.2 × 10⁻¹⁰. This confirms that P_R was variable over the entire time period.

In order to determine the timescale of variability of P_R, we have searched for all pairs of measurements between which the values of P_R differed by a factor ≥2, in order to calculate the halving or doubling timescale τ₂. We restrict the analysis to the well-sampled observations between September 14 and 21 (MJD: 58740–58747), with Δt < 50 hr. We then identify the shortest timescale of variability under the constraint that the measurements meet the same uncertainty criterion as we used for the flux density: $| {P}_{i}-{P}_{j}| \gt (3/2)({\sigma }_{{P}_{i}}+{\sigma }_{{P}_{j}})$. The date restriction reduces the number of measurements to 223, yielding ∼25,000 data pairs. Of these, only 341 pairs meet the minimum Δt and uncertainty criteria. For each of these data pairs, we assumed a linear change with time to determine τ₂.

Figure 19 shows the distribution of τ₂ obtained from this fitting method. The minimum value is τ₂ = 5 hr. The peak of the distribution lies in the 20–25 hr time bin, longer than that of the optical flux variations but shorter than than that of the X-ray variations.

Zoom In Zoom Out Reset image size

The position angle of the polarization χ_R was aligned with the direction of the 43 GHz parsec-scale radio jet (Jorstad et al. 2017) at the beginning of the monitoring period. Swings away from this polarization angle are seen in September during the high optical flux state.

At the beginning of the week of concurrent observations at all wavelengths, P_R rose from ∼5% to ∼9% near in time to the X-ray brightening event P1. This was accompanied by a swing in χ_R from parallel to the radio jet to ∼45° away. While χ_R returned to nearly parallel to the jet direction shortly thereafter, P_R remained high for two days, slowly decreasing from 9% to 6%.

Complicated polarization behavior is seen during the periods P1 and P2. A large EVPA swing and increase in P_R accompanied P1 when only the high-energy light curves showed a pronounced peak. However, a relatively stable EVPA ∼15° from the jet direction and varying P_R occurred during P2 when pronounced changes in flux occurred at all wavelengths. Such behavior, with a stable EVPA despite changing flux and degree of polarization, has been observed in BL Lac in the past (Gaur et al. 2014), which is able to be incorporated into a shock-in-jet model if the shock Lorentz factor is low.

Figure 20 shows the polarization changes in the normalized Stokes parameter q versus u plane for the time period September 14–20 based on their measured values. The polarization of BL Lac during this time was characterized by high positive values of q. A smooth connection can be made between the polarization states of BL Lac on different nights, appearing as swings or rotations in the q–u plane, although this trend does not change χ_R significantly (see Figure 11(e)).

Zoom In Zoom Out Reset image size

In order to characterize the origin of the rotation in the q–u plane, we display in Figure 21 the change in q and u as a function of the R-band flux density at flux peak P2 seen across all wavelengths (September 18 18:00:00–20 00:00:00 UT). The relationship between q and u is markedly nonlinear. A linear relation would be expected if a single variable source were responsible for the variability of both the flux and polarization (Hagen-Thorn & Marchenko 1999). The behavior could instead result from the superposition of emission from a number of components with different flux and polarization parameters, as we discuss in Section 8.2 below.

Zoom In Zoom Out Reset image size

7.1. Short-timescale Correlations

We perform a correlation analysis between each of the light curves presented above and the TESS light curve. We use the z-transformed discrete correlation function (ZDCF; Alexander 1997), which provides a properly normalized (between −1 and 1) interval of results, along with uncertainties. The latter are sampling errors based on the noise of the original data, which we calculate using the recommended 100 Monte Carlo draws.

We use the Peak Likelihood algorithm (PLIKE; Alexander 2013) to estimate the maximum of the cross-correlation function (CCF) and the 1σ fiducial distribution of the lag interval. The PLIKE algorithm provides an estimate of the CCF peak and the uncertainty of the time lag without any a priori assumptions about either the shape of the CCF peak or models of the light curves.

We verify the significance of the correlations by comparing the ZDCF results with the statistics of correlations on 3000 pairs of bootstrapped artificial light curves (ALCs). We have removed points with excessively large errors per the recommendation of Alexander (1997), thus eliminating any need to weigh selected points. We generate the 3000 ALCs by randomly selecting (with replacement) and placing flux points and associated uncertainties on the preserved observational dates. Once the ALCs have been built for each light curve to be compared, we randomly pair the ALCs (without replacement) and compute the ZDCF. Results of this analysis are used to derive a 2σ probability of obtaining a given coefficient of correlation by chance. In all cases where the ZDCF values of the actual light-curve correlation are greater than 0.8, the bootstrap analysis generally gives a lag time consistent with the peak of the ZDCF. When the ZDCF is weaker, the bootstrap analysis still generally agrees, but frequently with less than 2σ significance. It is important to note that no result smaller than the bin size of the data is meaningful. In the case of the γ-ray light curve, the shortest meaningful time delay is 0.25 day. To reduce the impact of the upper limits, we use flux values of half of the 2σ upper-limit values of the γ-ray data in all of the correlations.

Table 13 lists values of the ZDCF peaks and their PLIKE lag fiducial intervals. A graphical representation of the ZDCF of several light curves with the TESS light curve is shown in Figure 22. Due to the time binning of the γ-ray light curve, the results of the correlation analysis are consistent with zero lag between the γ-ray and optical (TESS) light curves. However, there is a significant lag between the optical and X-ray variations. In particular, the TESS light curve leads the X-ray variations by ${0.38}_{-0.11}^{+0.08}$, ${0.29}_{-0.09}^{+0.04},$ and ${0.051}_{-0.005}^{+0.007}$ days for the Swift 0.3–3 keV, NuSTAR 3–10 keV, and NuSTAR 10–22 keV light curves, respectively. This trend of optical leading X-ray variations is seen with the WEBT B- and R-band and Swift UV light curves as well.

Zoom In Zoom Out Reset image size

Table 13.  Results of Multiwavelength Correlation Analysis

Band 1	Band 2	PLIKE interval (days)			ZDCF value of peak	Bootstrap significance
		Lower bound	Peak	Upper bound		of ZDCF peak
Gamma-Ray	TESS	−0.158	−0.105	+0.016	${0.314}_{-0.100}^{+0.097}$	3.3
NuSTAR 10–22 keV	TESS	+0.046	+0.051	+0.058	${0.893}_{-0.063}^{+0.049}$	4.2
NuSTAR 3–10 keV	TESS	+0.196	+0.286	+0.326	${0.874}_{-0.029}^{+0.026}$	11.0
Swift XRT 0.3–3 keV	TESS	+0.262	+0.376	+0.457	${0.789}_{-0.065}^{+0.057}$	6.8
Swift UVW2	TESS	−0.018	+0.063	+0.089	${0.941}_{-0.020}^{+0.017}$	8.8
Swift UVM2	TESS	−0.021	−0.001	+0.035	${0.952}_{-0.019}^{+0.016}$	7.3
Swift UVW1	TESS	−0.020	+0.014	+0.053	${0.962}_{-0.013}^{+0.011}$	8.7
Swift U	TESS	−0.036	−0.013	+0.017	${0.965}_{-0.012}^{+0.010}$	8.8
WEBT B	TESS	−0.046	−0.017	−0.007	${0.970}_{-0.005}^{+0.004}$	18.6
WEBT R	TESS	−0.028	−0.015	−0.003	${0.963}_{-0.003}^{+0.003}$	32.7

Note. Positive lags indicate that band 2 leads band 1.

Download table as:  ASCIITypeset image

7.2. Long-timescale Correlations

It has been noticed in several FSRQs that the long-timescale optical and γ-ray variations are well correlated, with lags between 0 and ∼3 days, with the optical leading the γ-ray variations in some cases and the opposite in others (e.g., in the well-studied case of 3C454.3; Bonning et al. 2009; Gaur et al. 2012; Gupta et al. 2017; Kushwaha et al. 2017). We now use the 3 month long WEBT R-band and Fermi γ-ray light curves to investigate whether such a correlation occurred in BL Lac during our monitoring period. We have binned the WEBT R-band observations to match the binning of the γ-ray light curve, taking the average time of the R-band observations for the correlation.

The ZDCF of the data is shown in Figure 23. Three major peaks are seen in the correlation, corresponding to the WEBT R band leading the γ-ray variations by 0.2, 3.4, and 10.0 days. We have examined each peak more closely, and the PLIKE maximum-likelihood interval, value of the ZDCF peak, and bootstrap significance of each peak are given in Table 14. Again, no result smaller than the bin size of the data is meaningful, so the lag of ∼0.2 day is essentially the same as zero lag.

Zoom In Zoom Out Reset image size

Table 14.  Results of WEBT R-band and γ-ray Correlation Analysis

PLIKE interval (days)			ZDCF value	Bootstrap significance
Lower bound	Peak	Upper bound	of peak	of ZDCF peak
−0.21	+0.23	+2.04	${0.46}_{-0.06}^{+0.06}$	8.27
+0.25	+3.41	+5.09	${0.46}_{-0.64}^{+0.06}$	7.75
+9.95	+9.95	+12.94	${0.34}_{-0.07}^{+0.07}$	5.44

Note. Positive lags indicate that the R band leads the γ-ray band.

Download table as:  ASCIITypeset image

To visualize how such local maxima in the ZDCF arise, Figure 24 shows the γ-ray and WEBT R-band light curves, with the γ-ray light curve shifted by the lags indicated in Table 14. We have normalized the γ-ray light curve by its median value and the R-band light curve by half of its median value to more closely inspect the variations. Both light curves have been vertically shifted to not overlap.

Zoom In Zoom Out Reset image size

Shifting the γ-ray light curve in time relative to the optical light curve reveals the cause of the peaks in the correlation. While the majority of the light curves feature low-amplitude flux variations, two large-amplitude variations dominate the correlation: matching the γ-ray peak near P2 with the optical peaks near P1 and P2 causes the 0.2 and 3.4 days correlations, while matching the large-amplitude γ-ray peak, G, with optical peak P2 produces the local maximum in the ZDCF near 10 days. In fact, the matching of G with P1 is also seen in the ZDCF as a second, slightly less statistically significant bump at a lag of ∼13 days in Figure 23.

The analysis of the local maxima of the ZDCF given in Table 14 finds that the statistical significance of the local maximum of the ZDCF at 0.2 day lag is higher than that at the other lags. Therefore, the correlation between the γ-ray and optical flux variations that is prominent during the 5 day period of intensive monitoring is also present in the 90 day data. The local maxima at the longer lags are less significant, resulting from offsets of peaks in the light curves that may be physically unrelated.

8.1. Summary of Variability of BL Lac

8.2. Interpretation

The nonlinear relation of q and u suggests the superposition of a number of emission components with different flux and polarization parameters. This can be interpreted as evidence for turbulence in the jet (e.g., Marscher 2014), which can be roughly modeled as N_turb cells, each with a uniform but randomly oriented magnetic field (Burn 1966). Under such a model, the degree of linear polarization is $\langle P\rangle \approx {P}_{\max }{N}_{\mathrm{turb}}^{-1/2}$ (Burn 1966). Here P_max corresponds to the degree of polarization of incoherent synchrotron emission in a uniform magnetic field, related to the spectral index α as: ${P}_{\max }\,=(| \alpha | +1)/(| \alpha | +5/3)\,\ast \,100$ (%). Adopting the value $| \alpha | \,=1.17$ corresponding to the optical spectrum when the flux density was at its median value, F_R = 40 mJy (Table 12), we obtain P_max = 76%. The standard deviation of the polarization about this mean is predicted to be ${\sigma }_{P}\approx 0.5{P}_{\max }{N}_{\mathrm{turb}}^{-1/2}$ (Jones 1988). For the observed mean polarization of 6.7%, the required number of turbulent cells is N_turb ≈ 130. The measured standard deviation of 2.1% is somewhat less than the value of 3.4% expected for the turbulence model. This can be reconciled by the existence of partial ordering of the magnetic field. For example, the partial ordering along the jet direction inferred from the mean EVPA can be the consequence of compression of the turbulent plasma by a shock (e.g., Hughes et al. 1985; Cawthorne 2006; Marscher 2014) or the superposition of turbulence and a helical magnetic field (e.g., Raiteri et al. 2010; Gabuzda 2018).

The analysis presented in Section 7 reveals correlations between the TESS light curve and all others. The X-ray variations lag behind the TESS light curve by ∼0.38 ^+ 0.08_−0.11 days at 0.3–3 keV, ${0.29}_{-0.09}^{+0.04}$ days at 3–10 keV, and ${0.051}_{-0.005}^{+0.007}$ days at 10–22 keV, respectively. The cross-frequency correlations and time lags imply that the emission regions in the various optical, X-ray, and γ-ray bands are partially co-spatial. Such correlations and lags are expected if the variable X-ray emission in BL Lac is mainly caused by inverse Compton scattering of synchrotron seed photons in the frequency range of ∼10¹² to ∼10¹⁵ Hz by relativistic electrons that are energized at a front in the jet, e.g., a shock (Marscher & Gear 1985). In this case, higher-energy electrons maintain their energy over a shorter distance beyond the shock than do their lower-energy counterparts. This leads to progressively longer lags of the synchrotron or inverse Compton variations toward lower frequencies that can be derived as follows.

The time lag between acceleration and energy loss is given in the rest frame of the emitting plasma by

$$\begin{eqnarray}&&{t}_{\mathrm{loss}}^{{\prime} }\sim 7.7\times {10}^{8}{\left[\left({B}^{2}+8\pi {u}_{\mathrm{ph}}\right)\gamma \right]}^{-1}\,{\rm{s}},\end{eqnarray} \tag{ 3 }$$

where B is the magnetic field strength in Gauss, u_ph is the energy density of seed photons for inverse Compton scattering in erg cm⁻³, and γ is the electron energy in rest-mass units. In the observer's frame, ${t}_{\mathrm{loss}}={t}_{\mathrm{loss}}^{{\prime} }(1+z)/\delta$, where δ is the Doppler factor. The SED displayed in Figure 17 indicates that the synchrotron IR and inverse Compton γ-ray luminosities are comparable. This implies that B² ∼ 8πu_ph, which we assume to be the case. The value of γ relates to the frequency of observation ν (in Hz) as

$$\begin{eqnarray}&&\gamma \sim 6.0\times {10}^{-4}{\left[\nu (1+z)/(B\delta )\right]}^{1/2}\end{eqnarray} \tag{ 4 }$$

for synchrotron radiation, where δ ≈ 8 (Jorstad et al. 2017) is the Doppler factor, and z = 0.069 is the host-galaxy redshift. The mean value of γ for inverse Compton scattering also depends on frequency as ν^1/2. We then obtain

$$\begin{eqnarray}&&{t}_{\mathrm{loss}}\sim 0.8{\left[\displaystyle \frac{B(1+z)}{\delta }\right]}^{1/2}{\left(\displaystyle \frac{\nu }{{\nu }_{\mathrm{TESS}}}\right)}^{-1/2}{\left({B}^{2}+8\pi {u}_{\mathrm{ph}}\right)}^{-1}\,\end{eqnarray} \tag{ 5 }$$

where ν_TESS = 4 × 10¹⁴ Hz is the median frequency of the TESS band, and t_loss is in days. Figure 25 presents the cross-frequency lag data relative to the TESS light curve, along with the best fit to the equivalent frequency dependence from Equation (5), i.e., t_loss(ν) − t_loss(ν_TESS). In the fit, there is zero delay between the TESS light curve and that at 300 keV (a free parameter, since we did not observe at an X-ray energy where there was zero lag). The magnetic field value of the fit is ∼3 G. Despite the uncertainties in the lags between TESS and the optical-UV light curves, the magnetic field strength is well specified by only the X-ray lags, which reflect the radiative energy losses. The ν^−1/2 relation provides a good fit to the lag data, with a reduced χ² = 0.9.

Zoom In Zoom Out Reset image size

We can also equate the radiative energy loss timescale to the minimum timescale of variability in the TESS band, 0.5 hr, and solve Equation (5) for B. This independent calculation also results in a value of ∼3 G.

Although the X-ray lags relative to the TESS light curve agree with the above model, the SEDs displayed in Figure 17 pose a problem. Although the hard X-ray spectrum has a spectral index $| \alpha | \lt 1$ expected for inverse Compton scattering, the soft X-ray slope is steep during the high-flux states, α ≈ −1.4 (see Figure 1). This suggests that the soft (0.3–3 keV) X-rays arise from a combination of synchrotron radiation by the highest-energy electrons whose energy distribution is steepened by radiative energy losses, and inverse Compton scattering by electrons with lower energies for which the radiative losses are modest. The above explanation of the X-ray time lag requires inverse Compton scattering to be the dominant emission process of the variable component of the flux. The X-ray spectral slope during the first maximum of the light curves, P1, is quite steep, with α ≈ −1.5, which strongly suggests that the soft X-ray flare represents the high-frequency end of the synchrotron spectrum. This conclusion, which agrees with previous studies (e.g., Raiteri et al. 2010), is supported by the pronounced variability in the 0.3–3 keV energy band near the peak (see Figure 10), since the highest-energy electrons that emit synchrotron radiation at X-ray frequencies have the shortest timescales of energy loss. If this is true, then the soft X-ray variations should lead the TESS light curve during the flare. Indeed, Figures 10 and 11 show that a local maximum in the TESS light curve lags behind the soft X-ray maximum. The correlation that gives a delay between the TESS and 0.3–3 keV light curve then must arise from the long lower-flux period, when the soft X-ray spectrum was flatter, and therefore the contribution of the inverse Compton component was more important. This is corroborated by the broken power-law fits to the combined Swift and NuSTAR spectra: the break energy shifts from ∼6 to ∼2 keV between the high and low X-ray flux states.

The magnetic field value that we infer from the time lags, ∼3 G, is ∼10 times higher than that derived via the "core shift" method (Lobanov 1998) by O'Sullivan & Gabuzda (2009) for the 43 GHz VLBI core located at a distance of ∼0.5 pc from the vertex of the jet. However, the core shift method applies to the ambient jet rather than a shock that energizes electrons as in the above scenario. Compression of the magnetic field in such a shock is expected to increase the field strength by a factor roughly equal to the Lorentz factor, which is estimated to be ∼6 in BL Lac (Jorstad et al. 2017). This implies that the variable emission reported here occurs ∼0.3 pc from the jet vertex, upstream of the 43 GHz "core."

We have carried out a high time-resolution, multiwavelength observing campaign of BL Lacertae, including monitoring at 2 minute cadence with TESS, in order to investigate the short-timescale variability of the blazar. Our data set includes: (1) three months of observations with the Fermi-LAT and ground-based WEBT-affiliated telescopes, (2) 25 days of monitoring with TESS, and (3) five days of densely sampled NuSTAR and Swift measurements.

All of the optical, UV, and X-ray light curves exhibit a similar trend during the five days of concurrent monitoring. Two high-flux states are separated by a low-flux plateau. The fractional amplitude of the variations is proportional to frequency up to at least the NuSTAR hard X-ray band. The minimum timescale at optical wavelengths is very short, ∼30 minutes, while the average is 15 hr, very similar to the minimum observed X-ray timescale of 14.5 hr.

Our analysis of the observations confirms statistically significant correlations among the light curves at all frequencies. Frequency-dependent time lags relative to the TESS variations can be explained by a model involving energization of the radiating electrons at a front, such as a shock, beyond which radiative energy losses restrict the emission to smaller volumes at higher frequencies (Marscher & Gear 1985). Both the minimum timescale of variability in the TESS band and the values of the time lags agree with such a model if the magnetic field is ∼3 G.

Consistent patterns of light curves, SEDs, and polarization versus time have proven elusive to find in blazar data. This is a consequence of both complexity in the physical processes in blazar jets and gaps in time and frequency coverage of monitoring programs. As our study demonstrates, the latter deficiency can be overcome by organizing intensive monitoring programs with current space- and ground-based facilities. Of particular importance to such efforts are instruments capable of essentially continuous monitoring, such as TESS. Future similar campaigns with even longer durations are likely to provide further valuable insights into the time-variable phenomena that occur in relativistic jets.

We gratefully acknowledge the comments and suggestions provided by the anonymous referee that have improved this work. The data collected by the WEBT Collaboration are stored in the WEBT archive at the Osservatorio Astrofisico di Torino—INAF (https://www.oato.inaf.it/blazars/webt/); for questions regarding their availability, please contact the WEBT President Massimo Villata (massimo.villata@inaf.it). The research at Boston University was supported by NASA grants 80NSSC19K1731 (TESS Guest Investigator Program), 80NSSC20K0080 (NuSTAR Guest Investigator Program), 80NSSC17K0649 (Fermi Guest Investigator Program), and Massachusetts Space grant 316080, as well as by Boston University Hariri Institute Research Incubation Award 2019-03-007. This research was partially supported by the Bulgarian National Science Fund of the Ministry of Education and Science under grants DN 18-10/2017, DN 18-13/2017, KP-06-H28/3 (2018), and KP-06-PN38/4 (2019). The Skinakas Observatory is a collaborative project of the University of Crete, the Foundation for Research and Technology—Hellas, and the Max-Planck-Institut für Extraterrestrische Physik. The work at the University of Sofia was supported by Bulgarian National Science Fund under grant DN18-10/2017 and National RI Roadmap Projects DO1-277/16.12.2019 and DO1-268/16.12.2019 of the Ministry of Education and Science of the Republic of Bulgaria. E.S. and S.I. were supported by funds of the project RD-08-122/2020 of the University of Shumen, Bulgaria. K.M. acknowledges JSPS KAKENHI grant No. JP19K03930. G.D. gratefully acknowledges the observing grant support from the Institute of Astronomy and NAO Rozhen, BAS, via the bilateral joint research project "Gaia Celestial Reference Frame (CRF) and Fast Variable Astronomical Objects." This work is a part of the following projects: No. 176011 "Dynamics and Kinematics of Celestial Bodies and Systems," No. 176004 "Stellar Physics," and No. 176021 "Visible and Invisible Matter in Nearby Galaxies: Theory and Observations," supported by the Ministry of Education, Science, and Technological Development of the Republic of Serbia. S.O.K. acknowledges financial support by Shota Rustaveli National Science Foundation of Georgia under contract PHDF-18-354. This work is partly based upon observations carried out at the Observatorio Astronómico Nacional on the Sierra San Pedro Martir (OAN-SPM), Baja California, Mexico. E.B. acknowledges financial support from UNAM-DGAPA-PAPIIT through grant IN113320. Support for K.L.S. was provided by NASA through Einstein Postdoctoral Fellowship Award No. PF7-180168. This article is partly based on observations made with the LCOGT Telescopes, one of whose nodes is located at the Observatorios de Canarias del IAC on the island of La Tenerife in the Observatorio del Teide. The work at Colgate University was supported by the Colgate University Research Council. We gratefully acknowledge the contribution of data from D. Blinov and the Robopol program. This study includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA Explorer Program. The TESS data analysis benefited from conversations and developmental work that took place at the Expanding the Science of TESS meeting at the University of Sydney. This study has used Digitized Sky Survey images based on photographic data obtained with the Oschin Schmidt Telescope on Mount Palomar to aid in the reduction and visualization of the TESS data. The Palomar Observatory Sky Survey was funded by the National Geographic Society. The Oschin Schmidt Telescope is operated by the California Institute for Technology and Palomar Observatory. The plates were processed into the present compressed digital format with their permission. The Digitized Sky Survey was produced at the Space Telescope Science Institute (STScI) under U.S. Government grant NAG W-2166. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of Astropy,³⁶ a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018), and astroquery, an astronomical web-querying package in Python (Ginsburg et al. 2019).

Facilities: Fermi - Fermi Gamma-Ray Space Telescope (formerly GLAST), NuSTAR - , Swift - , TESS - , WEBT - , Gaia. -

Software: Astropy (Astropy Collaboration et al. 2013, 2018), astroquery (Ginsburg et al. 2019), eleanor (Feinstein et al. 2019), HEAsoft (NASA High Energy Astrophysics Science Archive Research Center (HEASARC), 2014), Fermi Science Tools (Fermi Science Support Development Team 2019), pyrallaxes (Luri et al. 2018).