VHE gamma-ray detection of FSRQ QSO B1420+326 and modeling of its enhanced broadband state in 2020

MAGIC is a system of two IACTs with a mirror dish diameter of 17 m each (Aleksić et al., 2016a). The telescopes are located in the Canary Islands, on La Palma ($28.7^{\circ}$ N, $17.9^{\circ}$ W), at a height of 2200 m above sea level.

QSO B1420+326 has been observed by MAGIC on a few occasions following enhanced states at lower energies. We report on observations between 2019 December 31 (MJD=58848) and 2020 February 6 (MJD=58885). The observations consisted of several triggers of MAGIC Target of Opportunity program and are therefore, also due to bad weather periods, not continuous. The data selection was based on the atmospheric transmission and rates of background events. The total amount of good quality data is 14.0 hr, out of which 2.9 hr, taken between 2019 December 31 (MJD=58848) and 2020 January 4 (MJD=58852), has been taken with a special low-energy analogue trigger: SUM-Trigger-II (García et al., 2014). The data are analyzed using MARS, the standard analysis package of MAGIC (Zanin et al., 2013; Aleksić et al., 2016b). The data selection is based mainly on the atmospheric transmission measured with LIDAR (Fruck & Gaug, 2015) and on hadronic background rates. The effect of atmospheric absorption is corrected using LIDAR information. The Sum-Trigger-II part of the dataset is analyzed with dedicated low-energy analysis procedures including a special image cleaning (the so-called MaTaJu cleaning) with the cleaning thresholds tuned to the extragalactic field of view of QSO B1420+326 (Shayduk, 2013; Ceribella et al., 2019). For the part of the dataset during which no signal is detected, upper limits on the flux are computed following Rolke et al. (2005) at 95% confidence level.

The Large Area Telescope (LAT) is a pair-conversion telescope, launched on 2008 June 11 (MJD=54628) as one of the two scientific instruments on board the Fermi Gamma-ray Space Telescope (Atwood et al., 2009). Its energy range extends down to $\sim$ 30 MeV and up to $\sim 300$ GeV, with peak sensitivity at $\sim 1$ GeV. In the 4FGL (Abdollahi et al., 2020), QSO B1420+326 is associated to the gamma-ray source 4FGL J1422.3+3223, which has a $>100$ MeV flux of $(9.1\pm 1.3)\times 10^{-9}$  cm^-2 s^-1 and a power-law spectrum with photon index $2.38\pm 0.07$, obtained from data between 2008 June 11 (MJD=54682.65) and 2016 August 2 (MJD=57602.24). QSO B1420+326 does not appear in any of the Fermi-LAT hard-spectrum source catalogs (Ajello et al., 2017, see, e.g., the 3FHL catalog,), which is consistent with its relatively steep spectrum.

We use the Python package Fermipy (Wood et al., 2017) to analyze the Fermi-LAT data. We use Pass8 event data (Atwood et al., 2013) and select photons of the SOURCE class, in a square region of interest (ROI) of 10${}^{\circ}\times 10^{\circ}$, centered at the position of the target source. We perform a binned analysis with 10 bins per decade in energy and 0.1^∘ binning in RA and Dec, in the energy range 0.1–300 GeV, adopting the instrument response functions P8R3_SOURCE_V2. A correction for energy dispersion is included for all sources in the model except for the isotropic diffuse components. We apply a cut to include only the gamma-rays with zenith angle $<90^{\circ}$ to limit contamination from the Earth’s limb. We include in the model of the region all sources listed in the 4FGL within 15^∘ from the ROI center, along with the Galactic (Acero et al., 2016) and isotropic diffuse emission models (gll_iem_v06.fits and iso_P8R3_SOURCE_V2_v1.txt, respectively).

We first perform a likelihood fit using the full Fermi-LAT data set available at the time of the analysis, including events in the time range 2008 August 4 to 2020 February 11 (MJD 54682.66–58890.00). We fit the full spectrum of the target source, the diffuse models, and the normalization of catalog sources within 5^∘ as free parameters. We also optimize the target source localization, taking advantage of the $\sim 3.5$ extra years of data with respect to 4FGL. The detection significance is estimated with the Test Statistic (Mattox et al., 1996, TS,). We search for new sources by performing a TS map of the ROI. No significant (TS$>$25) new gamma-ray source is found in this analysis. Although mild excesses with TS$\sim$10 are seen in the residuals, such fluctuations are to be expected when periods over long time ranges such as this are considered, and therefore we choose not to add these excesses as point sources in the ROI model. The model resulting from this initial fit is used as an input for the computation of the HE gamma-ray light curves.

