EDP Sciences
Press Release
Free Access
Issue
A&A
Volume 595, November 2016
Article Number A98
Number of page(s) 11
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201629461
Published online 04 November 2016

© ESO, 2016

1. Introduction

Even though there are already over 60 blazars detected in the very high energy (VHE 100 GeV) range, most of them are relatively close-by sources with redshift z ≲ 0.5. Until mid-2014, the farthest sources observed in this energy range were 3C 279 (z = 0.536, Albert et al. 2008), KUV 00311-1938 (z > 0.506, Becherini et al. 2012) and PKS1424+240 (z = 0.601, Acciari et al. 2010). In the last two years the MAGIC (Major Atmospheric Gamma Imaging Cherenkov) telescopes have discovered VHE gamma-ray emission from QSO B0218+357 at z = 0.944 (Mirzoyan 2014) and afterwards PKS1441+25 at z = 0.940 (Ahnen et al. 2015), almost doubling the boundaries of the known gamma-ray universe. Observations of distant sources in VHE gamma-rays are difficult owing to strong absorption in the extragalactic background light (EBL, see e.g. Gould & Schréder 1966). At a redshift of ~1 it results in a cut-off at energies1~100 GeV. Such energies are at the lower edge of the sensitivity range of the current generation of Imaging Atmospheric Cherenkov Telescopes (IACTs), making such observations challenging. To maximize the chance of detection, the observations are often triggered by a high state observed in lower energy ranges. In particular, Fermi Large Area Telescope (Fermi-LAT) scans the whole sky every 3 h and provides alerts on sources with high fluxes and information about the spectral shape of the emission in the GeV range.

QSO B0218+357, also known as S3 0218+35, is classified as a flat spectrum radio quasar (FSRQ, Ackermann et al. 2011). The classification is based on the optical spectrum (Cohen et al. 2003). It is located at a redshift of zs = 0.944 ± 0.002 (Cohen et al. 2003). One of the five features from which Cohen et al. (2003) derived the redshift was confirmed by (Landoni et al. 2015). The object is gravitationally lensed by the face-on spiral galaxy B0218+357 G located at a redshift of zl = 0.68466 ± 0.00004 (Carilli et al. 1993). Strong gravitational lensing forms multiple images of the source (see e.g. Kochanek et al. 2004). The flux magnification of an image is the ratio of the number of photons gravitationally deflected into a small solid angle centered on the observer to the number of photons emitted by the source in such a solid angle. The 22.4 GHz VLA radio image shows two distinct components with an angular separation of only 335 mas and an Einstein ring of a similar size (O’Dea et al. 1992). Observations of variability of the two radio components led to a measurement of a delay of 10–12 days between the leading and trailing images (Corbett et al. 1996; Biggs et al. 1999; Cohen et al. 2000; Eulaers & Magain 2011). In the radio image, the leading component (also called “image A” in the literature) is located to the west of the trailing component (image B). The delayed component had a 3.57–3.73 times weaker flux (Biggs et al. 1999). However, the observed ratio of magnification varies with the radio frequency (Mittal et al. 2006), presumably owing to free-free absorption in the lensing galaxy (Mittal et al. 2007). In the optical range the leading image is strongly absorbed (Falco et al. 1999).

In 2012 QSO B0218+357 went through a series of outbursts registered by the Fermi-LAT (Cheung et al. 2014). Even though Fermi-LAT does not have the necessary angular resolution to disentangle the two emission components, the statistical analysis of the light curve auto-correlation function led to a measurement of a time delay of 11.46 ± 0.16 days. Interestingly, the average magnification factor, contrary to radio measurements, was estimated to be ~1. Changes in the observed GeV magnification ratio were interpreted as microlensing effects on individual stars in the lensing galaxy (Vovk & Neronov 2016). Microlensing on larger scale structures has been considered as well (Sitarek & Bednarek 2016). The radio follow-up observations of QSO B0218+357 after the 2012 gamma-ray flare did not reveal any correlation between the two bands (Spingola et al. 2016).

Another flaring state of QSO B0218+357 was observed by Fermi-LAT on 2014 July 13 and 14 (Buson & Cheung 2014). Contrary to the results for the 2012 flaring period, in this case the ratio of the leading to delayed GeV emission was at least 4 (Buson et al. 2015). The 2014 flare triggered follow-up observations by the MAGIC telescopes, which in turn led to the discovery of VHE gamma-ray emission from QSO B0218+357 (Mirzoyan 2014).

In this work we present the results of the observations by the MAGIC telescopes and supporting multiwavelength instruments of the QSO B0218+357 during the flaring state in July 2014. In Sect. 2 we describe the instruments taking part in those observations and the data reduction. The effect of the lensing galaxy on the observed emission is discussed in Sect. 3. Section 4 is devoted to the results of the observations. In Sect. 5 we model the spectral energy distribution (SED) of the source. We use the Fermi-LAT and MAGIC observations to discuss constraints on the EBL in Sect. 6.

2. Instruments, observations and analysis

The VHE gamma-ray observations of QSO B0218+357 during the flaring state in July 2014 were performed with the MAGIC telescopes. The source was also monitored in GeV energies by Fermi-LAT, in X-ray by Swift-XRT and in optical by KVA.

2.1. MAGIC

MAGIC is a system of two 17 m Cherenkov telescopes located in the Canary Island of La Palma at a height of 2200 m a.s.l. The MAGIC telescopes combine large mirror area, allowing us to observe gamma rays with energies as low as ~50 GeV, with the stereoscopic technique providing strong hadronic background rejection, and hence good sensitivity at low energies. This makes them an excellent instrument for observations of distant FSRQs. In summer 2012 MAGIC finished a major upgrade (Aleksić et al. 2016a) greatly enhancing the performance of the instrument (Aleksić et al. 2016b). The sensitivity2 of the MAGIC telescopes achieved in the energy range 100 GeV is at the level of 1.45% of Crab Nebula flux in 50 h of observations. The angular resolution of MAGIC is of the order of 0.09°, i.e. insufficient for spatially resolving the emission from the two lensed image components of QSO B0218+357.

