Multiwavelength study of the gravitationally lensed blazar QSO B0218+357 between 2016 and 2020

ABSTRACT

1 INTRODUCTION

QSO B0218+357, also known as S3 0218+35, is one of only a handful of flat spectrum radio quasars (FSRQs) detected in very high energy (VHE; ≳100 GeV) gamma-ray emission. It has a redshift of z_s = 0.944 ± 0.002 (Cohen, Lawrence & Blandford 2003; Paiano et al. 2017). The source showed strong variability in the GeV range in 2012 (Cheung et al. 2014) when a series of flares was observed by Fermi Large Area Telescope (LAT). Another flare was observed by Fermi-LAT in 2014, and during the follow-up the source was discovered in VHE gamma-rays by MAGIC telescopes (Buson et al. 2015a, b; Ahnen et al. 2016). Similarly to QSO B0218+357, GeV emission from the second gravitationally lensed source detected by Fermi-LAT, PKS 1830−211, also shows evidence of a measured delay between different lens images (Barnacka, Glicenstein & Moudden 2011). Observations of PKS 1830−211 by the H.E.S.S. telescopes following a GeV flare did not show any significant gamma-ray emission (H. E. S. S. Collaboration 2019).

QSO B0218+357 is gravitationally lensed by B0218+357 G, a spiral galaxy seen face-on at a redshift of z_l = 0.68466 ± 0.00004 (Carilli, Rupen & Yanny 1993). Strong gravitational lensing is observed when the lens is a galaxy or a cluster of galaxies. Such a massive lens can produce multiple images of the source separated by ∼arcseconds. Thus, the images of strongly lensed sources can be well resolved at wavelengths from radio to X-rays with existing instruments.

Stars can cause further lensing effects within a lensing galaxy. In such cases, the deflection angle of lensed images is of the order of microarcseconds. Thus, the effect is called microlensing. The change in position of microlensed images cannot be observed with existing instruments. The microlensing effect is observed as changes in the flux of the strongly lensed image.

The relative flux densities observed for lensed images depend on the geometry of the source-lens-observer system, and can be affected by microlensing. Further, different geometrical paths and gravitational delays cause the emission to arrive at different times in various images. In the case of QSO B0218+357, the image is composed of two images A and B separated by only 335 mas and an Einstein ring of a similar size (O’Dea et al. 1992). The A component (located westwards) is brighter and this signal precedes that from component B.

Variable radio emission observed in 1992/1993 and 1996/1997 with the Very Large Array at 5, 8.4, and 15 GHz yields time delays between these two components of 12 ± 3 d (Corbett et al. 1996), 10.5 ± 0.4 d (Biggs et al. 1999), |$10.1^{+1.5}_{-1.6}$| d, (Cohen et al. 2000), 11.8 ± 2.3 d (Eulaers & Magain 2011), and 11.3 ± 0.2 d (Biggs & Browne 2018). The statistical analysis of the 2012 high state Fermi-LAT >0.1 GeV light-curve autocorrelation function led to a similar value of the time delay (11.46 ± 0.16 d; Cheung et al. 2014). These values are consistent with the delay between the two components of the 2014 Fermi-LAT flare (Ahnen et al. 2016).

The gamma-ray emission of QSO B0218+357 comprises many flares with time-scales as short as a few hours. The short time-scales of gamma-ray flares combined with the ability of the Fermi-LAT observatory to monitor the sky continuously allow one to search for delayed counterparts of flares and put constraints on the magnification ratio. For example, the two-night-long sub-TeV flare was observed contemporaneously with the detection of the image B flare in Fermi-LAT (Ahnen et al. 2016)

Unfortunately, since the MAGIC observations in 2014 only covered the time of the B image of the flare, no measurement of the magnification ratio or delay could be obtained. Monitoring of QSO B0218+357 with Cherenkov telescopes during flaring events could enable the capture of multiple flares and constrain models of the magnification ratio and time delays.

At radio frequencies, the B component is 3.57–3.73 times fainter than the A component (Biggs et al. 1999). However, the observed ratio of magnification varies with the radio frequency (Mittal et al. 2006), possibly due to free–free absorption in the lens (Mittal, Porcas & Wucknitz 2007). In the optical range, the leading image is strongly absorbed, inverting the brightness ratio of the two images (Falco et al. 1999). It has been suggested that the optical absorption occurs in the host galaxy rather than the lens (Falomo et al. 2017). Interestingly, the magnification ratio observed at GeV energies shows variability. The average GeV magnification factor during 2012 high state was estimated to be ∼1 (Cheung et al. 2014), while during the 2014 flare it was comparable to or even larger than that measured at radio frequencies (Ahnen et al. 2016). Changes in the observed GeV magnification ratio can be interpreted as microlensing effects either due to individual stars (Vovk & Neronov 2016) or due to larger scale structures in the lens (Sitarek & Bednarek 2016).

Because it takes about 1–2 d for Fermi-LAT data to be collected, downlinked, and processed, it is difficult to trigger observations for phenomena with similar durations, like the two-night 2014 flare. Therefore, the shortness of the VHE gamma-ray emission significantly hinders the possibility of Target of Opportunity observations of a flare in both images. In addition observational visibility constraints further limit the possibility of a follow up of the delayed emission with ground-based instruments. Thus, since 2016, we have taken advantage of the 11 d delay to trigger MAGIC. Observational windows that allow visibility under favourable zenith angle conditions in moonless nights 11 d after each slot have been identified. During these time slots, MAGIC observations were performed, and contemporaneous multiwavelength (MWL) coverage from radio to GeV gamma-rays was obtained. In this paper, the results of these observations are reported. Additionally, an MWL campaign on QSO B0218+357, organized in 2020 August in response to a hint of enhanced activity in the source, is also reported.

In Section 2, the instruments that took part in the MWL campaign, the data taken, and the analysis methods are described. The results are reported in Section 3. In Section 4, the broad-band emission of the low state of the source is modelled. The paper concludes with a summary of the results in Section 5.

2 OBSERVATIONS AND DATA ANALYSIS

QSO B0218+357 was observed over a broad energy range: radio (The Owens Valley Radio Observatory, OVRO), radio interferometry (a joint VLBI array of KVN (Korean VLBI Network) and VERA (VLBI Exploration of Radio Astrometry), KaVA), optical [Kunliga Vetenskapsakademi, KVA and Nordic Optical telescope, NOT; Neil Gehrels Swift observatory (Swift) Ultraviolet/Optical Telescope (Swift-UVOT) and X-ray Multi-Mirror Optical Monitor (XMM–OM)), X-ray (X-ray Telescope (Swift-XRT) and XMM–Newton], GeV gamma-rays (Fermi-LAT), and VHE gamma-rays (MAGIC). During the 2020 August MWL campaign, dedicated observations by Himalayan Chandra Telescope (HCT), Joan Oró Telescope (TJO), and Metsähovi were taken.

The historical data obtained via the Space Science Data Center¹ from the following catalogues are also used: CLASS (Myers et al. 2003), JVASPOL (Jackson et al. 2007), KUEHR (Kuehr et al. 1981), NIEPPOCAT (Nieppola et al. 2007), NVSS (Condon et al. 1998), Planck (Planck Collaboration VII, XXVIII, XXVI 2011, 2014, 2016), GB6 (Gregory et al. 1996), GB87CAT (Gregory & Condon 1991), WMAP5 (Wright et al. 2009), WISE (Wright et al. 2010), 1SWXRT (D’Elia et al. 2013), 1SXPS (Evans et al. 2014), and FGL (Abdo et al. 2010; Nolan et al. 2012; Acero et al. 2015).

2.1 MAGIC

MAGIC is a system of two imaging atmospheric Cherenkov telescopes with a mirror dish diameter of 17 m each. The telescopes are located in the Canary Islands, on La Palma (28.7° N, 17.9° W), at a height of 2200 m above sea level (Aleksić et al. 2016a). The data were analysed using mars, the standard analysis package of MAGIC (Zanin et al. 2013; Aleksić et al. 2016b). Wherever applicable, upper limits on the flux were computed following the approach of Rolke, López & Conrad (2005) using a 95 per cent confidence level and assuming a 30 per cent total systematic uncertainty on the collection area.