The light curves are calculated by performing a likelihood fit in each time bin. The fitting strategy is designed to adjust the number of free parameters to the photon statistics available in each bin. In the Fermi-LAT source catalogs, a source is considered detected if TS is at least 25. In each light curve bin, we fit the full spectrum of the target source and the normalizations of the sources within 3^∘ from the ROI center. If the target source has TS $<$ 25, we progressively restrict the free parameters in the fit, reloading the average model at each step. First, we reduce the sources with free normalization to a radius of 1^∘, then we freeze all sources except the target, and finally, if the target is not significantly detected, we only fit its normalization, leaving the spectral parameters fixed to the average value from the initial model. We consider the target source to be detected in a given time bin if TS $>9$, and the signal-to-noise ratio (that is, flux divided by its error, or $F/\Delta F$) in that bin is larger than two. If this is not the case, we reported a 95% confidence upper limit on the flux.

We report light curves with fixed binning of one day and 30 days, and one with adaptive binning (Lott et al., 2012), setting a constant relative flux uncertainty of 15%. The latter method provides an estimate of the shortest time scale over which it is possible to obtain a statistically significant detection and a robust determination of the target’s spectral parameters.

We calculate 0.1–300 GeV spectral energy distributions (SEDs) in the time intervals listed in Section 3, by performing a likelihood fit in several energy bins. The number of bins is optimized as a trade-off between energy resolution and photon statistics. We also perform an analysis including data from 2008 August 4 (MJD=54682) to 2017 January 1 (MJD=57754), to compute a quiescent-state Fermi-LAT SED to which the high-state ones could be compared. This time range was chosen based on the monthly light curve, which shows some early signs of flaring activity between the second half of 2017 and the first half of 2018 (see Appendix A.3). This time range is quite similar to the one corresponding to the 4FGL catalog, but provides a small increase in photon statistics due to the later end time.

Finally, we calculate the probability for each single gamma ray recorded by the Fermi-LAT of being associated with QSO B1420+326, using the gtsrcprob tool, in order to investigate the highest energy photons in the 0.1-300 GeV energy range.

The Neil Gehrels Swift Observatory (Gehrels et al., 2004) carried out 26 observations of the source between 2018 January 20 (MJD=58138) and 2020 February 10 (MJD=58889), in particular on 14 individual days between 2020 January 2 and February 10 (MJD=58850 – 58889). The pointed observations were performed with both the X-ray Telescope (XRT; Burrows et al., 2005, 0.2–10.0 keV) and the Ultraviolet/Optical Telescope (UVOT; Roming et al., 2005, 170–600nm). The hard X-ray flux of this source is below the sensitivity of the BAT instrument for the short exposures of these observations; therefore, the data from this instrument are not used. Moreover, the source is not present in the Swift BAT 105-month catalog (Oh et al., 2018).

All XRT observations were performed in photon counting (PC) mode. The XRT spectra are generated with the Swift XRT data products generator tool at the UK Swift Science Data Centre⁴⁴ 4 http://www.swift.ac.uk/user_objects (Evans et al., 2009, for details see). Ancillary response files are generated with xrtmkarf, and account for different extraction regions, vignetting and point-spread function corrections. We use the spectral redistribution matrices v014 in the Calibration data base maintained by HEASARC. Some of the spectra have very few photons, and so we are not able to use $\chi^{2}$ statistics. To maintain the homogeneity in our analysis, we grouped the obtained spectra using the task grppha to have at least one count per bin and we perform the fit with the Cash statistics (Cash, 1979). The data collected during 2019 June 27 and 29 (MJD=58661 and 58663) are summed in order to have enough statistics to obtain a good spectral fit. We fit the spectra with an absorbed power-law using the photoelectric absorption model tbabs (Wilms et al., 2000), with a neutral hydrogen column density fixed to its Galactic value (Ben Bekhti et al., 2016, $1.08\times 10^{20}$ cm^-2;). We apply also a log-parabola model to the XRT data, testing if this model is preferred over a single power-law model on a statistical basis by applying an F-test. The log-parabola model is preferred over a simple power-law model only on 2020 January 28 and at 95% confidence level. However, this can be due to the low statistics of the single XRT observations.