The telescopes could not immediately follow the flare alert published by Fermi-LAT in mid July 2014 from QSO B0218+357 as it occured during the full Moon time (the Cherenkov light from low energy gamma-ray showers would be hidden in the much larger noise from the scattered moonlight). The MAGIC observations started 10 days later, with the aim of studing the possible emission during the delayed flare component. The observations were performed during 14 consecutive nights from 2014 July 23 (MJD = 56 861, two nights before the expected delayed emission) to 2014 August 5 (MJD = 56 874). The total exposure time was 12.8 h and the source was observed at an intermediate zenith angle (20°–43°). The data reduction (stereo reconstruction, gamma/hadron separation and estimation of the energy and arrival direction of the primary particle) was performed using the standard analysis chain of MAGIC (Zanin et al. 2013; Aleksić et al. 2016b). The sky position of QSO B0218+357, contrary to the Crab Nebula used to estimate the MAGIC telescopes performance in Aleksić et al. (2016b), is not projected against the Milky Way optical background. A 30% smaller night sky background registered by the MAGIC telescopes for QSO B0218+357 allowed us to apply image cleaning thresholds lower by 15% with respect to the ones used in the standard analysis presented in Aleksić et al. (2016b). For the zenith angle range in which the observations were performed this resulted in the analysis energy threshold of about 85 GeV (measured as the peak of the Monte Carlo (MC) energy distribution3 for a source with the spectral shape of QSO B0218+357). The lower image cleaning thresholds were validated by applying the same procedure to the so-called pedestal events, i.e. events which contain only the light of the night sky and electronic noise. An acceptable fraction of about 10% of such images survived the image cleaning. The analysis was performed using a dedicated set of MC simulations of gamma rays with the night sky background and the trigger parameters tuned to reproduce as accurately as possible the actual observation conditions.

2.2. Fermi-LAT

Fermi-LAT is a pair-conversion telescope optimized for energies from 20 MeV to greater than 300 GeV (Atwood et al. 2009). Generally, Fermi-LAT is operated in scanning mode, providing coverage of the full sky every three hours. Starting on December 2013 and until December 2014, a new observing strategy that emphasized coverage of the Galactic center region was adopted. QSO B0218+357 data presented in this paper were obtained during this time interval. As a consequence, the coverage on the blazar position was on average a factor of 0.6 of the maximum one. Additionally, at the time of the expected delayed emission, Fermi performed a ToO (Target of Opportunity) observation on QSO B0218+357 to enhance exposure toward the source position. The ToO lasted approximately 2.7 days (2014 July 24, 00:30:01 UTC to 2014 July 26 18:24:00 UTC, MJD = 56 862.0256 864.77).

