© The Authors 2025

1. Introduction

Blazars are active galactic nuclei (AGNs) with relativistic jets pointing toward the observer. Their central engine is likely to be a supermassive black hole of 10⁷−10⁹ solar masses, fed by the infall of matter from a surrounding accretion disk. Although they constitute only a tiny fraction of the astrophysical sources observed in the optical band, blazars are by far the most common type of objects detected at γ-ray energies (Abdollahi et al. 2020).

Radiative emission from blazars is mostly nonthermal radiation and it ranges from radio to very high-energy (VHE; 100 GeV<E<100 TeV) γ rays. The broadband spectral energy distribution (SED) of blazars is characterized by two distinct bumps (e.g., Ghisellini et al. 2017). The low-energy one peaks in the infrared-to-X-ray energy range and it is commonly attributed to synchrotron radiation emitted by relativistic electrons accelerated in the jet. The high-energy bump peaks above mega-electronvolt energies and it is most likely due to inverse Compton (IC) scattering. The photon seeds for the IC scattering can be the synchrotron ones from the same electron population (synchrotron self-Compton, SSC; Konigl 1981; Maraschi et al. 1992) or can be external to the jet, such as radiation from the broad-line region (BLR), accretion disk, and dusty molecular torus (Shakura & Sunyaev 1973; Dermer & Schlickeiser 1994; Finke 2016). The presence of a subdominant hadronic component is also possible, as has been discussed in Aharonian (2000) and Murase et al. (2012) and suggested by the evidence for neutrino emission from active galaxies (Aartsen et al. 2018a, b; Ansoldi et al. 2018; Abbasi et al. 2022).

According to the features in their optical/UV spectra, blazars can be classified as flat spectrum radio quasars (FSRQs), which are characterized by the presence of strong emission lines in their optical/UV spectra, and BL Lacertae objects (BL Lacs), having no or weak emission lines. BL Lacs can be further divided into three subclasses based on their synchrotron peak frequency, ν_s; that is, low-frequency-peaked (LBLs, ν_s<10¹⁴ Hz), intermediate-frequency-peaked (IBLs, 10¹⁴ Hz<ν_s<10¹⁵ Hz), and high-frequency-peaked (HBLs, ν_s>10¹⁵ Hz) BL Lacs (e.g., Padovani & Giommi 1995). Most of the blazars detected up to VHE γ rays are HBLs. At the time of writing (March 2025), according to TeVCat¹, VHE γ-ray emission has been detected from 57 HBLs. Conversely, only ten IBLs and two LBLs have been detected at VHE γ rays, usually during flaring episodes (e.g., Ahnen et al. 2018). In high-energy (HE; 100 MeV<E<100 GeV) γ rays, the majority of the BL Lac objects detected by the Large Area Telescope (LAT) on board the Fermi Gamma-Ray Space Telescope are LBLs and IBLs (e.g., Ajello et al. 2022). The lack of low-frequency peaking sources in the VHE γ-ray band is mainly due to their high-energy bump being located at lower energies than for HBLs.

Blazars are characterized by high variability over very different timescales. Lightcurves can show long-term trends over periods lasting from several years to months and short-term activity on timescales from weeks to days. For several blazars, intra-night and intra-day variability on timescales of minutes to hours has been assessed (e.g., Sillanpaa et al. 1988; Goyal et al. 2018). The fastest variations are observed especially during flaring episodes. The amplitudes of the variations in time of the radiative emissions of blazars are energy-dependent. In the case of the HBL Mrk 421, the largest variability amplitudes have been detected in the energy bands corresponding to the falling tails of the two bumps of the SED, in hard X-rays (from a few kilo-electronvolts up to several tens of kilo-electronvolts) and in the VHE γ-ray band (e.g., Acciari et al. 2020a).

In order to perform detailed studies of the radiative emission mechanisms acting in blazars, it is fundamental to have a complete energy coverage of their emission, with multiwavelength (MWL) observations covering from radio to VHE. As blazars may show extremely fast variable behavior during flares, it is of crucial importance that these observations be performed (quasi-) simultaneously, in order to provide a proper characterization of the source emission state. However, especially for weak blazars such as the one studied in this work, complete coverage may be difficult to obtain in all energy ranges during low emission states.

The blazar B2 1811+31, located at RA and Dec (J2000) 18h13m35.2028s, +31d44m17.621s (Petrov 2011), is classified as a BL Lac object in the Fermi-LAT Fourth Source Catalog (4FGL, Abdollahi et al. 2020) and as an IBL in Laurent-Muehleisen et al. (1999). The blazar was listed as one of the most promising VHE γ-ray candidates in Fallah Ramazani et al. (2017). Optical spectroscopic observations reported in Giommi et al. (1991) settled its redshift to z = 0.117.

Following the detection by the Fermi-LAT of a high state from the source in the E>100 MeV energy range on October 1, 2020 (MJD 59123) (Angioni et al. 2020), a MWL observational campaign on B2 1811+31 was organized. The observations performed during this high-state period with the Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) telescopes led to the first-time detection of VHE γ-ray emission from the source (Blanch 2020). The telescopes on board the Neil Gehrels Swift Observatory, which are sensitive in the optical-to-X-ray range, joined the follow-up campaign, as well as several optical and radio ground-based telescopes. These observations allowed us to characterize in detail the properties of the source high state from radio to VHE γ rays.

In addition, we include in this paper a comprehensive MWL dataset from 2005 up to 2024 extending from radio to HE γ rays. B2 1811+31 was a frequent target of Swift observations since 2005. It has been monitored since 2009 at 15 GHz by the Owens Valley Radio Observatory (OVRO) and it was included over the years in several monitoring programs in the optical band. Ultimately, B2 1811+31 has been observed by Fermi-LAT for its entire mission thanks to the continuous sky-survey operations. The source has been significantly detected in HE γ rays already in its first months of operations and was already included in the Fermi-LAT First Source Catalog (1FGL, Abdo et al. 2010a). We present the analyses of these long-term MWL data to contextualize the high-state episode and to compare the spectral and temporal features of the flaring and steady low-state emissions.

The paper is structured as follows. Section 2 is dedicated to present the MWL dataset. The analysis results of data from MAGIC, Fermi-LAT, Swift-XRT and Swift-UVOT are presented, along with the reduction of the optical and radio data. Section 3 is dedicated to the variability, intraband and multiband correlation analyses carried out to infer insights about the emission regions responsible for the MWL emission. In Section 4, we report on the classification of the source high and low states according to the partition based on the synchrotron peak frequency. A discussion on the literature on the source redshift, along with the redshift indirect estimation from simultaneous HE and VHE observations, is presented in Section 5. In Section 6, a leptonic interpretation model for the SED reconstructed in quasi-simultaneous observations is presented and its self-consistency and physical implications are discussed. We summarize the main results and present our conclusions in Section 7.

2. Instruments and analysis

In this section we present the MWL datasets and analyses performed in each energy band of B2 1811+31. Table 1 reports the list of the instruments whose data from the 2020 γ-ray high state were included in the analysis, along with the time ranges of the observations. From the same instruments, we collected the available long-term data. In addition, we included long-term optical data from the Katzman Automatic Imaging Telescope (KAIT, Filippenko et al. 2001), the Catalina Real-Time Transient Survey (CRTS, Drake et al. 2009) and from the Tuorla blazar monitoring program² (Nilsson et al. 2018). Further details on each dataset are given in the following sections.

The long-term MWL lightcurve of B2 1811+31 collected from 2005 to 2023 is shown in Figure 1. Figure 2 shows a close-up view of the Fermi-LAT and optical R-band lightcurves over a 250-day period which approximately corresponds to the source HE γ-ray high state in 2020 (Section 2.2). Additionally, Figure 3 provides a zoom-in on the MWL lightcurve over a 70-day period covering the MAGIC observations and the MWL observational campaign following the high-state detection by Fermi-LAT on October 1, 2020 (MJD 59123).

Fig. 1. B2 1811+31 long-term lightcurve collected from 2005 to 2023. From top to bottom panels: HE γ-ray flux above 100 MeV from Fermi-LAT monthly binned data, X-ray flux in the 0.3−10 keV energy range from Swift-XRT, optical R-band data and radio data. The dashed red line marks the Fermi-LAT high-state detection on October 1, 2020 (MJD 59123). The shaded light orange, light blue, and green bands in the top panel indicate the “Pre-flare”, “Flare” and “Post-flare” periods, respectively (Table 2).

Fig. 2. MWL lightcurve of B2 1811+31 in HE γ-rays (weekly binned Fermi-LAT data, upper panel) and optical R-band (instruments are reported in the legend, lower panel) in a period of around 250 days surrounding the October 2020 high state. This period corresponds approximately to the ‘Flare’ period (Table 2). In the upper panel, the solid red lines show the double-exponential fit (Eq. (5)) of the rising and falling trends in the Fermi-LAT lightcurve (Section 3.1.2).

Fig. 3. MWL lightcurve of B2 1811+31 in a period of approximately 70 days surrounding the Fermi-LAT high-state detection on MJD 59123, marked by the dashed red line. From top to bottom panels: VHE γ-ray flux above 135 GeV from MAGIC, HE γ-ray flux above 100 MeV from Fermi-LAT in daily bins and using the ‘adaptive-binning’ method, X-ray flux in the 0.3−10 keV range from Swift-XRT, optical/UV data in the photometric filters of Swift-UVOT, optical R-band data and radio data. The shaded yellow and light violet bands indicate periods of 48 h and 72 h surrounding the two sets of MAGIC consecutive observation nights, respectively, defined as periods A and B in Section 6.

2.1. MAGIC

The MAGIC telescopes (Aleksić et al. 2016) constitute a stereoscopic system of two 17 m diameter imaging atmospheric Cherenkov telescopes located at the height of about 2200 m a.s.l. at the Observatorio del Roque de los Muchachos (La Palma, Spain). Following the high-state detection from Fermi-LAT in the HE γ-ray band, the MAGIC telescopes performed observations of B2 1811+31 from October 5, 2020 (MJD 59127), up to and including October 11, 2020 (MJD 59133). A total of about 5 hours of good quality data in good weather conditions, dark time and wide zenith range from 20° up to 65° was collected over five observation nights. Table C.1 reports the night-wise time intervals and zenith ranges of the MAGIC observations. Data were analyzed using the MAGIC analysis and reconstruction software MARS (Zanin et al. 2013). We used the standard variable θ², which is defined as the squared angular distance of the reconstructed shower direction with respect to the source direction, to look for any significant VHE γ-ray excess with respect to the background. Primarily, a unique dataset was obtained by combining the data collected during the 5 observation nights. On this dataset, the presence of VHE γ-ray emission from B2 1811+31 was established with a statistical significance of 5.3σ. The statistical significance was estimated using the Li&Ma formula reported in Li & Ma (1983).

We then derived the night-wise VHE γ-ray flux for energies above 135 GeV. The energy threshold was optimized in order to achieve a proper flux estimation for each night regardless of the observational conditions, such as zenith range, weather, night sky background level. The MAGIC lightcurve at E>135 GeV is reported in Table C.1, along with the night-wise significances of the VHE γ-ray signal from B2 1811+31. The lightcurve is shown in the top panel of Figure 3, where the 95% confidence level upper limits are indicated as downward arrows in VHE γ rays when the significance is below 3σ. Since the flux levels in each observation nights are compatible within the 1σ statistical uncertainty band, we conclude that no significant variability is seen in the VHE lightcurve. The weakness of the signal prevented any further quest for intra-night variability. For reasons that are discussed in Section 6.1, the MAGIC dataset was divided into two periods, one including the observations carried out on October 5 and October 6 (period A), and the second one from October 9 up to and including October 11 (period B). For each period, we evaluated the overall spectrum combining the data taken in the observations within the same period and fitted it with a power-law (PL) function

$$\frac{\mathrm{d}N}{\mathrm{d}E}=N_0{\left(\frac{E}{E_0}\right)}^{-{\mathrm{\Gamma }}_{\text{PL}}},$$(1)

with photon index Γ_PL, normalization constant N₀ and decorrelation energy E₀. In order to reconstruct the intrinsic spectra of the source in the two periods, the observed spectra were unfolded by the energy dispersion using the Tikhonov method (Albert et al. 2007) and then corrected for γ-ray absorption by the interaction with the extra-galactic background light (EBL) using the Domínguez et al. (2011) model. The intrinsic spectra are soft with Γ_PL = 4.16±0.63_stat, E₀ = 130.95 GeV and N₀=(4.21±1.51_stat)×10⁻¹⁰ TeV⁻¹ cm⁻² s⁻¹ for period A and Γ_PL = 3.75±0.40_stat, E₀ = 125.16 GeV and N₀=(7.36±1.99_stat)×10⁻¹⁰ TeV⁻¹ cm⁻² s⁻¹ for period B (Table C.4).

2.2. Fermi-LAT

The Large Area Telescope (LAT) on board the Fermi satellite is an imaging, wide field-of-view (∼2.4sr at about 1 GeV), pair conversion γ-ray instrument sensitive to photons from 30 MeV to about 1 TeV (Atwood et al. 2009). The Fermi-LAT observed a hard-spectrum flare at giga-electronvolt energies from the source on October 1, 2020 (MJD 59123). The daily γ-ray flux at E>100 MeV was found to be higher by a factor 11 than the average value reported in the 4FGL, with a daily photon index of 1.4 ± 0.2, harder than the 2.14 ± 0.6 value reported in the 4FGL. In addition, Fermi-LAT detected several E>10 GeV photons, the most energetic one having energy of 61 GeV, with probability >99% of association with B2 1811+31 (Angioni et al. 2020).