XMM-Newton (Jansen et al., 2001) observed the source on 2020 January 24 between 04:44:07 and 11:49:07 (MJD 58872.3 – 58872.5) for a total duration of 25 ks (observation ID 0850180201). The 3 EPIC cameras were operated in the large-window mode with medium filter. The data are reduced using the XMM-Newton Science Analysis System (SAS v16.0.0), applying standard event selection and filtering. Inspection of the background light curves show that no strong flares were present during the observation, with good exposure times of 20, 24 and 24 ks for the pn, MOS1 and MOS2, respectively. For each of the detectors the source spectrum is extracted from a circular region of radius 30 arcsec centered on the source, and the background spectrum from a nearby region of radius 30 arcsec on the same chip. All the spectra are binned to contain at least 20 counts per bin to allow for $\chi^{2}$-based spectral fitting. All spectral fits are performed over the 0.3–10 keV energy range using XSPEC v.12.10.1. The energies of spectral features are quoted in the source rest frame. All errors are given at the 90% confidence level. The data from the three EPIC cameras are initially fitted separately, but since good agreement is found ($<5\%$) we proceed to fit them together. Galactic absorption is included in all fits using the tbabs model.

Three different models are applied: a simple power-law, a broken power-law, and a log-parabola model. The results of the fits are presented in Table 1. The F-test shows an improvement of the fit using both a broken power-law and a log-parabola model with respect to the simple power-law, with the log-parabola model providing the best fit.

In order to check for the presence of intrinsic absorption, a neutral absorber at the redshift of the source is added to this model, but it does not improve the fit quality and thus is not required. Moreover, no Fe line was detected in the spectrum. The 90% upper limit on the equivalent width (EW) of a narrow emission line at 6.4 keV is EW$<10$ eV.

During the Swift pointings, the UVOT instrument observed QSO B1420+326 in all its optical ($v$, $b$ and $u$) and UV ($w1$, $m2$ and $w2$) photometric bands (Poole et al., 2008; Breeveld et al., 2010). For each epoch, possible multiple exposures in the same filter are first summed with the task uvotimsum and then analyzed using the uvotsource task included in the HEAsoft package (v6.28) with the 20201026 release of the Swift/UVOTA CALDB. We check if the observations are affected by small-scale sensitivity problems⁵⁵ 5 https://swift.gsfc.nasa.gov/analysis/uvot_digest/sss_check.html. Source counts are extracted from a circular region of 5 arcsec radius centered on the source, while background counts are derived from a circular region of 20 arcsec radius in a nearby source-free region.

The Optical Monitor (OM) on board of the XMM-Newton satellite observed the source in the $u$, $w1$, $m2$, and $w2$ filters in imaging mode together with a fast readout window. The total exposure times of the imaging observations are: 3500 s ($u$), 3500 s ($w1$), 4400 s ($m2$) and 4400 s ($w2$). The data are processed using the SAS tasks omichain and omfchain.

The UVOT and OM flux densities are corrected for extinction using the E(B–V) value of 0.010 from Schlafly & Finkbeiner (2011) and the extinction laws from Cardelli et al. (1989).

The REM telescope (Zerbi et al., 2001; Covino et al., 2004), a robotic telescope located at the ESO Cerro La Silla observatory (Chile), performed optical photometric observations of QSO B1420+326 in the period 2020 January 24–February 6 (MJD=58872 – 58885). Observations were carried out with the Optical Slitless Spectrograph (ROSS2) obtaining three 240 s integration images in the optical g^′, r^′, i^′ bands. The REM data presented here were obtained as ToO observations triggered by the high gamma-ray flux observed by Fermi-LAT after the MAGIC detection. Instrumental magnitudes are obtained via aperture photometry and absolute calibration is performed by means of secondary standard stars in the field reported by the AAVSO Photometric All-Sky Survey (APASS) catalog⁶⁶ 6 https://www.aavso.org/apass. Transformation between the u^′ g^′ r^′ i^′ z^′ and UBVRI photometric systems are performed using the equations reported in Jester et al. (2005)⁷⁷ 7 https://www.sdss.org/dr16/algorithms/sdssUBVRITransform/.

Optical photometric (BVRI) and polarimetric (R band) observations were carried out at the 1.83 m Perkins telescope (Flagstaff, AZ, USA), 40 cm LX-200 telescope (St.Petersburg, Russia) and 70 cm AZT-8 telescope (Crimea) from 2020 January 23 to April 21 (MJD 58871 -- 58960). The photometric data are reduced using differential aperture photometry with respect to comparison stars in the quasar field⁸⁸ 8 see https://vo.astro.spbu.ru/vlar/opt_thumbs/b21420_1.png. The polarimetric observations obtained at the Perkins telescope were performed and reduced in the same manner as described in Jorstad et al. (2010). The details of polarimetric observations carried out at the LX-200 and AZT-8 telescopes can be found in Larionov et al. (2008).