Fermi-LAT data were extracted from a circular region of interest (ROI) of 15° radius centered at the QSO B0218+357 radio position, , (J2000; Patnaik et al. 1992). The analysis was done in the energy range 0.1–300 GeV using the standard Fermi Science Tools (version v9r34p1) in combination with the P7REP_SOURCE_V15 LAT Instrument Response Functions. For obtaining the light curve, data collected between MJD = 56 849–56 875 (2014 July 11–2014 August 6) were used. For the spectral analysis only data spanning the two days during the flare observed by MAGIC (MJD = 56 863.125–56 864.5) were used. We applied the gtmktime filter (#3) cuts to the LAT data following the FSSC recommendations4. According to this prescription, time intervals when the LAT boresight was rocked with respect to the local zenith by more than 52° (usually for calibration purposes or to point at specific sources) and events with zenith angle >100° were excluded to limit the contamination from Earth limb photons.

The spectral model of the region included all sources located within the ROI with the spectral shapes and the initial parameters for the modeling set to those reported in the third Fermi-LAT source catalog (3FGL, Acero et al. 2015) as well as the isotropic (iso_source_v05.txt) and Galactic diffuse (gll_iem_v05.fit) components5. For generating the light curve the source of interest was modeled with a power-law spectral shape with normalization and index free to vary. To access the detection significance we used the Test Statistic (TS) value. It is defined as TS = −2log (L0/L), where L0 is the maximum likelihood value for a model without an additional source (the “null hypothesis”) and L is the maximum likelihood value for a model with the additional source at the specified location. The TS quantifies the probability of having a point gamma-ray source at the location specified and corresponds roughly to the square of the standard deviation assuming one degree of freedom (Mattox et al. 1996). As in our analysis the second model had two more degrees of freedom (i.e. normalization and index were left free), therefore TS = 9 (25) corresponds to significance of ~2.5 (4.6)σ, respectively. During the analyzed period, QSO B0218+357 was not always significantly detected. Flux upper limits at the 95% confidence level were calculated for each interval where the source TS was <9.

2.3. Swift

QSO B0218+357 was observed by the Swift satellite during 10 epochs, each with an exposure of about 4.5 ks. The observations did first follow the original alert of enhanced activity in GeV gamma rays, and then were resumed at the expected time of arrival of the delayed component. The data were reduced with the HEASoft package version 6.17. The Swift X-ray Telescope (XRT, Burrows et al. 2005) is a CCD imaging spectrometer, sensitive in the 0.2–10 keV band. We reduced the data using the calibration files available in the version 20140709 of the Swift-XRT CALDB. We run the task xrtpipeline with standard screening criteria on the observations performed in pointing mode. Observations were done in Photon Counting (PC) mode with count rates about 0.02 counts/s. The weak X-ray emission compelled us to merge different epochs to create a good quality spectrum (see Sect. 4.3). We combined different event files with the task xselect summing the corresponding exposure maps with the task image. The merged source and background counts were extracted with the task xrtproducts from a circular region of 35 for the source and 120 for the background. We grouped each spectrum with the corresponding background, redistribution matrix (rmf), and ancillary response files (arf) with the task grppha, setting a binning of at least 20 counts for each spectral channel in order to use the χ2 statistics. The spectra were analyzed with Xspec version 12.8.1. We adopted a Galactic absorption of NH = 5.6 × 1020 cm-2 from the LAB survey (Kalberla et al. 2005).

Simultaneous observations by the Ultraviolet Optical Telescope (UVOT, Roming et al. 2005), on board of Swift, did not result in a significant detection of the emission from the source in the UV range.

2.4. KVA

The optical R-band observations were done using the 35 cm Celestron telescope attached to the KVA 60 cm telescope (La Palma, Canary islands, Spain). The observations started on 2014 July 24 (MJD = 56 862.2) and continued on almost nightly basis until 2014 August 5 (MJD = 56 874.2). Further follow-up observations were performed in August and September. The data have been analyzed using the semi-automatic pipeline developed at the Tuorla Observatory (Nilsson et al., in prep.). The magnitudes are measured using differential photometry. We performed absolute calibration of the optical fluxes using stars with known magnitudes present in the field of view of the instrument during observations of all targets of a given night (see Table 3 of Nilsson et al. 2007, and references therein). QSO B0218+357 is rather faint in the optical range (about 19 mag) and the telescope is relatively small, therefore several images from the same night were combined for the measurement of the average flux. For the spectral analysis the optical flux was deabsorbed using a galactic extinction of AR = 0.15 (Schlafly & Finkbeiner 2011).

3. Influence of the lensing galaxy

The interpretation of the QSO B0218+357 observations is not trivial due to the influence of the lensing galaxy. The lensing introduces a (de-)magnification factor to the observed flux with respect to the intrinsic emission, due to both lensing by the galaxy itself, which is additionally causing the time delay between the images, and microlensing by individual stars in the lens galaxy (Vovk & Neronov 2016). Variability of the flux magnification caused by microlensing is larger (and thus can more significantly affect the measured light curve) for smaller emission regions in the source. Using a simple Singular Isothermal Sphere model (SIS, see e.g. Kochanek et al. 2004) we roughly estimate the absolute magnification of the leading and trailing images. A rigorous lens modeling performed by Barnacka et al. (2016) yielded a model consistent with SIS.

The ratio between the observed angular distances to the lens of the leading and trailing radio images of the source has been measured to be ~4 (York et al. 2005). In the framework of SIS model this results in the individual magnifications of the two images to be μleading ≈ 2.7, μtrailing ≈ 0.67.

Using the flux ratio between the two images measured in the radio frequency range, the same model allows us to also compute the absolute magnifications independently. With the value of μleading/μtrailing ≈ 3.6 (Biggs et al. 1999) we obtain very similar results μleading ≈ 2.8, μtrailing ≈ 0.77. Averaging both methods we assume μleading ≈ 2.7, μtrailing ≈ 0.7 in the further calculations. The radio emission in blazars is believed to originate from regions much larger than the ones involved in gamma-ray production. Therefore, the values given above for the individual magnifications of images are not affected by possible microlensing on individual stars of B0218+357 G.

On the other hand the microlensing can significantly modify the fluxes observed in the HE and VHE energy ranges (Neronov et al. 2015; Vovk & Neronov 2016). The flux magnification due to microlensing depends on the size of the emission region, which might vary with the energy e.g. due to cooling effects. Thus, it might modify the observed spectrum and this, in principle, can affect the EBL constraints and source modeling. However, during the 2014 flaring period the magnification ratio observed in Fermi-LAT was comparable to, or larger than, the radio one. This suggests that the microlensing, if present, might have a bigger effect on the leading rather than on the trailing image, which was observed by MAGIC. Namely, if a microlensing event amplified the observed emission during the delayed flare with a given magnification of μstar,trailing, the leading flare must have been also amplified with even larger magnification μstar,leadingμstar,trailing by an independent microlensing event to keep the observed ratio of fluxes. Assuming that a probability that the flux of the trailing image is magnified with a factor of μstar,trailing is ptrailing, the probability that both images are independently magnified resulting in the observed flux ratio is much smaller, roughly .

Absorption in the lensing galaxy can also affect the observed fluxes at different energies. Falco et al. (1999) interpreted the different reddening of the two images of QSO B0218+357 as an additional absorption of the leading image with the differential extinction ΔE(BV) = 0.90 ± 0.14. In fact the absorption is so strong that it inverts the brightness ratio of the two images in the optical range, making the trailing image brighter. Also, in the leading image, the H2 column density was estimated at the level of 0.5–5 × 1022 [cm-2] by an observation of a molecular absorption (Menten & Reid 1996). In addition the dependence of the radio flux ratio on the frequency could also stem from free-free absorption (Mittal et al. 2007). No absorption has been measured for the trailing image. Observations of 21 cm absorption feature in B0218+357 points to an HI column density of 1021(Ts/ 100 K)/(f/ 0.4) [cm-2], where TS is the spin temperature and f is the fraction of the flux density obscured by HI (Carilli et al. 1993). The absorption of sub-TeV emission is expected to be negligible in the lensing galaxies (Barnacka et al. 2014).

4. Results

In this section we discuss the spectral and temporal characteristics of the QSO B0218+357 emission obtained in different energy bands.

4.1. MAGIC

The VHE gamma-ray emission was detected on the nights of 2014 July 25 and 2014 July 26 (MJD = 56 863.2 and 56 864.2 respectively), during the expected arrival time of the delayed component of the flare registered by Fermi-LAT. The detection cuts were optimized to provide the best sensitivity in the 60–100 GeV estimated energy range (see Aleksić et al. 2016b). The total observation time during those 2 nights of 2.11 h yielded a statistical significance, computed according to Li & Ma (1983), Eq. (17), of 5.7σ (see Fig. 1).

thumbnail Fig. 1

Distribution of the squared angular distance, θ2, between the reconstructed source position and the nominal source position (points) or the background estimation position (shaded area). The vertical dashed line shows the value of θ2 up to which the number of excess events and significance are computed.

Open with DEXTER

The light curve above 100 GeV is shown in Fig. 2.

thumbnail Fig. 2

Light curve of QSO B0218+357 during the flaring state in July/August 2014. Top panel: MAGIC (points) above 100 GeV and a Gaussian fit to the peak position (thick solid line). Second panel from the top: Fermi-LAT above 0.3 GeV with the average flux from the 3rd Fermi Catalog (Acero et al. 2015) marked with a dashed line. Notice that, during the days where the trailing emission was expected Fermi-LAT was in pointing mode allowing the significant detection of lower flux levels. Third panel from the top: Swift-XRT count rate in the 0.3–10 keV range. Bottom panel: KVA in R band (not corrected for the contribution of host/lens galaxies and the Galactic extinction). The two shaded regions are separated by 11.46 days.

Open with DEXTER

A fit with a Gaussian function gives the peak position at MJD = 56 863.86 ± 0.30stat and a standard deviation of 0.75 ± 0.34stat days. The corresponding fit probability6 is 21%. The two flaring nights give a mean flux above 100 GeV of (5.8 ± 1.6stat ± 2.4syst) × 10-11 cm-2 s-1. The relatively large systematic error is mainly due to the 15% uncertainty in the energy scale.

The SED obtained from the two nights 2014 July 25 and 26 (MJD = 56 863.2 and 56 864.2) is presented in Fig. 3.

thumbnail Fig. 3

Gamma-ray SED of QSO B0218+357 as observed during the two flaring nights, 2014 July 25 and 26, by MAGIC (red filled circles) and after deabsorption in EBL according to the Domínguez et al. (2011) model (blue open squares). The shaded regions show the 1 standard deviation of the power-law fit to the MAGIC data. Black diamonds show the Fermi-LAT spectrum from the same time period. Black points show the average emission of QSO B0218+357 in the 3FGL catalog (Acero et al. 2015).

Open with DEXTER

The reconstructed spectrum spans the energy range 65–175 GeV and can be described as a power law (1)with the fit probability of 47%. The parameters obtained are: f0 = (2.0 ± 0.4stat ± 0.9syst) × 10-9 cm-2 s-1 TeV-1 and γ = 3.80 ± 0.61stat ± 0.20syst. The quoted systematic uncertainty on the spectral index takes into account also the small background estimation uncertainty for a weak low-energy source (see Eq. (3) of Aleksić et al. 2016b). As the redshift of the source is close to 1 the spectrum is severely affected by the absorption of VHE gamma-rays in the EBL. Correcting the observed spectrum for such absorption modeled according to Domínguez et al. (2011), we obtain an intrinsic spectral index of 2.35 ± 0.75stat ± 0.20syst. The corresponding normalization of the emission at 100 GeV is (4.6 ± 0.8stat ± 2.1stat) × 10-9 cm-2s-1TeV-1. The spectral points are obtained using the Bertero unfolding method, while the fit parameters are obtained using the so-called forward unfolding (Albert et al. 2007).

4.2. Fermi-LAT

The GeV light curve of QSO B0218+357 is shown in the second panel of Fig. 2. We used a minimum energy 0.3 GeV in the light curve (instead of 0.1 GeV) in order to increase the signal to noise ratio in the flux measurements: the spectrum of this source during this flaring episode is very hard (see below), while the diffuse backgrounds fall with energy with an index of >2.4, and the PSF of LAT at 0.1 GeV is about twice larger than that at 0.3 GeV. Significant GeV gamma-ray emission was detected by Fermi-LAT both during the leading flare and during the expected arrival time of the delayed emission (TS of 615 and 129 respectively). The spectrum contemporaneous to the MAGIC detection, derived between MJD = 56 863.07 and 56 864.85, can be described by a power-law function with slope γ = 1.6 ± 0.1. The corresponding flux above 0.1 GeV is F > 0.1 GeV = (1.7 ± 0.4) × 10-7 cm-2 s-1. For comparison, the leading flare was marginally harder, γ = 1.35 ± 0.09, with ~4 times higher flux F > 0.1 GeV = (6.7 ± 1.0) × 10-7 cm-2 s-1 The spectral index measured by Fermi-LAT during the 2014 outburst, is much harder than γ ~ 2.3 during both the flaring period of 2012 (Cheung et al. 2014) and the average state of this source reported in Third Fermi-LAT Catalog (Acero et al. 2015).

4.3. Swift-XRT

In the third panel of Fig. 2 we present the X-ray light curve of QSO B0218+357. The whole observed light curve shows only a small hint of variability. A constant fit gives χ2/Nd.o.f. = 21.3/9, corresponding to the probability of 1.1%. The source did not show an enhanced flux in the X-ray range during the trailing gamma-ray flare. The average count rate from the two observations during the enhanced gamma-ray flux results in (79.2 ± 7.7)% of the rate averaged from the remaining 8 pointings. The rate obtained in the 0.3–10 keV energy range is similar to the one obtained during the 2012 flaring period (0.027 ± 0.003, Donato et al. 2012).

As the source is a weak X-ray emitter and the observed variability is not very strong we have combined all the pointings for the spectral modeling of the source. Moreover, the lack of strong variability also implies that the observed emission is the sum of the two images of the source, with at least one of them affected by the hydrogen absorption. In order to provide higher accuracy per spectral point, we rebin the spectrum to 50 events per bin.

We model the X-ray spectrum as a sum of two power-law components, with the same intrinsic normalization and spectral slope, but magnifications fixed to 2.7 and 0.7 respectively (see Sect. 3). Following the detection by Menten & Reid (1996) of the molecular absorption line in the brighter image, we include hydrogen absorption at the redshift of the lens in the first (brighter) component. However, due to large uncertainty in the hydrogen column density, we leave it as a free parameter of the model. With such assumptions, X-ray intrinsic spectrum can be well described (χ2/Nd.o.f. = 42.2/34) by a simple power law (see Fig. 4).

thumbnail Fig. 4

Energy-binned counts observed by Swift-XRT from the direction of QSO B0218+357 (data points). The emission is modeled as a sum (solid red line) of two power-law components with the same spectral index. The first component (A image) is magnified by a factor 2.7 with an additional strong hydrogen absorption at the lens (dotted blue line). The second component (B image) is intrinsically weaker (magnification factor 0.7), but not absorbed at the lens (dashed green line).

Open with DEXTER

We obtained the following spectral fit parameters (2)where the reported uncertainties are statistical only. The corresponding column density (2.4 ± 0.5) × 1022 at. cm-2 is within the bounds given by the radio measurement of Menten & Reid (1996).

The X-ray spectrum can be alternatively described by a simpler model, considering only absorption of the total emission (i.e. same absorption is affecting both images). The resulting spectrum is then slightly harder, with an index of 1.59 ± 0.10. The corresponding effective hydrogen column density is smaller, (0.57 ± 0.17) × 1022 at. cm-2. The fit probability is however worse in this case, with χ2/Nd.o.f. = 54.7/34. Therefore, for the SED modeling (see Sect. 5) we use the spectrum obtained using the assumption that the absorption affects only the leading image.

4.4. KVA

The bottom panel of Fig. 2 shows the optical light curve of QSO B0218+357 in the R band. In all of our observations the source was fainter than 19 mag. The resulting error bars for the flux points were therefore relatively large and no significant variability was detected. We estimated the observed galaxy flux within our measurement aperture (5.0 arcsec radius) and taking into account that the calibration is made through an aperture of the same size. Using the data in Lehár et al. (2000) a lens galaxy flux of Fgalaxy = 13 μJy was derived. The resulting flux (corrected for both the Galactic absorption and the galaxy contribution) for the observation during the flare is then 70 ± 20 μJy.

5. Modeling of the broadband emission

In order to model the broadband emission spectrum of QSO B0218+357 we need to determine the magnification factors affecting different energy ranges and correct for them. As no strong variability was seen in either the optical or the X-ray range we can assume that the observed emission in those energy ranges is the sum of both lensed images. However, the optical leading image is strongly absorbed (Falco et al. 1999), thus the total magnification in the optical range is close to μtrailing. On the other hand in the X-ray range the emission 2 keV is not strongly absorbed in either of the two images. We correct the absorption of softer X-rays in the analysis (see Sect. 4.3). Therefore we assume that the magnification in the X-ray energy range is μleading + μtrailing ≈ 3.4. The strong variability in the GeV and sub-TeV gamma-ray range and the much harder GeV spectrum during the MAGIC observations point to the magnification in the GeV energy range at this time to be close to μtrailing. The broadband SED of QSO B0218+357 demagnified according to the numbers derived above and corrected for the X-ray and gamma-ray absorption is shown in Fig. 5. In green we report historical data, obtained from ASI Science Data Center (ASDC)7, tracking particularly well the low energy component. These historical data are the sum of the emission of the source passing through both of its images, however especially in optical and UV range are affected by strong absorption in the leading component. In order to derive the intrinsic flux of the source we apply the following correction factor to the flux 1/(μtrailing + μleading × TA(f)), where TA(f) is the frequency-dependent fraction of the leading image flux surviving the attenuation. In order to estimate TA(f) we use the differential extinction of the leading image, ΔE(BV) = 0.90 ± 0.14 (Falco et al. 1999), to scale the dust extinction curve of the Milky Way (Pei 1992; Xue et al. 2016), taking also into account the redshift of the lens. In the X-ray range the TA(f) shape was determined from the hydrogen column density obtained in Sect. 4.3.

thumbnail Fig. 5

Broadband SED of QSO B0218+357 modeled with a two-zone model. The reconstructed fluxes (red squares) are corrected for different magnifications in different energy ranges (see the text). Historical measurements (ASDC7) are shown with green circles and triangles (flux upper limit). Gray curves depict the emission from the region located within the BLR, while orange curves refer to the region located beyond the BLR. Long dashed curves show the synchrotron radiation, dotted the SSC emission and short dashed the external Compton emission. Dashed-dotted light blue line represents the accretion disk emission and its X-ray corona. The solid black line shows the sum of the non-thermal emission from both regions.

Open with DEXTER

The SED is dominated by the emission at GeV – sub-TeV energies, which is relatively common in flaring FSRQs (see e.g. Aleksić et al. 2014; Pacciani et al. 2014). Although the corrections for lensing are uncertain, the intrinsic GeV spectrum appears to be hard for this flaring state.

Table 1

Input parameters for the emission model shown in Fig. 5.

Interestingly, the gamma-ray flare seen by MAGIC was not accompanied by a similar increase in either optical or in X-ray flux. This is unusual for FSRQs, where a correlation is often seen. Comparison of the optical data to the archival measurements, shows that there was no large change in the position of the low energy peak during the high energy flare. The high energy peak position however moved from the sub-GeV range in low state to tens of GeV during the flare.

5.1. One-zone leptonic models

As the two peaks have a large separation in frequency the resulting spectrum cannot be easily explained in the framework of one-zone models. As detailed in Tavecchio & Ghisellini (2014) in such a case one-zone models inevitably require quite large Doppler factors and very low magnetic fields. For the specific case of QSO B0218+357  in order to explain the high energy peak photons with frequency νC = 1025νC,25Hz ≃ 1025 Hz, the Lorentz factors of the electrons emitting at the peak would need to reach (or exceed in case the scattering is in the Thomson regime) 8 × 104νC,25δ-1(1 + zs), where δ is the Doppler factor of the emitting region. Since the same electrons are responsible also for the synchrotron radiation of the low energy peak, a very low value of the magnetic field is required: , where νs,12 is the synchrotron peak frequency in units of 1012 Hz. If the high energy component is interpreted as synchrotron-self-Compton (SSC) emission, the ratio of the high energy peak luminosity to the synchrotron luminosity has to be equal to the ratio between the synchrotron photon energy density and the magnetic field energy density. This condition, coupled to the value of the magnetic field derived above, allows us to derive the required Doppler factor (see Tavecchio & Ghisellini 2014, for details): , where Ls,46 is the synchrotron luminosity (measured in the units of 1046 erg s-1) and Δt1d is the variability timescale in units of days. With such a large value of the Doppler factor it is rather unlikely that the radiation energy density in the jet frame is dominated by the synchrotron one. Instead, as usually assumed in the modeling of a FSRQ (e.g. Sikora 1994), it is most plausible that the high energy component is produced by the scattering of external photons (from the broad line region, BLR, the disk or the molecular torus).

In the case of such an external Compton (EC) scenario, some constraints can be derived considering that the SSC emission, expected now to peak in the X-ray band, cannot have a flux exceeding the value fixed by the XRT data. Similarly to the discussion above for the SSC case, we can thus derive a constraint on the Doppler factor: . We conclude that the extremely large values of the Doppler factors, plus the lack of simultaneous optical and X-ray variability to the GeV and sub-TeV flare, strongly disfavor one-zone models, pointing instead to a two-zone model, as discussed in the case of the flaring phase of the FSRQ PKS 1222+216 (Tavecchio et al. 2011).

Another important element to consider for the modeling is the huge opacity for gamma rays characterizing the innermost regions of a FSRQ. In particular, gamma rays with energies exceeding a few tens of GeV produced within the radius of the BLR would be strongly attenuated. Therefore, the highest energy part of the spectrum, observed by MAGIC, should have been emitted close to or beyond the BLR radius (see e.g. Pacciani et al. 2014, and references therein).

5.2. Two-zone external Compton model

Considering the conclusions above, we reproduced the broadband emission of the source with a two-zone model, inspired by the scenario c) of Tavecchio et al. (2011). The two emission regions are moving with the same Doppler factor along the jet. We take the simplest assumption that the first emission region is located, as in the case of other FSRQs, inside the BLR. The opacity condition forces however the second emission region, the production site of the VHE gamma rays, to be outside of the BLR. The gamma-ray emission is the sum of the SSC and EC components on the radiation field of BLR and dust torus. Both radiation fields are included in the calculations of both emission zones, however the BLR radiation field dominates the EC in the zone closer to the black hole, and the radiation field of the torus dominates the farther zone. The luminosity of the accretion disk is taken to be Ld = 6 × 1044 erg s-1 (Ghisellini et al. 2010). This value is quite low, if compared to a typical FSRQ. The radius of the BLR and that of the torus, calculated accordingly to the scaling laws of Ghisellini & Tavecchio (2009), are RBLR = 7.7 × 1016 cm and Rtorus = 2 × 1018 cm.