The analysis of Fermi-LAT data was performed using the ScienceTools v2.0.8 and the fermipy³ v1.0.1 Python package. Photon events from the Pass 8 P8R3 (Atwood et al. 2013; Bruel et al. 2018) SOURCE class, with reconstructed energy in the 100 MeV−1 TeV range and direction within 15° from the nominal position of the source were selected, applying standard quality cuts (‘DATA_QUAL>0 && LAT_CONFIG= = 1’). The P8R3_SOURCE_V3 instrument response functions were adopted. To reduce the contamination from the Earth limb, a zenith angle cut of 90° was used. All the localized sources included in the Third Data Release of the 4FGL (4FGL-DR3, Abdollahi et al. 2022) with position within 20° from the source were included in the model, along with the Galactic interstellar diffuse and residual background isotropic emission, as modeled, respectively, in gll_iem_v07.fits and iso_P8R3_SOURCE_V3_v1.txt⁴.

We evaluated the significance of the γ-ray signal from the source employing a maximum-likelihood test statistic (TS) defined as $\text{TS}=2\left(log{L}_1-log{L}_0\right)$, where L is the likelihood of the data given the model with (L₁) or without (L₀) a point-like source at the position of B2 1811+31.

Dedicated spectral analyses were performed in time windows corresponding to the low and high states of the source. Table 2 provides the definition of the “Flare” period, corresponding to the 2020 γ-ray high state, as well as the definitions of the “Pre-flare” and “Post-flare” time intervals. The Pre-flare period lasts from August 2008, the beginning of Fermi-LAT data taking, up to the beginning of the 2020 outburst. The Post-flare period extends from the end of the Flare period up to January 2023. The temporal boundaries of these time intervals were evaluated from the long-term Fermi-LAT lightcurve as discussed later in this section. In the top panel of Figure 1, the three periods are indicated by the shaded areas. In each of these time intervals, a binned maximum-likelihood analysis was performed dividing the energy range in 8 bins per decade and the angular coordinates’ space in pixels with size of 0.1°. In the fit, the normalizations of the diffuse components and of the significantly detected (TS > 25) sources within 5° from B2 1811+31 were left free. The spectral parameters of the other 4FGL sources were fixed to the published 4FGL-DR3 values. The HE γ-ray spectrum of the source was primarily modeled with a PL function (Eq. (1)) with free normalization and spectral index. The resulting SEDs in the Pre-flare, Flare and Post-flare periods are shown in Figure 4, whereas the PL best-fit parameters are reported in Table 2.

Fig. 4. SED in the 100 MeV−1 TeV energy range resulting from Fermi-LAT data from the “Pre-flare” (red), “Flare” (blue) and “Post-flare” (green) periods, as defined in Table 2. The average spectral indices and fluxes in the same energy range for the three periods are reported in Table 2.

During the high state, on average the source experienced strong flux enhancement in the HE γ-ray band, by a factor of ∼6 compared to the average radiative states before and after the 2020 flare. In addition, strong spectral hardening occurs, as the photon index changes from 2.11±0.03 in low state to 1.83±0.02 in high state. These evidence indicate that the spectral break of the HE bump shifted to higher energies with respect to the quiescent state. We provide a qualitative interpretation of this evidence in Section 4, in view of the classification of the source broadband emission during the high and low states.

To test the curvature of the HE γ-ray spectra during the three periods, maximum-likelihood fits were performed assuming log-parabola (LP), power-law with exponential cutoff (PLEC) and broken power-law (BPL) models⁵. The TS_curv test statistic, defined as ${\text{TS}}_{\text{curv}}=2\left(log{L}_{\text{curved model}}-log{L}_{\text{PL}}\right)$ as in 4FGL, was employed. Results are summarized in Table C.2. No significant improvement was found up to the 3σ significance level⁶. The spectral break of the HE bump of 2020 flare SED is reasonably lying around tens of giga-electronvolts, where the uncertainties become relevant given the Fermi-LAT sensitivity in this energy range. Nevertheless, the MAGIC data presented in Section 2.1 are able to characterize the SED above the spectral break, therefore we postpone a detailed investigation of the break point of the HE SED bump to Section 6, where the modeling of the high-state broadband emission is presented.

To investigate the long-term variability of B2 1811+31 in the HE γ-ray band, the lightcurve from August 2008 to January 2023 was evaluated in bins of 30 days by performing dedicated maximum-likelihood fits in each time bin. In each time bin, the fit was performed leaving free the normalization and photon index of the B2 1811+31 PL spectral model, whereas all the spectral parameters of the other sources included in the model were fixed to the best-fit values obtained from the analysis of Fermi-LAT data integrated from August 2008 to January 2023. The lightcurve is shown in Figure 1. Upper limits at a 95% confidence level are shown as downward arrows for each time bin in which the TS of the source was found to be smaller than 4, which corresponds approximately to the 2σ significance level. As a sanity check, we verified that the flux to flux error ratio is approximately proportional to the $\sqrt{\mathit{TS}}$ for all time bins, as prescribed in Abdollahi et al. (2023). We repeated the same sanity check also for the other lightcurves derived from Fermi-LAT data throughout this work. The long-term lightcurve shows that the source experienced a high-state period starting several months before the Fermi-LAT high-state detection on MJD 59123. In Table 2, we report the beginning of the Flare period to be MJD 58940. This date approximately corresponds to the first bin of the monthly Fermi-LAT lightcurve showing at least a two-fold flux increase with respect to the average flux in the previous bins. Analogously, the end of the Flare period was set to MJD 59190, since the flux drops significantly after that time. No such high states of B2 1811+31 were recorded at γ rays either between 2008 and the 2020 outburst or after that event. In view of the discussion on the multiband correlations (Section 3.3), it is worth to highlight the following. While most of the monthly bins during the 2010−2012 period resulted in upper limits, as the time of the 2020 γ-ray flare got closer, the upper limits become less frequent and the average reconstructed flux slowly increases, up to the 2020 high state. This indication of evolution on timescales of years can be traced in connection with the evolution of the emission in the other wavebands (Section 2.4).

The fast variability of the HE γ-ray flux during the 2020 flare was investigated by computing the weekly binned lightcurve during the Flare period (Figure 2). The variability analyses performed on this lightcurve, described in Section 3.1.2, were used to identify a period of approximately 70 days surrounding the MAGIC observations, corresponding approximately to the double-peak structure at MJD 59110−59160 in Figure 2. In this period, the daily Fermi-LAT lightcurve was computed. The time bins of the daily lightcurve were determined by requiring that the bins in which the MAGIC observations fell were approximately centered on the MAGIC observing times, reported in Table C.1. Additionally, we constructed lightcurves with adjustable-width time bins using the ‘adaptive-binning’ method (Lott et al. 2012). This method naturally adjusts the bin widths by requiring a constant relative flux uncertainty, σ_F/F, so that low-activity periods are encapsulated within longer bins and high-activity periods are divided into shorter bins, without favoring any a priori arbitrary timescale. Both the daily lightcurve and the adaptive-binning with σ_F/F≤0.3 lightcurve are shown in Figure 3. With respect to the average flux above 100 MeV in the quiescent state, (1.7 ± 0.1) ×10⁻⁸ cm⁻² s⁻¹ (Table 2), the daily flux on the day of the Fermi-LAT high-state observation is (3.0 ± 1.0) ×10⁻⁷ cm⁻² s⁻¹, hence indicating a daily flux increase by a factor 18±6 in the HE γ-ray band. The fast variability in the daily Fermi-LAT lightcurve can be exploited to constrain the size of the emission region responsible for the γ-ray flare, as discussed in Section 3.1.1.

2.3. Swift

The Neil Gehrels Swift satellite (Gehrels et al. 2004) is a rapidly slewing, multiwavelength observatory primarily designed to observe gamma-ray bursts and their X-ray and optical/UV afterglows with its three instruments on board: the Burst Alert Telescope (BAT, Barthelmy 2005) sensitive in the range 15−350 keV, the X-Ray Telescope (XRT, Burrows et al. 2005), sensitive to 0.2−10 keV photons, and the UV/Optical Telescope (UVOT, Roming et al. 2005) covering the 190−600 nm range. Data from Swift-BAT were not included in this work as the hard X-ray fluxes from blazars are generally below the BAT sensitivity for the typical exposure times on the order of several ks. Standard Swift data-taking consists in simultaneous observations with its three telescopes.

Besides being part of a long-term monitoring by Swift, B2 1811+31 was observed 8 times in 2020 as part of the MWL campaign, gathering total exposure times of 14.1 ks for both Swift-XRT and Swift-UVOT. Additionally, the source lies in the field-of-view of the SAFE POINTING 3 observing mode. The Swift safe pointings are predetermined locations on the sky that are observationally safe for UVOT and to which the spacecraft points when observing constraints prevent observations of automated or pre-planned targets (Gehrels et al. 2004). When the spacecraft operates in safe pointing mode, only data from XRT and BAT are taken, whereas the exposure of UVOT is zero. The first observations date back to 2005, whereas the latest ones we included were acquired in February 2024, allowing us to characterize in detail the long-term optical-to-X-ray SED evolution.

2.3.1. Swift-XRT

The reduced XRT data products used for this study were downloaded from the UK Swift Science Center⁷ (Evans et al. 2009). Calibrated source and background files, along with the appropriate response files, were used as inputs to the spectral fitting package XSPEC v12.13.1 (Arnaud 1996) through PyXspec v2.1.2.

The spectra in the 0.3−10 keV band were fitted with the tbabs photon absorption model implemented in XSPEC folded with a PL (Eq. (1)) and a BPL model of the form