The regular monitoring observations were performed between MJD 57397 and 58875 in dark night conditions. The monitoring time slots were scheduled so as to allow for possible observations in ∼11 d if enhanced emission was seen. This results in possible observation slots (up to two per moon period) lasting between 1 and 6 d. Motivated by the 2-d duration of the 2014 VHE flare, in such slots observations every second night were scheduled (on some occasions this scheme was modified due to weather conditions or competing sources). After the data selection, based mainly on the atmospheric transmission measured with LIDAR (Fruck & Gaug 2015) and on hadronic background rates, the data set consists of 72.7 h, spread over 73 nights.

Since MJD 58122 the data have been taken with the novel Sum-Trigger-II (Dazzi et al. 2021). The Sum-Trigger-II part of the data set consists of 38.4 h and was analysed with dedicated low-energy analysis procedures including a special image cleaning procedure (Shayduk 2013; Ceribella et al. 2019).

Additionally, during the 2020 August campaign, 2 h of good quality data were taken on MJD 59081 and 59082. Due to a forest fire observations on MJD 59083 could not be used.

2.2 Fermi-LAT

In order to estimate weekly and daily fluxes of QSO B0218+357, the number of free spectral parameters is limited. The sources located within 7° from QSO B0218+357 had their normalization set as a free parameter if their variability index was higher than 18.483,⁵ while all sources within 5° from QSO B0218+357 had their normalization free to vary if their TS was higher than 40 integrated over the full time period. The spectral indices of all the sources with free normalization were left as a free parameter if the source showed a TS value higher than 25 over the integration time (weekly or daily), in all the other cases the indexes were frozen to the value obtained in the overall fit.

The spectral analysis was also performed in two different smaller intervals corresponding to the optical/GeV flare (MJD 57600-57700) and to the X-ray flare (MJD 58860.7–58866.7). The source was significantly detected in both time intervals, with a significance of 45.9σ and 6.6σ, respectively.

2.3 XMM–Newton

XMM–Newton (Jansen et al. 2001) observed the source four times between 2019 August and 2020 January. The integration times of the observations were in the range of 11.3–19.8 ks. The EPIC pn CCD camera (Strüder et al. 2001) operated in full-frame mode with the medium filter applied during all the observations. The data were processed using the XMM-Newton Science Analysis System (sas v.18.0.0; Gabriel et al. 2004) following standard settings and using the calibration files available at the time of the data reduction. The EPIC pn Observation Data Files (ODFs) were processed with the sas-task epproc in order to generate the event files. Event files were cleaned of bad pixels, and events spread at most in two contiguous pixels (PATTERN≤4) were selected. Periods of high background levels were removed by analysing the light curves of the count rate at energies higher than 10 keV. The resulting net-exposure times are reported in Table 4. In order to include all of the source counts and simultaneously minimize the background contribution, source counts were extracted from a circular region of radius between 30 and 35 arcsec. The background counts were extracted from a circular region of radius 50–65 arcsec located on a blank area of the detector close to the source. Response matrices for spectral fitting were obtained using the sas-task rmfgen and arfgen. All the spectra were binned in order to have no less than 20 counts in each background-subtracted spectral channel, and the instrumental energy resolution was not oversampled by a factor larger than 3.

X-ray spectral analyses were carried out with xspec v.12.9.1 (Arnaud 1996). No variability in the spectra of the XMM–Newton observations performed at MJD 58697, 58721, and 58724 (obs. ID 0850400301, 0850400401, 0850400501) was observed, thus all the observations were combined with the sas-task epicspeccombine for the spectral modelling of the source. In contrast, the observation performed at MJD 58863.7 (obs. ID 0850400601) indicated an increase of the X-ray flux by a factor of ∼1.4 with respect to previous observations, thus this spectrum was fitted separately.

2.4 Swift-XRT

The X-ray Telescope (XRT; Burrows et al. 2004) onboard the Neil Gehrels Swift observatory (Swift) observed the source four times between 2016 January and 2020 January. Additionally seven pointings were taken around the time of the 2020 August campaign. Due to the source faintness, all of these observations were performed in photon counting mode. The event lists for the period of interest were downloaded from the publicly available swiftxrlog (Swift-XRT Instrument Log).⁶ The data were processed using the standard data analysis procedure (Evans et al. 2009) and the configuration described by Fallah Ramazani, Lindfors & Nilsson (2017) for blazars. The spectra of each observation were binned in a way that each bin contains one count. Therefore, the maximum likelihood-based statistic for Poisson data (Cash statistics; Cash 1979) method was used in the spectral fitting procedure and flux measurements of individual observations.

No spectral variability was observed within the Swift-XRT data. In order to evaluate the average state of the source during the monitoring, two combined spectra were produced using the observational data taken during 2016–2017 (OBSIDs 00032533003, 00032533005, 00032533006, and 00032533007) and 2020 August (OBSIDs 00032533008, 00032533009, 00032533010, 00032533011, 00032533012, 00032533013, 00032533014, 00032533015). These spectra are binned in a way that each bin contains 20 counts. Therefore, the maximum likelihood-based statistic for Gaussian data method is the spectral fitting procedure of these two spectra.

2.5 UV

During the four monitoring Swift pointings, the UVOT instrument observed the source in the u optical photometric band (Poole et al. 2008; Breeveld et al. 2010). The data were analysed using the uvotimsum and uvotsource tasks included in the heasoft package (v6.28) with the 20201026 release of the Swift/UVOTA CALDB. Source counts were extracted from a circular region of 5 arcsec radius centred on the source, while background counts were derived from a circular region of 20 arcsec radius in a nearby source-free region. The source was not detected with a significance higher than 3σ in the single observations, therefore the four UVOT images were summed using the uvotimsum task and analysed the summed image with the uvotsource task.

The Optical Monitor (OM) observed the source four times in u filters in imaging mode. The total exposure times of the imaging observations were 16 400, 17 700, 11 300, and 19 800 s. The data were processed using the SAS task omichain. For the count rate to flux conversion the conversion factors given in the SAS watchout dedicated page⁷ were used.

During the 2020 August campaign five observations with Swift-UVOT were taken in u band. None of these pointings resulted in the detection of a signal above 2σ significance. Similarly to 2016–2020 monitoring data, a stacked analysis was performed to evaluate the average emission in this period.

The UVOT and OM flux densities were corrected for Galactic extinction using a value A_U = 0.299 mag (Schlafly & Finkbeiner 2011).

2.6 Optical

The source was monitored with the NOT and 35 cm Celestron telescope attached to the KVA telescope. Both telescopes are located at the same site as the MAGIC telescopes and the NOT observations were carried out quasi-simultaneously with MAGIC, while KVA performed additional monitoring more frequently. Starting in 2020 July, the source was also monitored with the TJO at the Montsec Astronomical Observatory.⁸ In addition during the August MWL campaign the source was observed with the HCT.⁹

NOT observations were carried out using ALFOSC in B, V, R, and I bands, while the KVA and TJO observations operated in R band only. The data were analysed using the semi-automatic pipeline and standard procedures of differential photometry (Nilsson et al. 2018). The same comparison and control stars were used as in Ahnen et al. (2016)¹⁰. The r-, g-, and i-band magnitudes of the stars were available in the PANSTARRS database.¹¹ The magnitudes in B, V, and I band were calculated from i-, r-, and g-band magnitudes using the formulae of Lupton (2005).¹² Using those formula, consistency of the R-band values derived in Ahnen et al. (2016) was checked. The I-band filter used at the NOT differs from the standard I-band filter enough for the colour correction to become significant for a very red input spectrum as in this case (Falomo et al. 2017). The spectrum obtained by Falomo et al. (2017) was downloaded from the ZBLLLAC repository¹³ and synthetic photometry was performed through the standard I band and the NOT I band. This showed that the NOT I-band magnitudes needed to be corrected by +0.08 mag to transform to the standard system. The aperture used was 3 arcsec, slightly smaller than the 4 arcsec used for KVA and TJO data.

The source was observed with HCT in four epochs (MJD 59081–59085). The observations were carried out in the Bessell U, B, V, R, and I bands available with HFOSC. The data were reduced in a standard manner using various tasks available in iraf. Aperture photometry was performed on the source and nearby stars. The standard magnitude of the source was obtained using differential photometry with the same comparison and control stars as used by Ahnen et al. (2016).