Photometric optical observations of QSO B1420+326 were carried out with the MFOSC-P instrument (Srivastava et al., 2018), used in imaging mode, mounted on the 1.2 m telescope of Mount Abu IR Observatory (MIRO⁹⁹ 9 https://www.prl.res.in/~miro/). The source was observed in B, V, R, and I bands (Johnson-Cousins filters) on 2020 February 2 and 6 (MJD=58880 and 58885). MIRO is located at Gurushikhar peak in Mount Abu, India, at altitude of 1680 m and is operated by Physical Research Laboratory (PRL), Ahmedabad, India.

The optical data from National Astronomical Observatory (NAO) Rozhen, Bulgaria were obtained between 2020 January 24 and 26 (MJD=58872 – 58874). We used the 2-m RCC telescope with Andor iKon-L CCD camera (2048x2048 px, 13.5 $\mu$m/pixel) and the 50/70 cm Schmidt telescope with FLI PL-16803 CCD camera (4096x4096 px, 9 $\mu$m/pixel). Additional observations were carried out 2020 January 31 to February 2 (MJD=58879 - 58881) at Student Astronomical Observatory (SAO) Plana (Ovcharov et al., 2014) with 35 cm Newton telescope and SBIG STL-11000M CCD Camera (4008x2672 px, 9um/pixel). All cameras are equipped with standard photometric UBVRI Johnson-Cousins filters.

The data are reduced (including bias subtraction, flat-fielding, and cosmic-ray correction) and analyzed using standard photometry packages from IRAF¹⁰¹⁰ 10 http://iraf.noao.edu/. For each image the PSF value is measured and aperture photometry is applied. Standard stars from the SDSS DR12 and VizieR catalogs are used for photometric calibration after applying transformation equations¹¹¹¹ 11 http://www.sdss3.org/dr8/algorithms/sdssUBVRITransform.php##Lupton2005.

The Astronomical Observatory of the University of Siena observed QSO B1420+326 in its program devoted to optical photometry of blazars in support of MAGIC. The observatory runs a remotely operated 30 cm Maksutov-Cassegrain telescope installed on a Comec 10micron GM2000-QCI equatorial mount. The focal plane hosts a Sbig STL-6303 camera equipped with a 3072 x 2048 pixels KAF-6303E sensor; Johnson-Cousins BVRI filters are available. Multiple 300 s images of QSO B1420+326 were acquired at each visit. After standard dark current subtraction and flat-fielding, images are averaged and aperture photometry is performed on the average frame by means of the MaximDL software package. Reference and check stars in the field of view are selected from the APASS9 (Henden, et al., 2016) catalog. The reference R magnitudes are derived from those reported in the same APASS9 catalog after conversion between the two different photometric systems, following a formula taken from Munari (2012).

Additionally we use publicly available data in V-band and g-band of ASAS-SN (Shappee et al., 2014; Kochanek et al., 2017).

Conversion of magnitudes to energy fluxes is done using zero points of Bessell et al. (1998). Correction for the Galactic extinction is applied using the E(B–V) value of 0.011 from Schlafly & Finkbeiner (2011) and the extinction laws from Cardelli et al. (1989).

We have performed observations of optical spectra of the quasar QSO B1420+326 using the 4.3 m Lowell Discovery Telescope (LDT; Lowell Obs., Flagstaff, AZ) equipped with the DeVeny spectrograph and the Large Monolithic Imager (LMI), in response to the detection of the source at VHE by MAGIC on 2020 January 21 (MJD=58869). We employed a grating setting of 300 grooves per millimeter, which provides spectra from 3300Å to 7500Å with a dispersion of 2.17 Å per pixel, a blaze wavelength of 5000Å  and a slit width of 2.5${}^{\prime\,\prime}$. The spectroscopic observations were performed on 2020 January 28 (MJD 58876.578), February 8 (MJD 58887.497) and 25 (MJD 58904.368). During this month the brightness of the quasar fell from 14.6 mag to 16.25 mag in R band. Each observation of the quasar consisted of 3 exposures of 600 s (900 s on February 25). Observations of a comparison star HD126944 of A type, located $\sim$86^′ from the quasar, were performed before and after target observations, with two 30 s exposures for each observation. Bias and flat-field images were obtained regularly. The LDT allows a switch between different instruments in 2-3 min. Therefore, photometry of the quasar using the LMI in V and R filters were performed just before or after spectral observations to calibrate the spectra. The observations are reduced using programs written in IDL (v.8.6.1) that implement the technique described in Vacca et al. (2002) developed for reduction of spectra obtained with SpeX at the NASA Infrared Telescope Facility on Mauna Kea.

The NIR observations were carried out with the camera CANICA (Carrasco et al., 2017), along with the Guillermo Haro 2.1m telescope (OAGH), located at Cananea Sonora Mexico. The camera is based on a Hawaii detector of 1024 by 1024 pixels with a plate scale of 0.32 arcsec per pixel. The data are part of the monitoring program "NIR photometry of AGNs with Gamma Ray emission detected by Fermi-LAT" . Relative photometry is obtained with respect to the 2MASS point sources included in the field of view (5.5 arcmin). Absolute fluxes are obtained adopting the zero point values of 2MASS derived by Cohen et al. (2003). The host galaxy is not detected in IR by the 2MASS survey in any of the three bands, resulting in an upper limit of H $\sim$ 17.7 mag, much weaker than 13–11 mag observed by CANICA during the investigated period. Therefore the effect of the host galaxy is negligible.

The 37 GHz observations were made with the 13.7 m diameter Metsähovi radio telescope. The detection limit of the telescope at 37 GHz is on the order of 0.2 Jy under optimal conditions. Data points with a signal-to-noise ratio < 4 are handled as non-detections. The flux density scale is set by observations of DR 21. Sources NGC 7027, 3C 274 and 3C 84 are used as secondary calibrators. A detailed description of the data reduction and analysis is given in Teräsranta et al. (1998). The error estimate in the flux density includes the contribution from the measurement rms and the uncertainty of the absolute calibration.

We also use publicly available 15 GHz OVRO data (Richards et al., 2011) and 8.63 GHz data from Badary RT32 reported in (Kharinov, 2020).

We requested Director’s Discretionary Time (DDT) with the Very Long Baseline Array (VLBA) following the gamma-ray activity and VHE detection of the quasar QSO B1420+326  and were granted 6 epochs of observations of the source separated by approximately 1 month intervals (ID BD227), with 8 hrs per epoch. Thus far, we have obtained 3 epochs of observations under the project, on 2020 March 8, May 10, and June 6 (MJD=58916, 58979, 59006). The observations were performed with all 10 antennas in continuum, dual circular polarization mode at 43 GHz using 4 intermediate frequency bands (IFs), each of 64 MHz width. The data were correlated at the National Radio Astronomy Observatory (NRAO, Soccoro, NM) using the VLBA DiFX software correlator. Five sources were observed at each epoch (QSO B1420+326  3C 279, 3C 345, PKS 1055+18, and B2 1308+326), with 60% of the time devoted to the main target, QSO B1420+326. The sources 3C 279, 3C 345, and PKS 1055+018 are used for fringe finding during the correlation. The data are reduced using the Astronomical Image Processing System software (van Moorsel, Kemball, & Greisen, 1996, AIPS,)) and Difmap package (Shepherd, 1997) in the same manner as described in Jorstad et al. (2017), except without averaging of the final calibrated data over IFs. We use the sources observed along with QSO B1420+326 to perform absolute calibration of the electric vector position angle (EVPA), since these sources are monitored in the VLBA-BU-BLAZAR program¹²¹² 12 www.bu.edu/blazars/VLBAproject.html, so that their polarization properties at 43 GHz are known.