According to the scenario presented here, the GeV and sub-TeV emission is mostly produced in the EC and SSC process in the farther region (see orange curves in Fig. 5). This allows it to escape strong absorption of sub-TeV emission in the BLR radiation field. The size of the emission region is sufficiently small to account for the one-day variability timescales observed in this energy band (see Fig. 2). On the other hand the optical and X-ray emission comes mostly from the inner region. Lack of strong variability in those energy bands seen in Fig. 2 points to the stability of the emission from this region on the timescales of at least a fortnight. It is also self-consistent with the procedure of demagnification of the flux described in Sect. 4.3. We recall that blazar emission models, reproducing the innermost regions of the jet (distance from the black hole below 1 pc), cannot account for the radio emission (frequencies at which the region is optically thick) which, instead, is produced by farther, optically thin regions of the jet. The spatial separation of “Jet in” and “Jet out” might in principle introduce a delay between emission observed from them. If the same population of electrons, traversing along the jet, encounters first “Jet in” and afterwards “Jet out” we can expect to observe a delay of: ~(1 + zsRdist/ (cΓδ), where ΔRdist is the distance between the two regions. Using the modeling paramaters reported in Table 1 one would obtain a time delay of only ~6.9 h, which is significantly shorter than the duration of the flare, and very small on the temporal scale of Fig. 2. Moreover the delay would not be observable if, as assumed above, the emission from the “Jet in” region is quasi-stable.

The parameters for the region inside the BLR (Table 1) are in the range of those typically derived for a FSRQ with leptonic models (e.g. Ghisellini & Tavecchio 2015). For the outer regions there is a strong constraint on the luminosity of the synchrotron component, which – given the large Lorentz factors of the electrons required to produce the high energy component – peaks in the UV-soft X-ray band. To keep the synchrotron component below the limits and, at the same time, reproduce the powerful high energy IC component, the magnetic field must be kept to quite low values. This is similar to the case of PKS 1222+216 discussed in Tavecchio et al. (2011). As in that case, a possibility to explain such low values could be to assume that this emission region is the product of processes involving magnetic reconnection, in which magnetic energy is efficiently converted to electron energies (e.g. Sironi et al. 2015).

6. EBL constraints

The VHE gamma-ray observations of distant sources can be used to constrain the level of EBL. A wide range of methods have been applied in the past, starting from comparing spectral shape in the unabsorbed GeV range with the one in the TeV range (e.g. Dwek & Slavin 1994), and progressing to more elaborate methods, such as testing a grid of generic EBL spectral shapes and excluding the ones resulting in a pile up (i.e. convex spectrum) or a too hard intrinsic spectrum (Mazin & Raue 2007). More recently population studies have been performed. In particular Abramowski et al. (2013) used the specific shape of the EBL-induced feature in the spectrum. They performed a joint fit of multiple TeV spectra of sources with redshift z< 0.2 assuming smoothness of the intrinsic spectrum and putting constraints on the normalization factor of an EBL model with ~15% precision. A somewhat different approach was applied by Ackermann et al. (2012a) to Fermi-LAT data of about 150 BL Lacs. In this case low photon statistics in the sub-TeV regime does not provide a good handle on the particular spectral shape of the feature. However, as the observations below a few tens of GeV are not affected by the EBL absorption, the authors could obtain a direct measure of the absorption at the sub-TeV energies, assuming that the intrinsic spectral shape can be described by a log parabola. In the last class of methods (see e.g. Domínguez et al. 2013) instead of generic function shapes in TeV, the spectral form based on modeling of broadband emission with synchrotron-self-Compton scenario is used.

thumbnail Fig. 6

Probability of a SED fit as a function of the EBL scaling parameter. Different styles and colors of the lines represent different spectral shapes: power law (solid, black), power law with an exponential cut-off (dotted, green), log parabola (long dashed, blue), log parabola with exponential cut-off (dot-dashed, red). The vertical lines show the scaling for which the best probability is obtained (solid) and + 1 change in χ2 of the fit from this maximum (dotted). Nominal light scale of MAGIC is assumed in the middle panel. The light scale is decreased (increased) by 15% in left (right) panel.

Open with DEXTER

Even though the above methods applied to nearby (z ≲ 0.2) sources improved greatly our EBL knowledge in this redshift range, and the detection of 3C 279 (Albert et al. 2008) led to a major revision of the EBL models, the measurements at higher redshift are still sparse and are burdened with large uncertainties. The 1σ error band of the Ackermann et al. (2012a) measurement for the sources with redshift 0.5 <z< 1.6 allows for about a factor of two uncertainty in optical depth for EBL absorption. More recently, PKS1441+25 observations with MAGIC and VERITAS resulted in constraints on the scaling factor of optical depth predicted by the current EBL models to be 1.5–1.7 (Ahnen et al. 2015; Abeysekara et al. 2015).

Here we use the Fermi-LAT and MAGIC data collected from QSO B0218+357 during the flare to perform an independent measurement of EBL absorption at z = 0.944. Since the two instruments measured a similar timescale of the flare it is plausible to assume that the GeV and sub-TeV emission originates from the same emission region. This assumption is further supported by the SED modeling presented in Sect. 5. The dependence of the size of the emission region on the energy, combined with microlensing, might affect the observed GeV spectrum (see Sect. 3), introducing additional systematic uncertainties in the derived constraints on EBL. The spectrum observed by MAGIC from QSO B0218+357 gives us a chance to probe the EBL at wavelengths of ~0.3–1.1 μm. We use a method adapted from Abramowski et al. (2013). We perform a joint spectral fit combining Fermi-LAT and MAGIC points using a set of possible spectral shapes. To cover better the energy range of the EBL induced cut-off for this study we use a finer binning of the MAGIC data than presented in Fig. 3, resulting in 5 bins. The intrinsic spectral shapes are attenuated by EBL according to optical depths presented in Domínguez et al. (2011); however we allow an additional scaling parameter α of the optical depth. The following spectral models (power law, power law with a cut-off, log parabola, log parabola with a cut-off) are used where we apply additional source physics-driven conditions: Ecut> 0, b> 0. For each spectral shape we compute the χ2 value of the fit as a function of α. We determine the best fit and the best estimation of α from the minimum of such a curve. The 1σ statistical uncertainty bounds of the α parameter can be obtained as the range of α in which the χ2 increases by 1 from the minimum value.

The fit probability as a function of the EBL scaling parameter is shown in the middle panel of Fig. 6. Out of the phenomenological function shapes (Eqs. (3)–(6)) the highest fit probability is obtained with the simple power-law spectral model. Using this spectral model we obtain an estimation of the EBL scaling parameter of α = 1.19 ± 0.42stat at a redshift of 0.944. Such an assumption of a single power-law between 3 and 200 GeV, even though slightly preferred by the best fit probability, might be at odds with the FSRQ emission models. The spectral models allowing for an intrinsic curvature/cutoff exhibit a slower dependence of χ2 on the EBL scaling for the low values of α resulting in less constraining bounds. Notably, all the tested spectral shapes provide a 1σ upper bound below the value for a simple power-law spectral shape. Spectral shapes with an additional intrinsic cut-off result in only a small increase of χ2 from the scaling factor of 1 (nominal EBL) to 0 (no EBL). Therefore no strict lower bound can be derived on α.

Systematic uncertainties can affect the obtained results. In particular, since we use a combination of MAGIC and Fermi-LAT data, a shift in energy scale or flux normalization between the spectra obtained from the two experiments could affect the result. Due to a very steep spectrum in the sub-TeV range the dominant systematic effect is the 15% uncertainty of the energy scale of MAGIC (Aleksić et al. 2016b). It can shift the reconstructed spectrum in energy causing a large, up to 40%, shift in estimated flux. This effect is much larger than the pure flux normalization uncertainty reported in Aleksić et al. (2016b). The uncertainty of the spectral slope, due to the limited energy range of the spectrum, has a negligible effect. Finally the systematic uncertainty of Fermi-LAT (see e.g. Ackermann et al. 2012b, where 2% accuracy on the energy scale is reported) is also negligible compared to the ones of MAGIC in the case of this source. Therefore in order to investigate the systematic uncertainty on the EBL scaling parameter we performed full analysis using telescope response with a modified light scale by ± 15% following the approach in Aleksić et al. (2016b) and Ahnen et al. (2016).

In all three cases (see Fig. 6) the best probability of the fit out of the assumed phenomenological function shapes is obtained with a simple power-law fit. However, the corresponding EBL scaling parameter shifts by 0.25 for a power-law case. Also the statistical error for an increased light scale in this case is slightly larger. Therefore, allowing for an intrinsic curvature (log-parabola spectral shape, and/or an exponential cut-off) we obtain a 95% C.L. upper limit of α< 2.7. This limit is less constraining than the one obtained with PKS 1441+25 (Ahnen et al. 2015).

We repeated the analysis substituting the EBL model of Domínguez et al. (2011) by other currently considered models: Franceschini et al. (2008), Finke et al. (2010), Gilmore et al. (2012), Inoue et al. (2013). In the case of all the models the highest fit probability was obtained with a power law spectrum. The results of the best scaling parameter of the optical depth of these models are summarized in Table 2.

Table 2

Limits on the scaling parameter α of the optical depths in various EBL models.

As in the case of the model of Domínguez et al. (2011) we report the limits on the optical depth scaling factor for a power-law intrinsic spectral shape and a more conservative 95% C.L. upper limit for intrinsic spectral shapes allowing an arbitrary steepening or a cut-off. The combined Fermi-LAT and MAGIC spectrum is consistent with all five EBL models considered here.

7. Conclusions

MAGIC has detected VHE gamma-ray emission from QSO B0218+357 during the trailing component of a flare in July 2014. It is currently the most distant source detected with a ground-based gamma-ray telescope, and the only gravitationally lensed source detected in VHE gamma-rays. The VHE gamma-ray emission lasted for two nights achieving the observed flux of ~30% of the Crab Nebula at 100 GeV. Using the EBL model from Domínguez et al. (2011), the intrinsic spectral index in this energy range was found to be 2.35 ± 0.75stat ± 0.20syst. The VHE gamma-ray flare was not accompanied by a simultaneous flux increase in the optical or X-ray energy range. We have modeled the X-ray emission as a sum of two components with different magnifications, the weaker one absorbed with column density of (2.4 ± 0.5) × 1022 at. cm-2. The combined Fermi-LAT and MAGIC energy spectrum is consistent with the current EBL models; however, these constraints are not very strong. The EBL density scaling parameter is less than 2.1–2.8 of the one predicted by the tested models. The broadband emission of QSO B0218+357 is modeled in a framework of a two-zone external Compton model. According to this scenario, the quasi-stable optical and X-ray emission originates mostly in the inner zone. The enhanced gamma-ray emission during the flare is produced in the second zone, located outside of the BLR.


1

Unless specified otherwise, the energies are given in the Earth’s frame of reference.

2

Defined as the flux of the source with a Crab-like spectral shape that gives a gamma-ray excess with a significance of 5σ.

3

Note that it is also possible to reconstruct the flux slightly below such a defined threshold.

6

The fit probability throughout the paper is defined as the probability that the observed χ2 of points distributed along the used model shape exceeds by chance the value of χ2 obtained from fitting the data points with this model.

Acknowledgments

We would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF and MPG, the Italian INFN and INAF, the Swiss National Fund SNF, the ERDF under the Spanish MINECO (FPA2012-39502), and the Japanese JSPS and MEXT is gratefully acknowledged. This work was also supported by the Centro de Excelencia Severo Ochoa SEV-2012-0234, CPAN CSD2007-00042, and MultiDark CSD2009-00064 projects of the Spanish Consolider-Ingenio 2010 programme, by grant 268740 of the Academy of Finland, by the Croatian Science Foundation (HrZZ) Project 09/176 and the University of Rijeka Project 13.12.1.3.02, by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3, and by the Polish MNiSzW grant 745/N-HESS-MAGIC/2010/0. 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’Énergie 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.

References

All Tables

Table 1

Input parameters for the emission model shown in Fig. 5.

Table 2

Limits on the scaling parameter α of the optical depths in various EBL models.

All Figures

thumbnail Fig. 1

Distribution of the squared angular distance, θ2, between the reconstructed source position and the nominal source position (points) or the background estimation position (shaded area). The vertical dashed line shows the value of θ2 up to which the number of excess events and significance are computed.

Open with DEXTER
In the text
thumbnail Fig. 2

Light curve of QSO B0218+357 during the flaring state in July/August 2014. Top panel: MAGIC (points) above 100 GeV and a Gaussian fit to the peak position (thick solid line). Second panel from the top: Fermi-LAT above 0.3 GeV with the average flux from the 3rd Fermi Catalog (Acero et al. 2015) marked with a dashed line. Notice that, during the days where the trailing emission was expected Fermi-LAT was in pointing mode allowing the significant detection of lower flux levels. Third panel from the top: Swift-XRT count rate in the 0.3–10 keV range. Bottom panel: KVA in R band (not corrected for the contribution of host/lens galaxies and the Galactic extinction). The two shaded regions are separated by 11.46 days.

Open with DEXTER
In the text
thumbnail Fig. 3

Gamma-ray SED of QSO B0218+357 as observed during the two flaring nights, 2014 July 25 and 26, by MAGIC (red filled circles) and after deabsorption in EBL according to the Domínguez et al. (2011) model (blue open squares). The shaded regions show the 1 standard deviation of the power-law fit to the MAGIC data. Black diamonds show the Fermi-LAT spectrum from the same time period. Black points show the average emission of QSO B0218+357 in the 3FGL catalog (Acero et al. 2015).

Open with DEXTER
In the text
thumbnail Fig. 4

Energy-binned counts observed by Swift-XRT from the direction of QSO B0218+357 (data points). The emission is modeled as a sum (solid red line) of two power-law components with the same spectral index. The first component (A image) is magnified by a factor 2.7 with an additional strong hydrogen absorption at the lens (dotted blue line). The second component (B image) is intrinsically weaker (magnification factor 0.7), but not absorbed at the lens (dashed green line).

Open with DEXTER
In the text
thumbnail Fig. 5

Broadband SED of QSO B0218+357 modeled with a two-zone model. The reconstructed fluxes (red squares) are corrected for different magnifications in different energy ranges (see the text). Historical measurements (ASDC7) are shown with green circles and triangles (flux upper limit). Gray curves depict the emission from the region located within the BLR, while orange curves refer to the region located beyond the BLR. Long dashed curves show the synchrotron radiation, dotted the SSC emission and short dashed the external Compton emission. Dashed-dotted light blue line represents the accretion disk emission and its X-ray corona. The solid black line shows the sum of the non-thermal emission from both regions.

Open with DEXTER
In the text
thumbnail Fig. 6

Probability of a SED fit as a function of the EBL scaling parameter. Different styles and colors of the lines represent different spectral shapes: power law (solid, black), power law with an exponential cut-off (dotted, green), log parabola (long dashed, blue), log parabola with exponential cut-off (dot-dashed, red). The vertical lines show the scaling for which the best probability is obtained (solid) and + 1 change in χ2 of the fit from this maximum (dotted). Nominal light scale of MAGIC is assumed in the middle panel. The light scale is decreased (increased) by 15% in left (right) panel.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.