The observed magnitudes were corrected for the galactic extinction using values obtained from the NED¹⁴ (Schlafly & Finkbeiner 2011) and listed in Table 1. The magnitudes were converted to flux densities using formula F = F₀ × 10^−mag/2.5 with F₀ = 4260 Jy in B, F₀ = 3640 Jy in V, F₀ = 3080 Jy in R, and F₀ = 2550 Jy in I.

The flux densities needed to be corrected for the contribution from the host galaxy of QSO B0218+357 at (z_s = 0.944) and the lens galaxy at z_l = 0.684. Jackson, Xanthopoulos & Browne (2000) imaged this target with the HST through the F555W, F814W, and F160W filters and measured the flux density from the face-on spiral galaxy to be (6 ± 2), (13 ± 2), and (15 ± 2) (10⁻¹⁸ erg s⁻¹ cm ⁻² Å⁻¹), respectively. Since QSO B0218+357 is classified as a FSRQ (Paiano et al. 2017; Abdollahi et al. 2020), the host galaxy is likely to be a luminous (M_K ∼ −26.5) bulge-dominated galaxy (e.g. Olguín-Iglesias et al. 2016). The R-band magnitude of such a galaxy at the redshift of 0.944 would be ∼22.5 mag. An early-type galaxy template was taken from Mannucci et al. (2001), redshifted to z = 0.944 and integrated over the R-band filter transmission. Then the template was scaled to match the integrated flux density to R = 22.5. The scaled spectrum corresponds to ∼40 per cent of the flux densities observed by Jackson et al. (2000), i.e. a significant part of the ‘spiral galaxy’ surrounding component B could actually be the host galaxy. This is what Falomo et al. (2017) propose based on a high signal-to-noise ratio spectrum of B0218+357. Their spectrum shows gaseous absorption lines at the lens redshift, but no stellar photospheric lines are detected, which led them to propose that the spiral structure belongs to the host galaxy, not the lens. It is impossible to determine the relative contributions of the lens galaxy and the host galaxy from the present data, especially since the latter may also be lensed and absorbed by the former. A simple assumption, that 100 per cent of the flux densities determined by Jackson et al. (2000) arise from the lens is used. Thus a fit of a late type (Sa) galaxy template from Mannucci et al. (2001), redshifted to 0.684, to the Jackson et al. (2000) flux densities was carried out. Then synthetic photometry was performed through the BVRI bands to the fitted template to obtain flux densities within the aperture for each filter, and these values are reported in Table 1. These values were then subtracted from total flux densities.

The 2020 August observations with NOT were interrupted by the forest fire on MJD 59084. The data from MJD 59083 have lower signal to noise ratio and gradients in the background, resulting in larger than usual reported uncertainties in our analysis.

2.7 Radio

Between 2017 January and 2019 January, QSO B0218+357 was frequently observed with KaVA at 22 and 43 GHz. A total of 16 sessions were performed during this period. In most cases each session lasted two consecutive nights, with a 5–8-h track at 22 GHz on the first day and a similar track at 43 GHz on the following day. By default seven stations (three from KVN and four from VERA) joined each session. However, occasionally VERA-Mizusawa or VERA-Ishigaki was missing due to local issues. In addition, triggered by the 2020 August campaign, KaVA performed follow-up observations at 43 GHz for a total of nine sessions between 2020 August 22 and October 8. The observing time of each follow-up session was 3.5–4 h and on average five to six stations joined. All the data were recorded at 1 Gbps (a total bandwidth of 256 MHz with eight 32-MHz subbands) with left-hand circular polarization and correlated by the Daejeon hardware correlator (Lee et al. 2015). The initial data calibration (amplitude, phase, bandpass) was performed using the National Radio Astronomy Observatory (NRAO) Astronomical Image Processing System (AIPS; Greisen 2003) based on the standard KaVA/VLBI data reduction procedures (Niinuma et al. 2014; Hada et al. 2017). Imaging was performed using difmap software (Shepherd, Pearson & Taylor 1994) with the standard CLEAN and self-calibration procedures.

During the 4 KaVA 43 GHz sessions made on 2017 January 15th (MJD 57768), 2017 October 17th (MJD 58043), 2017 November 12th (MJD 58069), and 2018 January 5th (MJD 58123), additional simultaneous observations of the source with the KVN-only array with a 43 GHz/86 GHz dual-frequency recording mode were carried out. A wideband 4 Gbps mode was used where each frequency band was recorded at 2 Gbps (a bandwidth of 512 MHz for each band). 86 GHz fringes were detected by transferring the solutions derived at 43 GHz using the frequency-phase transfer (FPT) technique (e.g. Algaba et al. 2015; Zhao et al. 2019). Imaging was carried out in difmap.

Some of the KaVA/KVN results are presented in Hada et al. (2020) together with detailed radio images and analyses at each frequency. See Hada et al. (2020) for full details of the KaVA/KVN data reduction and imaging procedures. The typical angular resolution of KaVA (a maximum baseline length D = 2300 km) is 1.2 mas (22 GHz) and 0.6 mas (43 GHz), that of KVN (D = 560 km) is 1 mas at 86 GHz. Here we report on the whole data set, and investigate the kinematics of the jet.

The OVRO 40-m telescope uses off-axis dual-beam optics and a cryogenic receiver with 2 GHz equivalent noise bandwidth centred at 15 GHz. The double switching technique (Readhead et al. 1989), where the observations are conducted in an ON–ON fashion so that one of the beams is always pointed on the source, was used to remove gain fluctuations and atmospheric and ground contributions. Until 2014 May a Dicke switch was used to alternate rapidly between the two beams. Since 2014 May a 180 degree phase switch has been used, with a new pseudo-correlation receiver. Gain drifts were compensated with a calibration relative to a temperature-stabilized noise diode. The primary flux density calibrator was 3C 286 with an assumed value of 3.44 Jy (Baars et al. 1977). DR21 was used as secondary calibrator. Richards et al. (2011) describe the observations and data reductions in detail. Since the telescope is a single dish, it measures total flux densities integrated over the whole lensed structure: A, B and the Einstein ring.

The 37 GHz observations taken during the 2020 August campaign were made with the 13.7 m diameter Metsähovi radio telescope. The observations were ON–ON observations, alternating the source and the sky in each feed horn. A typical integration time to obtain one flux density data point was between 1200 and 1400 s. The detection limit of the telescope at 37 GHz was of the order of 0.2 Jy under optimal conditions. Data points with a signal-to-noise ratio <4 were treated as non-detections. The flux density scale was set by observations of DR 21. Sources NGC 7027, 3C 274, and 3C 84 were 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.

3 RESULTS

The MWL light curves measured during the monitoring campaign are presented in Fig. 1.

Figure 1.

MWL light curve of QSO B0218+357 between 2016 January and 2020 January. From top to bottom: MAGIC flux above 100 GeV, Fermi-LAT flux above 0.1 GeV, Fermi-LAT spectral index, X-ray flux in 0.3–10 keV range (corrected for Galactic absorption) measured with Swift-XRT and XMM–Newton, U-band observations from Swift-XRT-UVOT and XMM-OM, B-band observation from NOT, optical observations in R band with KVA and NOT. KaVA VLBI observations at 22 GHz (filled symbols show A image, empty ones B image) and 86 GHz (sum of A and B images shown with stars). OVRO monitoring results at 15 GHz. Flux upper limits are shown with downward triangles. Optical data are corrected for the host/lens galaxy contribution and galactic absorption. The points in red are contemporaneous (within 24 h slot) with MAGIC observations. The grey filled regions mark the enhanced emission periods F1 and F2.

Open in new tabDownload slide

3.1 Search for VHE emission

No significant VHE gamma-ray emission was found in the total data set of MAGIC monitoring data (see the left-hand panel of Fig. 2). Due to expected variability of the emission an additional analysis separating the data set into individual nights was performed. The distribution of the significances of the measured excess is shown in the right-hand panel of Fig. 2, and the upper limits on the flux above 100 GeV are reported in the top panel of Fig. 1. As the source is a known VHE gamma-ray emitter, we also report the nominal flux values on each observation night, however, none of them is significant, comparing to the corresponding uncertainty bar. The distribution is consistent with the lack of a measurable gamma-ray excess. By using the Fermi-LAT data, an additional study was performed to evaluate the expected VHE gamma-ray flux on individual nights (see Appendix  A), however, no clear hard GeV states could be identified.