The first detection of the VHE gamma-ray emission from QSO B1420+326 was achieved on 2020 January 20 (MJD=58868). A highly significant detection of $14.3\sigma$ is obtained in 1.6 hr of effective time (see Fig. 3).

In the subsequent period of 2020 January 26 to February 1 (MJD=58874 – 58880) further hints of signal are obtained, with the highest significance of the excess ($6.6\sigma$) on the night of 2020 January 31 (MJD=58879), which has also the longest exposure of 2.5 hr.

During the flare (period C), the flux observed by MAGIC above $100\,$GeV reached $(7.8\pm 1.3_{\rm stat})\times 10^{-11}\mathrm{cm^{-2}s^{-1}}$. The observed spectrum in this period can be described by a power-law: $dN/dE=(2.49\pm 0.31_{\rm stat})\times 10^{-9}(E/100\,\mathrm{GeV})^{-4.22\pm 0.24_{\rm stat}}\mathrm{[TeV^{-1}cm^{-2}s^{-1}]}$. Correcting for the EBL absorption according to Domínguez et al. (2011), the unattenuated spectrum can be described as $dN/dE\,e^{\tau_{\rm EBL}(E)}=(4.04\pm 0.54_{\rm stat})\times 10^{-9}(E/100\,\mathrm{GeV})^{-3.57\pm 0.29_{\rm stat}}\mathrm{[TeV^{-1}cm^{-2}s^{-1}]}$.