$$\frac{\mathrm{d}N}{\mathrm{d}E}=\{\begin{array}{cc}K{\left(\frac{E}{E_0}\right)}^{-{\mathrm{\Gamma }}_{1,\text{BPL}}}\hfill & E<E_{\text{break}}\hfill \\ K{E}_{\text{break}}^{{\mathrm{\Gamma }}_{2,\text{BPL}}-{\mathrm{\Gamma }}_{1,\text{BPL}}}{\left(\frac{E}{E_0}\right)}^{-{\mathrm{\Gamma }}_{2,\text{BPL}}}\hfill & E\ge E_{\text{break}}\hfill \end{array}.$$(2)

In the fit, the neutral hydrogen (HI) column density N_H was fixed to the value reported in Willingale et al. (2013), N_H = 6.09×10²⁰ atoms/cm², and we adopted the molecular abundancies reported in Wilms et al. (2000). The Cash-statistic (Cash 1979) was employed in the fit procedures, the free parameters being the normalization and the spectral index (indices, in case of BPL model) of the source model. The XSPEC implementation of the F-test⁸ was used to test whether the BPL model was preferred over the PL one on a statistical basis, the cut being fixed to 3σ significance level.

The results of the spectral analyses of the B2 1811+31 Swift-XRT observations taken during the 2020 γ-ray high-state period are reported in Table C.6. In 2 out of the 8 observations during the 2020 high state, the BPL model was preferred over the PL model for describing the X-ray spectrum of the source, with a spectral break at around 2−3 keV. The X-ray spectral points were corrected for the ISM extinction effects, in order to reconstruct the intrinsic X-ray SED. The bin-by-bin correction factor was constructed by integrating the e^{−σ_ISM(E)N_H} absorption term over each energy bin, using the HI column density value previously reported, where σ_ISM(E) is the total photo-ionization cross-section of the ISM, normalized to the total hydrogen number density N_H along the line-of-sight, in atoms cm⁻².

During the MAGIC observation period, two B2 1811+31 observations were carried out by Swift-XRT, in MJD 59128 and MJD 59132. The spectra acquired in these observations yield X-ray fluxes of (11.0 ± 0.6) and (11.4 ± 0.6), in units of 10⁻¹² erg cm⁻² s⁻¹, and PL spectral indices of 2.5±0.1 and 2.6±0.1, respectively. In Section 6, we employed these measurements to gather insights into the particle population dominating the X-ray flux from the source.

The long-term evolution of the B2 1811+31 SED in the 0.3−10 keV energy range, corresponding to the 0.7×10¹⁷−2.4×10¹⁸ Hz frequency range, is shown in Figure 5. A variability in the X-ray flux by more than 2 orders of magnitude at 10¹⁸ Hz is observed. Most of the highest-flux SEDs in Figure 5 correspond to observations taken during the 2020 high state. In addition, some high-state X-ray SEDs, shown in red in Figure 5, were computed from observations carried out during an X-ray flare in 2005. The available optical observations in 2005 from CRTS do not indicate an analogous high state in the optical R-band (Figure 1). Figure 5 shows that during the 2020 flare the source exhibited a harder and brighter X-ray state compared to its low state. A similar trend is observed in HE γ rays when comparing the average spectra between high and low states, as shown in Figure 4 and discussed in Section 2.2. The interpretation of this behavior, along with the classification of the high and low states of the source, is discussed in Section 4.

Fig. 5. Long-term evolution of the optical-to-X-ray SED of B2 1811+31 reconstructed from the Swift-UVOT and Swift-XRT observations.

2.3.2. Swift-UVOT

The UVOT data were processed and analyzed using the uvotimsum and uvotsource tasks integrated in the HEASOFT 6.32 software package, along with the Swift/UVOTA Calibration Database (CALDB) files v20211108. UVOT observations were performed using the optical v, b, u and UV w1, m2, w2 photometric broadband filters (Poole et al. 2008; Breeveld et al. 2010). The uvotsource tool performs aperture photometry on a single source and returns its count rates, magnitude and flux density accounting for instrumental effects embedded in the CALDB. Source counts were extracted from a circular region of 5 arcseconds radius centered on the source, while background counts were derived from a circular region of 20 arcseconds radius in a nearby region free of contaminating sources. To correct for Galactic interstellar extinction effects, the standard procedure described in Fitzpatrick (1999) was employed, using E(B−V) = 0.0431±0.0017 (Schlafly & Finkbeiner 2011). The reconstructed magnitudes from the observations performed during the 2020 high state are reported in Table C.7. The evolution of the long-term optical/UV SED is in Figure 5 and it shows that during the 2020 high state the optical/UV flux experienced a two-fold increase with respect to the quiescent state. In Section 4, we discuss the classification of the low and high states based on the selection of simultaneous UVOT and XRT data.

2.4. Optical

In this project, we include optical data from the following surveys, observatories and programs: Zwicky Transient Facility (ZTF, Bellm et al. 2019), Katzman Automatic Imaging Telescope (KAIT, Filippenko et al. 2001) at the Lick Observatory, Catalina Real-Time Transient Survey (CRTS, Drake et al. 2009), Astronomical Observatory of the University of Siena⁹ (0.3 m), Würzburg Observatory¹⁰ (0.5 m) and the Tuorla blazar monitoring program¹¹ (Nilsson et al. 2018). The Tuorla blazar monitoring program employed data from the 0.35 m Kungliga Vetenskapsakademien (KVA) telescope until the end of 2019 and since 2022 it uses data from the Jean Orò Telescope (TJO)¹².

The survey data have been extracted from the databases, while the dedicated observations of this source from Siena Observatory, Würzburg Observatory and Tuorla blazar monitoring were analyzed using standard procedures of differential photometry (e.g., Nilsson et al. 2018). An aperture of 5 arcseconds was used for the photometry.

In order to construct the long-term lightcurve, we combined the survey data and dedicated observations, following the prescriptions of Kouch et al. (2024). The dedicated observations were performed in the R-band, while the survey data were measured using many different filters. In addition, the data from different telescopes can be affected by other instrumental systematic discrepancies (e.g., due to different apertures and different comparison stars). In order to account for these effects, we used offsets of a few mJy, computed employing Tuorla blazar monitoring data as reference. The shifts were determined using data from simultaneous or quasi-simultaneous nights.

The fluxes were corrected for host galaxy contribution, which was estimated to 0.015 mJy within aperture of 5 arcseconds, using the radial profile reported in Nilsson et al. (2003). Galactic extinction is 0.094 mag in R-band and the observed magnitudes were corrected accordingly (Schlafly & Finkbeiner 2011).

The 2020 flare corresponds to a coherent high state in optical, X-rays and γ rays. Moreover, the long-term optical lightcurve shows that the 2020 high state occurred at the apex of the long-term average increasing trend in the optical band lasting for several years, since 2015−2016. After the high-state episode, the optical flux dropped with a decay time significantly shorter with respect to the timescale of the increasing trend leading to the 2020 high state. The γ-ray long-term evolution indicates an analogous long-term rising-peaking-falling trend.

2.5. Radio

We included 15 GHz radio data from the AGN monitoring program of the 40 m Telescope of the Owens Valley Radio Observatory (OVRO) and from the TELAMON program conducted with the Effelsberg 100 m telescope at 14−17 GHz and 19−25 GHz. The TELAMON flux densities are averages over multiple sub-frequencies within each band. For the 14−17 GHz band, the two sub-frequencies are 14.25 and 16.75 GHz, whereas the four sub-frequencies within the 19−25 GHz band are 19.25, 21.15, 22.85 and 24.75 GHz. High-level data were provided by the respective instrument teams. The data analysis chains for the two radio telescopes are described in Richards et al. (2011) and Eppel et al. (2024), respectively, for OVRO and TELAMON.

B2 1811+31 was included in the OVRO monitoring program in 2009, thus providing a long-term radio coverage that allowed us to cross-correlate the long-term radio emission with those in the optical and γ-ray bands as reported in Section 3. Qualitatively, the long-term radio lightcurve follows a different trend from the one displayed by the lightcurves at higher frequencies. In particular, the two radio flares in 2017 and 2021, visible in the bottom panel of Figure 1, have no trivial correspondence at higher frequencies.

3. Variability and correlations

3.1. Variability analysis

The amplitude of the variations in time of the radiative emissions from AGNs can be orders of magnitude greater than in astrophysical sources such as stars and non-active galaxies. The timescales of the fast coherent variations in the blazars emission in a given waveband are usually employed, through causality arguments, to constrain the size R_b of the emission region dominating the radiative output in the same waveband. The observed timescales t_var for variation in the radiation are longer than the light-crossing time, which results in

$$R_{\text{b}}\lesssim \frac{c{\delta }_{\text{D}}{t}_{\text{var}}}{1+z}$$(3)

after accounting for Doppler boosting between the observer frame and the reference frame comoving with the emission region and for cosmological effects through the source redshift z. The relativistic Doppler factor is defined as

$${\delta }_{\text{D}}={\left[\mathrm{\Gamma }(1-\beta cos\theta )\right]}^{-1}$$(4)

and governs the transformation of photon energies between the frames with relative Lorentz factor $\mathrm{\Gamma }={\left(1-{\beta }^2\right)}^{-1/2}$, where βc denotes the velocity of the moving emission region and θ is the angle formed by βc with the line of sight from the observer.

3.1.1. Short-timescale variability

In the literature, constraints on the emission region size have been estimated from the fast variability in many different wavebands, such as X-rays, HE γ rays, and VHE γ rays (e.g., MAGIC Collaboration 2024; Foschini et al. 2013; Aharonian et al. 2007, respectively). In this project, we use the method described in Foschini et al. (2011) to infer the short-timescale variability of the HE γ-ray flux during the γ-ray high-state period and to quantify the significance of the estimate. The method consists in scanning the Fermi-LAT daily lightcurves in Figure 3 to find the minimum doubling/halving time τ, defined from F(t) = F(t₀) 2^{−(t−t₀)/τ}, where t and t₀ are the centers of two consecutive time bins in which the source fluxes F(t) and F(t₀) were reconstructed without ending up into upper limits. The significance of the reconstructed τ was computed as $S=|F(t)-F(t_0)|/\sqrt{{\sigma }_{F(t)}^2+{\sigma }_{F(t_0)}^2}$, where σ_F(t) and ${\sigma }_{F(t_0)}$ are the uncertainties on F(t) and F(t₀), respectively.

The resulting variability timescale with the highest significance levels are reported in Table 3. The MWL follow-up campaign on B2 1811+31 was organized after the Fermi-LAT observation of elevated γ-ray emission from the source. Indeed, the days in which the daily flux yielded the most significant variability timescales closely correspond to those surrounding the Fermi-LAT detection of high state. The resulting variability timescale with the highest significances are $t_{\text{var}}\approx \left(3-6\right)\text{h}$, so the size of the emission region dominating the γ-ray flux with relativistic Doppler factor δ_D has to be smaller than ${\text{R}}_{\text{max}}\approx \left(3-6\right)\times {10}^{14}{\delta }_{\text{D}}\text{cm}$. Variability timescales on the order of few hours are compatible with those found for other VHE blazars in flaring state (e.g., Acciari et al. 2021; Abe et al. 2024a).

3.1.2. Variability in timescales of several days

In this section, we present the analysis of the HE γ-ray flux variability in timescales of several days in a period of approximately 250 days surrounding the Fermi-LAT high-state observation. Indeed, as shown in the long-term MWL lightcurve in Figure 1, B2 1811+31 persisted in a relatively high state (with respect to the average flux several years prior) from April 2020 (MJD≈58940) to December 2020 (MJD≈59190), thus from several months before the Fermi-LAT high-state observation. This period was denoted as “Flare” in Table 1.

Figure 2 shows the Fermi-LAT lightcurve with weekly bins and the optical R-band lightcurve in the Flare period. The Fermi-LAT lightcurve with bins of 7 days was computed using the method described in Section 2.2. Visual inspection of the two lightcurves qualitatively indicates a common trend between the flux evolution in the two energy bands. The cross-correlation between the optical and HE γ-ray fluxes is quantitatively evaluated in Section 3.3. The lightcurves show repeating rising and falling trends. We employed a double-exponential function,

$$F(t)=A·\{\begin{array}{cc}{2}^{\left(\frac{t-t_{\text{peak}}}{t_{\text{rise}}}\right)}\text{if}t\le t_{\text{peak}}\hfill & \hfill \\ 2^{\left(\frac{t_{\text{peak}}-t}{t_{\text{decay}}}\right)}\text{if}t>t_{\text{peak}}\hfill & \hfill \end{array},$$(5)

to fit the peaks of the Fermi-LAT weekly binned lightcurve around the three highest-flux bins, for which the optical R-band lightcurve shows similar trends. In the time interval MJD 58940 − 58970, the fit yielded t_rise = (5.6 ± 1.6) d and t_decay = (4.5 ± 1.6) d for the bump peaking at t_peak = (58962 ± 2) MJD. Analogously, the fit of the lightcurve in the MJD 59020 − 59050 time interval yielded t_rise = (7.4 ± 2.1) d and t_decay = (6.2 ± 2.0) d, with t_peak = (59042 ± 2) MJD. The double-peak structure in the MJD 59110−59160 time interval was fitted with the sum of two double-exponential functions (Eq. (5)). We retrieved estimates of t_rise = (18 ± 4) d for the first peak and t_decay = (11 ± 3) d for the second one, while the other temporal parameters of the two double-exponential function were poorly constrained. Figure 3 shows the MWL lightcurve in the MJD 59110−59160 period. In Section 6 we discuss how these estimates can act as guide for modeling the broadband SED during the 2020 high state.

3.2. Intraband correlations

3.2.1. X-rays (Swift-XRT)

In this section, we discuss the intraband correlation between the flux, F, and photon index, Γ, in the X-ray and HE γ-ray bands. Figure 6 shows the scatter plot of the photon index, Γ_X, and flux, F_X, in the 0.3−10 keV energy range from the Swift-XRT data presented in Section 2.3.1. We primarily tested the correlation between Γ_X and F_X separately during the Pre-flare, Flare and Post-flare periods using the Pearson coefficient r, a statistical measure of linear correlation. The statistical uncertainties were included in the computation of both the correlation coefficient and significance through simulations, as described in Appendix A. We find that during the Pre-flare (Flare) period the Pearson coefficient between Γ_X and log₁₀F_X resulted to be r=−0.75 (r=−0.71), with p−value = 4×10⁻³ (p−value = 2×10⁻²), corresponding to a significance of 2.9σ (2.3σ) for the null hypothesis of non-correlation. Therefore, log₁₀F_X is linearly anti-correlated with Γ_X with a confidence level of 95% in both the Pre-flare and Flare periods. However, the limited number of Swift-XRT observations in the Post-flare period prevents a reliable statistical analysis of the source behavior during this phase. Over the full dataset, which includes data from all periods combined, we find r=−0.77 with p−value = 4×10⁻⁶ (significance of 4.6σ). The trend was fitted with a linear function between Γ and log₁₀F

$$\mathrm{\Gamma }=p_0+p_1{log}_{10}\left(F/F_0\right),$$(6)

with F₀ = 1 erg cm⁻² s⁻¹, yielding best-fit values p₀=−1.1±0.6 and p₁=−0.33±0.05.

Fig. 6. Scatter plot of the intrinsic power-law spectral index ΓX and flux FX in the 0.3−10 keV energy range from the long-term Swift-XRT observations of B2 1811+31. The observations carried out the “Pre-flare”, “Flare” and “Post-flare” periods are marked in red, blue and green, respectively. The orange line marks the fit of the full dataset with Eq. (6).

The analysis reported here indicates that the long-term evolution of the X-ray spectrum of the source follows a harder-when-brighter trend; in other words, higher X-ray fluxes correspond to harder spectra. This trend has been observed in the evolution of the X-ray emission of several IBLs and HBLs (e.g., Mrk 421, Aleksić et al. 2015) and it can be explained in terms of the dynamics of the electron population that is dominating the emission in the X-ray band. Indeed, for HBLs and IBLs with ν_s close to 10¹⁵ Hz (Section 4), the 0.3−10 keV energy band typically hosts the falling edge of the synchrotron SED bump. Therefore, this radiation is emitted by the highest-energy electrons in the jet, above the break of the particle spectrum.

Acceleration mechanisms result in the injection of particles in the emission regions with a hard spectrum, e.g., $\frac{\mathrm{d}N}{\mathrm{d}E}\propto E^{-2}$ from first-order Fermi acceleration. As a consequence, after an increase in the particle injection rate, the X-ray spectrum hardens and the X-ray flux increases. In other words, high X-ray fluxes correspond to hard X-ray spectra.

The highest-energy electrons primarily cool through synchrotron radiation emission. Moreover, synchrotron cooling proceeds faster with increasing particle energy (Section 6). Therefore, after the injection of accelerated particles with a hard spectrum, synchrotron cooling simultaneously softens the high-energy tail of the particle spectrum and lowers the density of particles with sufficient energy to radiate photons in the 0.3−10 keV band (e.g., Tramacere et al. 2022). This implies that low X-ray fluxes correspond to soft X-ray spectra.

3.2.2. High-energy γ rays (Fermi-LAT)

In Section 2.2, we pointed out that during the 2020 high state the source exhibited on average a harder and brighter HE γ-ray state compared to the low states, as shown in Figure 4. Conversely, in this section we provide evidence that the source spectral properties at HE γ rays on daily timescales during the 2020 high state and on monthly timescales during the low states show hints of softer-when-brighter trends instead.

Figure 7 shows the scatter plot of the photon index, Γ_γ, and flux, F_γ, in the 100 MeV−1 TeV energy range from the Fermi-LAT lightcurves with bins of 30 days in the top panel and with daily bins in the bottom panel. Only the temporal bins in which the TS of the source was found to be above 25 were included. The top panel shows that the 2020 high state is well separated in the Γ_γ−F_γ plane from the low states, that is, the Pre-flare and Post-flare periods. Nonetheless, in two temporal bins within the Flare period with flux F_γ≈5.5×10⁻⁸ ph cm⁻² s⁻¹, the source likely transitioned between the high and low states. By employing the Pearson coefficient r, we analyzed the linear correlation between Γ_γ and log₁₀F_γ during the low states using the lightcurve with 30-day bins from the Pre-flare and Post-flare periods (top panel). Conversely, the correlation during the 2020 high state (bottom panel) was investigated using the daily lightcurve in Figure 3. The statistical uncertainties were included in the analysis as described in Appendix A. In addition, in each state the trends were fitted using Eq. (6), with F₀ = 1 ph cm⁻² s⁻¹. The results are reported in Table C.3.

Fig. 7. Scatter plot of the PL spectral index Γγ and flux Fγ in the HE γ-ray data. In the top (bottom) panel, Γγ and Fγ derive from the Fermi-LAT lightcurve with 30-day bins (daily bins) shown in Figure 1 (Figure 3). The dotted lines mark the fits using Eq. (6).

A positive linear correlation with r = 0.41 is found with p−value = 0.033 (significance of 2.1σ) during the Pre-flare, thus indicating with a ∼95% confidence level a softer-when-brighter trend; in other words, higher fluxes correspond to softer spectra. In the Flare and Post-flare periods, respectively, we find positive r values equal to r = 0.32 and r = 0.25, with p-values of 0.10 and 0.32 (1.7σ and 1.0σ). The lower significance values compared to the Pre-flare period are likely due to the lower number of bins and to the higher influence of statistical uncertainties.

Although the correlation is statistically significant above the 2σ level only during the Pre-flare period, it is noteworthy that in all three periods we find hints of positive correlation between the flux F_γ and the photon index Γ_γ. This suggests that the softer-when-brighter trend could be an intrinsic feature of the source spectral evolution at HE γ rays within a given activity state, high or low. However, due to the limited statistics, a definitive conclusion cannot be drawn from the current dataset.

In the following, we report two key aspects regarding the energy range and main cooling mechanism of the electrons responsible for the source HE γ-ray emission. We then discuss the interpretation of the softer-when-brighter trend at HE γ rays, in the context of the processes occurring within the jet.

The energy range of the electron distribution dominating the emission in the HE γ-ray band is different compared to the one responsible for the X-ray emission (Section 3.2.1). Figure 4 shows that the peak energy of the SED high-energy component is likely below 1 GeV during the low states and it shifts at energies above tens of giga-electronvolts during the high state. In both activity states, the emission at energies between 100 MeV and a few giga-electronvolts is mainly produced by particles with energies around or below the spectral break of the particle distribution (Abdo et al. 2010b).

Moreover, for electrons in this energy range, the timescales for the escape from the emission region are shorter than those for cooling through synchrotron radiation emission (Section 6). This implies that synchrotron emission is not the dominant cooling mechanism for electrons in this energy range. Therefore, the spectral properties of the source emission at HE γ rays are primarily determined by the processes of acceleration and escape, rather than from synchrotron cooling (Tramacere et al. 2022).

Although this is not the first indication of softer-when-brighter trends in the spectral evolution at HE γ rays of a blazar (e.g., Prince et al. 2022; Diwan et al. 2023), the interpretation of this effect is still unclear. Theoretical studies (e.g., Böttcher & Baring 2019) indicate that acceleration processes tend to lead to harder-when-brighter trends in HE γ rays, rather than softer-when-brighter ones. One possibility is that softer-when-brighter trends may result when the total HE γ-ray flux of the source is the sum of the radiation from different emission components, e.g., SSC and external Compton (Khatoon et al. 2024), or from multiple emission regions. In the case of B2 1811+31, in Section 6 we investigate the influence of multiple emission regions to the broadband emission of the source during the 2020 high state. However, a quantitative investigation of the softer-when-brighter trends through a time-dependent model accounting for multiple emission regions (e.g., Böttcher & Baring 2019; Khatoon et al. 2024) is beyond the scope of this work.

3.3. Multiband correlations

Correlations between the observed fluxes at different wavelengths have been long employed to infer insights on the dynamical processes acting in AGN, for instance yielding evidence for BLR stratification (Clavel et al. 1991; Peterson et al. 1991). In recent years, multiband cross-correlations of blazar emissions from radio to VHE γ rays were investigated to distinguish the signatures of possible multiple emission components in the jet. These indications are then interpreted as guide for the SED modeling (Lindfors et al. 2016).

For this task, we employed the z-transformed Discrete Correlation Function (zDCF, Alexander 1997), which is optimized for sparse and unevenly sampled astronomical time-series. The zDCF provides a more robust method for evaluating the delay-dependent multiband cross-correlations than the Discrete Correlation Function (DCF, Edelson & Krolik 1988), due to the utilization of equal-population binning and Fisher's z-transform.

The significance of the zDCF is estimated by computing the zDCF of populations of simulated uncorrelated lightcurves sharing the same power spectral density (PSD) as the observed ones. Lightcurves were generated with a Python 3 implementation¹³ of the DELCgen package (Connolly 2016), which allows for the generation of lightcurves using both Timmer & Koenig (1995) (T&K) and Emmanoulopoulos et al. (2013) (EM) methods. We decided to use the latter, as in the EM method the fluxes of the simulated lightcurves could be generated from the probability density function of the fluxes of the observed lightcurve. In the original T&K method, no constraints on the simulated flux distribution are imposed, thus allowing for the generation of lightcurves with bins of negative fluxes. For this reason, in order to produce lightcurves that are both more physical and more resembling the observed ones, we employed the EM method.

We computed the cross-correlations of the long-term HE γ-ray, optical R-band and 15 GHz radio lightcurves. The PSDs of these long-term lightcurves were fitted with power-law functions in the form PSD(ω) = Aω^−β. The best-fit index for the γ-ray lightcurve is β_γ = 1.1, in agreement with Tarnopolski et al. (2020) and references therein, which report typical values between 1.1 and 1.6 for the PSD indices for the Fermi-LAT lightcurves of a sample of 104 blazars. The optical R-band lightcurve yields a best-fit power index of β_R = 1.5, whereas the 15 GHz radio lightcurve yielded β_15 GHz = 2.3. These PSD power-law indices for the optical and radio AGN lightcurves are compatible with those reported in Lindfors et al. (2016) for a sample of 32 VHE γ-ray emitting BL Lac objects. For each long-term cross-correlation to be tested, 10³ pairs of uncorrelated lightcurves were generated with PSD modeled as power-law with the corresponding indices and with the same flux distributions as the observed ones. The results of the zDCF analysis are shown in Figure 8.

Fig. 8. Summary of the zDCF analyses on the long-term B2 1811+31 lightcurves: (a) zDCF of the long-term HE γ-ray and R-band lightcurves, (c) zDCF of the long-term HE γ-ray and 15 GHz lightcurves, (d) zDCF of the long-term R-band and 15 GHz lightcurves. The 1σ and 3σ bands are indicated in dark and light gray, respectively. Panel (b) shows the distribution of the reconstructed zDCF at zero time-lag for 103 simulated uncorrelated HE γ-ray and optical lightcurves, whereas the red line marks the value reconstructed with the corresponding observed B2 1811+31 lightcurves in Figure 1.

Positive correlation with no delay between the long-term γ-ray and R-band lightcurves with a confidence level of 95% is established, as Figures 8a and 8b show. This trend is generally seen for LBLs, IBLs and FSRQs (Lindfors et al. 2016) and it is related to the peak frequency ranges of the two main bumps characterizing the broadband SED. This evidence can be interpreted in terms of the same emission regions contributing significantly to the fluxes in the two energy bands. On the other hand, the radio emission shows neither correlated flares nor correlated long-term trends with the emissions in the optical band and at HE γ rays, as shown in Figures 8c and 8d. The absence of correlated trends between the long-term radio and optical R-band lightcurves is not usual among BL Lacs, as shown in Lindfors et al. (2016). This evidence suggests that the bulk of the radio emission may originate from different zones from those dominating the radiative output in the optical band and at HE γ rays. The zones responsible for the most of the radio emission are expected to be larger, less relativistic and located further along the jet than those responsible for the bulk of the emission at higher frequencies.

4. Source classification

BL Lacs can be empirically classified according to the frequency of the SED low-energy bump into LBLs, IBLs and HBLs, as described in Section 1. During high-state periods, the synchrotron peak frequency can shift significantly (Pian et al. 1998) and the high-state classification can be different from the low-state one (Böttcher 1999). In this section, we deal with the classification of the source during both the low and high radiative states.

B2 1811+31 is classified as IBL in Laurent-Muehleisen et al. (1999). Figure 5 shows the evolution of the optical-to-X-ray SED from 2005 to 2024 using Swift-UVOT and Swift-XRT data. As the variability of the optical-to-X-ray SED is remarkable, we focused on simultaneous observations in order to select source states to be classified. The Swift-UVOT and Swift-XRT observations chosen as characteristic for the B2 1811+31 low state are those carried out on May 13−24, 2011 (MJD 55694−55705). In order to increase the statistics and have more precise flux estimation, we included also the Swift-UVOT and Swift-XRT observations taken in January 2015, as the spectral properties in these observations were fully compatible with the 2011 observations. For the low state, in order to have coverage also in infrared, we included temporally close data taken from the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010) at 3.4, 4.6, 12 and 22 μm wavelengths (W1, W2, W3, W4 filters). The considered observations were performed in March 2010 and the high-level data were downloaded from the SED Builder tool¹⁴ provided by the ASI Space Science Data Center (SSDC). In addition, they connect smoothly to the optical/UV spectral points derived from the Swift-UVOT observations characterizing the source low state, so we decided to include them as characteristic of the low state. For the high state, we selected data from simultaneous Swift-XRT and Swift-UVOT data from the 2020 high state (Tables C.6 and C.7, respectively). The synchrotron peak frequency was estimated by fitting the infrared-to-X-ray SED with a quadratic function in the log₁₀ν−log₁₀νF_ν plane

$$\nu F\left(\nu \right)=f_0{10}^{-b({log}_{10}\left(\nu /{\nu }_{\text{s}}\right){)}^2},$$(7)

where ν_s is the synchrotron peak frequency, f₀ is a normalization constant and b is the coefficient of the second-order term in log₁₀ν. The resulting ν_s values for the low and high states are

1. Low state: log₁₀(ν_s/Hz) = 14.71 ± 0.03.

2. High state: log₁₀(ν_s/Hz) = 15.21 ± 0.23.

The fits of the selected data are shown in Figure 9. The quiescent behavior is compatible with the one reported in Laurent-Muehleisen et al. (1999) for this source. Instead, the high-state classification is borderline between IBL and HBL, as the reconstructed ν_s is slightly above 10¹⁵ Hz. Thus, a significant shift of the synchrotron peak frequency during the 2020 high state with respect to the quiescent state is established.

Fig. 9. Fit of the infrared-to-X-ray SEDs characteristic of the B2 1811+31 low and high states with Eq. (7), in blue and red, respectively. The spectral points in the optical/UV range were reconstructed from Swift-UVOT observations, whereas the X-ray ones derive from Swift-XRT data. Infrared points are WISE data. The shaded areas correspond to the 1σ band fit uncertainties, while the best-fit synchrotron peak frequencies and their uncertainties are marked by the vertical solid and dotted lines.

The blazar 1ES 1215+303, discovered at VHE by MAGIC (Aleksić et al. 2012), has shown several MWL similarities with B2 1811+31, such as the long-term optical-gamma-ray flux increase, presented in Section 2.4, and the shift in the synchrotron peak frequency during high states (Valverde et al. 2020). As B2 1811+31, during low states 1ES 1215+303 is classified as an IBL, whereas it shows HBL behavior during high-state episodes. These similarities indicate common physical conditions in the jets of the two blazars. Moreover, they provide evidence for sources crossing the boundaries of the standard partition of BL Lacs based on the synchrotron peak frequency.

In the following, we provide an interpretation for the shift at higher frequencies of the synchrotron peak during the 2020 high state, as well as for the average harder and brighter X-ray and HE γ-ray spectra of B2 1811+31 compared to the low states. These trends result mainly from the average hardening, brightening and higher maximum energies of the particle spectrum injected in the regions dominating the source emission during the high state. Further contributions can include enhanced magnetic field strengths, higher bulk Lorentz factors and closer alignment with our line of sight (e.g., Raiteri et al. 2017). These would result in enhancing the beaming and Doppler-boosting effects.

5. Redshift estimation

The estimation of the redshift of BL Lac objects is typically a difficult task, as the BL Lac category is characterized by definition by weak spectral lines, if any, in their optical/UV spectra. Nonetheless, a precise measurement of the redshift is a crucial information for γ-ray blazars. Since γ rays above a few hundreds of giga-electronvolts propagating across cosmological distances suffer from extinction by pair-production with ambient background photons composing the EBL. This leads to a redshift-dependent suppression of the observed VHE γ-ray spectrum which becomes more severe for higher energies and redshift values (Dwek & Krennrich 2013). In addition, in BL Lac objects the photon fields external to the jet are expected to be too weak to contribute significantly to the absorption of VHE γ rays produced by the emission regions in the jet.

The determination of the B2 1811+31 redshift, z = 0.117, dates back to Giommi et al. (1991). In Nilsson et al. (2003), it is suggested that the redshift of B2 1811+31 could actually be higher than the one quoted in literature. The argument supported by the authors derives from the analyses of imaging observations carried out with the Nordic Optical Telescope (NOT). The optical images of a sample of BL Lac objects were fitted with a two-component spatial model accounting for the emissions from both the BL Lac core and the host galaxy. For B2 1811+31, an effective half-light angular radius of the host galaxy of θ_eff = 2.9±0.3 arcsec was reported, which corresponds to r_eff = 7.0±0.7 kpc if the source is located at z = 0.117. The corresponding surface brightness of the host galaxy at θ_eff was found equal to μ_eff = 23.26 ± 0.22 mag arcsec⁻². This value is incompatible with the general trend followed by BL Lac objects in the μ_eff−r_eff space, which is the Kormendy relation (Kormendy & Djorgovski 1989; Kormendy et al. 2009) for elliptical galaxies, the typical host galaxies of BL Lac objects. As the authors suggest, there is no clear reason (bad PSF, bright companions nearby...) why the 2D fit of the B2 1811+31 image should be significantly worse. Therefore, by requiring that the radial profile of the source follows the Kormendy relation, an estimate of z = 0.28±0.03 for the source imaging redshift can be derived, based on the photometric detection of the host galaxy reported in Nilsson et al. (2003).

Given the importance of redshift estimates in evaluating properly the blazar dynamics, we applied the method presented in Prandini et al. (2010a) to constrain the redshift value of the source from the simultaneous spectra from MAGIC and Fermi-LAT. The method is based on the assumption that the slope of the EBL-corrected VHE spectrum should not be harder than the one measured by Fermi-LAT at HE γ rays. An iterative procedure in fine steps of redshift from 0.01 to 2.0 was implemented. For each redshift value, the observed MAGIC spectral points were corrected for EBL extinction effects using the Domínguez et al. (2011) model and fitted with a power-law. The higher the redshift, the harder the photon index of the EBL-corrected VHE spectrum. The range of redshift values at which the best-fit power-law index coincides within the statistical uncertainties with the Fermi-LAT data yields an estimate of the upper limit z^* on the source redshift. The empirical linear relation between z^* and z_rec, in other words the estimate of the source true redshift, with the updated parameters in Prandini et al. (2010b), was employed to estimate z_rec. This method has a reported systematic uncertainty of 0.05 on the reconstructed redshift z_rec (Prandini et al. 2010a). These systematic effects are mainly related to the uncertainty on the location of the intrinsic break in the spectrum of the source, as well as to a spread in the intrinsic spectral indices in the VHE band and to the uncertainties derived from the choice of different EBL models. The application of this method to B2 1811+31 simultaneous Fermi-LAT and MAGIC data in MJD 59124−59126 yielded an estimate of z_rec = 0.22±0.14_stat±0.05_syst, which is compatible with both the literature redshift z = 0.117 (Giommi et al. 1991) and the value that can be derived from the Kormendy relation. The major relevance of the statistical uncertainty on this estimate reflects the poor characterization of the spectrum due to the low significance of the detection. Nonetheless, as the inferred redshift value is coherent with the literature one, we used the latter as reference value for the broadband SED modeling.

6. SED modeling

The broadband radiative emission from blazars is commonly interpreted as originated from nonthermal particles accelerated at relativistic speeds in the jet. In the simplest scenario, the emission region is assumed to be a spherical blob populated by electrons spiraling in the comoving magnetic field B. The blob can be imagined as an idealized emission zone which can represent a variety of relativistic plasma configurations within the jet, such as superluminal knots, recollimation shocks, or standing shocks. The bulk relativistic motion of the emission zone is supported by the observation of superluminal motion in radio, optical, and at X-rays along the jet of bright AGN, respectively, at parsec scale in radio and kilo-parsec scale in optical and X-rays (Whitney et al. 1971; Cohen et al. 1971; Biretta et al. 1999; Snios et al. 2019). The radiation observed in the Earth reference frame is highly beamed and Doppler-boosted by the emission zone bulk relativistic motion and by the Earth lying at small view angles, typically up to a few degrees, from the jet axis.

Efficient particle acceleration occurs in blazar jets. Several mechanisms can contribute to the acceleration of particles up to the highest energies, such as magnetic reconnection (Spruit et al. 2001; Giannios & Spruit 2006), first-order Fermi acceleration at relativistic shock fronts such as recollimation or internal shocks (Kirk & Schneider 1987; Heavens & Drury 1988), stochastic (second-order Fermi) acceleration (e.g., Schlickeiser & Dermer 2000), and shear flow acceleration (e.g., Rieger & Mannheim 2002). Several cooling mechanisms, radiative and non-radiative, contribute in parallel to the overall particle cooling (Kardashev 1962). Electrons can radiatively lose energy via the emission of synchrotron photons and IC scattering of various seeds of low-energy photons.

In BL Lacs, the intensity of photon fields external to the jet, such as from the accretion disk or BLR, is thought to be rather minimal with respect to the nonthermal radiation originated from the jet, as manifested by the weakness or absence of lines in the optical/UV spectra. SSC models, in which the electrons upscatter via IC the same synchrotron photons that they have radiated, are typically employed for describing the SEDs of BL Lacs (Böttcher et al. 2013). In the specific context of HBLs and IBLs, SSC models are generally preferred. However, it has been shown in Acciari et al. (2020b) that one-zone SSC models can lead to far from energy equipartition solutions requiring extreme model parameters. Two-zone SSC models are often successful in reproducing the broadband SEDs of IBLs while providing more physical solutions than one-zone models, e.g., in the cases of S5 0716+714 (MAGIC Collaboration 2018) and VER J0521+211 (MAGIC Collaboration 2025). For a few IBLs, such as W Comae (Acciari et al. 2009), 3C 66A (Abdo et al. 2010) and PKS 0903-57 (Shah et al. 2021), one-zone solutions including IC scattering of external photon fields have been preferred over one-zone SSC solutions. Conversely, for LBLs SSC models fail more frequently. Indeed, the SEDs of LBLs typically show high values of Compton dominance; that is, the ratio between the luminosities at the peak frequencies of the high-energy and low-energy bumps (e.g., Finke 2013). This feature is hardly described by pure SSC models and a substantial external Compton component is often required, as in the case of OT 081 (Abe et al. 2024b).

Alternatively, hadronic models have been employed for modeling the radiative states of several blazars (e.g., Zech et al. 2017; Cerruti et al. 2017). In recent years, hadronic models found increasing interest from the community thanks to the results from the IceCube Neutrino Observatory. In the jets of blazars, neutrinos can be produced from the decay of mesons created in the interactions of relativistic protons with low-energy photon fields or with nuclear targets. IceCube reported the indication of a link between a neutrino track event and the blazar TXS 0506+056 during a high state in γ rays (Aartsen et al. 2018b, a) and the unveiling of several hot spots in the neutrino sky coincident with the directions of active galaxies (Aartsen et al. 2020; Abbasi et al. 2021). The classification of the blazar TXS 0506+056 is not firmly established, although its synchrotron peak frequency would locate it into the IBL/HBL transition region (ν_synch≲10¹⁵ Hz, Aartsen et al. 2018a). In Padovani et al. (2019), it was pointed out that TXS 0506+056 belongs to the category of masquerading BL Lacs. These sources are characterized by luminosities as high as those observed for FSRQs, although showing synchrotron peak frequency values typical of IBLs and HBLs. In recent years, this category of blazars garnered increasing attention from the community due to their potential as neutrino source candidates (e.g., de Clairfontaine et al. 2023; Rodrigues et al. 2024). However, currently B2 1811+31 is not associated with IceCube neutrino events (e.g., Abbasi et al. 2023; Giommi et al. 2020), thus we investigated a purely leptonic scenario for modeling the radiative state of the source during the 2020 flare.

As the source showed fast variability, care was given to properly select data representative of compatible source states; in other words, observations that show no significant evolution in the spectrum. Section 6.1 is devoted to the details of the data selection for the broadband emission modeling.

6.1. Data selection for the broadband SED modeling

We divided the B2 1811+31 MAGIC observations into two datasets based on the results of dedicated MWL analyses in the periods surrounding the MAGIC observations reported in Table C.1. The first dataset is composed of the first two observations, carried out on approximately MJD 59126.9 and MJD 59127.9; that is, October 5–6, 2020. The following three MAGIC observations, which constitute the second dataset, were carried out from around MJD 591230.9 to MJD 59133.9; that is, October 9–11, 2020. In Figure 3, the yellow vertical band marks a region of 48 hours centered on the first two MAGIC observations times (period A), extending from MJD 59126.4 to MJD 59128.4, whereas the light violet band indicates a period of 72 hours centered on the following three MAGIC observations (period B), extending from MJD 59130.4 to MJD 59133.4.

In order to look for differences in the source states during the two periods, dedicated Fermi-LAT analyses were carried out using the same prescriptions presented in Section 2.2. The results are summarized in Table C.5 and shown in Figure 10. The TS for the detection of the source in the two periods were TS_A = 24 and TS_B = 37. The reconstructed HE γ-ray spectra during the A and B periods are significantly different, providing evidence for different source states in the two periods. The HE γ-ray spectrum is harder in state B than in state A, the resulting PL spectral indices are Γ_A = 2.5 ± 0.3 and Γ_B = 1.7 ± 0.2. The reconstructed fluxes in the energy range above 100 MeV during the two states are different as well: 21 ± 7 and 8 ± 4, in units of 10⁻⁸ ph cm⁻² s⁻¹. We also performed analyses of Fermi-LAT data over periods of 24 hours centered on each MAGIC observation night. Although these short-duration analyses employ datasets with poorer statistics with respect to analyses on longer integration times, they followed the general trend of showing harder spectra in period B than in A. It is worth pointing out that these results follow the hint of softer-when-brighter trend within the 2020 high state described in Section 3.2.2.

Fig. 10. Spectral energy distribution at HE and VHE γ rays of B2 1811+31 during periods A and B. For the two periods, spectral points are shown in brown and violet and the 1σ uncertainty bands are in shaded yellow and light violet, respectively. The points with full round markers are derived from Fermi-LAT data. For the two periods, the best-fit parameters of the spectrum reconstructed from Fermi-LAT observations are reported in Table C.5. The points with full square markers show the MAGIC EBL-corrected SEDs. The observed MAGIC spectral points are indicated with hollow square markers and dashed error bars.

In Figure 10, the EBL-corrected MAGIC spectra derived from the observations in periods A and B are shown (Section 2.1). While in period B there is a smooth connection between the Fermi-LAT and MAGIC spectra, in period A the two spectra do not seem to be easily compatible. As reported in Section 3, during the high state the source shows variability in timescales of a few hours. Thus, this effect might be due to fast variations in the γ-ray spectrum of the source that may have been observed by MAGIC and that were smoothed out in the analysis of Fermi-LAT data integrated over 24 and 48 hours. The B2 1811+31 SED in state B was constructed using data simultaneous or quasi-simultaneous with the MAGIC observations. We merged the data from the three MAGIC observations within state B and included in the SED to be modeled the resulting average EBL-corrected source spectrum presented in Section 2.1.

We did not perform the modeling of the source state in period A as it is poorly constrained in both the HE and VHE γ-ray bands. Therefore we performed the modeling of a single snapshot of the source broadband SED in flaring state, in particular the one in period B, rather than a time-dependent modeling.

6.2. Broadband SED modeling

We employed the agnpy (Nigro et al. 2022) open-source code for computing the synchrotron and SSC SEDs for modeling the emission of the source during the high state. The models were investigated using a ‘fit-by-eye’ strategy. As discussed previously, the broadband SED of BL Lacs is generally modeled with SSC models due to the low intensity of photon fields external to the jet.

6.2.1. One-zone modeling of the SED

We first adopted a one-zone SSC scenario, with a unique spherical emission region being responsible for most of the radiative output. The emission region in blazars is moving with relativistic bulk Lorentz factor Γ along an axis forming a viewing angle θ with the line of sight from Earth. It is assumed to be populated by relativistic electrons with electron energy density (EED), $n_{\text{e}}\left(\gamma \right)$, with γ being the Lorentz factors in the frame comoving with the emission region. We employed a broken power-law EED defined as

$$n_{\text{e}}\left(\gamma \right)=k_{\text{e}}\left[{\left(\frac{\gamma }{{\gamma }_{\text{b}}}\right)}^{-p_1}H\left(\gamma ;{\gamma }_{\text{min}},{\gamma }_{\text{b}}\right)+{\left(\frac{\gamma }{{\gamma }_{\text{b}}}\right)}^{-p_2}H\left(\gamma ;{\gamma }_{\text{b}},{\gamma }_{\text{max}}\right)\right]$$(8)

between γ_min and γ_max, with break at γ_b, spectral indices p₁ and p₂ below and above the break, respectively. The H function is the Heaviside function, defined such that $H\left(x;a,b\right)=1$ if a≤x≤b and $H\left(x;a,b\right)=0$ otherwise. k_e is the EED value at γ_b, in cm⁻³, and the electron density is defined as $K_{\text{e}}={\int }_{{\gamma }_{\text{min}}}^{{\gamma }_{\text{max}}}\mathrm{d}\gamma n_e\left(\gamma \right)$. The magnetic field tangled to the region is B. The energy densities in electrons and magnetic field of the emission region are $u_{\text{e}}=m_{\text{e}}{c}^2{\int }_{{\gamma }_{\text{min}}}^{{\gamma }_{\text{max}}}\mathrm{d}\gamma \gamma n_e\left(\gamma \right)$ and $u_B=\frac{B^2}{8\pi }$ in Gaussian cgs units, respectively. The energy equipartition parameter, k_eq, is defined by $k_{\text{eq}}=\frac{u_e}{u_B}$.

The MWL spectral analyses presented in Section 2 can be employed to estimate the EED spectral indices p₁ and p₂. Synchrotron radiation theory (e.g., Dermer & Menon 2010) predicts that a power-law spectrum of electrons, with EED $n_{\text{e}}\left(\gamma \right)\propto {\gamma }^{-p}$, radiates a synchrotron SED $\nu F_{\nu ,\text{synch}}\propto {\nu }^{\frac{3-p}{2}}$. As shown in Section 4, the X-ray spectrum of B2 1811+31 lies deep in the falling edge of the synchrotron bump and it is produced by the electrons with γ≫γ_b, where $n_{\text{e}}\left(\gamma \right)\propto {\gamma }^{-p_2}$. In period B, the X-ray flux is well described with a PL (Eq. 1) with spectral index Γ_X = 2.58±0.10 (Table C.6). Hence, we derive p₂ = 2Γ_X−1≈4.2±0.2. Conversely, the Swift-UVOT 2020 high-state data yield that νF_ν∝ν^0.25±0.10 in the optical/UV band. This would imply a spectral index of 2.5±0.2 for the electron population dominating the radiative output in this band. However, as shown in Section 4, the high-state synchrotron peak lies close to the optical/UV band, hence the electrons contributing to the flux in this band have Lorentz factors close to γ_b, where the EED spectral index changes. Thus, for the B2 1811+31 high state, the argument to estimate p₁ from the optical/UV spectral index cannot be applied safely.

The Doppler factor (Eq. (4)) is a crucial parameter for the SSC models, since it strongly affects the properties of the observed radiation ($\nu F_{\nu }\propto {\delta }_{\text{D}}^4$). In literature, standard values commonly used for blazars single-zone models are δ_D≈10−20 (e.g., Hovatta & Valtaoja 2009). In the construction of the one-zone model solution, we tried several values for the Doppler factor δ_D, up to values around 30. While keeping fixed the δ_D value, the other model parameters were tuned to account for the observed SED and for the constraint deriving from the observed variability at HE γ rays, presented in Section 3.1.1. A single-zone SSC solution with δ_D = 20 is shown in the left panel of Figure 11a. The parameters of this solution are reported in Table 4. In Section 3.1.1, we derived an upper limit on the emission region size of ${\text{R}}_{\text{max}}\approx \left(6-12\right)\times {10}^{15}\text{cm}$ for δ_D = 20. The solution in Figure 11a has δ_D = 20, R_b = 10.6×10¹⁵ cm and magnetic field B = 0.13 G. In the one-zone solutions obtained with δ_D fixed to values lower than 20, higher values of both the region size and the magnetic field were required to reproduce the spectra in all the energy bands. These solutions violated the constraint imposed by the observed variability at HE γ rays and showed a worse agreement with the observed SED with respect to the solution with δ_D = 20 shown in Figure 11a.

Fig. 11. Leptonic modeling of the broadband SED of B2 1811+31 during period B within the 2020 high state in VHE γ rays, marked by the red points. Archival data are indicated by the hollow gray points. The shaded blue and orange areas indicate the 1σ uncertainty bands on the intrinsic spectrum of the source derived from Fermi-LAT and MAGIC observations, respectively. The two insets show a zoom-in on the optical/UV and X-ray bands, respectively. For the two-zone SSC model, the contributions from the blob and core are indicated, respectively, by the orange and green lines, whereas the total is computed as the sum of the two. (a) One-zone SSC modeling. (b) Two-zone SSC modeling.

As is discussed in Tavecchio & Ghisellini (2016), the magnetic field strengths for single-zone models tend to be significantly lower than the values required for equipartition. The one-zone SSC scenario presented in Figure 11a has k_eq = 31.3. In literature, values on the order of 100 and higher are usually found for single-zone models (e.g., Acciari et al. 2020b). Although the reported solution has lower k_eq, it is still far from equipartition (k_eq = 1). To reproduce the optical-to-X-ray spectrum, the break between the spectral indices p₁ and p₂ is required to be quite extreme, with p₁ = 2.0 and p₂ = 4.4. The resulting p₂ value is compatible within the uncertainties with the estimate p₂≈4.2±0.2 derived from the X-ray spectral index. Such an extreme EED spectral break can be a consequence of describing the broadband SED with a model composed by a single emission region. The presented one-zone SSC solution underestimates the radio flux. Single zones employed for modeling high-state SEDs are known to hardly manage to account for the radio emission (e.g., Acciari et al. 2020b). This is due to the dominance of synchrotron-self-absorption effects at radio frequencies for compact and energetic emission regions.

6.2.2. Two-zone modeling of the SED

We investigated a two-zone SSC model in line with the one presented in Tavecchio et al. (2011), with a few modifications. Due to the BL Lac nature of the source, we neglect the IC scattering contribution from photon seeds external to the jet. Given the MWL properties presented in the previous sections, we propose a physical scenario as follows. A small and energetic region, denoted as ‘blob’, is required to dominate the emissions in the X-ray and VHE γ-ray bands. The second zone, denoted as ‘core’, is required to be larger and less energetic (lower bulk Lorentz factor and EED at lower γ values) than the blob. This larger region dominates the optical/UV emission and yields a significant contribution to the HE γ-ray band, as suggested by the observed correlation between the long-term fluxes in these energy bands, presented in Section 3.3. The blob region is considered to be located closer to the central engine than the core region, which would account for the blob region being more energetic than the core. Moreover, as the magnetic field strength along the relativistic jet scales with the distance from the central engine, the magnetic field tangled to the core is required to be smaller than the blob one. As the core region is assumed to be located further along the jet than the blob one, the interplay between the two regions is expected to be negligible: the synchrotron photons of one zone do not serve significantly as seed photons for the IC of the electrons from the other region (Tavecchio et al. 2011). The construction of the two-zone model solution was performed following a similar procedure to that employed for the one-zone model solution. Several values for the Doppler factors δ_D of two regions were tried. For each value of δ_D, the parameter space was explored accounting for the constraints resulting from the MWL analyses presented in the Section 3. As is described later in this section, in the parameter space scan, we also investigated the possibility to account for the observed broadband SED in a self-consistent scenario; that is, with emission regions in which the particle distributions are compatible with those expected at equilibrium under the effects of escape from the emission region and radiative cooling via the emission of synchrotron and SSC radiation (Inoue & Takahara 1996; Cerruti 2020). The resulting two-zone SSC solution is presented in Figure 11b, the parameters are reported in Table 4. This solution characterizes the optical-to-X-ray SED more naturally than the single-zone one (Appendix B), without invoking hard spectral breaks in the distributions of the electrons in both regions. The p_{1, core} and p_{2, core} parameters of the core region are 1.8 and 2.7, respectively, with a spectral break close to 1. The spectral break of the blob region EED also mitigates, with p_{1, blob} = 2.0 and p_{2, blob} = 3.8. As discussed above, the p₂ parameter of the region dominating the X-ray flux is required to be soft by the observed X-ray spectral index. The relativistic Doppler factors of the two regions are δ_{D, blob} = 10 and δ_{D, core} = 4. For δ_D = 10, the variability analysis reported in Section 3.1.1 yields an upper limit of ${\text{R}}_{\text{max}}\approx \left(3-6\right)\times {10}^{15}\text{cm}$ on the size of the zone dominating the γ-ray flux during the high state. The size of the blob region in the presented solution is R_blob = 4.1×10¹⁵ cm, thus it satisfies the constraint set by the causality argument.

In addition, δ_{D, blob} = 10 is consistent with the range of Doppler factor values allowed by the constraint for the region to be transparent to γ rays up to the highest-energy ones observed by MAGIC; that is, E_{ph, max} = 400 GeV in the reference frame of the observer (Figure 10). The blob region dominates the source emission in the X-ray and VHE γ-ray bands. VHE γ-ray photons can interact with soft X-rays to produce e⁺e⁻ pairs, hence suppressing part of the VHE γ-ray flux emitted by the region (e.g., Dermer & Menon 2010). A lower limit to the Doppler factor of the blob can be derived by requiring that the emission region is not opaque to the observed highest-energy γ rays. This lower limit can be estimated as (e.g., Dondi & Ghisellini 1995)

$${\delta }_{\mathrm{D},\hspace{0.5em}\mathrm{blob}}>\sqrt[6]{\frac{{\sigma }_{\text{T}}{d}_{\text{L}}^2{\left(1+z\right)}^2{f}_{\text{X}}{ϵ}_{\text{max}}}{4m_{\text{e}}{c}^4{t}_{\text{var}}}},$$

σ_T being the Thomson cross section, d_L the luminosity distance, f_X the observed X-ray flux, t_var the fast variability timescale and ${ϵ}_{\text{max}}=\frac{E_{\mathrm{ph},\hspace{0.5em}\max }}{m_{\text{e}}{c}^2}$. Using z = 0.117 and f_X≈10⁻¹¹ erg cm⁻²s⁻¹ (Table C.6, Figure 11b), we obtain for t_var = 3−6 h (Table 3) a lower limit of 8−9 to δ_{D, blob}. Hence, the value δ_{D, blob} = 10 would allow the blob region to be significantly transparent to γ rays up to the ones at the highest energies observed by MAGIC.

As far as the variability analyses presented in Section 3.1.2, it is interesting to notice that the light-crossing time argument on t_rise = (18 ± 4) d, observed in correspondence to the VHE detection by MAGIC, yields an estimated upper limit on the emission region size R_max = (1.8 ± 0.4) ×10¹⁷ cm for a Doppler factor δ_D = 4, using Eq. (3). This value is fully compatible with the value R_core = 2.1×10¹⁷ cm employed for the core region. The value of t_rise = (18 ± 4) d was estimated from the variations in the HE γ-ray flux using the Fermi-LAT weekly binned lightcurve. In the two-zone SSC solution shown in Figure 11b, both the core and blob regions contribute significantly to the HE γ-ray total flux. However, the two regions produce variability on very different timescales; ∼10 days for the core and ∼5  hours for the blob. Therefore, the core is likely responsible for the slower variations in the HE γ-ray flux at the timescales from days to weeks, shown in Figure 2, whereas the blob is likely producing the short-timescale variations on top of the long-term ones, shown in Figure 3.

To account properly for both the X-ray and VHE γ-ray spectra, the energy range of the electrons populating the blob region extends up to E_max=m_eγ_maxc²≈150 GeV. Moreover, the core component is able to account more naturally for the radio flux, which was severely underestimated by the one-zone solution.

In the presented solution, the core region is at energy equipartition, whereas the blob equipartition coefficient is improved (k_eq≈18) with respect to the one-zone scenario. This is consistent with Acciari et al. (2020b), in which the authors presented a systematic study on two-component scenarios employed for modeling blazar broadband SEDs and found that two-zone models are systematically closer to equipartition with respect to single-zone ones, but still particle-dominated.

As previously mentioned, in the construction of the one-zone and two-zone model solutions, we investigated the possibility to account for the observed broadband SED with self-consistent particle populations in the emission regions. The temporal evolution of the EED is typically studied through a continuity equation under the effects of particle acceleration/injection, escape, adiabatic expansion and radiative cooling (e.g., Kardashev 1962). Figure 12 shows, for the two regions presented in Figure 11b, the competition among the light-crossing timescale, t_lc, and the main radiative cooling timescales. The light-crossing timescale can be interpreted as the characteristic time for the particle escape from the emission region. In the two-zone solution in Table 4, all particles are ultra-relativistic, hence the light-crossing time is energy-independent, with t_lc∼R_b/c. The timescales for synchrotron energy losses were computed as (Dermer & Menon 2010)

$$t_{\text{synch}}\left(\gamma \right)=\gamma {\left(\frac{d\gamma }{\mathit{dt}}\right)}_{\text{synch}}^{-1}=\frac{3m_{e}c}{4{\sigma }_T{u}_B\gamma }\propto {\gamma }^{-1}.$$(9)

Fig. 12. Competition among the cooling timescales for the two-zone model in Figure 11b, as a function of the electron Lorentz factor γ. The times are in the observer frame; that is, corrected by the corresponding $\frac{1+z}{{\delta }_{\text{D}}}$ with respect to the frames comoving with the two regions. Dashed (dotted) curves refer to the blob (core) region. Blue and red curves indicate the synchrotron and SSC cooling timescales, respectively, whereas the green lines represent the light-crossing times. The orange markers indicate the γb values of the EED within the two emission regions (Table 4).

For each emission region, the SSC cooling timescale corresponds to that of the IC scattering off the (isotropic) synchrotron photon number density $n_{\text{ph}}(ϵ)=\frac{\mathrm{d}N}{\mathrm{d}V\mathrm{d}ϵ}$, with $ϵ=\frac{h\nu }{m_{\text{e}}{c}^2}$ being the dimensionless photon energy and n_ph(ϵ) dV dϵ the number of photons in a volume dV having dimensionless energy between ϵ and ϵ+dϵ. Following Nigro et al. (2022), n_ph(ϵ) of the synchrotron fields in the two regions were computed as

$$n_{\text{ph}}(ϵ)=\frac{1}{m_{\text{e}}{c}^2ϵ}{u}_{\text{ph}}(ϵ)=\frac{3}{4}\frac{1}{m_{\text{e}}{c}^2ϵ}\frac{R_{\text{b}}}{c}\frac{\sqrt{3}{e}^{3}B}{h}{\displaystyle {\int }_1^{\infty }}d\gamma n_{\text{e}}(\gamma )R(x),$$(10)

where u_ph(ϵ) is the photon energy density, the factor 3/4 accounts for averaging the radiation in a sphere and R(x) is the single-electron synchrotron spectral power averaged over the pitch angles (e.g., Eq. 7.45 in Dermer & Menon 2010), expressed as a function of the variable $x=\frac{4\pi ϵm_{\text{e}}^2{c}^3}{3\mathit{eBh}\gamma }$. We computed the cooling timescales of relativistic electrons due to IC scattering with the isotropic synchrotron radiation field, $t_{\text{SSC}}\left(\gamma \right)=\gamma {\left(\frac{d\gamma }{\mathit{dt}}\right)}_{\text{SSC}}^{-1}$, as in Moderski et al. (2005), appendix C. Therefore we accounted for IC scattering energy losses both in Thomson and Klein-Nishina (KN) regimes, in which t_SSC∝γ⁻¹ and t_SSC∝γ/log(γ), respectively. The steady-state solution of the EED continuity equation when synchrotron and IC scattering (Thomson limit only) are the main radiative cooling mechanisms have EED spectral break p₂−p₁ = 1 and γ_b equal to the value at which the escape and radiative cooling timescales coincide (e.g., Inoue & Takahara 1996; Cerruti 2020). We find that, for both regions in Figure 11b, the γ_b values in our solution (Table 4) appear to be equal those obtained comparing the radiative cooling timescales with t_lc (Figure 12). The spectral break in the core region is p_{2, core}−p_{1, core}≈0.9, thus it is consistent with the one expected at equilibrium when t_cool∝γ⁻¹. Conversely, for the blob of the two-zone solution, a stronger break is necessary to account for the observed spectra. This suggests that further effects may be relevant, such as the particle distribution being hardly at equilibrium or possible effects related to cooling in the KN regime. In contrast, in the investigation of the one-zone model presented in Section 6.2.1, we find it is not possible to properly account for the observed SED in all energy bands while requiring that the γ_break is the one expected at equilibrium or requiring that the spectral break is equal to 1. Hence, in Figure 11a, we present the solution which fits better the broadband SED, although the value of γ_break for which the synchrotron and light-crossing timescales coincide is γ_break = 1.3×10⁵. Although the two-zone model solution presented here is not time-dependent, this analysis supports the self-consistency of the proposed physical two-zone model scenario for the MWL emissions of B2 1811+31 during the 2020 high state.

In SSC models, the predicted νF_ν from an emission region depends on the bulk Lorentz factor, Γ, and viewing angle, θ, only through δ_D; that is, the predicted flux does not depend independently on Γ and θ. For an emission region with given δ_D, Eq. (4) implies that $\mathrm{\Gamma }\ge {\mathrm{\Gamma }}_{\text{min}}=\frac{1+{\delta }_{\text{D}}^2}{2{\delta }_{\text{D}}}$ and $\theta \le {\theta }_{\text{max}}=arccos\sqrt{1-{\delta }_{\text{D}}^{-2}}$, with θ_max corresponding to Γ=δ_D. In the presented two-zone solution, the blob region has δ_D = 10, which results in an upper limit of θ_max≈5.7° on the jet viewing angle. Conversely, the jet power depends only on Γ and is independent of θ. For each emission region, from the energy densities u_e and u_B in electrons and magnetic field, the corresponding luminosities can be computed as $L_{\text{e}}=2\pi R_{\text{b}}^2\beta c{\mathrm{\Gamma }}^2{u}_{\text{e}}$ and $L_{\text{B}}=2\pi R_{\text{b}}^2\beta c{\mathrm{\Gamma }}^2{u}_{\text{B}}$, respectively. Thus, for viewing angles between 1° and 5° and for Γ<δ_D, the jet luminosity in particles of the blob (core) region is L_e≈1×10⁴³ erg s⁻¹ (4×10⁴³ erg s⁻¹), whereas the corresponding jet luminosity in magnetic field is scaled down by a factor k_eq = 18.3 (k_eq = 1).

Long-term MWL campaigns for several blazars have shown radio emission correlated with a significant delay, from weeks to years, with respect to the optical and γ-ray emission (e.g., Pushkarev et al. 2010). This evidence is interpreted as related to synchrotron-self-absorption effects and to the emission region becoming transparent in radio later than at higher frequencies (e.g., Blandford & Königl 1979). In the case of B2 1811+31, the 15 GHz long-term lightcurve shown in the bottom panel of Figure 1 reports a high state in 2021, nearly a year after the 2020 VHE outburst. However, the cross-correlation analysis presented in Section 3.3 showed that the radio flux is not significantly correlated with the emissions in the optical and HE γ-ray bands, even when delays up to a few years are taken into account. This suggests that other regions, possibly located further along the jet, may dominate the radio flux from the source. Because of this, in the two-component solution shown in Figure 11b, we left the core region able to account for most of the radio emission (up to ∼50%), but not entirely. For instance, the core would account entirely for the observed radio data if γ_min was lowered from 0.5×10³ to 0.35×10³ and the electron density, K_e, was increased to 1.0 cm⁻³, keeping the other parameters in Table 4 fixed. With such modifications, the agreement between the data and the flux predicted from the two-zone model at higher frequencies would not be compromised and the energy equipartition coefficient of the core region would be k_eq≈1.1. One possibility for the interpretation of the 2021 radio flare is adiabatic expansion of the core region during its propagation along the jet (Tramacere et al. 2022). However, we consider the available data not sufficient to confirm nor exclude that the origin of the 2021 radio flare is related to the 2020 MWL outburst.

The parameter space in SSC models is quite vast and SSC solutions are usually degenerate. The solution shown in Fig. 11b represents one possible set of parameters, but all observational constraints have been taken into account, the model well reproduces the observed SED and the consistency of the proposed physical scenario for the two regions has been investigated.

7. Summary and conclusions

In this paper, we have presented a long-term MWL view of B2 1811+31, with the main emphasis given to the 2020 γ-ray flare of the source. During this high-state period, the MAGIC telescopes reported the first-time detection of VHE γ-ray emission from the source. The MWL observational campaign during the flaring state was organized after the detection of a high state in HE γ rays by the Fermi-LAT and it provided coverage from radio to the VHE band.

The long-term MWL behavior of B2 1811+31 indicates that the optical and γ-ray emissions correlate, whereas the evolution of the radio flux is significantly different from that at higher frequencies. This provides evidence for the same emission regions dominating the radiative output in the optical and HE γ-ray bands. The 2020 flaring state occurred at the apex of an increasing trend in the optical band that persisted for approximately eight years. In addition, after the flare, the source showed a falling trend of the optical flux that led the source to a state compatible to the one before the rising trend. The duration of the falling trend was approximately two years, much shorter than the rising trend. Indications of this long-term rising-peaking-falling evolution can be found in X-rays and HE γ rays as well. The long-term X-ray flux around 1 keV exhibits the harder-when-brighter trend and it shows variations by more than 2 orders of magnitude.

We have confirmed the IBL classification of the source during the low state making use of archival IR data and of dedicated optical/UV and X-ray analyses. During the 2020 high state, the average HE γ-ray spectrum became harder compared to the low states. A similar behavior has been observed in X-rays. Conversely, during different activity periods, we find harder-when-brighter trends in X-rays and a hint of softer-when-brighter trends at HE γ rays. In addition, we claim a significant shift of the synchrotron peak frequency during the 2020 flare, which led the source to a borderline state between the IBL and HBL classes.

Analysis of the fast variability in the HE γ-ray band during the 2020 flare was used to estimate an upper limit on the size of the emission region dominating the γ-ray output. We tested one-zone and two-zone SSC scenarios for modeling the broadband SED probed with the MAGIC telescopes. The widely used one-zone SSC model provides a good estimation of the broadband emission, with some known limitations, such as the underestimation of the radio flux and a hard break in the electrons distribution. In addition, this scenario is far from energy equipartition. For the two-zone model, we investigated a scenario with a small and fast blob dominating the flux at X-rays and VHE γ rays, and a larger and slower blob, denoted as the core, which dominates the optical flux and contributes significantly at HE γ rays. The two-zone model manages to describe the SED in the optical/UV and X-ray bands without invoking extreme parameters for the electrons distribution. The proposed two-zone SSC solution is closer to energy equipartition and accounts more naturally for the radio flux than the one-zone scenario. In addition, the spectral breaks of the EEDs within the two emission regions are comparable to the ones expected from radiative cooling, supporting self-consistency for the proposed model.

Data availability

All data shown in Figures 1, 2, 3, 4, 5, 9, 10 and 11 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/697/A172

Acknowledgments

We would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF, MPG and HGF; the Italian INFN and INAF; the Swiss National Fund SNF; the grants PID2019-104114RB-C31, PID2019-104114RB-C32, PID2019-104114RB-C33, PID2019-105510GB-C31, PID2019-107847RB-C41, PID2019-107847RB-C42, PID2019-107847RB-C44, PID2019-107988GB-C22, PID2022-136828NB-C41, PID2022-137810NB-C22, PID2022-138172NB-C41, PID2022-138172NB-C42, PID2022-138172NB-C43, PID2022-139117NB-C41, PID2022-139117NB-C42, PID2022-139117NB-C43, PID2022-139117NB-C44 funded by the Spanish MCIN/AEI/ 10.13039/501100011033 and “ERDF A way of making Europe”; the Indian Department of Atomic Energy; the Japanese ICRR, the University of Tokyo, JSPS, and MEXT; the Bulgarian Ministry of Education and Science, National RI Roadmap Project DO1-400/18.12.2020 and the Academy of Finland grant nr. 320045 is gratefully acknowledged. This work was also been supported by Centros de Excelencia “Severo Ochoa” y Unidades “María de Maeztu” program of the Spanish MCIN/AEI/ 10.13039/501100011033 (CEX2019-000920-S, CEX2019-000918-M, CEX2021-001131-S) and by the CERCA institution and grants 2021SGR00426 and 2021SGR00773 of the Generalitat de Catalunya; by the Croatian Science Foundation (HrZZ) Project IP-2022-10-4595 and the University of Rijeka Project uniri-prirod-18-48; by the Deutsche Forschungsgemeinschaft (SFB1491) and by the Lamarr-Institute for Machine Learning and Artificial Intelligence; by the Polish Ministry Of Education and Science grant No. 2021/WK/08; and by the Brazilian MCTIC, CNPq and FAPERJ. The Fermi-LAT Collaboration acknowledges support for LAT development, operation and data analysis from NASA and DOE (United States), CEA/Irfu and IN2P3/CNRS (France), ASI and INFN (Italy), MEXT, KEK, and JAXA (Japan), and the K.A. Wallenberg Foundation, the Swedish Research Council and the National Space Board (Sweden). Science analysis support in the operations phase from INAF (Italy) and CNES (France) is also gratefully acknowledged. This work performed in part under DOE Contract DE-AC02-76SF00515. We acknowledge the use of public data from the Swift data archive. This work has made use of data from the Joan Oró Telescope (TJO) of the Montsec Observatory (OdM), which is owned by the Catalan Government and operated by the Institute for Space Studies of Catalonia (IEEC). This research has made use of data from the OVRO 40-m monitoring program (Richards et al. 2011), supported by private funding from the California Institute of Technology and the Max Planck Institute for Radio Astronomy, and by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911. S.K. was funded by the European Union ERC-2022-STG – BOOTES – 101076343. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. This work is based on observations with the 100-m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg. Author Contributions. D. Cerasole: project leadership, Fermi-LAT data analysis, Swift data analysis, multiwavelength data analysis, theoretical modeling and interpretation, paper drafting and editing; L. Di Venere: supervising the Fermi-LAT analysis, paper editing; T. Hovatta: OVRO data analysis; E. Lindfors: optical data analysis and collection, theoretical interpretation, paper drafting and editing; S. Loporchio: MAGIC data analysis, theoretical interpretation, paper drafting and editing; S. Kiehlmann: OVRO data analysis; P. Kouch: optical data analysis and collection; L. Pavletic: MAGIC analysis cross-check; A. C. S. Readhead: OVRO data analysis. The rest of the authors have contributed in one or several of the following ways: design, construction, maintenance and operation of the instrument(s); preparation and/or evaluation of the observation proposals; data acquisition, processing, calibration and/or reduction; production of analysis tools and/or related Monte Carlo simulations; discussion and approval of the contents of the draft.

References

Appendix A: Details of the intraband correlation analysis including statistical uncertainties

In this section, we report details of the intraband correlation analysis between the flux, F, and photon index, Γ, in the X-ray and HE γ-ray bands. The results have been presented in Section 3.2.

Consider a lightcurve of a source in a given energy range. Let this be constituted of n temporal bins, F=(F₁, …, F_n) and Γ=(Γ₁, …, Γ_n) being the arrays of the flux F and photon index Γ in the n bins. Due to the measurement uncertainties, each F_i (Γ_i) is a random variable whose mean value ${\overline{F}}_i$ (${\overline{\mathrm{\Gamma }}}_i$) and standard deviation ${\overline{\sigma }}_{F_i}$ (${\overline{\sigma }}_{{\mathrm{\Gamma }}_i}$) have been estimated from the data. $\overline{F}$ and ${\overline{\sigma }}_F$ ($\overline{\Gamma }$ and ${\overline{\sigma }}_{\mathrm{\Gamma }}$) denote the arrays of the mean values and standard deviations of F (Γ). Let x be defined as $x=\left(F,\Gamma \right)$, with analogous definitions for $\overline{x}$ and ${\overline{\sigma }}_x$. Since the Pearson coefficient r of the correlation between F (or log₁₀F) and Γ is a function of random variables, r is a random variable itself; that is, $r=r\left(x\right)$. The probability density function of r depends on the one of x and on the functional form of r. If the function $r\left(x\right)$ can be well approximated by its first-order Taylor expansion in a region of size ${\overline{\sigma }}_x$ around the mean $\overline{x}$, then the mean value $\overline{r}$ of the distribution of r can be approximated by $r_0\equiv r\left(\overline{x}\right)$; that is, the Pearson coefficient r computed using the mean values $\overline{x}$ (Cowan 1997). If the uncertainties σ_F and σ_Γ are significantly lower than $\overline{F}$ and $\overline{\Gamma }$, up to factors of a few per cent, this condition is satisfied. Conversely, when the statistical uncertainties are not negligible and the level of (anti-)correlation indicated by r₀ is non-zero, the functional form of the Pearson coefficient shows significant non-linearity in a region around $\overline{x}$ of a size comparable to ${\overline{\sigma }}_x$, thus preventing us to employ r₀ to approximate $\overline{r}$ safely. In this case, statistical uncertainties lower the level of (anti-)correlation compared to r₀ and Monte-Carlo simulations have to be used to estimate $\overline{r}$ properly (Cowan 1997). An example of this effect can be shown using the scatter plots of the photon index and flux at HE γ rays in Figure 7. The best-fit values $\overline{F}$ and $\overline{\Gamma }$ would indicate high levels of correlation in the Pre-flare, Flare and Post-flare periods, the r₀ values between log₁₀F and Γ being equal to 0.89, 0.71 and 0.93, respectively. However, if the statistical uncertainties around $\overline{F}$ ($\overline{\Gamma }$) are included by extracting flux (index) values from distributions centered at $\overline{F}$ ($\overline{\Gamma }$) and with standard deviation ${\overline{\sigma }}_F$ (${\overline{\sigma }}_{\mathrm{\Gamma }}$), the simulated flux-index pairs tend to show, on average, higher scattering around the general trend compared to $\overline{F}$ and $\overline{\Gamma }$, therefore indicating a lower correlation compared to r₀.

In the simulations, we assumed that the joint probability density function of the fluxes F and indices Γ can be approximated by independent Gaussian functions centered at $\overline{F}$ and $\overline{\Gamma }$, respectively, with standard deviations given by ${\overline{\sigma }}_F$ and ${\overline{\sigma }}_{\mathrm{\Gamma }}$. Each lightcurve was simulated by extracting the F and Γ arrays from the assumed joint probability density function and the corresponding Pearson coefficient $r=r\left({log}_{10}F,\Gamma \right)$ was computed. A large number N of lightcurves was simulated and these were employed to estimate $\overline{r}$. The guidelines for choosing N are reported below. In the case of the correlation between log₁₀F_γ and Γ_γ at HE γ rays in Figure 7, the $\overline{r}$ values in the Pre-flare, Flare and Post-flare periods were found to be equal to 0.41, 0.32 and 0.25 (Table C.3), respectively, thus indicating significantly lower correlation than the corresponding r₀ values reported above. We checked that assuming log-normal distributions with mode $\overline{F}$ and ${\overline{\sigma }}_F$ rather than Gaussian functions for the flux F leads to estimates of $\overline{r}$ compatible within a few per cent with those in Table C.3. Conversely, in the case of the correlation between log₁₀F_X and Γ_X in the X-ray band (Figure 6), for each period the differences between $\overline{r}$ and r₀ were found to be below 0.01, due to the low relative uncertainties in X-rays (Table C.6).

To evaluate the significance of the correlation quantified by $\overline{r}$ under the null hypothesis of non-correlation, we considered the N lightcurves simulated as described above and, for each lightcurve, the F and Γ arrays were shuffled randomly and independently. In this way, we obtained a set of N lightcurves in which the flux and index are uncorrelated. We used these to construct the distribution of Pearson coefficients under the null hypothesis, from which we derived the p-value and significance of $\overline{r}$. To probe a correlation with a given p-value, we used a number of simulations N≳100/p-value. It is worth mentioning that including the statistical uncertainties has no major effect on the distributions of the Pearson coefficient under the hypothesis of non-correlation.

Appendix B: Comparison between the one-zone and two-zone model solutions

The estimation of the parameters of the one-zone and two-zone model solutions presented in Figures 11a and 11b, respectively, was performed using a ‘fit-by-eye’ strategy. Nonetheless, we report here quantitative indications that the two-zone model is more adequate to describe the broadband SED than the one-zone model. Table B.1 reports the χ² values for the one-zone and two-zone solutions, computed using the spectral points from Swift-UVOT, Swift-XRT, Fermi-LAT and MAGIC. The radio band was not included in Table B.1 since the one-zone model does not reproduce the spectrum at the lowest frequencies, as discussed in Section 6. The χ² values computed in the optical/UV band and in X-rays are higher than those in γ rays. This is a consequence of the small uncertainties of the spectral points in the optical/UV band and in X-rays with respect to those at HE and VHE γ rays. The χ² values computed for the two-zone solution are lower than those obtained for the one-zone model, indicating that the former provides a better agreement with the data. However, a more sophisticated procedure should be employed to statistically claim which model provides a more accurate description of the MWL data. Such a procedure should take into account the number of free parameters involved in the optimization and the systematic uncertainties on the measurements from each instrument. Given the high dimensionality of the model parameter space, such a procedure is outside the scope of this work.

Appendix C: Details of the multiwavelength analysis

In this section, we provide more detail on the observations and the analyses performed in the different energy ranges.

All Tables

All Figures

Fig. 1. B2 1811+31 long-term lightcurve collected from 2005 to 2023. From top to bottom panels: HE γ-ray flux above 100 MeV from Fermi-LAT monthly binned data, X-ray flux in the 0.3−10 keV energy range from Swift-XRT, optical R-band data and radio data. The dashed red line marks the Fermi-LAT high-state detection on October 1, 2020 (MJD 59123). The shaded light orange, light blue, and green bands in the top panel indicate the “Pre-flare”, “Flare” and “Post-flare” periods, respectively (Table 2).

In the text

Fig. 2. MWL lightcurve of B2 1811+31 in HE γ-rays (weekly binned Fermi-LAT data, upper panel) and optical R-band (instruments are reported in the legend, lower panel) in a period of around 250 days surrounding the October 2020 high state. This period corresponds approximately to the ‘Flare’ period (Table 2). In the upper panel, the solid red lines show the double-exponential fit (Eq. (5)) of the rising and falling trends in the Fermi-LAT lightcurve (Section 3.1.2).

In the text

Fig. 3. MWL lightcurve of B2 1811+31 in a period of approximately 70 days surrounding the Fermi-LAT high-state detection on MJD 59123, marked by the dashed red line. From top to bottom panels: VHE γ-ray flux above 135 GeV from MAGIC, HE γ-ray flux above 100 MeV from Fermi-LAT in daily bins and using the ‘adaptive-binning’ method, X-ray flux in the 0.3−10 keV range from Swift-XRT, optical/UV data in the photometric filters of Swift-UVOT, optical R-band data and radio data. The shaded yellow and light violet bands indicate periods of 48 h and 72 h surrounding the two sets of MAGIC consecutive observation nights, respectively, defined as periods A and B in Section 6.

In the text

	Fig. 4. SED in the 100 MeV−1 TeV energy range resulting from Fermi-LAT data from the “Pre-flare” (red), “Flare” (blue) and “Post-flare” (green) periods, as defined in Table 2. The average spectral indices and fluxes in the same energy range for the three periods are reported in Table 2.
In the text

	Fig. 5. Long-term evolution of the optical-to-X-ray SED of B2 1811+31 reconstructed from the Swift-UVOT and Swift-XRT observations.
In the text

	Fig. 6. Scatter plot of the intrinsic power-law spectral index ΓX and flux FX in the 0.3−10 keV energy range from the long-term Swift-XRT observations of B2 1811+31. The observations carried out the “Pre-flare”, “Flare” and “Post-flare” periods are marked in red, blue and green, respectively. The orange line marks the fit of the full dataset with Eq. (6).
In the text

	Fig. 7. Scatter plot of the PL spectral index Γγ and flux Fγ in the HE γ-ray data. In the top (bottom) panel, Γγ and Fγ derive from the Fermi-LAT lightcurve with 30-day bins (daily bins) shown in Figure 1 (Figure 3). The dotted lines mark the fits using Eq. (6).
In the text

Fig. 8. Summary of the zDCF analyses on the long-term B2 1811+31 lightcurves: (a) zDCF of the long-term HE γ-ray and R-band lightcurves, (c) zDCF of the long-term HE γ-ray and 15 GHz lightcurves, (d) zDCF of the long-term R-band and 15 GHz lightcurves. The 1σ and 3σ bands are indicated in dark and light gray, respectively. Panel (b) shows the distribution of the reconstructed zDCF at zero time-lag for 103 simulated uncorrelated HE γ-ray and optical lightcurves, whereas the red line marks the value reconstructed with the corresponding observed B2 1811+31 lightcurves in Figure 1.

In the text

Fig. 9. Fit of the infrared-to-X-ray SEDs characteristic of the B2 1811+31 low and high states with Eq. (7), in blue and red, respectively. The spectral points in the optical/UV range were reconstructed from Swift-UVOT observations, whereas the X-ray ones derive from Swift-XRT data. Infrared points are WISE data. The shaded areas correspond to the 1σ band fit uncertainties, while the best-fit synchrotron peak frequencies and their uncertainties are marked by the vertical solid and dotted lines.

In the text

Fig. 10. Spectral energy distribution at HE and VHE γ rays of B2 1811+31 during periods A and B. For the two periods, spectral points are shown in brown and violet and the 1σ uncertainty bands are in shaded yellow and light violet, respectively. The points with full round markers are derived from Fermi-LAT data. For the two periods, the best-fit parameters of the spectrum reconstructed from Fermi-LAT observations are reported in Table C.5. The points with full square markers show the MAGIC EBL-corrected SEDs. The observed MAGIC spectral points are indicated with hollow square markers and dashed error bars.

In the text

Fig. 11. Leptonic modeling of the broadband SED of B2 1811+31 during period B within the 2020 high state in VHE γ rays, marked by the red points. Archival data are indicated by the hollow gray points. The shaded blue and orange areas indicate the 1σ uncertainty bands on the intrinsic spectrum of the source derived from Fermi-LAT and MAGIC observations, respectively. The two insets show a zoom-in on the optical/UV and X-ray bands, respectively. For the two-zone SSC model, the contributions from the blob and core are indicated, respectively, by the orange and green lines, whereas the total is computed as the sum of the two. (a) One-zone SSC modeling. (b) Two-zone SSC modeling.

In the text

Fig. 12. Competition among the cooling timescales for the two-zone model in Figure 11b, as a function of the electron Lorentz factor γ. The times are in the observer frame; that is, corrected by the corresponding $\frac{1+z}{{\delta }_{\text{D}}}$ with respect to the frames comoving with the two regions. Dashed (dotted) curves refer to the blob (core) region. Blue and red curves indicate the synchrotron and SSC cooling timescales, respectively, whereas the green lines represent the light-crossing times. The orange markers indicate the γb values of the EED within the two emission regions (Table 4).

In the text