Figure 2.

Left-hand panel: Distribution of the squared angular distance between the nominal and reconstructed source position of the total data set of MAGIC data (points) and corresponding background estimation (shaded region). Right-hand panel: Distribution of significances of individual nights of MAGIC observations, the red line shows the expected distribution for lack of significant emission, i.e. a Gaussian distribution with a mean of 0 and standard deviation of 1.

Open in new tabDownload slide

The SED upper limits were computed from the total monitoring sample of the MAGIC observations and compared with the extrapolation of the Fermi-LAT SED (see Fig. 3). The Fermi-LAT data for this comparison are quasi-simultaneous, i.e. 24 h-long time windows centred on each MAGIC observations are stacked together. The MAGIC upper limits are an order of magnitude below the flux observed during the flare in 2014 (Ahnen et al. 2016). However, within the uncertainties of the Fermi-LAT extrapolation the upper limits are in agreement with a power-law SED from GeV to sub-TeV range.

Figure 3.

Comparison of the MAGIC SED upper limits (black empty squares) with the extrapolation (black line and shaded region) of the Fermi-LAT spectrum (black filled circles). The extrapolation assumes a power-law behaviour and an extragalactic background light (EBL) absorption following Domínguez et al. (2011) model. For comparison, the Fermi-LAT and MAGIC SED during the 2014 flare (Ahnen et al. 2016) are shown with grey markers.

Open in new tabDownload slide

In order to constrain the VHE gamma-ray duty cycle of the source, the nights with optimal exposure were selected. The data set contains 37 nights with exposure >100 GeV of at least |$1.4\times 10^{12}\, \mathrm{cm^{2}s}$| (corresponding to about 1 h of observation with a typical effective area of |$4\times 10^{8}\, \mathrm{cm^{2}}$|⁠). All but one of those nights provide 95 per cent C.L. flux upper limits stronger than the VHE gamma-ray flux observed during the 2014 flare (⁠|$5.8\times 10^{-11} \, \mathrm{cm^{-2} s^{-1}}$|⁠). Following the 2014 event, a duration of an individual flare of at least 2 d was assumed. Using Monte Carlo simulations, we related the assumed rate of flares with the corresponding probability of at least one of them being caught in the observation slots of MAGIC. We found that the VHE duty cycle is consistent with less than 16 flares per year at 95 per cent C.L. with a flux >100 GeV of at least |$5.8\times 10^{-11} \, \mathrm{cm^{-2} s^{-1}}$|⁠.

3.2 Enhanced emission periods

Flux variations were detected across different energy bands during the 4-yr-long multi-instrument observations of QSO B0218+357 (see Fig. 1 and Table 2). Enhanced GeV emission was observed by Fermi-LAT around MJD ∼ 57650. The rise of the GeV emission was gradual. In our study, the period from MJD 57600 to MJD 57700 (dubbed as F1) was selected, which covers a time interval where the GeV flux increases and decreases (see Fig 1). Based on the obtained light curve, the resulting spectrum reported in this manuscript is not expected to depend strongly on the exact definition of the start and end of this time interval, and that similar results would have been obtained by modifying this time interval by a few days. During this time interval, three optical measurements (MJD 57627.2, 57639.2, and 57640.2) yielded a flux nearly an order of magnitude larger than that of the low state of the source. Comparing to the 2014 flare discussed in Ahnen et al. (2016), the GeV emission is at a similar level (however, with a softer spectrum), but the optical emission is nearly an order of magnitude higher. Interestingly the first two points are separated by 12 d, similar to the time delay between the two lensed images of the source. The optical flux density between those two measurements returned to the low state level. However, due to poor optical sampling of the source, the hypothesis that the two optical flares are indeed the two images of the same flare cannot be validated. Two of the nights of the enhanced optical activity had simultaneous MAGIC data. No significant emission was observed (significances of −0.21σ and −0.54σ). The flux upper limit above 100 GeV is ∼3 × 10⁻¹¹ cm⁻² s⁻¹, i.e. constrained to be at least two times below the level observed during the 2014 flare. Previously, gravitational lensing was used to predict possible sites of gamma-ray flares detected in 2012 and 2014 (Barnacka et al. 2016). They entertained a hypothesis that if both flares were caused by the same plasmoid created in the vicinity of a supermassive black hole and travelling towards the radio core then the interaction of the plasmoid with the radio core should be observed around 2016 July. While OVRO monitoring at 15 GHz does not show a significant increase in flux density (similarly to Spingola et al 2016), the predicted interaction coincides with the beginning of F1 in gamma-rays, and available observations in R band show a significant increase in flux density during the predicted interaction of the plasmoid with the radio core. This example illustrates the complexity of studying emissions from these sources but also points to a unique potential and the importance of long-term monitoring of QSO B0218+357 to elucidate the MWL origin of emission.

A hint of enhanced activity was observed in the X-ray band by XMM–Newton on MJD 58863.7 (dubbed as F2). The flux density increased by (44 ± 19) per cent with respect to the previous XMM–Newton measurements. The contemporaneous MAGIC observations did not yield any significant detection (excess at the significance of 2.1σ). These observations were used to derive 95 per cent C.L. upper limit on the flux of ∼2.8 × 10⁻¹¹cm⁻²s⁻¹ above 100 GeV, which is two times below the VHE flux measured during the flare in 2014. No significant excess is observed in the MAGIC observations in the two neighbouring nights further suggesting that the marginal excess in MAGIC data during the X-ray flare is a background fluctuation. No excess of GeV flux was observed during the X-ray flare. Interestingly, while the optical flux density did not change considerably during F2, a hint of increase in the UV flux density by (70 ± 41) per cent was found.

3.3 August 2020 campaign

Besides the multi-instrument observations described above, QSO B0218+357 is one of the sources that are regularly checked for GeV flares in the Fermi-LAT data, and additional observations are organized if flares or hints of flares are found. On MJD 59071.502 a photon with estimated energy of 59.4 GeV was observed from the vicinity of the source. Additionally, quasi-simultaneous Swift-XRT observations performed on MJD 59069.566 as well as TJO on MJD 59070.993 showed hints (at ∼2.5σ level) of enhanced flux density. Immediate follow up with the MAGIC telescopes was not feasible, due to the presence of bright moonlight, which would have substantially increased the energy threshold of the observations. Instead, similarly to the 2014 event, an MWL campaign was organized at the expected arrival of the B image of the flare, at the assumption that the observed one was the A image.

The observations are summarized in Fig. 4. A second HE photon with energy of 20 GeV was observed on MJD 59082.826. The association probability of each photon above 10 GeV was assessed with QSO B0218+357 using the standard Fermi-LAT tool (gtsrcprob). The whole data set was divided in the four point spread functions (PSF) classes.¹⁵ This is needed since there are relevant PSFs differences between the four classes and those differences have an important effect on the association probability. During the August observations, only these two photons above 10 GeV were detected with an association probability higher than 80 per cent. The probability of the two photons being associated with the source is 99.91 per cent and 86.88 per cent for the first and the second photon, respectively. The time difference of these two photons is 11.324 d, which is curiously consistent with the previously measured delay between the two images.

Figure 4.

MWL light curve of QSO B0218+357 during the 2020 August MWL campaign. Vertical lines: MAGIC observation nights (red) and Fermi-LAT >10 GeV photons (blue).

Open in new tabDownload slide

In order to evaluate statistical chance probability of occurrence of HE photons close in time with the emission of the source in a broader time-scale is investigated. 120 such photons spanning the total observations of the source by Fermi-LAT (MJD = 54683–59299) are obtained. Conservatively assuming the time window for the second photon of ±1 d (motivated by the spread of radio delays of ∼10–12 d) a chance probability of 5.2 per cent is obtained. The time delay of 11.324 d is also within the 1σ uncertainty of that measured from 2012 Fermi-LAT high state (11.46 ± 0.16 d; Cheung et al. 2014). For such a narrow window the corresponding chance probability is |$0.83{{\ \rm per\ cent}}$|⁠.