After the flare (period D), significant gamma-ray emission is detected again with $5.6\sigma$, but at about half the flare level: $F_{D}(>100\,\mathrm{GeV})=(3.9\pm 0.6_{\rm stat})\times 10^{-11}\mathrm{cm^{-2}s^{-1}}$. The observed spectrum in this period can be described: $dN/dE=(0.91\pm 0.13_{\rm stat})\times 10^{-9}(E/100\,\mathrm{GeV})^{-3.90\pm 0.25_{\rm stat}}\mathrm{[TeV^{-1}cm^{-2}s^{-1}]}$. Correcting for the EBL absorption according to Domínguez et al. (2011), the spectrum can be described as $dN/dE\,e^{\tau_{\rm EBL}(E)}=(1.64\pm 0.22_{\rm stat})\times 10^{-9}(E/100\,\mathrm{GeV})^{-2.87\pm 0.36_{\rm stat}}\mathrm{[TeV^{-1}cm^{-2}s^{-1}]}$. Despite enhancement of the VHE gamma-ray flux by a factor of two, the spectral indices during and after the flare are consistent within $1.5\sigma$. It should be noted however that, in particular in period D, the uncertainties of the photon index are large.

Before the flare (period A), possibly due to shorter observations under less favorable zenith angle, no significant emission is detected and only a 95% C.L. limit of $<4.1\times 10^{-11}\mathrm{cm^{-2}s^{-1}}$ can be placed for the flux above 100 GeV. The limit is comparable to the emission detected from the source in period D. The SED of QSO B1420+326 observed by MAGIC in different periods is shown in Fig. 4.

The first detection of a HE outburst from B2 1420+326 was reported in 2018 December (Ciprini, 2018), with a flux increase by more than two orders of magnitude with respect to the average 4FGL value and significant spectral hardening. A similar spectral hardening has been reported in 2019 July (Angioni, 2019), where the first evidence of $>10$ GeV photons was provided, and again in 2020 January (Ciprini & Cheung, 2020).

The daily Fermi-LAT light curve is shown in Fig. 1, including the flux (second panel) and photon index (third panel). Both the flux and the photon index are significantly variable in this time interval, based on a simple $\chi^{2}$ test. The Fermi-LAT recorded a peak daily flux (E $>$ 100 MeV) of $(2.6\pm 0.2)\times 10^{-6}$ cm^-2 s^-1 on 2020 January 19 (MJD=58867), corresponding to about $300$ times the average value reported in the 4FGL catalog. The highest-energy photon observed by the Fermi-LAT was recorded two days prior (2020 January 17, MJD=58865), with an energy of $\sim 150$ GeV, providing the first indication of VHE emission from QSO B1420+326¹³¹³ 13 The probability of this photon being associated with the target source is $>99.9$%, as obtained from the gtsrcprob tool.. Accordingly, the Fermi-LAT recorded the hardest daily spectrum on the same date, with a photon index of 1.72 $\pm$ 0.08, a flux of $(1.2\pm 0.2)\times 10^{-6}$ cm^-2 s^-1 and a TS of 496. We also note that the photon index was consistently harder than the 4FGL catalog value ($2.38\pm 0.07$) during most of the time range shown in Fig 1.

The adaptively binned light curve is shown in Fig. 5. The shortest adaptive bin is centered on 2020 January 19, 00:17:45 (MJD=58867.012), and has a width of $\sim 6$ hours. The highest flux is recorded in the same bin, reaching a value of $(3.6\pm 0.5)\times 10^{-6}$ cm^-2 s^-1, i.e. 400 times higher than the 4FGL value.

We note that a test for spectral curvature was performed in all time bins of all light curves, and a power-law spectrum was found to be the best representation in all time intervals.

Additionally, we perform a likelihood fit over the full time interval included in Fig. 1, to characterize the average source properties in the HE band in this flaring state. For this time interval, the LogParabola model is preferred ($TS_{curv}=36$ ¹⁴¹⁴ 14 $TS_{curv}=2(\ln{L_{LP}}-\ln{L_{PL}})$ is used to check if a statistically significant curvature is detected using a LogParabola model compared with the PowerLaw model, where $\ln{L}$ indicates the logarithm of the likelihood for each model. A source is considered to have a statistically significant curvature if $TS_{curv}>25$.) with respect to a simple power law to describe the gamma-ray spectrum of the source. It yields a photon flux $(1.09\pm 0.02)\times 10^{-6}$ cm^-2 s^-1, and spectral parameters $\alpha=1.97\pm 0.02$, and $\beta=0.07\pm 0.01$ ¹⁵¹⁵ 15 The functional form of the LogParabola spectral model is $dN/dE=N_{0}(E/E_{b})^{-[\alpha+\beta\log(E/E_{b})]}$, where $N_{0}$ is the normalization, $E_{b}$ is the reference energy at which $N_{0}$ is measured, $\alpha$ is the slope and $\beta$ is the curvature parameter.. We also verify that there are no new significant point sources in addition to the initial 4FGL model during this period.