In the optical range a hint of increase of the R-band emission (2.5σ difference to the previous point and 4.2σ to the next) occurred close to the arrival time of the first Fermi-LAT photon. The observations during the planned monitoring at MJD = 59081–59085 were performed with higher cadence with additional instruments (NOT and HCT). The NOT measurement at MJD = 59083.1 is 2.5σ above the previous HCT observation, 2.4σ above the following HCT measurement. Similarly, comparing the enhanced NOT point to the previous NOT measurement at MJD = 59082.2, the difference shows a similar weak hint of 2.4σ. The B-band flux density increase was not accompanied by a similar R-band increase at the expected time of arrival of the second component of the flare. This would favour the interpretation of the B-band increase as a statistical fluctuation.

While some small hint of variability (2.3σ difference) can be seen between the first two Swift-XRT points, no variability can be seen during the expected time of arrival of the delayed component. This is understandable since the delayed component is expected to have a lower flux density at the peak, hence any small variability present in the leading emission can be easily missed in the trailing one.

In Fig. 5, radio follow-up light curves of the 2020 August campaign obtained with OVRO at 15 GHz and KaVA at 43 GHz is shown. For KaVA data in which A and B are spatially separated, the radio flux density for the core (the brightest component) was measured in each of A and B images. The core of A in 2020 August is found to be significantly brighter than the average flux density level in 2017–2018, and subsequently it shows continuous decrease in flux density at least until the end of our follow-up period (October 8), where the 43 GHz core flux density reaches a lowest level since the start of our KaVA monitoring from 2017. In contrast, we caution that the observed light curve of the weak component B is less defined because the observing conditions of KaVA follow-up sessions in 2020 were generally severe compared to the regular sessions in 2017–2018 (shorter integration time and smaller number of stations). In Fig. 5, there may be a significant amount of missing flux density in the light curve for B. This prevents us from cross-correlating the light curves of A and B. The 15 GHz OVRO monitoring did not cover the period in which the 43 GHz flux density decreased. The OVRO data 3 weeks before and 2 weeks after the KaVA flux density minimum show consistent flux densities. Metsähovi data cover the period in which the flux density measured by KaVA starts to decay, however, higher uncertainties and a larger integration region prevents sensitive probing of variability.

Figure 5.

Radio follow-up light curve of the 2020 August campaign of QSO B0218+357 with OVRO (black points) and KaVA (red filled and blue empty circles for core image A and B, respectively). Black squares show the significant 37 GHz Metsähovi flux densities (empty squares show the times of observations that resulted in signal-to-noise ratio below 4). The dotted lines show the average flux density from the monitoring period between 2017 January and 2018 December. The vertical blue lines show the times of arrival of the two HE Fermi-LAT photons during the 2020 August campaign.

Open in new tabDownload slide

3.4 Radio jet image

In Fig. 6, the evolution of the radio images of the source between 2017 January and 2018 December is presented. Here KaVA 43 GHz images are shown since the highest angular resolution is available. In addition to the core, a strong jet component is clearly detected (at ∼1.5/1.3 mas from the core in A/B, respectively, which corresponds to the projected distance between the two components of the order of 10 pc) in all the images. No additional knots have been observed throughout the observations. The jet direction is different in Image A and Image B, however, this is a geometrical effect caused by the lensing.

Figure 6.

Selected KaVA images at 43 GHz spanning 2 yr period. Top panel and bottom panels are images A and B of the source. For all images, contours start from −1 (dotted lines), 1, 2,... times 1.8 mJy beam⁻¹ (approximately 3σ) and increase by factors of 2^1/2. Blue and red dots represent the average (over all the images) position of the core and jet components. Two cyan points in image A represent the average positions of the sideway components, which were obtained based on five epochs where the sideway extension was clearly detected. For all the points, the position uncertainties were estimated based on the scatter from multiple epochs. Grey circle represents the smoothing kernel of the image.

Open in new tabDownload slide

In order to examine the kinematics of the jet component, the core and jet in the KaVA 43 GHz images were fitted by two 2D elliptical Gaussians, and the core-jet separation was measured as function of time. As can be seen in Fig. 6, the core-jet angular separations for both A and B images are quite stable over ∼2 yr of our observing period. Assuming that the same jet component is traced over 2 yr, a simple linear fit to the data is performed. Best-fitting values of 0.06 ± 0.03 and 0.04 ± 0.04 mas yr⁻¹ are obtained for Jet-A and Jet-B, respectively (corresponding to 3.1 ± 1.5c and 2 ± 2c).

In addition to the bright core and jet, the KaVA images of A also reveal diffuse extension in the direction perpendicular to the jet (see also Biggs et al. 2003; Hada et al. 2020). Two Gaussian models were additionally fitted to each KaVA 43 GHz image of A to characterize the positions of these extended structures. Since the sideways structures are generally diffuse and weak, reasonable fitting results were obtained on these features only for five epochs (MJD = 58069, 58123, 58434, 58450, 58476; 2017 November 12, 2018 January 5, November 12, November 28 and December 24) where the image quality was relatively high and SNR ∼ 4–10 were obtained for the fitted sideways components. In Fig. 6, the positions of the these additional features (averaged over the five epochs) are plotted in cyan. With respect to the core, the apparent ‘opening angle’ of this perpendicular extension (the angle between two vectors from the red point to each cyan point in Fig. 6) is estimated to be ∼62° (if only the extension of the main jet is considered, the opening angle is roughly a half of this). Possible proper motions in these components were also searched for but both of the components are essentially stationary, similar to the main jet component.

Regarding the mas-scale jet morphology during the 2020 August campaign, higher noise levels in the KaVA images (due to shorter integration time, smaller number of stations, higher humidity in the summer season) than in the 2017–2018 sessions made our image analysis more challenging (especially for the image B). Nevertheless, the overall radio morphology was quite similar to that in 2017–2018 (Fig. 6). While the core of A in 2020 August was significantly brighter than the average flux density level in 2017–2018 (see Fig. 5), no clear ejection of new components from the core during our KaVA follow-up period was found.

3.5 Lens geometry model

The observations of QSO B0218+357 with the Hubble Space Telescope (HST) show that the lensing galaxy is isolated pointing to a simple gravitational lens potential. The mass distribution of the lens has been shown to be well represented by a Singular Isothermal Sphere (SIS) model (Wucknitz, Biggs & Browne 2004; York et al. 2005; Larchenkova, Lutovinov & Lyskova 2011; Barnacka et al. 2016). The SIS model of QSO B0218+357 predicts a time delay of ∼10 d and a magnification ratio of ∼3.6. The predicted magnification ratio fits the observed ratio between radio images. However, the time delay is ∼1 d shorter as compared to the time delay measured at gamma-rays (Cheung et al. 2014; Barnacka et al. 2016).

Hada et al. (2020) used a Singular Elliptical Power-law (SEP) model and parameters fixed to the values obtained by Wucknitz et al. (2004). The SEP model provides a higher rate of change in the lens potential. As a result, the model with the same image positions predicts longer time delays and a larger magnification ratio between the images. The SEP model adapted by Hada et al. (2020) predicts a time delay of 11.6 d for the core, which matches better the observed time delay at gamma-rays. However, the predicted magnification ratio is |$\sim \,$|4.8, which significantly deviates from a reported average value of 3.5 from a broad range of radio observations (Patnaik, Porcas & Browne 1995).

At first parameters of the lens model presented in Barnacka et al. (2016) are adopted, which reconstructed the positions of the lensed images with an accuracy of |$1\,$|mas. The previous model was based on 15 GHz radio observations (Patnaik et al. 1995). Here, the reconstruction of the lens model is further improved by using high-resolution KaVA observations of the radio core and jet listed in Table 3.

A softened power-law potential (Keeton 2001; Barnacka et al. 2016) is used, which includes an Einstein radius, a scale radius of a flat core (s in mas), a projected axial ratio (q), and a power-law exponent (α). Similarly, like in Barnacka et al. (2016), the best lens model of softened power-law potential resulted in s ∼ 0, q ∼ 1, and α ∼ 1. Thus, the model reduced to the SIS potential.

Fig. 7 shows the radio observations with the prediction of the SIS model including the Fermat potential, predicted positions of the lensed images, as well as reconstructed position of the jet and core. The image is centred at the position of the lens located at (244.55,100.85) with respect to image A.

Figure 7.

Compilation of observations and lens model predictions. The colour map shows the Fermat potential of the lens model for the best reconstructed position of the radio core. The images form at the extreme of the Fermat surface (blue). The grey contours show the lensed images of the core and the Einstein ring observed at 1.687 GHz (Wucknitz et al. 2004). The image is centred at the reconstructed position of the lens indicated as the blue asterisk. The red open circles correspond to the reconstructed position of the images of the jet (A on the right, and B on the left). The red dash–dotted line connects the images. For the SIS lens model, the source (green asterisk) is located at half distance between images. The observed and reconstructed images of the radio core are also shown, however, they are superimposed with the red open circles.

Open in new tabDownload slide

Only positions of the radio images were used to find the best fit, as the observed time delays and magnification ratios are a subject of discussion. Fig. 8 shows the position of reconstructed images as red and green open circles. The model reconstructs the positions of the lensed images with average accuracy of 0.03 mas for the radio core and 0.15 mas for the jet.

Figure 8.

Comparison of the KaVA 43 GHz image taken on MJD = 54470 (2018 January 5), A in left-hand panel, B in right-hand panel compared with the reconstructed positions of the lens model images (stars).

Open in new tabDownload slide

The projected distance between reconstructed positions of the core and jet is |$1.2\pm 0.1\,$|mas, which correspond to |$9.8\,$| pc in the source plane, consistent with the result obtained by Hada et al. (2020). The SIS model predicts a time delay of 10.36 d for the core and 10.30 d for the jet. The shorter time delay for the jet indicates that the radio jet is positioned towards the lens centre as argued by Barnacka et al. (2016) based on observations of the Einstein radius at 1.687 GHz (Wucknitz et al. 2004). The Einstein ring forms from emission of the kpc-scale jet aligned with the lens centre. The Einstein ring observed at 1.687 GHz (Wucknitz et al. 2004) is also added to Fig. 7.

The predicted time delays reported here are calculated using |$H_0 = 67.3\pm \mathrm{1.2\, km\, s^{-1}\, Mpc^{-1}}$| (Planck Collaboration 2014). Note that the time delay is inversely proportional to H₀. Thus, larger values of |$H_0=73.3^{+1.7}_{-1.8}\mathrm{\, km\, s^{-1}\, Mpc^{-1}}$| as reported by Wong et al. (2020) would result in an even shorter time delay, and as such would further increase the discrepancy between the lens models and observations. The Hubble parameter estimation based on gravitationally induced time delays of variable sources with relativistic jets is in particular prone to biases as the emission can originate from multiple sites, which can introduce systematical errors (Barnacka et al. 2015).

Both the SEP (Hada et al. 2020) and improved SIS reconstruction based on KaVA observations, provide accurate reconstruction of the positions of the lensed images and consistent distance separation of the radio core and jet components. However, the predictions of these two models differ in terms of time delay and magnification ratios. The SIS scenario predicts magnification ratios of 3.81 and 3.67 for the core and jet, respectively. The flux densities of the lensed images at 43 GHz reported in Table 3 (Hada et al. 2020) indicate magnification ratios of 3.7 and 3.2 for the core and jet, respectively. For comparison, the SEP model predicts magnification ratios of 5 and 4.84 for the core and jet, respectively.

In principle, magnification ratio could be affected by microlensing or ‘substructure lensing’ such as compact clumps in the lens galaxy (Sitarek & Bednarek 2016). To match the prediction of the magnification ratio of ∼5 of the SEP model with the observed ∼3.5, either the brighter image A would have to be magnified by a factor of 1.35, or image B would have to be demagnified by the same factor. Moreover, the lensed images of both jet and core would have to be magnified/demagnified simultaneously by a similar factor. As a result, the diameter of the perturbing mass needs to be |$\gg 10\,$|pc, precluding microlensing. In principle, possible clumps of tens of pc of a giant molecular cloud (GMC; see e.g. Chevance et al. 2020) could result in additional moderate amplifications by a factor of 1.5 (see Sitarek & Bednarek 2016).

Interestingly, the observed magnification ratio fits the prediction of the SIS model well. The SIS model predicts correctly not only the values but also the degree to which the magnification ratio of the core is greater than the magnification ratio of the jet. However, the time delay of 10.3 d predicted by the SIS model is one day shorter than the time delay measured at gamma-rays. Such a discrepancy in the time delay could be explained if there is an offset of |$\sim 50\, \mathrm{ pc}$| between sites of radio and gamma-ray emission (for a review, see Barnacka 2018).

The degeneracy between the SIS and SEP lens models could be broken by precise measurement of the radio time delay. The time delay obtained in the SIS model depends only on H₀ and the image angular separation. Thus, the SIS model can be ruled out if the radio time delay is not 10.3 d.

However, measuring time delay at radio with an accuracy of hours is difficult as radio observation of B2 0218+35 shows almost no variability, in addition to gaps between the observations. Further high-resolution KaVA observations combined with detailed modelling of the lens could elucidate the true potential of the lens and test the clump scenario. KaVA observations provide well-resolved images of both the jet and core, thus providing multiple tests of the lens potential. One of the predictions of the SIS model is that the diameter of the Einstein ring is equal to the distance between the lensed images of the source. The current KaVA observations show that the distance between the lensed images of both the core and jet are equally within the range of uncertainty of the observations – as such, consistent with the SIS model. More precise observations, or detailed predictions of the SEP model on the expected difference between the lensed images of the core and jet could help exclude the SEP model. Moreover, elliptical lens models should result in formation of the odd number of images (Gottlieb 1994; Zhang et al. 2007; Petters & Werner 2010). Thus, detection of the third image in the vicinity of the predicted lens centre could provide further constraints on the model of the lens.

Here, we focused on the two most general lens models, namely SIS and SEP, and on the possibility to break the degeneracy between them. However, a potentially more general class of models might be required to reconstruct both the time delay and the magnification ratio.

The in-depth observations of the object B2 0218+35, combined with a precise model of the lens, have the potential to provide unique insights on the origin and site of the gamma-ray emission, the Hubble constant, or substructures in the lensing galaxy.

3.6 Modelling of dust in the lens

A Galactic absorption of N_H = 5.56 × 10²⁰ cm⁻² was adopted from the Leiden/Argentine/Bonn (LAB) survey (Kalberla et al. 2005). The X-ray flux density, corrected for the Galaxy absorption, in the (0.3–10) keV band is f = (1.53 ± 0.11)× 10⁻¹² erg cm⁻² s⁻¹ for the low flux density state, and f = (2.25 ± 0.24)× 10⁻¹² erg cm⁻² s⁻¹ for the high-flux density state from XMM–Newton observations.

In order to evaluate and correct the effect of additional absorption in the host or lens galaxies an approach similar to Ahnen et al. (2016) is applied. The higher sensitivity of XMM–Newton compared to the Swift-XRT telescope allows us to study a few alternative models of absorption and intrinsic source spectrum and select between them. For the investigations of the low state three cases are considered: (i) no additional absorption, (ii) absorption at the host, and (iii) absorption at the lens. In the case (ii) the absorption will affect the total emission observed from the source (in both images). In the case (iii), since the two images cross different parts of the lens galaxy, the absorption would be different for them. There are reasons to believe that in such a situation the absorption would mainly affect the brighter, A, image (see the discussion in Ahnen et al. 2016). Therefore in the case (iii) the observed emission is assumed to be composed of two ‘virtual’ sources located at the nominal location of QSO B0218+357 in the sky. The first component is affected only by Galactic absorption, and the second one is additionally absorbed by a hydrogen column density (N_{H, z}) at the redshift of the lens (z = 0.68).

Two spectral models are investigated: a simple power law and a log parabola (defined as F(E) = k(E/E₀)^−Γ and |$F(E)=k(E/E_0)^{-(\alpha +\beta \log (E/E_0))}$|⁠, where E₀  = 1 keV in all the models). For the absorption, the Tuebingen–Boulder interstellar medium absorption model (Wilms, Allen & McCray 2000) available in xspec was adopted. The column density N_{H, z} and the parameters of the log-parabola or power-law components were fitted by imposing that the values of α, β, or Γ are the same for the two components, and that the normalization of the component with only Galactic absorption is a factor 0.71/2.72 = 0.261 (corresponding to magnification ratio, see Section 3.5) lower than the normalization of the component with also internal absorption.