Finally, as mentioned in Section 2, we perform separate likelihood fits corresponding to the periods used to build time-resolved SEDs. The results of these fits are summarized in Table 3, together with the ones for the quiescent state.

There is no strong variability of the X-ray flux in the investigated period with an increasing trend from 2020 January 05 (MJD=58853) to 2020 January 25 (MJD=58873), where a peak flux a factor of 2 higher is observed. However there is a clear variability of the X-ray spectral index from hard values before the flare to much softer at the time around the optical and gamma-ray flares, back to hard values. This indicates a shift of the synchrotron peak of the source and connected with it, the shift of the crossing point between the synchrotron and the IC component. This is clearly visible in the X-ray spectra during the flare and during XMM-Newton observations (see Fig. 2). Results of the spectral fits to individual days of Swift-XRT observations are given in Appendix A.2.

Compared to the historical measurements, the optical emission is $\sim 1.5$ orders of magnitude higher throughout the investigated period. Moreover, during that period a strong optical flare is observed on 2020 January 19 (MJD=58867) with a variability time scale of the order of a few days. V-band observations performed during one of the Fermi-LAT peaks show an increase by nearly an order of magnitude with respect to observations at the beginning of Period A. Similar variability pattern is seen also in IR and UV ranges. However, the spectral shape in the IR-UV range varies during the investigated period, with the spectrum becoming bluer (harder) during the optical flare (see Fig. 2). This is consistent with the X-ray behavior of the source that also suggests a shift of the synchrotron peak position. Recently other occurrences of comparable flaring activity in the optical had been observed (Marchini et al., 2019, e.g., in July 2019,), reporting a $R=13.7$ magnitude even slightly brighter than the brightest $R$ point in Fig 1.

Also optical polarization shows interesting behavior with a dip of the polarization percentage at a few per cent level and a concurrent rotation by $\sim 150^{\circ}$. Similar EVPA rotation during low level of polarization has been also seen contemporaneous with VHE emission in, e.g., PKS 1510–089 (Aleksić et al., 2014). Large variability of EVPA (down to a time scale of 3 hours) has been also seen in PKS 0736$+$017, another FSRQ during the period of VHE gamma-ray detection (H. E. S. S. Collaboration et al., 2020).

Fig. 6 presents the optical spectra of QSO B1420+326 obtained with LDT at three epochs of different activity states (taken mostly after the gamma-ray flaring activity). The spectra show the presence of MgII emission line at $\lambda=$4706 Å (rest wavelength of 2798 Å) and a bump between 3800 Å and 4100 Å (rest frame $\sim$2260$-$2450 Å). The latter appears to be part of an FeII emission complex, whose strongest UV lines fall in the 2200-2600 Å range (Baldwin et al., 2004, e.g., ). The spectra also include three prominent absorption lines at $\lambda\sim$4001Å  4337Å  and 4477Å that intensify as the quasar brightens. The absorption lines are most likely intervening MgII systems at redshifts of $z\approx 0.43$, 0.55, and 0.60, respectively.

We analyze the characteristics of the MgII emission line and FeII bump as a function of continuum brightness. Fig. 7 shows an approximate Gaussian fit to the MgII line, while Fig. 8 plots similar modeling of the FeII bump at all three epochs. The parameters of the Gaussian and velocities of clouds, as well as the flux of the lines, are given in Table 4, although we are unable to estimate the velocity of gas producing FeII lines, since the bump consists of $>$100 lines. There is a significant difference between the MgII line and FeII bump behavior: i) the flux of the MgII line remains constant within 1$\sigma$ uncertainty independent of the continuum brightness, while the FeII bump increases in flux with the continuum level (see Fig. 9); ii) the central wavelength of the MgII line fits does not show a shift with respect to the rest wavelength, while $\lambda_{\circ}$ of the FeII bump shifts toward the blue side as the time after the VHE event passes. Unfortunately, it is not possible to distinguish whether the shift is due to a relative change of the brightness of lines that form the FeII bump, or due to gas motion toward the observer; and iii) the FWHM of the FeII bump is very stable despite the correlation of its flux with the continuum level, while the velocity of gas where the MgII line originates increases with time after the VHE event. We note also a significant change of the equivalent width (EW) of the MgII line with the continuum, with EW decreasing as the continuum rises. This questions the identification of QSO B1420+326 as a FSRQ, however, Table 4 shows that at the lower levels of activity EW$>$5Å for the MgII line¹⁶¹⁶ 16 EW of 5Å is the classical threshold between BL Lac and FSRQ objects, see, e.g., Sambruna et al. (1996).