The results of the fits are summarized in Table 4. The simple log-parabola model, without an additional absorption is not sufficient to explain the XMM–Newton data (χ²/N_dof = 413.4/374). Including an additional absorption at the lens for the brighter image the χ² improves by 38.4 at the cost of one additional degree of freedom (hydrogen column density at the lensed image A). The model is compared with observed XMM–Newton rates in Fig. 9.

Figure 9.

Differential energy flux of QSO B0218+357 folded with the response of XMM–Newton observed from low-state pointings (top panel) and from the X-ray flare (bottom panel). Points show the observed rate, while lines show the LP model for image A (blue dashed), image B (green dot–dashed), and total emission (red solid).

Open in new tabDownload slide

Since not all the investigated models are nested, to compare them Akaike Information Criterion (Akaike 1974) is used. For the case of χ² statistics, the relative difference of AIC parameter of two models can be computed as ΔAIC = 2Δn_p + Δχ². The relative likelihood of the models can be computed (see e.g. Burnham, Anderson & Huyvaert 2011) as p = exp (ΔAIC/2). Including the absorption at z = 0.68 the intrinsic curvature of the emission is preferred. The difference of χ² values is 11.6 for one additional parameter (describing the curvature) corresponds to p = 0.008 relative likelihood of the power-law model to the log parabola model. On the other hand, comparing models with absorption at z = 0.68 and z = 0.94, with the same number of free parameters in both models, the former has a lower χ² value by 23.9. Therefore, the model with absorption at the source is only 1.8 × 10⁻⁵ as likely as the model with the absorption at the lens. Summarizing, for the low state of the source, the preferred model of the emission involves an intrinsic log-parabola spectrum and absorption by a column density of (8.10 ± 0.93_stat) × 10²¹cm⁻² at the lens.

Comparing to the result obtained in Ahnen et al. (2016) 24 ± 5_stat × 10²¹cm⁻², the absorption obtained in this work is more precise, but also lower, and both values are consistent in the broad range derived by Menten & Reid (1996) (5–50 × 10²¹cm⁻²). The actual absorption in the lens might have changed if the region emitting X-rays has moved along the jet, or changed its size compared to the observations in 2014. However it is equally likely that additional systematics (assumption of the intrinsic spectral model, flux density and spectral variability, absorption in the other image of the lens and in the host galaxy) affected one or the other measurement.

The proper correction for the absorption and lensing during the MJD 58863.7 high X-ray state is more uncertain. Since the observations are separated by 140 d from the previous X-ray flux measurement, is not clear if the X-ray high state was a short duration flare, or a longer time-scale high state. In a case of a short flare the observations might have happened when the corresponding image A or image B reached the observer, resulting in a different absorption. On the other hand, if the enhanced state was significantly longer than the ∼11 d delay between the two images the observed emission should be the average from both images, similar to the case of the low state. This is further supported by the fact that the data collected by Swift-XRT for MJD 59069–59108 show also higher X-ray fluxes, similar to the last XMM–Newton measurement. To analyse MJD 58863.7 data of XMM–Newton the assumption that the observed increase in the X-ray emission was over a longer time-scale is applied, therefore the log-parabola spectral model was used with absorption of the brighter image and fixed flux ratio between both components. Such a model describes sufficiently well the observations (χ²/N_dof = 90.5/94), and provides a somewhat harder and more curved X-ray spectrum than during the low state. The derived N_{H, z} is roughly consistent (at 1.3σ level) with the values obtained from the low state fit. An alternative model has been tested in which the high state emission is a sum of the low state emission, with the spectral shape and flux fixed to the low-state-fitted values and an additional flaring component with an absorption at z = 0.68 (i.e. the flare originating from the brighter image A). The fit results for such an additional flaring component are reported in ‘LP + low’ row of Table 4. The two investigated models for the flaring state are not nested, however, they have the same number of free parameters and result in nearly the same χ², therefore neither can be rejected. For the sake of simplicity the same absorption model is used for the flaring state as for the low state.

As discussed in Section 2.4, two combined spectra are used for the spectral study using Swift-XRT data during 2016–2017 and 2020 August. Due to the restriction aroused by Swift-XRT sensitivity, only two spectral models within the scenario of case (iii) are investigated. In addition to the assumption implemented for XMM–Newton observations, N_{H, z} is also assumed to be equal to 8.10 × 10²¹ cm⁻². The log-parabola model cannot describe the observed spectrum better than the power-law model. The results are presented in Table 4.

4 LOW STATE BROAD-BAND SED MODELLING

The MWL behaviour of the source in different states is summarized in Fig. 10. During the low state, the GeV emission was significantly lower than during the 2014 flare described in Ahnen et al. (2016), and even during the F1 enhanced state the GeV spectrum was softer than in 2014. Except for the strong optical flares in F1, the rest of the MWL SED does not vary strongly with respect to the 2014 flare and historical measurements. The difference in reported radio measurements to the historical ones is most likely caused by much smaller integration region (only the inner jet of the source) achieved in radio interferometry with KaVA.

Figure 10.

MWL SED of QSO B0218+357 in different periods: average state of the monitoring (excluding optical flare data) in red (for clarity Swift-XRT and MAGIC data are shown with empty symbols, while XMM–Newton and Fermi-LAT and the rest of MWL data are shown with full symbols), X-ray flare in cyan, and optical/GeV flare in green. For the case of optical flare and minimum and maximum optical flux density during investigated period are plotted. For comparison historical data (obtained from SSDC service) are plotted in grey, while the 2014 flare (Ahnen et al. 2016) is in black.

Open in new tabDownload slide

For the average state a detailed modelling is performed. The emission is modelled in the framework of an External Compton (EC) model, which is a common scenario for FSRQs. There is growing evidence that the main target for FSRQs EC process is the reprocessed Dust Torus (DT; see e.g. Costamante et al. 2018; van den Berg et al. 2019). Moreover, even while no VHE gamma-rays has been detected from the QSO B0218+357 during the monitoring period, the source is a known emitter in this band (Ahnen et al. 2016), again supporting EC-DT scenario.

The recent measurement of the accretion disc luminosity of |$L_d=4.3\times 10^{44}\mathrm{erg\, s^{-1}}$| (Paliya et al. 2021) is used. The value is close to the |$L_d=6\times 10^{44}\mathrm{erg\, s^{-1}}$| estimated by Ghisellini et al. (2010) and applied in Ahnen et al. (2016). Using the updated value of L_d the sizes of the BLR and DT are computed using the scaling laws of Ghisellini & Tavecchio (2009), resulting in R_BLR = 6.6 × 10¹⁶ cm, and R_DT = 1.6 × 10¹⁸ cm. The temperature of the DT is set to 1000 K and its luminosity to 0.6L_d. A conical jet geometry is used, with half-opening angle of 1/Γ, where Γ is the Lorentz factor of the jet.

For the modelling the Doppler factor of the jet D = Γ = 15 is assumed. The electron energy distribution (EED) is assumed to follow a power law with an index of p₁ up to γ_b, where γ_b is the Lorentz factor of the electrons for which the time-scale for the dominating energy loss process is equal to the dynamic scale (see Acciari et al. 2021 for details). Above such a cooling break the EED steepens by 1 up to γ_max, which is determined from balancing the acceleration gain with the dominating energy loss process. The radiation processes are calculated using the agnpy¹⁶ code (Nigro et al. 2020), which implements the synchrotron and Compton processes following the prescriptions described in Dermer & Menon (2009) and Finke (2016). While the γ_b and γ_max are calculated assuming Thompson regime of the inverse Compton scattering, the actual spectra are computed using the full Klein–Nishina cross-section formula.

The emission region (hereafter ‘Close’ region) responsible for the high-energy bump is assumed to be located at the distance of d₁ = 2 × 10¹⁷ cm, i.e. a factor of a few more distant than the size of the BLR but deep in the DT radiation field. The model is confronted with the observations in Fig. 11, taking into account the magnification induced by the lensing (using the strong lensing magnifications derived in Section 3.5), and the absorption of emission from one of the images in the lens. The possible effect of microlensing is not corrected for, however, we expect it to have a minor influence on the long-term average spectrum. The free and derived parameters are summarized in Table 5.

Figure 11.