The increase in the flux of the FeII bump with the continuum and a possible motion of gas producing FeII lines toward the observer are quite interesting. The short time lag between the continuum and line variability suggests that the FeII emission-line region is much smaller than the region producing the MgII line and/or lies close to the line of sight. It is possible that the non-thermal jet is interacting with an FeII-emitting cloud, while the MgII emission is excited by the underlying thermal accretion disk continuum, which varies on a much longer time scale.

We adopt an approach suggested by Ghisellini et al. (2014) (see also references therein), who used an estimate of the accretion disk (AD) luminosity based on the luminosity of the BLR, L${}_{AD}\approx$10 L_BLR. The known flux density of the MgII line, combined with the BLR template constructed by Vanden Berk et al. (2001) for a composite emission spectrum of a quasar using SDSS spectra, allow us to estimate the total luminosity of the BLR in QSO B1420+326  $L_{BLR}$ = (1.8$\pm$0.2)$\times$10⁴⁵ erg s^-1. This, for luminosity distance of 4256.4 Mpctranslates to $L_{AD}\approx 2\times 10^{46}$ erg s^-1, with an uncertainty of a factor of $\sim 2$, as suggested by Ghisellini et al. (2014). Interestingly, the obtained luminosity of the accretion disk is rather high, in the high-end part of values shown for other sources Ghisellini et al. (2014).

Moderate variability is seen in radio observations (see also Appendix A.3). OVRO data during the investigated period show a gradual increase of the flux. No monotonous behavior is seen in flux at 37 GHz by Metsähovi, but a constant fit can be rejected at $\chi^{2}/N_{\rm dof}=51/21$. The amplitude of the variability is $\sim 10\%$.

The total and polarized intensity VLBA maps of QSO B1420+326 are presented in Fig. 10.

The images exhibit the bright VLBI core located at the north-western end of the jet and a weak extended jet at position angle $\sim 130^{\circ}$.

The total and polarized intensity images are modeled by circular components with Gaussian brightness distributions. Two bright features are apparent at each epoch: the core A0, and a knot K20. We assume that the core is a stationary feature of the jet and calculate parameters of K20 with respect to the core. Parameters of the modeling are given in Table 5.

According to Fig. 10 and Table 5, knot K20 is the most polarized feature of the jet. Fig. 11 shows the separation of K20 from the core as a function of time. According to a linear approximation, K20 moves with a proper motion $\mu$ = 0.30$\pm$0.04 mas/yr, which translates into superluminal motion down the jet with apparent speed $\beta_{\text{app}}=12.0\pm 1.7$ in units of $c$. Such motion suggests that the ejection¹⁷¹⁷ 17 passage of the centroid of K20 through the centroid of A0 of K20 through the VLBI core occurred on MJD 58831$\pm$21 (2019 December 13). Fig. 11 also shows the light curve of K20, which reveals a very fast decrease of the knot flux density by a factor of $\sim$4 in three months, which corresponds to the timescale of variability of 0.16$\pm$0.03 yr, as defined according to Burbidge, Jones, & O’Dell (1974). Applying the formalism proposed by Jorstad et al. (2005), this timescale of variability and the average size of K20 (0.084$\pm$0.015 mas according to Table 5) give a value of the Doppler factor of K20 $\delta$=33$\pm$9. The latter, along with the apparent speed of the knot, allow us to estimate the Lorentz factor, $\Gamma_{\text{b}}$=19$\pm$9, and viewing angle $\Theta_{\circ}$=1.1${}^{\circ}\pm$0.6^∘ of K20.

Using the proper motion and size of K20, we can estimate that during the time of VHE emission the upstream edge of the knot was passing through the centroid of the core. The direction of polarization in K20 is close to the jet direction, which suggests that K20 is most likely a moving shock whose surface is oriented transverse to the jet axis. This is supported by a higher degree of polarization of the knot compared with the core, which implies ordering of the magnetic field as expected by such a shock.