MWL SED of QSO B0218+357 composed of contemporaneous data to the MAGIC observations (red points), historical data from SSDC (grey), and data from the 2014 flare (Ahnen et al. 2016, in black) compared with the broad-band model derived from a two-zone SSC + EC scenario (with parameters reported in Table 5). Optical, UV, and X-ray data are corrected for the Galactic absorption, optical data are in addition corrected for the host/lens galaxy contribution. The lensing magnification, absorption in the lens galaxy, and EBL attenuation are corrected for in the model curves. For the closer region, dotted curve is the synchrotron emission, dashed the SSC, dot–dashed EC on DT, dot–dot–dashed EC on BLR. For the farther region, long-dotted is the synchrotron emission. The total emission is shown with an orange line.

Open in new tabDownload slide

The gamma-ray emission is explained as EC process on DT photons (which is also the dominating energy loss process of the electrons). On the other hand, according to the model, the X-ray emission is mostly caused by SSC process. The synchrotron emission corresponding to the ‘Close’ region can (largely) explain the optical and the rapidly falling UV emission.

However, the region is too compact for explaining the low-frequency radio emission which is heavily absorbed in the ‘Close’ region by synchrotron-self-absorption. Such low-energy emission is expected to originate from a larger scale jet. A commonly applied solution is the assumption of two emission regions (see e.g. MAGIC Collaboration 2020 and references therein). Therefore, motivated also by the radio knot observed by KaVA, a second region (hereafter ‘Far’) is added, located at the distance of 100 pc. The distance of the emission region is motivated by (deprojected) distance of the jet component in the KaVA image. The low-energy slope in this region is set to 2, and equipartition (i.e. u_e = B²/(8π)) is applied. Then the magnetic field strength and the acceleration coefficient are fixed to the values explaining the broad-band synchrotron emission. The two emission regions are assumed not to be co-spatial (‘Far’ region is more distant in the jet) and thus contrary to e.g. MAGIC Collaboration (2020) are not interacting with each other. The ‘Far’ region is distant and large-enough such that the dominating energy loss process is the synchrotron cooling, which, again due to size of the region and low values of the acceleration efficiency, does not introduce a cooling breakup to the maximum reached energies. Inverse Compton emission of the ‘Far’ region is negligible in comparison to the ‘Close’ region (the region is beyond the DT radiation field for EC process to play any role, and the energy density of the electrons is too low for SSC to be effective).

Combination of both emission regions can describe remarkably well the whole broad-band emission of the source. In fact, the previous modelling of the source, presented in Ahnen et al. (2016), also suggested two-zone model for this source. That modelling was, however, used to explain the flaring episode of the source, and neglected the radio and microwave emission from the large-scale jet. It is therefore likely that the three regions contribute to the time-variable, broad-band emission of the source: large scale jet responsible for the radio and microwave emission, emission region within DT responsible for the broad-band, high-energy, low state emission of the source and a third region (or a sub-region of the second one) in which VHE and HE gamma-ray flares occur. As the low state modelling attributes most of the radio emission to the ‘Far’ region, it is curious to observe radio variability in 2020 August campaign KaVA data over time-scales of tens of days (see Fig. 4). Such variability might not be connected with the source itself, but rather with the absorption and milli-lensing effects of large scale structures in the lensing galaxy. In Appendix  B, a possible scenario is discussed for explaining the MWL flares seen from the source during the monitoring.

5 CONCLUSIONS

Broad-band (radio, optical–UV, X-ray, gamma-ray) monitoring of QSO B0218+357 has been performed. The monitoring was aimed at the detection of a VHE gamma-ray flare of the source in time periods selected to allow additional follow up at the expected time of arrival of the second image. The deep exposure of 72 h of data did not reveal low-state VHE gamma-ray emission and constrained it to be less than about of an order magnitude below the level observed during 2014 flare. VLBI radio images obtained with KaVA show clear core-jet structure in both lensed images. No significant movement of the VLBI radio features was seen. No significant variability has been seen in KaVA images during the 2016–2019 monitoring, however, the follow up of 2020 August campaign showed a clear decay of core flux density in image A. The radio data have been used to improve the lens modelling to evaluate image magnifications and time delays for the core and jet component of the source. Precise measurements of the X-ray spectrum with XMM–Newton instrument was used to evaluate the absorption in the lens, and fit the hydrogen column density of the lens in the image A to a value of (8.10 ± 0.93) × 10²¹cm⁻².

The low-state broad-band emission of the source can be described with a two-zone model, in which the EED shape is determined self-consistently from the cooling, acceleration, and dynamical time-scales. Most of the radio and far-infrared emission is explained to originate from a large region with size/location motivated by the radio jet component. UV (and partially optical) data are explained as the synchrotron emission of the smaller region, that is also responsible for the gamma-ray emission (produced in EC scenario) and X-ray emission (generated via SSC process).

No short VHE gamma-ray flares have been observed in the night-by-night analysis. Comparison with the Fermi-LAT state of the source shows that it is unlikely that the source has reached a comparable flare to the one of 2014 during the monitoring. The MAGIC data were used to place a 95 per cent C.L. limit on the VHE gamma-ray duty cycle of the source: below 16 flares per year.

Monitoring data have revealed, however, a few flares/hints of enhanced states in optical, X-ray, and gamma-rays, during which no VHE gamma-ray emission was detected. While the limited MWL data and variability during enhanced periods do not allow us to properly model the enhanced states, a plausible scenario explaining qualitatively the change of behaviour of the source during those states by change of the basic parameters of the model is presented.

Additional MWL campaign triggered by hints of enhanced emission in gamma-rays, X-rays, and optical has been also discussed. Unfortunately lack of MAGIC data on the predicted night of the flare prevents us from drawing a firm conclusion on possible hardening of the electron energy distribution during the campaign.

While the primary goal of the MWL monitoring of the source has not been achieved due to in general low gamma-ray activity of the source in the last years, the campaign resulted in multiple interesting results, and observations of a few interesting events. Since the achieved constraints on the low-state VHE gamma-ray emission approach the extrapolation of the GeV emission, it is expected that the future Cherenkov Telescope Array (Acharya et al. 2013) will allow us to study it in detail.

ACKNOWLEDGEMENTS

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 ERDF under the Spanish Ministerio de Ciencia e Innovación (MICINN) (FPA2017-87859-P, FPA2017-85668-P, FPA2017-82729-C6-5-R, FPA2017-90566-REDC, PID2019-104114RB-C31, PID2019-104114RB-C32, PID2019-105510GB-C31, PID2019-107847RB-C41, PID2019-107847RB-C42, PID2019-107988GB-C22); 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-268/16.12.2019 and the Academy of Finland grant nr. 320045 is gratefully acknowledged. This work was also supported by the Spanish Centro de Excelencia ‘Severo Ochoa’ SEV-2016-0588 and CEX2019-000920-S, and ‘María de Maeztu’ CEX2019-000918-M, the Unidad de Excelencia ‘María de Maeztu’ MDM-2015-0509-18-2 and the ‘la Caixa’ Foundation (fellowship LCF/BQ/PI18/11630012) and by the CERCA program of the Generalitat de Catalunya; by the Croatian Science Foundation (HrZZ) Project IP-2016-06-9782 and the University of Rijeka Project uniri-prirod-18-48; by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3; the Polish National Research Centre grant UMO-2016/22/M/ST9/00382; and by the Brazilian MCTIC, CNPq, and FAPERJ. The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515. We thank Director of Indian Institute of Astrophysics for allotting us observing time with HCT under DDT. We also thank the staff of IAO, Hanle and CREST, Hosakote, that made HCT observations possible. The facilities at IAO and CREST are operated by the Indian Institute of Astrophysics, Bangalore. The Joan Oró Telescope (TJO) of the Montsec Astronomical Observatory (OAdM) 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 which was supported in part by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G, and NSF grants AST-0808050 and AST-1109911, and private funding from Caltech and the MPIfR. SK acknowledges support from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme under grant agreement No. 771282. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. This publication makes use of data obtained at the Metsähovi Radio Observatory, operated by Aalto University in Finland. Part of this work is based on archival data, software or online services provided by the Space Science Data Center - ASI. We would like to thank the anonymous journal reviewer for his/her comments that helped to improve the manuscript.

DATA AVAILABILITY

The data used in this article were accessed from the MAGIC telescope http://vobs.magic.pic.es/fits/. The derived data generated in this research will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES