© H.E.S.S. and MAGIC Collaborations 2021

1. Introduction

Rapid flares in blazars – that is, in active galactic nuclei where the relativistic jet is aligned with the observer’s line of sight (Blandford & Rees 1974) – are among the most puzzling events in the Universe. These flares have been observed in all blazar types and they can exhibit flux variations on the order of a few minutes and (Gaidos et al. 1996; Aharonian et al. 2007; Tavecchio et al. 2011; Ackermann et al. 2015). The spectral energy distribution (SED) of blazars exhibits two broadband components, which are usually attributed to synchrotron and inverse-Compton emission by relativistic electrons (Dermer & Schlickeiser 1993; Sikora et al. 1994; Mastichiadis & Kirk 1997; Błazejowski et al. 2000; Böttcher et al. 2013), or hadronic interactions (Böttcher et al. 2013; Mannheim 1993; Mücke et al. 2003). Within the standard one-zone model, the bulk of the broadband radiation (from the far-infrared to VHE γ-ray energies) is produced in a single emission region, while the radio emission is produced in the extended jet. The single emission region supposedly fills the cross-section of the jet and, therefore, fast flux variations and their inferred small emission regions are often considered to have originated from a region close to the base of the jet. Observations of fast-variable VHE γ-ray emission have challenged this picture (see e.g., Aleksić et al. 2011).

PKS 1510–089 is a blazar at cosmological redshift of z = 0.361 (Tanner et al. 1996), which corresponds to a luminosity distance of 1.9 Gpc, assuming a cosmic expansion with Hubble parameter H₀ = 70 km s⁻¹ Mpc⁻¹, matter density Ω_m = 0.3, and dark energy density Ω_Λ = 0.7. PKS 1510–089 exhibits broad emission lines in its optical spectrum, leading to its classification as a flat-spectrum radio quasar (FSRQ). The broad lines are produced by gas in the so-called broad-line region (BLR), which is in close proximity to the black hole. As BLR photons absorb VHE γ-ray photons through pair production, any VHE γ-ray emission detected from FSRQs is most likely produced beyond the BLR (e.g., Costamante et al. 2018; Meyer et al. 2019).

Multiwavelength monitoring observations of PKS 1510–089 revealed its complex behaviour (Brown 2013; Foschini et al. 2013; Saito et al. 2013, 2015; Kushwaha et al. 2016), with changing correlation patterns between different energy bands. In turn, the SED modelling has resulted in different interpretations, such as two distinct photon fields for the inverse-Compton process (Nalewajko et al. 2012) or two spatially separated emission zones (Barnacka et al. 2014). While these models are not exclusive, they emphasise the complex nature of this blazar. In

addition, one of the highest recorded speeds of apparent superluminal motion – up to ∼46 c – has been detected in radio observations of the jet of PKS 1510–089 (Jorstad et al. 2005).

PKS 1510–089 was detected for the first time in VHE γ rays with H.E.S.S. during a bright high-energy (HE, E >  100 MeV) γ-ray flare in 2009 (H.E.S.S. Collaboration 2013). Subsequent observations were initially guided by multiwavelength flares leading to a detection with MAGIC during another high state of HE γ-ray emission in 2012 (Aleksić et al. 2014a). Systematic monitoring efforts of VHE γ rays began several years later (Zacharias et al. 2019a), leading to the detection of PKS 1510–089 in VHE γ-rays with MAGIC during low HE γ-ray states (MAGIC Collaboration 2018a).

In 2016, H.E.S.S. monitoring led to the detection of a strong flare at VHE γ rays that was not triggered by multiwavelength events, which was followed up with MAGIC observations and allowed for the first studies on sub-hour time scales in VHE γ rays for this source. The results of these observations are presented here, along with detailed observations in HE γ-rays and optical frequencies, as well as long-term radio very-long-baseline interferometry (VLBI) observations (Sect. 2). The flare behaviour is discussed in Sect. 3, including its variability, the absorption of VHE γ rays, evolution on large scales, and constraints on the emission region parameters. The summary and conclusions are given in Sect. 4.

2. Data analysis and results

The flare took place on JD 2457538 (May 29−30, 2016) and JD 2457539 (May 30−31, 2016; hereafter, ‘flare night’), and has been observed with H.E.S.S. and MAGIC in VHE γ rays, with Fermi-LAT in HE γ rays, and with ATOM in the optical R band. Radio VLBI monitoring data at 43 GHz and 86 GHz taken over several months complement the data set.

2.1. H.E.S.S.

H.E.S.S. is an array of five Imaging Atmospheric Cherenkov Telescopes (IACTs) located in the Khomas Highland in Namibia at an altitude of about 1800 m. Its first phase was initiated in 2004 with four telescopes (CT1-4) of a 107 m² mirror area each, laid out in a square of 120 m side length, giving an optimal energy threshold of ∼100 GeV. In 2012, a fifth telescope (CT5) with a 600 m² mirror area was added to the centre of the array, reducing the energy threshold to ∼50 GeV. Due to instrumental issues in the considered observation period (JD 2457536 – 2457546, May 27 – June 06 2016), the analysis only covers events recorded with CT2-4. The reduced number of telescopes, and observation conditions raise the energy threshold of this data set to ≳190 GeV.

A total number of 31 runs (1 run lasts for about 28 min) passed the standard quality selection (Aharonian et al. 2006), resulting in an acceptance corrected observation time of 13.6 h. The data set was analysed with the Model analysis chain using LOOSE CUTS (de Naurois & Rolland 2009). The results were cross-checked and verified using the independent reconstruction and analysis chain ImPACT (Parsons & Hinton 2014) with consistent results.

Over the entire observation period, the source is detected with a significance of 29.5σ above an energy threshold of E_thr = 196 GeV. A power-law fit,

$$\begin{array}{c}\hfill F(E)=N(E_0)\times {\left(\frac{E}{E_0}\right)}^{-\mathrm{\Gamma }},\end{array}$$(1)

to the measured spectrum results in a normalisation of $N(0.265\mathrm{TeV})=(1.12\pm 0.04{{}_{\text{stat}}}_{-0.4}^{+0.7}{}_{\text{sys}})\times {10}^{-10}$ ph cm⁻² s⁻¹ TeV⁻¹ with an index of Γ = 4.7 ± 0.1_stat ± 0.15_sys. The normalisation is derived at the decorrelation energy, E₀. The large systematic errors are a consequence of the energy uncertainty of 15% (Aharonian et al. 2006) folded with the very soft spectrum of the source. Using a log-parabola function,

$$\begin{array}{c}\hfill F(E)=N(E_0)\times {\left(\frac{E}{E_0}\right)}^{-\alpha -\beta log(E/E_0)},\end{array}$$(2)

with a normalisation, $N(0.232\mathrm{TeV})=(2.0\pm 0.1{{}_{\text{stat}}}_{-0.7}^{+1.2}{}_{\text{sys}})\times {10}^{-10}$ ph cm⁻² s⁻¹ TeV⁻¹, and parameters α = 3.8 ± 0.3_stat ± 0.15_sys and β = 1.9 ± 0.6_stat marginally improves the fit by 2.8σ for one additional free parameter. The measured spectrum of the entire observation period is shown in the left panel of Fig. 1. The average flux of the period above 200 GeV is $(2.1\pm 0.1{{}_{\text{stat}}}_{-0.7}^{+1.3}{}_{\text{sys}})\times {10}^{-11}$ ph cm⁻² s⁻¹.

Fig. 1. Measured VHE γ-ray spectra from H.E.S.S. with spectral type and time frames as indicated. Points have a significance of at least 3σ, while upper limits are at a 99% confidence level. The ‘full data set’ spectrum is the average spectrum of all observations from JD 2457536 – 2457546.

The light curve with nightly averages and energy integrated above 200 GeV, using the average spectrum of the full data set, is shown in Fig. 2a. A fit of a constant to the light curve is ruled out on a run-by-run binning with more than 10σ. However, it is obvious that the source is not significantly detected apart from the two nights on JD 2457538 and 2457539, respectively. Given the high statistics of these two nights, they have been analysed individually.

Fig. 2. Light curves of PKS 1510–089 during the observation period. a: nightly averaged VHE γ-ray light curve from H.E.S.S. (red points) and MAGIC (green open squares). The grey bars mark the systematic uncertainties. b: 24 h-average HE γ-ray light curve from Fermi-LAT. The dashed line marks the 11 years average with F̅ = 3.92 × 10−7 ph cm−2 s−1. Arrows mark 95% CL upper limits for bins with test statistics (TS) < 9.0. c: 24 h-binned HE γ-ray spectral index assuming a power-law spectrum. The dashed line marks the 11-year average with the index = 2.402. d: nightly averaged optical R-band light curve (de-reddened) from ATOM. The horizontal bars mark the observation time.

On JD 2457538 (May 29−30, 2016), H.E.S.S. observed PKS 1510–089 for three runs and an acceptance-corrected observation time of 1.3 h. The spectrum of this night is compatible with a power-law with a normalisation, $N(0.274\mathrm{TeV})=(1.9\pm 0.2{{}_{\text{stat}}}_{-0.7}^{+1.1}{}_{\text{sys}})\times {10}^{-10}$ ph cm⁻² s⁻¹ TeV⁻¹, and an index of Γ = 4.3 ± 0.3_stat ± 0.15_sys. A log-parabola, while not statistically preferred over the power-law, is described by a normalisation of $N(0.240\mathrm{TeV})=(2.5\pm 0.3{{}_{\text{stat}}}_{-0.9}^{+1.5}{}_{\text{sys}})\times {10}^{-10}$ ph cm⁻² s⁻¹ TeV⁻¹ and parameters, α = 3.4 ± 0.5_stat ± 0.15_sys and β = 3 ± 1_stat. The spectrum is shown in the middle panel of Fig. 1. Deriving a power-law spectrum that accounts for the absorption by the extragalactic background light (EBL) using the model of Franceschini et al. (2008), gives a normalisation, $N(0.274\mathrm{TeV})=(5.2\pm 0.4{{}_{\text{stat}}}_{-1.8}^{+3.1}{}_{\text{sys}})\times {10}^{-10}$ ph cm⁻² s⁻¹ TeV⁻¹, and an index of Γ = 2.3 ± 0.3_stat ± 0.15_sys. The average flux above 200 GeV is $(3.8\pm 0.3{{}_{\text{stat}}}_{-1.3}^{+2.3}{}_{\text{sys}})\times {10}^{-11}$ ph cm⁻² s⁻¹, with no significant deviation from this value in the individual runs.

Over the next night on JD 2457539 (May 30−31, 2016), the source was observed for four runs and an acceptance-corrected observation time of 1.8 h. A power-law fit to the spectrum results in a normalisation of $N(0.268\mathrm{TeV})=(7.0\pm 0.3{{}_{\text{stat}}}_{-2.5}^{+4.2}{}_{\text{sys}})\times {10}^{-10}$ ph cm⁻² s⁻¹ TeV⁻¹ and an index of Γ = 4.8 ± 0.1_stat ± 0.15_sys. A log-parabola function with a normalisation, $N(0.233\mathrm{TeV})=(14.3\pm 0.7{{}_{\text{stat}}}_{-5.0}^{+8.6}{}_{\text{sys}})\times {10}^{-10}$ ph cm⁻² s⁻¹ TeV⁻¹, and parameters α = 3.3 ± 0.4_stat ± 0.15_sys and β = 3.1 ± 0.7_stat improves the fit with respect to the power-law at a level of 4.4σ. This is the first detection of significant curvature in the VHE γ-ray spectrum of PKS 1510–089. The spectrum is shown in the right panel of Fig. 1. Apparently, the spectral parameters α and β do not change significantly between the two nights. They are also consistent within errors with the spectral parameters of the entire observation period. Analyzing the spectra for the individual runs does not give a significant deviation from the power-law parameters of the entire night. Applying the EBL model of Franceschini et al. (2008) to the spectrum of the whole night, the curvature is much reduced and the data are well-described by a power-law with a normalisation of $N(0.268\mathrm{TeV})=(19\pm 0.8{{}_{\text{stat}}}_{-7}^{+11}{}_{\text{sys}})\times {10}^{-10}$ ph cm⁻² s⁻¹ TeV⁻¹ and an index of Γ = 2.9 ± 0.2_stat ± 0.15_sys.

The light curve on JD 2457539 is shown in Fig. 3a. The average flux of the night above 200 GeV amounts to $(14.3\pm 0.6{{}_{\text{stat}}}_{-5.0}^{+8.6}{}_{\text{sys}})\times {10}^{-11}$ cm⁻² s⁻¹ or ∼56% of the Crab nebula flux (or Crab unit, CU) as measured with H.E.S.S. above the same energy threshold. However, the flux deviates with 4.5σ from this average value on run-wise time-scales, clearly indicating, for the first time, the VHE γ-ray intranight variability in PKS 1510–089. The peak flux of the night is equal to $(20\pm 1{{}_{\text{stat}}}_{-7.0}^{+12.0}{}_{\text{sys}})\times {10}^{-11}$ cm⁻² s⁻¹ or ∼80% CU above 200 GeV.

Fig. 3. Light curves of PKS 1510–089 of the flare night, JD 2457539. a: VHE γ-ray light curve from H.E.S.S. (red points) and MAGIC (green open squares). The binning is 28 min and 20 min for the H.E.S.S. and MAGIC light curves, respectively. The grey bars mark the systematic uncertainties. b: detected energies from Fermi-LAT above 1 GeV. The clustering of points in two bunches is a result of the orbital motion of the Fermi satellite. c: optical R-band light curve from ATOM showing individual exposures of about 8 min duration each.

2.2. MAGIC

MAGIC is a system made up of two IACTs with reflectors of 17 m diameter each. They are located on La Palma in the Canary Islands, at a height of 2200 m.a.s.l. (Aleksić et al. 2016a). As PKS 1510–089 is a southern source, only observable by MAGIC at zenith angle above 38°, the corresponding trigger threshold for a Crab-like spectrum is ≳90 GeV (Aleksić et al. 2016b), about 1.7 times larger than for the low zenith observations.

PKS 1510–089 has been regularly monitored by the MAGIC telescopes since 2013. The last observation before the flare night took place on JD 2457536.48, that is, three days prior. During the flare night, MAGIC observed PKS 1510–089 for 2.53 h. After the flare, nightly follow-up observations were performed until JD 2457543.50, resulting in a total observation time of 7.5 h. The data have been reduced and analyzed using MARS, the standard MAGIC analysis software (Zanin et al. 2013; Aleksić et al. 2016b).

During the period of investigation, the only night for which a significant signal was observed is the one on JD 2457539 (May 30−31, 2016). The average flux above 200 GeV in this night is $(7.36\pm 0.40{{}_{\text{stat}}}_{-2.6,\mathrm{sys}}^{+4.4})\times {10}^{-11}$ ph cm⁻² s⁻¹ corresponding to $(32.5\pm 1.8{{}_{\text{stat}}}_{-11,\mathrm{sys}}^{+19})\%$ CU as measured with MAGIC above the same energy threshold¹ (Aleksić et al. 2016b). The flux shows an increase by a factor of ≳5 with respect to the flux upper limits measured on neighboring nights, as can be seen in Fig. 2a. The large systematic uncertainty is dominated by a ∼15% uncertainty in the energy scale (Aleksić et al. 2016b) convolved with the very soft spectrum of PKS 1510–089. The effect of the systematics in the energy scale on the integral flux was estimated by performing dedicated Monte Carlo simulations in which the signals obtained in all the camera’s pixels were scaled by +15% or −15%.

The observed spectrum on the flare night was reconstructed between 60 and 700 GeV; shown in Fig. 4. It exhibits clear curvature and can be described by a log-parabola, Eq. (2), with normalisation, $N(0.175\mathrm{TeV})=(23.8\pm 1.4{{}_{\text{stat}}}_{-10}^{+11}{}_{\text{sys}})\times {10}^{-10}$ ph cm⁻² s⁻¹ TeV⁻¹, and parameters α = 4.31 ± 0.12_stat ± 0.15_sys and β = 1.60 ± 0.41_stat. Correcting for the EBL absorption according to Franceschini et al. (2008), the spectrum can be described by the simpler power-law form with normalisation, $N(0.175\mathrm{TeV})=(35\pm 2{{}_{\text{stat}}}_{-15}^{+17}{}_{\text{sys}})\times {10}^{-10}$ ph cm⁻² s⁻¹ TeV⁻¹, and an index of Γ = 3.37 ± 0.09_stat ± 0.15_sys. Similar to the case of the integral fluxes, the systematic uncertainties on the flux normalisation are dominated by the 15% uncertainty in the energy scale of MAGIC (Aleksić et al. 2016b).

Fig. 4. Measured VHE γ-ray spectrum from MAGIC for JD 2457539 in green open squares, and the VHE γ-ray low state spectrum from MAGIC Collaboration (2018a) in grey filled squares.

The intranight light curve binned in 20 min bins, as shown in Fig. 3a, shows a clear variability since a constant fit can be excluded with more than 10σ. The flux drops from ∼50% CU at the beginning of the observation down to ∼7.5% CU at the end of the night.

2.3. Fermi-LAT

Fermi-LAT monitors the HE γ-ray sky every three hours in the energy range from 20 MeV to beyond 300 GeV (Atwood et al. 2009). A binned analysis of the P8R3 SOURCE class events between energies of 100 MeV and 316 GeV was performed for a region of interest (ROI) of 15° ×15° centred at the position of PKS 1510–089². Sources within a region of 30° ×30° that is centred on the source listed in the 4FGL catalog (Ajello et al. 2020) have been accounted for in the likelihood analysis. In order to reduce contamination from the Earth Limb, a zenith angle cut of 90° was applied. The analysis was performed with FERMIPY³ v.0.17.4 and the FERMITOOLS⁴ v.1.0.0 software packages using the P8R3_SOURCE_V2⁵ instrument response function (Atwood et al. 2013) and the GLL_IEM_V07 and ISO_P8R3_SOURCE_V2_V01 models⁶ for the Galactic and isotropic diffuse emission (Acero et al. 2016), respectively.

In a first step, the average spectrum of PKS 1510–089 is derived within the time window JD 2457201 – 2457750, roughly 550 days centred around the VHE γ-ray flare. In the fit, all spectral parameters of sources within 3° from PKS 1510–089 are left free to vary along with the normalisation of sources up to 10° from the ROI centre and of the background templates. All other source parameters are fixed to their respective 4FGL values. The normalisations of the background templates are left as additional free parameters. Additionally, spectral indices are frozen if a hint of emission from the source is detected with a significance of 5 <  TS <  9⁷ and all parameters are fixed if TS <  5 or the predicted number of photons is less than 0.01. Neither the residual nor TS maps show any particular hot spots above (or below) a significance at the ∼2σ level. Therefore, the best-fit model describes the ROI well. The spectrum of PKS 1510–089 is described with a log-parabola function, Eq. (2), with normalisation, N(743.5 MeV) = (9.6 ± 0.2_stat) × 10⁻¹¹ ph cm⁻² s⁻¹ MeV⁻¹, and parameters, α = 2.33 ± 0.01_stat and β = 0.041 ± 0.008_stat.

The best-fit ROI model is then used to derive a daily light curve of PKS 1510–089. In each time bin, the spectrum of PKS 1510–089 is set to a simple power law and only the normalisation of the source 4FGL J1514.8−0949, associated with PMN J1514−0948, a blazar of an unknown type, located at a distance of only 0.9° of PKS 1510–089, is allowed to vary. This is the brightest source within 5° of PKS 1510–089, which corresponds approximately to the 68% confidence radius of the point spread function (PSF) of the LAT at 100 MeV. The light curve bins have been centred around the VHE γ-ray observation windows. The resulting light curve is shown in Fig. 2b. The average flux of PKS 1510–089 in the HE γ-ray energy band (in the ∼500 day window chosen above) is (7.04 ± 0.14) × 10⁻⁷ ph cm⁻² s⁻¹. The flux variation during the VHE γ-ray flare is not particularly strong and by no mean does it match the HE γ-ray outbursts that were already detected in this source (e.g., Marscher et al. 2010; Saito et al. 2013; Foschini et al. 2013; Orienti et al. 2013; Ahnen et al. 2017; Meyer et al. 2019). The spectral index of the power-law fit to each 24 h bin shows a significant hardening of the HE γ-ray spectrum during the VHE γ-ray flare (see Fig. 2c), reaching an index that is as hard as 1.6 ± 0.1. During the H.E.S.S. observation window, the spectral index further appears to harden to 1.4 ± 0.2, whereas during the MAGIC observations, the spectrum marginally softens to 1.7 ± 0.2 (compatible results are obtained with the unbinned analysis pipeline of the FERMITOOLS)⁸ This indicates that the flare primarily influenced the highest energies, while lower energy emission was not particularly enhanced. Hence, the peak frequency of the inverse Compton SED shifted significantly from ≲0.1 GeV on average to more than 10 GeV during the outburst. This is further supported by looking at the highest energy photons detected with Fermi-LAT during the H.E.S.S. and MAGIC observation windows (see Fig. 3b). During the peak VHE γ-ray flux, several photons above 10 GeV are detected that originate from PKS 1510–089 with a probability > 0.9⁹. The highest energy photon during the H.E.S.S. (MAGIC) observations has an energy of 20.6 GeV (6.0 GeV) and a probability of belonging to the source of 99.0% (99.5%).

2.4. ATOM

The Automatic Telescope for Optical Monitoring (ATOM) is an optical telescope with 75 cm aperture located on the H.E.S.S. site (Hauser et al. 2004). In operation since 2005, ATOM provides optical monitoring of γ-ray sources and observed PKS 1510–089 with a high cadence in the R band. The data were analysed using the fully automated ATOM Data Reduction and Analysis Software and their quality has been checked manually. The resulting flux was calculated via differential photometry using five custom-calibrated secondary standard stars in the same field of view.

After the detection of the beginning VHE γ-ray flare, the frequency of observations was increased to several exposures per night resulting in a detailed light curve. Figure 2c shows the nightly averaged optical light curve, while Fig. 3c displays individual exposures for the flare night, JD 2457539. It is apparent that the optical flux exhibited an outburst, even though the flux of the flare was several factors below the historical high state in July 2015 (Jankowsky et al. 2015; Zacharias et al. 2019a). Interestingly, the details between the VHE γ-ray and the optical light curve differ significantly during the flare night. While the VHE γ-ray flux shows a single peak, the optical flux exhibits a double peak structure. More detailed analyses regarding this fact are provided in Sect. 3.1.

2.5. VLBA and GMVA

PKS 1510–089 has been monitored monthly with the Very Long Baseline Array (VLBA) at 43 GHz (7 mm) under the VLBA-BU-Blazar program¹⁰ since June 2007. The resulting total and polarised intensity images of the FSRQ have been analyzed at 16 epochs from March 2016 to June 2017, the period relevant to the VHE γ-ray activity reported here. The data reduction for the VLBA analysis was performed using the Astronomical Image Process System (AIPS, version 31DEC17) and Difmap (version 2.4) software packages in the manner described by Jorstad et al. (2017). Table 1 gives information about the observations at each epoch: number of antennas, restoring beam dimensions and position angle, and flux density correction factor, f_amp. The flux density correction factor at a given epoch was derived by calculating the total flux densities in the images of sources that have weak emission outside the angular size range of the VLBA images at 43 GHz (PKS 0235+164, S5 0716+71, PKS 1055+01, B2 1156+29, and PKS 1749+096). These sources were observed alongside PKS 1510–089 in the program. The value for each source was compared with the total flux density measured at 37 GHz at the Metsähovi Radio Observatory (Aalto University, Finland) within one to two days of the VLBA observation. The final factor, f_amp, which is the average of correction factors over the sources, was applied to the PKS 1510–089 images to adjust the flux density scale.

The total intensity images were modelled using components with circular Gaussian brightness distributions, with each component (knot) characterised by the following parameters: flux density, S, distance from the core, r, position angle with respect to the core, Θ, and angular size of the component, a (FWHM of the Gaussian). The parameters used correspond to the best fit between the model and data with respect to the χ² value provided by Difmap. The uncertainties of the parameters were calculated using the formalism given in Jorstad et al. (2017), which relates the uncertainties to the brightness temperature of the respective knot. The analysis of polarised intensity images was carried out with the Interactive Data Language (IDL, version 8.6.1) software, using I, Q, and U Stokes parameter maps produced via Difmap to calculate the degree of polarisation and electric vector position angle (EVPA).

Figure 5 shows a sequence of the VLBA images representing the jet behaviour following the VHE γ-ray activity. There are five main components in the jet, designated in Fig. 5 as A0 (the core), A1, A2, K15, and K16. The core, A1, and A2 are presumably stationary features of the jet, while knots K15 and K16 move with respect to them. Table 2 gives the parameters of the knots. The appearance of knot K15 and its possible association with the previous VHE γ-ray outburst of PKS 1510–089 in May 2015, was discussed in Ahnen et al. (2017). The parameters of K15 listed in Table 2 are based on additional epochs to those in Ahnen et al. (2017); they confirm the properties of K15 that were reported previously.

Fig. 5. Sequence of total (contours) and polarised (colour-scale) intensity images of PKS 1510–089 at 43 GHz, convolved with a beam of full width at half maximum (FWHM) dimensions 0.38 × 0.15 mas2 along PA = −10°. The global total intensity peak is 2677 mJy beam−1 and the global polarised intensity peak is 147 mJy beam−1. Black line segments within each image show the direction of the polarisation electric vector; the length of each segment is proportional to the polarised intensity value; the black vertical line indicates the position of the core, A0; circles within the images mark positions of superluminal knots K16 (red) as well as K15 (blue), and quasi-stationary features A1 (yellow) and A2 (gray), according to the modelling.

PKS 1510–089 was also observed with the Global Millimeter-VLBI Array¹¹ (GMVA) at 3 mm (86 GHz) at two epochs close to the VHE γ-ray event (May 22 and October 2, 2016), with a resolution of ∼70 μas. The data reduction was performed in the same manner as described in Casadio et al. (2019), which resulted in total and polarised intensity images shown in Fig. 6. There are two contemporaneous images at different frequencies (3 mm on May 22, 2016, and 7 mm on June 11, 2016). They reveal the same main jet features: the core A0 and two knots positioned on the east and west sides of the core that are associated with K15 and K16, respectively. In addition, the 3 mm image on October 2, 2016, shows two knots besides A0 and K16, which can be associated with stationary knots A1 and A2 according to their positions.

Fig. 6. Sequence of total (contours) and polarised (colour-scale) intensity images of PKS 1510–089 at 86 GHz, with a beam of FWHM dimensions 0.10 × 0.07 mas2 along PA = −15°. The global total intensity peak is 1137 mJy beam−1 and the global polarised intensity peak is 74 mJy beam−1; the contours are 0.5%, 1%, 2%, ...64%, and 95% of the global intensity peak. Black line segments within each image show the direction of the polarisation electric vector; the length of the segment is proportional to the polarised intensity values; the black vertical line indicates the position of the core, A0; the blue, red, yellow, and grey circles mark positions of knots K15, K16, A2, and A1, respectively, according to the modelling.

Figure 7 left follows the separation of knot K16 from the core. The figure shows also quasi-stationary knots, A1 and A2, detected at some epochs, which are identified with quasi-stationary features A1 and A2, respectively, reported previously by Jorstad et al. (2001, 2005, 2017) based on similarity of their parameters. As can be seen in Fig. 7 left, K16 appears to accelerate beyond A2. In order to investigate whether the break in the motion of K16 is statistically significant, two mathematical representations have been fitted to the measured distance values of K16 from A0 as function of time: a straight line (with two free parameters) and a broken line (with four free parameters). While the χ²/N_d.o.f. = 17.1/12 obtained with the linear fit is statistically acceptable, a much lower value of χ²/N_d.o.f. = 1.4/10 is obtained for the broken linear fit. The χ² values are likely underestimated due to the involvement of a correlated systematic error connected with the amplitude calibration, which dominates the uncertainties in the brightness temperature calculations. Nevertheless, the large difference in the computed χ² points to a better description of the movement of knot K16 by the broken line. The break occurs on JD 2457716 ± 26 when K16 is at the same distance from the core as A2. The apparent speeds of K16 before and after the break are β_app = 10.8 ± 2.5 and β_app = 25.8 ± 2.9, respectively, in units of c. The ejection of K16 from the core took place on JD 2457361 ± 58 and the passage of K16 through A1 is estimated on JD 2457497 ± 46.

Fig. 7. Left: separation of knots K16 (red/brown circles at 7 mm/3 mm, respectively) from the core A0 (black dotted line) according to the 43 GHz and 86 GHz maps; the red lines approximate the motion of K16; dotted yellow and blue lines show the average positions of A1 (yellow/gray triangles at 7 mm/3 mm) and A2 (cyan/black triangles at 7 mm/3 mm); the red line segment at the position of A0 indicates the 1σ uncertainty of the ejection time of K16. Right: light curves of the core A0 (black), A1 (yellow), A2 (cyan), K15 (blue), K16 (red), the sum of all components (magenta) at 43 GHz, and from the entire source at 37 GHz (gray). The red vertical line indicates the ejection time (and its 1σ uncertainty) of K16 from A0. In both panels, the black vertical line marks the VHE γ-ray flare.

Figure 7 (left) also presents the positions of knot K16 with respect to the core in the 3 mm images, under the assumption that the locations of the 3 mm and 7 mm cores are the same. It is possible that the VLBI core at 3 mm is located closer to the central engine than the 7 mm core owing to frequency-dependent opacity (Königl 1982). However, positions of K16 at 3 mm agree very well with the initial (linear) motion derived from the 7 mm images, implying that any separation between the 3 mm and 7 mm cores is significantly less than the resolution of the 7 mm images.

Table 2 gives the physical parameters (Lorentz factor, Γ, Doppler factor, δ, and viewing angle, Θ_∘) for K15 and K16 using the approach suggested by Jorstad et al. (2005). The apparent speed β_app in the host galaxy frame depends on two parameters, Γ and Θ_∘, according to β_app = β sin (Θ_∘)/[1 − β cos (Θ_∘)], where $\beta =\sqrt{1-{\mathrm{\Gamma }}^{-2}}$ is the velocity in units of the speed of light. A change of the apparent speed can be caused by a variation in either of these parameters. The values of Γ ∼ 23 and Θ_∘ ∼ 0.6° of K16 listed in Table 2 are obtained near the core where radiative energy losses dominate over adiabatic losses (Jorstad et al. 2005). Unfortunately, the method cannot be used farther down the jet where adiabatic losses become important. However, the measurement of β_app ∼ 26 of K16 observed later implies that the Lorentz factor of K16 becomes Γ ≥ 26. If a minimum increase of Lorentz factor to Γ = 26 is assumed, then β_app ∼ 26 requires Θ_∘ ∼ 2.2°. This is a significant change of the viewing angle with respect to the value derived near the core. The change is comparable with the opening angle of the jet, ∼1° (Jorstad et al. 2017), and should be noticeable in projection; however, the projected position angle, Θ, of K16 does not change significantly. One possibility is that the change in jet direction is along the same position angle, as has been found to be the case between the parsec and kiloparsec scale jet of PKS 1510–089 (Homan et al. 2002). If, instead, the viewing angle of the feature were to remain stable during the epochs discussed here, K16 would need to have accelerated from Γ ∼ 23 to Γ ∼ 38. Although the MOJAVE survey (Homan et al. 2015) and VLBA-BU-BLAZAR program (Jorstad et al. 2017) have found that, statistically, relativistic jets of quasars accelerate (increase in Lorentz factor) gradually on scales from several to 100 pc, the change in the speed of K16 appears to occur specifically at the location of quasi-stationary feature A2. However, deceleration of a knot was found previously at the same position (knot B4 in Jorstad et al. 2017). In any case, a change of knot speeds after passage through the location of A2 supports the hypothesis that an interaction between the moving and quasi-stationary features had occurred.

Agreement between the timing of passage of K16 through A1 and the VHE γ-ray event makes the stationary feature especially interesting. The average distance of A1 is ∼0.18 mas downstream of A0, which is translated into a de-projected distance of ∼40 pc for an average viewing angle of the jet of 1.2° (Jorstad et al. 2017). The distance between the black hole and VLBI core at 43 GHz in several blazars is ∼10 pc, as reported by, for example, Fromm et al. (2015), Karamanavis et al. (2016). In the case of PKS 1510–089, this appears as a reasonable assumption, given that the 15 GHz VLBI core has been inferred to be located ∼18 pc from the black hole (Pushkarev et al. 2012). Therefore, the distance of A1 from the black hole is ∼50 pc. Uncertainties in the jet viewing angle, angular distance to A1 (cf. Table 2), as well as ∼50% relative uncertainty in the distance of the 43 GHz core, translate to a large uncertainty of ∼15−20 pc of this distance estimate, which is nonetheless incompatible with the typical BLR or dusty torus (DT) distances, which are ∼1−3 pc from the central engine (e.g., Kaspi et al. 2000).

The same VLBA data at 7 mm as discussed here have been analyzed by Park et al. (2019). The authors have also found two moving components, K15 and J15, which generally correspond to K15 and K16, respectively, of the present study. There are some differences in the identification of knots and fitting of polynomials yielding different ejection times of the knots. However, independent of these differences, according to Fig. 4 in Park et al. (2019) the time of passage of J15 through stationary feature A1 is JD ∼2457500 (April 21, 2016), which agrees within a 1σ uncertainty with that of K16 derived above. While Park et al. (2019) do not identify the quasi-stationary knots A1 and A2, it is a justifiable assumption to assume that they are present in their modelling as well: for example, their knot J15b can be identified with the stationary feature A2, and an unidentified knot near the core, seen in the last three epochs in their Fig. 2, has parameters similar to those of feature A1.

Figure 7 (right) presents the light curves of the moving knots, core, and stationary features of A1 and A2, along with the Metsähovi light curve of PKS 1510–089 at 37 GHz and the flux densities of the sum of all components detected at a given epoch. Both the core and 37 GHz light curves exhibited a significant increase in flux density at the time of ejection of K16. The light curves of K15 and K16 have been used to estimate the timescale of variability, τ, of the knots, which is ∼0.1 yr, similar for both knots. This value is applied to calculate the Doppler factor, δ, of K15 and K16 given in Table 2. Knot K16 possesses higher Doppler and Lorentz factors than K15. Although the knots have similar viewing angles, according to the VLBA and GMVA images, they should be located on opposite sides of the projected jet axis. The ratio of the maximum flux density of K16 to that of K15, ∼2, agrees with δ being the main factor for the modulation of the fluxes of knots, along with S_ν ∝ δ^3 + α, where −α is the slope of the optically thin flux density spectrum.

The magnetic field geometry of the jet can be inferred by analyzing the polarisation parameters of the core and K16. Figure 8 shows the degree, P, and electric-vector position angle, EVPA, of linear polarisation of the core. The core polarisation is moderate around the ejection time of K16, and the EVPA rotates as K16 separates from the core, with a large swing of ∼70° just after the VHE γ-ray event. The rotation of optical EVPAs during VHE γ-ray events have been reported for a number of blazars, for example, in S5 0716+71 (MAGIC Collaboration 2018b). In order to study the magnetic field of K16, profiles of the polarisation parameters have been constructed along a line drawn through the centre of the knot and perpendicular to the jet. This line moves away from the core with time, as K16 does. An example is shown in Fig. 9 (left). The profiles of P and the EVPA at some epochs are plotted in Fig. 9 (right). The profiles show a low degree of polarisation in the centre of K16, and a significant increase of P closer to the edges of the feature, where the EVPA is oblique or perpendicular to the jet direction. Such polarisation properties are typical for a spine-sheath structure of the jet, as discussed, for example, in MacDonald et al. (2015). A low polarisation degree in the centre of K16 indicates a very turbulent magnetic field in the spine.

Fig. 8. Top: degree of linear polarisation of the core, A0, computed from the VLBA data at 43 GHz. Bottom: position angle of polarisation of the core. The red vertical line marks the time of ejection of K16 with its 1σ uncertainty. The black vertical line marks the VHE γ-ray flare.

Fig. 9. Left: total (contours) and polarised (colour-scale) intensity images of PKS 1510–089 at 43 GHz. The total intensity peak is 2357 mJy beam−1 and the polarised intensity peak is 69 mJy beam−1; the contours are 0.5%, 1%, 2%...64%, and 95% of the total intensity peak. The red circle indicates the position of K16 according to modelling; the blue line across K16 and perpendicular to the jet direction shows the profile line for calculating polarisation parameters. Right: profiles of degree (top) and position angle (bottom) of polarisation at different epochs; the magenta horizontal line indicates the average jet direction.

3. Results and discussion

3.1. Variability study

From the VHE γ-ray light curves in Figs. 2a and 3a, the flare duration t_dur can be constrained to at most 48 h at these energies. The HE γ-ray spectral index (Fig. 2c) and the optical light curve (Fig. 2d) suggest that the flare may even have started a day earlier, implying a maximum duration of 72 h. The light curves also suggest a slow rise (up to 2 days) compared to a faster decay (potentially within a single night). The low cadence of observations in the VHE γ-ray and optical bands coupled with the low level of statistics in the HE γ-ray band do not allow us to draw firmer conclusions on the full extent of the flare.

It is not even clear whether the maximum in the VHE γ-ray light curve corresponds to the true peak of the flare or a secondary peak. While the HE γ-ray index indicates that the global peak probably took place within this one-day bin, which is centred on the H.E.S.S. observation window, any time between the H.E.S.S. observations on JD 2457538 and JD 2457539 is possible. The same is true for the optical light curve. The average fluxes of nights JD 2457538 and JD 2457539 are equal. However, the former measurement is only from a single exposure; whereas on the latter night, 15 exposures were conducted and display clear variability. Interestingly, the VHE γ-ray and optical light curves of JD 2457539 (Figs. 3a and c) show a different type of evolution. While the VHE γ-ray light curve is suggestive of a single peak and a subsequent decay, the optical light curve exhibits two peaks with the second one brighter than the first one. The decay after the second peak appears steeper than any other optical flux variation of the night.

The variability in PKS 1510–089 is further investigated using a phenomenological description assuming exponential variability (a test with a linear variability definition giving fully consistent results is provided in Appendix A) of the emission with the goal to determine the characteristic variability time scales and potential changes thereof. The exponential variability time scale between subsequent flux points F_i = F(t_i) is defined as

$$\begin{array}{c}\hfill t_{\exp }=\frac{t_{i+1}-t_i}{|ln{F}_{i+1}-ln{F}_i|}·\end{array}$$(3)

The inverse of this time scale is plotted in Fig. 10b for the VHE γ-ray light curve. The time scales are compatible with the harmonic mean of t̄_exp ∼ 1.8 h ∼ 108 min except for the next to last time step. Here, the exponential time scale t_exp = (16 ± 5_stat) min indicates a deviation from the harmonic mean of 2.4σ.

Fig. 10. Light curves and time scales of PKS 1510–089 of the flare night, JD 2457539. a: VHE γ-ray light curve from H.E.S.S. (red points) and MAGIC (green open squares). The binning is 28 min and 20 min for the H.E.S.S. and MAGIC light curves, respectively. The grey bars mark the systematic uncertainties. b: inverse exponential time scales between subsequent points of the VHE γ-ray light curve. Systematic uncertainties have been considered for the open symbol, as it marks the step from H.E.S.S. data points to MAGIC data points. The grey dash-dotted line marks the harmonic mean with ${\overline{t}}_{\exp }^{-1}\sim 0.56$ h−1. c: inverse exponential time scales between subsequent points of the R-band light curve. The grey dash-dotted line marks the harmonic mean with ${\overline{t}}_{\exp }^{-1}\sim 0.03$ h−1. d: optical R-band light curve from ATOM showing individual exposures of 8 min duration each.

3.1.1. Assessing the significance of the VHE γ-ray flux drop

In order to verify the steep decline of the flux in the penultimate time step, further tests have been conducted on the MAGIC light curve. The energy threshold of 200 GeV used for the light curve in Fig. 10a was selected to facilitate a comparison of the γ-ray emission above the same analysis threshold for both H.E.S.S. and MAGIC. However, due to the steepness of the source spectrum the choice of the common threshold is not optimal for the investigations of light curve features. By applying the new energy threshold of 100 GeV, the total excess of γ-ray events over the background obtained on the flare night becomes 5.5 times larger and 2.3 times more significant (4443 ± 86 compared to 805 ± 36 excess events). The observed drop in the flux is by more than a factor of two, therefore it cannot be explained by possible variations in the atmospheric transmission. Additionally, the coincident rate of γ-like events from background control regions does not show a similar drop.

Several fitting functions have been tested on the > 100 GeV MAGIC light curve (see Fig. 11). Due to the fast variability, for the fits the average value of the function inside the light curve bin is used rather than the value corresponding to the bin centre. The light curve fitted with an exponential decay in the time range covered by the MAGIC observations results in a flux halving time of (95 ± 7_stat) min. However the corresponding χ²/N_d.o.f. = 44.0/7 results in chance probability of only 2.2 × 10⁻⁷ falsifying the hypothesis of a simple exponential decay. Fitting the MAGIC light curve with a broken exponential, χ²/N_d.o.f. = 8.7/5 is obtained, indicating that the broken exponential is preferred over a simple exponential at the level of 5.5σ. The resulting fit shows the change of the exponential time scale from (180 ± 40_stat) min before the break to (21 ± 7_stat) min. It was also tested whether the light curve can be described with another function, in particular with a linear decline. However the resulting value of χ²/N_d.o.f. = 26.5/7 disfavours such a fit. More importantly, as the above fits are done from the light curves of a single instrument, and due to large amplitude of the observed variability, the uncertainties of the exponential time scales are dominated by the statistical rather then by systematic uncertainties.

Fig. 11. Various fits to the MAGIC light curve above 100 GeV on JD 2457539. Olive green line shows the exponential fit, magenta the broken exponential and cyan the linear fit. The exponential fit to the first seven points of the light curve is shown with red line, with the red shaded region showing the 68% C.L. uncertainty band and the green line shows its extrapolation to the end of the light curve.

An independent method was used to confirm the significance of the break in the light curve. The fit of the MAGIC light curve was conducted with a single exponential in a limited range (the first seven points). Such a fit gives a good description of the beginning of the light curve, with χ²/N_d.o.f. = 7.0/5. Then the fit is extrapolated to the last two points. The corresponding χ²/N_d.o.f. = 35.9/2 is calculated taking into account both the uncertainty of the fit extrapolation and of the last two flux measurements. Conservatively applying eight trials (motivated by the fact that the light curve could be broken between each two consecutive points), results in a chance probability of 1.3 × 10⁻⁷, corresponding to 5.1σ, compatible with the preference of the broken exponential over a single one.

As a final check of the significance of the light curve feature, a toy Monte Carlo (MC) study is performed based on the assumption that the source exhibits an exponential decay following the time scale and the flux normalisation of the single exponential fitted to the first seven points of the LC. For each bin, the integrated flux from this fit was calculated and using the collection area and effective time, a number of γ events has been computed. The latter is used as expected value of a Poisson distribution from which random numbers of γ events are drawn. The number of background events was also computed from Poissonian distribution following the measured background rates. The excess is then converted into the flux using the same collection area and effective time. The uncertainties of the flux take into account both the Poissonian fluctuations and the uncertainty of the collection area (the uncertainty of the effective time is negligible). Then, 2 × 10⁸ of such light curves were generated and the same analysis repeated as in the second method. A fraction of 4.7 × 10⁻⁷ light curves achieve the χ² value greater than the one obtained from the data. Correcting the chance probability conservatively for trials is 3.8 × 10⁻⁶, corresponding to 4.5σ.

In order to test if the variability in flux is accompanied by spectral variability the intrinsic (i.e. corrected for EBL absorption), the spectral index was computed for each 20 min long observation run of MAGIC. The results are shown in the left panel of Fig. 12. The first seven runs, even if a hint of spectral softening can be seen are still consistent with a constant spectral index (χ²/N_d.o.f. = 7.9/6). Including the points after the break into the fit, the fit chance probability gets much lower (χ²/N_d.o.f. = 16.9/8). Additionally, the average spectra from the seven runs of MAGIC before the break and from the two runs after the break have been constructed. Both spectra are shown in the right panel of Fig. 12. Consistently with the run-by-run analysis, the spectrum obtained from the end of observations with MAGIC appears softer than from the observations before the break in the light curve. The spectral index changes from 3.3 ± 0.1_stat to 4.4 ± 0.5_stat. A softening of the γ-ray spectrum points to a softening of the underlying electron distribution.

Fig. 12. Intrinsic spectra as measured by MAGIC on the night of the flare. Left: evolution of the spectral index with the line indicating the average value. Right: MAGIC spectra of the source obtained before (red points and magenta shaded region) and after (blue points and cyan shaded region) the break in the light curve. The dashed line shows the Crab nebula spectrum for a reference (Aleksić et al. 2015).

In conclusion, a simple exponential decay of the VHE γ-ray light curve is inconsistent with data observed by MAGIC above 100 GeV at a >4.5σ level. Instead, it shows evidence of a faster decline of the flux with a time scale of ∼20 min and a hint of accompanying spectral softening.

3.1.2. Variability in the optical light curve

For the optical light curve, the variability was assessed using the phenomenological function in Eq. (3), with the results presented in Fig. 10c. Additionally, the linear variability function in Appendix A provides consistent results, too. The harmonic mean of the optical exponential time scales is t̅_exp ∼ 35 h, about a factor 30 longer than the time scales at VHE γ rays during the flare. This slow variability time scale may be influenced by a significant background flux in the optical domain from parts of the jet not participating in the flare or from the accretion disk (D’Ammando et al. 2011), as the optical flux rises above the adjacent flux levels by only a factor of ≲2, see Fig. 2d.

There are two outliers in the optical timescales (see Fig. 10c), with the second one being more significant than the first one. The first outlier corresponds to a time scale of t_exp = (5 ± 1) h. It deviates by 2.8σ from the harmonic mean. This deviation is strongly influenced by the higher cadence of observations compared to the rest of the light curve. For the second outlier, t_exp = (5.8 ± 0.5) h is obtained. Here, the deviation is more than 9σ from the harmonic mean.

3.1.3. Comparison between the multiwavelength light curves

Both the VHE γ-ray and optical light curves exhibit intranight variability, with evidence for a significant steepening of the flux decay occurring in both bands almost simultaneously near the end of the flare. This implies a change in the physical conditions of the emission region, namely from the interplay of injection, acceleration, and cooling to cooling only (e.g., Sitarek & Bednarek 2010, their Fig. 10). The change in source’s behaviour is further supported by a hint of spectral softening in the VHE γ-ray band coincident with the drop in the observed light curve.

The total flux variation in the VHE γ-ray light curve is about an order of magnitude during the flare night, while the optical variation is only about 10%. This is also reflected in the variability time scales as mentioned above. However, both differences could be influenced by a significant steady background flux in the optical domain – such as from other regions of the jet, or the accretion disk – which would reduce the effect of the variable flux component on the total flux. Furthermore, the details in the light curves are different. The VHE γ-ray light curve exhibits a single peak followed by the decay extensively discussed above. On the contrary, the optical light curve displays two peaks with neither being coincident with the peak in the VHE γ-ray light curve. The simultaneous steep flux decay near the end of the observations may therefore point to a common termination of the flare processes.

On timescales of several days (see Fig. 2), the optical light curve matches the flux evolution of the HE γ-ray light curve within the errors. Both light curves exhibit similar flux variations of about a factor two. While this correlation can provide insights for modelling attempts, it is noteworthy that it is not a common feature in flares in PKS 1510–089 (e.g., Nalewajko et al. 2012; Ahnen et al. 2017; Prince et al. 2019).

While it is the first time that a fast flare is observed in VHE γ rays in PKS 1510–089, fast flares have been observed in other energy bands in this source before. Notably, 20−30 min flux doubling timescales – corresponding to a 30−40 min exponential variability timescale, as used in the present paper – were derived in HE γ rays during a very bright flare in 2011 (Foschini et al. 2013; Meyer et al. 2019) with peak fluxes exceeding 10⁻⁵ ph cm⁻² s⁻¹. As the HE γ-ray flux is about an order of magnitude lower during the event discussed here, the characteristics of the flares are different despite the similar variability timescale.

In previous cases of fast variability in the VHE γ-ray band in blazars, the observations revealed mainly either the smoothly raising or falling component of the flare (see Aleksić et al. 2011 for the FSRQ PKS 1222+216, and MAGIC Collaboration 2019 for the BL Lac object BL Lacertae), a combination of individual flares (see Aharonian et al. 2007 for the BL Lac object PKS 2155−304 and Aleksić et al. 2014b for the BL Lac object/radio galaxy IC 310), or nearly symmetric flares, without any evidence of a faster drop of the emission (see Albert et al. 2007 for the BL Lac object Mrk 501). Therefore, the sudden decrease of the flux in the light curve of PKS 1510–089 marks a rare occurrence of the cessation of a flare observed at VHE γ rays. A similar event to the 2016 flare of PKS 1510–089 was observed in the BL Lac object PKS 2155−304 on MJD 53945.97 (see Fig. 1 in Aharonian et al. 2009) and possibly also on MJD 53944 (see Fig. 8 in H.E.S.S. Collaboration 2010). The observation of such a feature in the FSRQ PKS 1510–089 is particularly interesting since the BLR absorption of VHE γ rays in this class of objects makes it impossible for the emission to originate in the inner, more compact parts of the jet.

3.2. The γ-ray spectrum

The EBL-deabsorbed HE and VHE γ-ray spectra of the peak night (JD 2457539) are plotted in Fig. 13. The apparent softening of the HE γ-ray spectrum from one observation window to the other, despite being not statistically significant, is probably influenced by the (non-)detection of photons with energies exceeding 10 GeV during the observation windows (see Fig. 3b). The peak of the γ-ray SED is located somewhere between 10 and 60 GeV. The grey data points display the low-state spectrum of PKS 1510–089 (MAGIC Collaboration 2018a) with a peak position at 100 MeV or less. During the flare, the peak position shifted by more than a factor 100 to higher energies.

Fig. 13. Observed HE and VHE γ-ray spectra of PKS 1510–089 during the flare night, JD 2457539. The Fermi-LAT butterflies are integrated over the precise H.E.S.S. (red) and MAGIC (green) observation windows and are plotted until the highest detected photon energy. The VHE γ-ray spectra are corrected for EBL absorption using the model of Franceschini et al. (2008), with the H.E.S.S. spectrum in red and the MAGIC spectrum in green. The grey data points are the low-state spectra from MAGIC Collaboration (2018a). For this spectrum, HE points above 1 GeV are omitted due to potential bias (for details, see MAGIC Collaboration 2018a).

The EBL-deabsorbed VHE γ-ray spectra are devoid of curvature. Interestingly, the spectral break between the HE and VHE γ-ray ranges does not change significantly between the observation windows. During the H.E.S.S. window, it is ΔΓ = 1.5 ± 0.3_stat, where as during the MAGIC window, it is ΔΓ = 1.7 ± 0.2_stat. Hence, the cause of the break appears steady during the evolution of the flare. Such a break is absent in the low-state data, with the connection between HE and VHE γ-ray data appearing smoother than during the flare.

There are a couple of possible explanations for the break. It could correspond to the maximum of the accelerated particle distribution, determined by the maximum efficiency of the acceleration process. A second possibility could be a break due to Klein-Nishina effects depending on the energies of the particles and the ambient photon field that is inverse Compton scattered to produce the γ rays. These possibilities are further discussed in Sect. 3.4. Here, the focus is on a third interpretation, namely absorption of the γ rays by a soft photon field located close to the central engine, like the broad-line region (BLR). By assessing the maximum absorption allowed by the data, the distance of the emission region from the black hole can be estimated. The following study is based on Meyer et al. (2019) and H.E.S.S. Collaboration (2019).

The detailed Fermi-LAT, H.E.S.S., and MAGIC spectra, covering more than three decades of energy, have enabled a simultaneous fit of a set of assumed intrinsic spectra folded with absorption patterns by the BLR. The latter is derived following the model of Finke (2016), which is motivated by reverberation mapping and assuming that the accretion disc radiation is absorbed by the BLR clouds and re-emitted as monochromatic lines at fixed distances from the black hole.

Two geometries of the BLR are implemented in the study. In the shell geometry, BLR photons are emitted in infinitesimal thin shells around the black hole; whereas in the ring geometry, the BLR photons originate from thin rings orthogonal to the jet axis. The model includes emission lines from Lyϵ to Hα but neglects any contribution from the thermal continuum. Motivated from reverberation mapping, each line has an associated luminosity and is emitted in a shell or a ring at a fixed distance (see Table 5 in Finke 2016). As input, the model requires the black hole mass, M_•, and the luminosity of the Hβ line, L(Hβ). For PKS 1510–089, the reverberation-mapping results obtained by Liu et al. (2006) are used, namely, log₁₀(M_•/M_⊙) = 8.2 with the solar mass, M_⊙, and L(Hβ) = 1.77 × 10⁴³ erg s⁻¹. Using the relations summarised in Finke (2016) between L(Hβ) and L(5100 Å), as well as between L(5100 Å) and the radius of the Hβ, the radius of the Lyα emitting shell or ring is determined, R_Lyα = 7.69 × 10¹⁶ cm. Its luminosity is the highest in the model (a factor of 12 higher than L(Hβ)) and, therefore, responsible for most of the absorption. The total BLR luminosity in this model is equal to 5.3 × 10⁴⁴ erg s⁻¹ and the Lyα emission accounts for 40 % of the total luminosity. Following Finke (2016), the optical depths τ_γγ(r, E′) are calculated as a function of the distance r of the emission region to the black hole and the γ-ray energy in the host galaxy frame E′. At E′ = 50 GeV (100 GeV), the ring geometry results in an optical depth of 1.1 (7.9) when the emission region is located at 0.5 R_Lyα and decreases to 0.4 (2.0) when the emission region is placed at R_Lyα.

The shell geometry results generally in higher values of the optical depth (compare also Fig. 14 in Finke 2016). In addition, a BLR model with a ring geometry is tested, which has a L(Hβ) luminosity 13 times larger than the value reported by Liu et al. (2006). This results in a BLR radius, R_BLR ∼ R_Lyα = 2.6 × 10¹⁷ cm, about three times as large as before, which corresponds to the value used in Aleksić et al. (2014a). This model is dubbed the ‘high luminosity BLR ring geometry’.

The distance, r, from the black hole is constrained by simultaneously fitting an intrinsic spectrum, F(E), modified by the BLR absorption exp(−τ_γγ) to the Fermi-LAT and IACT data, where the IACT data are corrected for EBL absorption following Franceschini et al. (2008) ¹². The Fermi-LAT data that were used are contemporaneous with the H.E.S.S. and MAGIC observation windows. Since these windows encompass less than three hours apiece and PKS 1510–089 was not in a flaring state in the LAT energy band, the likelihood is extracted as a function of the normalisation, N₀, and spectral index, Γ, and this likelihood surface, ℒ_Fermi(N₀, Γ, θ|D_Fermi), is used in the combined fit rather than flux points or likelihood curves for each energy bin as done by Meyer et al. (2019). In the likelihood function, θ denotes nuisance parameters (spectral parameters of the other sources in the ROI) and D_Fermi denotes the data. The likelihood surfaces for the H.E.S.S. and MAGIC observation windows are shown in Fig. 14.

Fig. 14. Two-dimensional likelihood surface of a power-law fit to the Fermi-LAT data contemporaneous to the H.E.S.S. observations (left) and MAGIC observations (right). The best-fit values are marked in red for H.E.S.S. and green for MAGIC. The likelihood surfaces are used in a simultaneous fit with IACT data to search for BLR absorption features. See main text for further details.

Different spectral shapes are tested since intrinsic curvature or a cut-off are degenerate with a cut-off induced by γ-ray absorption. Only spectral functions are used that can be described as power laws in the Fermi-LAT energy band in order to use the extracted likelihood surface, namely a simple power law (PL), a power law with sub-exponential cut-off (EPL), and a smoothly broken power law (BPL). It should be noted that for these functions, which are additionally multiplied by the BLR absorption, ℒ_Fermi is only an approximation of the true likelihood value, since ℒ_Fermi does not account for possible curvature at the highest Fermi-LAT energies. As BLR absorption only sets in above ∼20 GeV and the curvature of the chosen models also occurs at energies close to or beyond the highest energy photon detected with the LAT, the results are not expected to change if a complete likelihood formulation had been used instead.

From a physical point of view, a break in the spectrum is expected also for the intrinsic emission, as the SED γ-ray peak in blazars is usually located inside the observed energy range, except for only the most extreme high-synchrotron-peaked blazars (e.g., Singh et al. 2019; Biteau et al. 2020). Therefore, the intrinsic spectral shape should be modelled either with a BPL or EPL. However, the PL model is retained as a test case to evaluate whether the observed spectral break could in principle be explained with BLR attenuation only.

For each combination of intrinsic spectrum and assumed BLR geometry (ring, high luminosity ring, or shell), the parameters of the intrinsic spectrum and r are optimised. For the MAGIC data, the objective function χ² − 2lnℒ_Fermi is minimised, where χ² is the χ² value of the fit to the MAGIC data points which takes into account the full correlation matrix of the flux points. For the H.E.S.S. data, on the other hand, it is assumed that the likelihood for a flux value in the ith energy bin can be described with a Gaussian distribution, ℒ_{H.E.S.S., i}, as done in H.E.S.S. Collaboration (2019). Therefore, in this case, the objective function 2(lnℒ_Fermi + ∑_ilnℒ_{H.E.S.S., i}) is maximised. For flux upper limits, one sided Gaussian distributions are used. The sum over the likelihoods in each energy bin is then combined with ℒ_Fermi.

The best-fit spectra for all spectral functions and the BLR ring geometry are shown in Fig. 15. The best-fit index for the H.E.S.S. and Fermi-LAT spectrum for the PL spectrum deviates strongly from the best-fit index in the Fermi-LAT energy band, indicating the presence of an intrinsic spectral curvature. However, due to the broad likelihood profile (see Fig. 14), the change in ℒ_Fermi is only −3.8 for the PL spectrum.

Fig. 15. HE and VHE γ-ray spectra of PKS 1510–089 during the flare night, JD 2457539. The Fermi-LAT butterflies are integrated over the precise H.E.S.S. (red, left panel) and MAGIC (green, right panel) observation windows and are plotted up to the highest detected photon energy. The VHE γ-ray spectra are corrected for the EBL absorption with the model of Franceschini et al. (2008), with the H.E.S.S. spectrum in red and the MAGIC spectrum in green. The lines mark the best-fit spectra for combined Fermi-LAT and H.E.S.S. and MAGIC fits for the BLR ring geometry and different intrinsic spectra: a power law (PL), power law with sub-exponential cut-off (EPL), and a broken power law (BPL). In the PL case, curvature is only provided by the BLR absorption, whereas the EPL and BPL models include intrinsic spectral curvature. The best-fit values of the distance r from the black hole are provided in the legend.

In order to judge whether BLR absorption is significantly detected, the likelihood profiles of the distance r are derived for the different intrinsic spectra, BLR geometries, and also by taking into account the systematic uncertainties of the IACT measurements by shifting the IACT flux points to lower and higher flux values following the assessments in Sect. 2.

The profiles are shown in the panels of Fig. 16 for the H.E.S.S. (left panel) and MAGIC (right panel) combined fits; the spread of the profiles corresponds to the spread caused by the systematic uncertainties. For both MAGIC and H.E.S.S., when a simple PL is assumed as the intrinsic function, BLR absorption is required to provide the spectral break in order to describe the data. This results in a clear minimum in the likelihood curves for the PL model. The statistical preference of the fit with BLR absorption can be assessed by the TS and Δχ² values for the maximum values of r in Fig. 16, since for such large values of r the absorption is negligible in the considered energy bands. However, as described below, the PL fit is not significantly preferred over the intrinsic spectra with curvature. As a result, a detection of BLR absorption cannot be claimed. Furthermore, as discussed above, it is expected for an FSRQ like PKS 1510–089 that the SED peak of the high-energy emission falls within the γ-ray energy band, which is captured in intrinsic models with curvature. The preference of one intrinsic spectral model over another one is assessed with the Akaike information criterion (Akaike 1974):

$$\begin{array}{c}\hfill \mathrm{AIC}=2n_{\mathrm{par}}-2ln\mathcal{L},\end{array}$$(4)

where n_par is the number of fit parameters and ℒ is the total summed likelihood (2lnℒ = −χ² + 2lnℒ_Fermi for the combined Fermi-LAT and MAGIC fit). Models with an AIC larger than the minimum value, min(AIC), are less likely to minimise the loss function with a probability p(AIC) = exp[(min(AIC)−AIC)/2] (e.g., Burnham et al. 2011). The AIC and p(AIC) values are provided in Table 3. Even though for the MAGIC spectra the PL fits provide the minimum AIC value, the other spectral models are almost as likely to minimise the loss function, with p(AIC) between 0.2 and 0.5. For H.E.S.S. all spectral models seem to prefer a fit including BLR absorption (assuming the ring geometry). However, for the EPL and BPL models, the preference is not significant when the systematic uncertainties on the observed flux are taken into account. Judging from the AIC, models including curvature appear to be marginally more likely to minimise the cost function compared to a simple power law. However, this preference is not at a significant level.

Fig. 16. Likelihood and Δχ2 profiles for the distance r in the BLR ring geometry for the combined Fermi-LAT and H.E.S.S. (left) and Fermi-LAT and MAGIC fits (right). The spread of the profiles shows the effect of the IACT systematic uncertainties on the measured fluxes. For the H.E.S.S. and Fermi-LAT fits, the TS value is defined as twice the difference between the best-fit log-likelihood and the log-likelihood profiled over r. The dashed lines show where the TS or Δχ2 values are equal to 2.71, which corresponds to lower limits at 95% confidence.

The resulting best-fit spectra and observations are shown in Fig. 15. The best-fit values for the distance, r, are between 2 × 10¹⁷ cm and 8 × 10¹⁷ cm, which corresponds to 2.6 R_BLR ≲ r ≲ 10.4 R_BLR. The cut-off caused by the BLR is clearly seen for spectra without intrinsic curvature, which prefer a lower value of r, causing a steep and sharp drop in the spectra shown in Fig. 15.

Obviously, it is not possible to disentangle an intrinsic cut-off from a BLR-induced one, that is, it is not possible to put an upper limit on the distance of the emission region. Even though the broken power law seems to be the only function providing a satisfactory fit to the Fermi-LAT and H.E.S.S. data, the other spectral shapes cannot be ruled out. Therefore, only lower limits (at 95% confidence level) on the distance, r_lim, have been derived. Out of the tested BLR geometries and parameter space, the BLR ring geometry produces the least γ-ray attenuation and, hence, provides the most conservative (least constraining) limits. The limit values also demonstrate that the uncertainty introduced by the unknown exact BLR geometry dominates the systematic errors. For the combined H.E.S.S. and Fermi-LAT fit, the lowest lower limit is r >  2.0 × 10¹⁷ cm = 2.6 R_BLR. This limit is a factor of three stronger than previous limits (Meyer et al. 2019) derived from Fermi-LAT data alone. Due to the softer Fermi-LAT and MAGIC spectra, the limit is relaxed during the MAGIC observation window to r >  7.2 × 10¹⁶ cm = 0.9 R_BLR.

As the emission region is unlikely to move backwards (against the jet flow) between the H.E.S.S. and MAGIC observations, the lower limit during the H.E.S.S. observation window provides the minimum distance of the emission region from the black hole. Therefore, the emission region must be located outside of the BLR at a position where the jet has significantly widened.

It has also been investigated if absorption of the sub-TeV photons in DT radiation can be used to put further constraints on the location of the emission region. For specific parameters (DT temperature of 1000 K and emission region deep in the DT radiation field), the absorption can considerably affect the observed spectrum above ∼400 GeV. Nevertheless the large statistical and systematic uncertainties in the measured flux at the highest energies detected by H.E.S.S. and MAGIC, combined with poorly constrained DT parameters prevent deriving any robust limits on the location of the emission region with this method.

Finally, the fits with curved models can be used to estimate the possible position of the inverse Compton peak. Using the break energy in the BPL model and the BLR ring geometry, the peak is determined to be located at E_peak = (24.1 ± 23.2_stat) GeV for the H.E.S.S.observation window and (57.2 ± 47.5_stat) GeV for the MAGIC window. The large errors indicate that the peak can only be weakly constrained due to the large gap in energy between the Fermi-LAT and IACT observations. Within the uncertainties, the peak positions are consistent between the two measurements.

3.3. The large-scale jet structure

In the previous section, it is established that the emission region must have been located beyond the BLR. Interestingly, in coincidence with the flare a fast radio knot, K16, moved through the jet. Here, a possible connection between K16 and the flare is analysed.

The radio data suggest that K16 was ejected from the core between JD 2457303 (October 07, 2015) and JD 2457419 (January 31, 2016). A sketch showing the separation of K16 from A0, is presented in Fig. 17 (bottom panel), along with A1 and A2. The parameters of A1 and A2 (see Table 2) during the current observations are similar to those reported previously (Jorstad et al. 2017), and the projected position angle of A1 coincides with that of K16 at the time of the VHE γ-ray flare. Although the viewing angles of K15 and K16 are similar, according to the images the knots should be moving along trajectories on opposite sides of the jet cross-section resulting in different position angles. Using the proper motion of K16 determined from the VLBA images (Table 2), the knot should reach the distance of the A1 stationary feature in ∼4.5 month, which yields an epoch of possible interaction between the centroids of K16 and A1 of JD 2457497±46 (March 3–June 3, 2016), taking into account the uncertainties of the proper motion of K16 and the position of A1. Interestingly, the 1σ uncertainty of the interaction time includes the time of the VHE γ-ray flare, suggesting a possible connection between the events in question.

Fig. 17. Top: HE γ-ray light curve of PKS 1510–089 in 24 h bins. Grey arrows are upper limits (95% confidence level for TS < 9.0), and the grey dashed line indicates the 11-year average. Bottom: separation of knot K16 (red) from the core A0 (black line) according to the 43 GHz maps. The red line approximates the motion of K16 with the 1σ uncertainty given by the shaded region. The same is shown for A1 (blue) and A2 (magenta). In both panels the green dashed line indicates the time of the VHE γ-ray flare.

An inspection of Fig. 17, which also presents the HE γ-ray light curve over the entire time frame, suggests interesting associations between γ-ray events and activity in the parsec-scale jet: (1) the passage of K16 through the stationary feature A1 coincides with the VHE γ-ray flare and the beginning of some activity in the HE γ-ray band; and (2) prominent activity at HE γ-rays begins as K16 starts to interact with the upstream portion of A2. This is interesting as K16 may have changed direction or accelerated after the interaction with A2. Unfortunately, the large number of radio knots detected in PKS 1510–089, and the uncertainties in their crossing times make an unambiguous identification difficult. However, each of the four VHE γ-ray flares detected previously in PKS 1510–089 can be associated with a superluminal knot propagating in the jet (Marscher et al. 2010; H.E.S.S. Collaboration 2013; Aleksić et al. 2014a; Lindfors 2015; Ahnen et al. 2017).

The analysis of the polarisation parameters of the core and K16 (Figs. 8 and 9) reveals a highly variable, albeit low value of the polarisation P in K16, which varies from 3% to an undetectable level as K16 separates from the core. The core polarisation increases after the separation of K16 from the core, and its EVPA rotates by ∼70° just after the TeV detection (see Fig. 8). This implies that the knot moves down a turbulent spine that possesses a weakly ordered magnetic field component oriented perpendicular to the jet (Fig. 5). From this orientation, K16 may be interpreted as a transverse shock.

3.4. Modelling constraints

The most common interpretation of FSRQ γ-ray emission is inverse-Compton scattering of BLR or DT photons (see e.g., Böttcher et al. 2013; van den Berg et al. 2019). Here, the lack of absorption by the BLR photons places the emission region at least 2.0×10¹⁷ cm from the black hole outside of the BLR (for the assumed BLR model with R_BLR = 7.69 × 10¹⁶ cm). A common assumption in blazar modelling is that the emission region fills the entire width of the jet. In turn, the maximum distance r of the emission region from the black hole in a conical jet with opening angle (in units of radians) α = 0.26/Γ_b (cf. Pushkarev et al. 2009), can be calculated as

$$\begin{array}{c}\hfill r\approx \frac{t_{\mathrm{VHE}}^{\ast }c\delta }{(1+z)\alpha }=2.1\times {10}^{17}\mathrm{cm}\left(\frac{t_{\mathrm{VHE}}^{\ast }}{20\min }\right)\left(\frac{\delta }{45}\right)\left(\frac{{\mathrm{\Gamma }}_{\mathrm{b}}}{45}\right),\end{array}$$(5)

employing the observed minimum variability time scale $t_{\mathrm{VHE}}^{*}\sim 20$ min. Ignoring for the moment the possible association of the VHE γ-ray flare with the interaction of the knot K16 with the standing feature A1, the usual approximation δ = Γ_b was chosen. Even with the rather extreme assumptions on the opening angle and the Doppler factor, the distance is barely compatible with the minimum distance allowed by the absorption study. If the emission region is located beyond this distance, it cannot fill the width of the jet (Lister et al. 2013; Wehrle et al. 2016).

While an absorption limit cannot be obtained for the DT, the lower limit derived from the BLR absorption is compatible with the distance of ∼50 pc from the black hole, where the interaction between the knot K16 and the standing feature A1 took place. If the flare was indeed triggered by the interaction of K16 with A1, the emission region is not immersed in either the BLR or the DT photon fields, which are located at most a few parsecs from the black hole. This assumption is used in the following. In turn, going forward, the kinematic values derived in Sect. 2.5 are used, namely δ = 43 and Γ_b = 23.

During most observed states, the peak position of the inverse-Compton component in the SED of PKS 1510–089 is located at energies < 100 MeV (Saito et al. 2015; Brown 2013; Barnacka et al. 2014; H.E.S.S. Collaboration 2013; MAGIC Collaboration 2018a). During this flare, the peak position has shifted by more than a factor 100 to ∼50 GeV (Sect. 3.2), while the spectral fluxes at 100 MeV barely change (cf. Figs. 2 and 13). This implies that the distribution of the flare electrons is either hard or narrow with a high minimum electron Lorentz factor γ_min.

In case of a hard electron distribution, the break in the γ-ray spectrum could represent the Klein-Nishina break. The scattered photon energy in the Thomson regime is E_γ = 4γ²E_ph, where γ is the electron Lorentz factor and E_ph the soft photon energy. In the Klein-Nishina regime, nearly all of the electron energy is transferred to the γ ray, that is, E_γ ≈ γm_ec². At the Klein-Nishina break, the quantity 4γE_ph = m_ec², resulting in E_ph = m_ec²/(4γ) = (m_ec²)²/(4E_γ). The energy E_ph of the soft photon in the comoving frame then is

$$\begin{array}{c}\hfill E_{\mathrm{ph}}=\frac{\delta {(m_{\mathrm{e}}{c}^2)}^2}{4(1+z)E_{\gamma ,\mathrm{br}}^{\ast }}=41\mathrm{eV}\left(\frac{\delta }{43}\right){\left(\frac{E_{\gamma ,\mathrm{br}}^{\ast }}{50\mathrm{GeV}}\right)}^{-1},\end{array}$$(6)

where $E_{\gamma ,\mathrm{br}}^{\ast }$ is the γ-ray energy at the Klein-Nishina break in the observer’s frame. The electron Lorentz factor at the break is ${\gamma }_{\mathrm{br}}=m_{\mathrm{e}}{c}^2/(4E_{\mathrm{ph}})=3200{\left(\frac{\delta }{43}\right)}^{-1}\left(\frac{E_{\gamma ,\mathrm{br}}^{*}}{50\mathrm{GeV}}\right)$. If the photon is produced inside the jet, it should be a synchrotron photon. In this case, the synchrotron component would have extended well into the soft X-ray domain, which would be unusual for FSRQs, where the X-ray domain typically belongs to the inverse-Compton component. Unfortunately, no data are available to probe this possibility.

If the soft photon is provided by external, isotropic photon sources, another transformation to the frame of the photon production region needs to be applied, giving:

$$\begin{array}{c}\hfill E^{\prime }\sim E_{\mathrm{ph}}/{\mathrm{\Gamma }}_{\mathrm{b}}=1.8\mathrm{eV}\left(\frac{\delta }{43}\right){\left(\frac{E_{\gamma ,\mathrm{br}}^{\ast }}{50\mathrm{GeV}}\right)}^{-1}{\left(\frac{{\mathrm{\Gamma }}_{\mathrm{b}}}{23}\right)}^{-1}\end{array}$$(7)

in the case when this frame is that of the host galaxy. This is in the (red) optical domain. Such photons could be synchrotron emission from the radio feature A1 within a ring-of-fire scenario (MacDonald et al. 2015).

As the emission region is smaller than the full diameter of the jet, the surrounding jet material could also provide synchrotron seed photons within a spine-sheath scenario (Tavecchio & Ghisellini 2008). In this case, the relative Lorentz factor between the emission region and the surrounding jet material is less than between the emission region and a stationary (in the host-galaxy frame) source such as A1. Assuming a relative Lorentz factor of Γ_rel ∼ 5, the required synchrotron photon energy in the sheath is ∼8 eV.

The observed variability is a combination of particle acceleration and cooling effects, as well as the light-crossing time through the source. Which one is dominating at any given time is difficult to say without proper modelling. Nonetheless, some important clues can be derived if cooling dominates the cessation of the VHE γ-ray and optical flare.

Assuming again that, first, only internal processes dominate the cooling (resulting in an SSC situation), the Compton dominance – the ratio of the γ-ray peak flux to the synchrotron peak flux – is probably larger than unity¹³ and thereby dictates that the energy density in the synchrotron photons dominates over the magnetic field energy density. In such a situation, the cooling becomes nonlinear (Zacharias & Schlickeiser 2012), which could help explaining the fast cooling during the flare cessation. While simple estimates indicate that SSC is a viable option, proper constraints cannot be derived as the synchrotron spectrum during the flare is unknown.

If the cooling is dominated by inverse-Compton scattering of external photon sources, constraints on that photon source can be derived. The inverse-Compton cooling time scale is given by:

$$\begin{array}{c}\hfill t_{\mathrm{cool}}=\frac{3m_{\mathrm{e}}{c}^2}{4c{\sigma }_{\mathrm{T}}}{\left(u_{\mathrm{tot}}{F}_{\mathrm{KN}}\gamma \right)}^{-1},\end{array}$$(8)

with σ_T the Thomson cross section, u_tot the total energy density, and

$$\begin{array}{c}\hfill F_{\mathrm{KN}}\approx \{\begin{array}{cc}1,\hfill & 4E_{\mathrm{ph}}\gamma \ll m_{\mathrm{e}}{c}^2\hfill \\ {(1+\frac{4E_{\mathrm{ph}}\gamma }{m_{\mathrm{e}}{c}^2})}^{-1},\hfill & 4E_{\mathrm{ph}}\gamma \gg m_{\mathrm{e}}{c}^2\hfill \end{array}\end{array}$$(9)

as a Klein-Nishina correction factor (Moderski et al. 2005). Assuming that the VHE γ-ray photons with energies > 200 GeV are emitted in the (deep) Klein-Nishina domain, the total energy density can be approximated as

$$\begin{array}{c}\hfill u_{\mathrm{tot}}\approx \frac{3E_{\mathrm{ph}}}{c{\sigma }_{\mathrm{T}}{t}_{\mathrm{cool}}}=0.26\frac{\mathrm{erg}}{{\mathrm{cm}}^3}{\left(\frac{t_{\mathrm{VHE}}^{\ast }}{20\min }\right)}^{-1}{\left(\frac{E_{\gamma ,\mathrm{br}}^{\ast }}{50\mathrm{GeV}}\right)}^{-1},\end{array}$$(10)

where Eq. (6) is used, and the cooling time is estimated from the observed time scale of the flux drop $t_{\mathrm{VHE}}^{\ast }=t_{\mathrm{cool}}(1+z)/\delta$. Equation (10) assumes a narrow soft photon distribution peaking at E_ph, as no other handle is available for the nature of these photons. While it may be a bad representation of the true soft photon distribution and may influence the results, it should not invalidate the main conclusion.

In case of the ring-of-fire scenario, the total energy density that must be provided by A1 is

$$\begin{array}{c}\hfill u_{\mathrm{rf}}^{\prime }\sim u_{\mathrm{tot}}/{\mathrm{\Gamma }}_{\mathrm{b}}^2=5\times {10}^{-4}\frac{\mathrm{erg}}{{\mathrm{cm}}^3}{\left(\frac{t_{\mathrm{VHE}}^{\ast }}{20\min }\right)}^{-1}{\left(\frac{E_{\gamma ,\mathrm{br}}^{\ast }}{50\mathrm{GeV}}\right)}^{-1}{\left(\frac{{\mathrm{\Gamma }}_{\mathrm{b}}}{23}\right)}^{-2},\end{array}$$(11)

while in the spine-sheath scenario the sheath must provide:

$$\begin{array}{c}\hfill u_{\mathrm{ss}}^{\prime }\sim u_{\mathrm{tot}}/{\mathrm{\Gamma }}_{\mathrm{rel}}^2=9\times {10}^{-3}\frac{\mathrm{erg}}{{\mathrm{cm}}^3}{\left(\frac{t_{\mathrm{VHE}}^{\ast }}{20\min }\right)}^{-1}{\left(\frac{E_{\gamma ,\mathrm{br}}^{\ast }}{50\mathrm{GeV}}\right)}^{-1}{\left(\frac{{\mathrm{\Gamma }}_{\mathrm{rel}}}{5}\right)}^{-2}·\end{array}$$(12)

These are high energy densities for a jet region 50 pc from the black hole (cf., the modelling of a pc-scale jet in Zacharias & Wagner 2016, which gives u ≲ 10⁻⁶ erg cm⁻³). Such energy densities may leave observable signatures of the jet in the optical or UV domain.

The energy densities can be transformed into fluxes at Earth and compared to the total observed optical fluxes from the light curve in Fig. 2. The maximum optical flux during the flare is ∼1.3×10⁻¹¹ erg cm⁻² s⁻¹. The synchrotron energy densities are transformed to fluxes according to

$$\begin{array}{c}\hfill F_i^{\ast }=u_i^{\prime }{\delta }_i^{4}c{\left(\frac{a_{\mathrm{pc}}}{2D_{\mathrm{L}}}\right)}^2,\end{array}$$(13)

where $u_i^{\mathrm{\prime }}$ is the result from either Eq. (11) or Eq. (12), δ_i the respective Doppler factor, a_pc is the diameter of A1 in pc, which can be calculated from the angular size of A1 in Table 2, and D_L the luminosity distance of PKS 1510–089. Inserting $t_{\mathrm{VHE}}^{*}\approx 20$ min and $E_{\gamma ,\mathrm{br}}^{\ast }\approx 50$ GeV, the total flux in the ring-of-fire scenario is $F_{\mathrm{rf}}^{\ast }=8\times {10}^{-14}$ erg cm⁻² s⁻¹ with δ_rf = 1 for a stationary source in the host galactic frame. In the spine-sheath scenario, the result is $F_{\mathrm{ss}}{}^{*}=6\times {10}^{-10}{\left(\frac{{\delta }_{\mathrm{ss}}}{4.5}\right)}^4$ erg cm⁻² s⁻¹, where δ_ss is the Doppler factor of the sheath which has been estimated using the observation angle of the jet, 1.2° (Jorstad et al. 2017), and Γ_rel = Γ_bΓ_ss(1 − β_bβ_ss) (Ghisellini et al. 2005), where β_b and β_ss are the normalised speeds of K16 and the sheath, respectively, and Γ_b and Γ_ss are the corresponding Lorentz factors.

The expected flux in the ring-of-fire scenario is more than two orders of magnitude below the optical flux limit. This scenario is also in line with the fact that the variations of the γ-ray flux on the flare night by over an order of magnitude are accompanied with only much smaller variations in the optical flux. On the other hand, the flux in the spine-sheath scenario exceeds the optical flux limit by a factor 60. This disfavours the spine-sheath scenario as a viable option under the assumption that A1 (and thus the sheath) exhibits a homogeneous radiation density. The flux could be substantially reduced if the sheath Doppler factor, δ_ss, is lower or the relative Lorentz factor, Γ_rel, is higher. Both imply a lower sheath Lorentz factor Γ_ss. However, in order to drop the flux by the factor 60 requires Γ_ss ∼ 1.5. This almost resembles the ring-of-fire scenario.

The estimates in the ring-of-fire scenario can also be used to constrain a model in which the radiation field comes from luminous stars. Such models could naturally explain fast variations of the emission (see e.g., Banasiński et al. 2016). The estimated value of $F_{\mathrm{rf}}^{\ast }$ however corresponds to the isotropic luminosity of 4 × 10⁴³ erg s⁻¹, which is a few orders of magnitude larger than the luminosity of bright stars. Therefore, even close to the surface of the star, the radiation density would not be enough to cause the cooling break. An alternative argument against such a scenario, operating in the case of the flare discussed here, is the connection with the radio component at the distance of 50 pc from the base of the jet. The size of the jet at such a distance is several orders of magnitude larger than the radius of stars, therefore, only a small fraction of the blob filled with relativistic particles would be immersed in the strong radiation field from the star.

Combining the ring-of-fire scenario with models for fast variability – such as the Turbulent Extreme Multi-Zone (TEMZ) model of Marscher (2014) or the magnetic reconnection model of Giannios et al. (2009) – could explain the observations. Both models describe the radiation of small emission regions within a larger turbulent zone. In the TEMZ model, efficient acceleration of very high energy electrons that cause rapid VHE γ-ray flares results from a temporary alignment of the turbulent magnetic field with a direction relative to the shock normal that is conducive to such acceleration (Baring et al. 2017). In the magnetic reconnection model of Giannios et al. (2009), regions of plasma containing oppositely directed magnetic fields come into contact with each other, perhaps as a result of turbulent motions. The magnetic fields reconnect, which creates ‘mini-jets’, containing extremely energetic electrons that stream with high bulk Lorentz factors relative to the already highly relativistic ambient jet flow. Both models can produce sporadic VHE γ-ray flares with very short timescales of variability.

In summary, if the association between the VHE γ-ray flare and the interaction of K16 with A1 is true, the most reasonable options seem to be turbulently-variable SSC or the ring-of-fire scenario. In the former case, photons are provided by the flaring component itself, while in the latter case, the standing feature, A1, provides the target photons for the inverse-Compton scattering process.

4. Summary and conclusions

The dense monitoring of the FSRQ PKS 1510–089 in the VHE γ-ray band with H.E.S.S. (Zacharias et al. 2019a) and MAGIC (MAGIC Collaboration 2018a) led to the detection of a bright and short flare at this energy range, which lasted for two nights. The observed flux surpassed the previous flares by a factor of ten.

On the first night (JD 2457538), H.E.S.S. observed a significant flux increase compared to previous nights, but no variability within the 1.5 h of observation. On the second night (JD 2457539), an even brighter flux was detected with H.E.S.S. including significant variability. This is the first time that VHE γ-ray intranight variability was detected in PKS 1510–089. From a peak flux of about 80% of the Crab Nebula flux above an energy of 200 GeV, the flux decayed throughout the rest of the night, as shown through observations with MAGIC. Variability analyses have revealed that the VHE γ-ray light curve exhibited a common trend with a variability timescale of 1.5 h. Close to the end of the observations, a deviation from this common trend occurred revealing a variability time scale of about 20 min. A detailed analysis indicates that this sudden cessation is significant with more than 4.5σ.

Optical R-band observations with ATOM show a complex light curve on JD 2457539, which deviates from the behaviour at VHE γ-rays. Unlike the VHE γ-ray light curve, the optical one exhibits a double-peaked structure. The first peak occurs after the peak in the VHE γ-ray light curve, while the second optical peak has no correspondence in the VHE γ rays. The common variability time scale in the optical light curve is ∼15 h. Interestingly, at the same time as the cessation in the VHE γ-ray light curve, a similar steepening of the decay also takes place in the optical domain. This may imply that a common process is responsible for the faster drops of the respective light curves.

Despite the strong flare at VHE γ-ray energies, the HE γ-ray light curve (E >  100 MeV) was less affected by the outburst exhibiting only a mild flux variation. On the other hand, the HE γ-ray spectrum significantly hardened, reaching a value of the photon index of 1.6, compared to the usual value of ∼2.4. This implies a shift of the Compton peak for this flare towards the VHE range. Hence, the flare mainly influenced the higher energies of the Fermi-LAT energy range.

The γ-ray spectrum also shows remarkable features. For the first time in PKS 1510–089, the observed VHE γ-ray spectrum is significantly curved. The curvature is fully explained through absorption by the EBL. By combining the VHE with the HE γ-ray spectrum, additional absorption features by the BLR were searched for. The spectra are compatible with negligible absorption by the BLR, implying that the emission region is located at a distance from the black hole of at least r >  2.6 R_BLR. This is in line with results from other FSRQs (Aleksić et al. 2011, 2014a; Costamante et al. 2018; Zacharias et al. 2019b; Meyer et al. 2019; H.E.S.S. Collaboration 2019), indicating that the γ-ray emission region is located beyond the BLR.

The VLBI observations at 43 GHz over the course of a few months around the flare revealed a fast moving knot, K16 (δ ∼ 43), which crossed the standing jet feature, A1, around the time of the flare. A1 is located at a de-projected distance of about 50 pc from the black hole. Hence, if the association of the VHE γ-ray flare with the crossing of a stationary emission feature 50 pc from the black hole by a fast, bright knot is true, the production of VHE γ rays can occur far down the jet where turbulent plasma crosses a standing shock (Böttcher & Dermer 1998; MacDonald et al. 2015).

The change in the variability time scale during the cessation of the VHE γ-ray and optical flare indicates that the injection and acceleration of particles had stopped. In turn, the particle spectrum was only influenced by cooling processes, for instance, possibly inverse-Compton cooling on the optical photons of the feature A1, producing the steep decline in the light curves (Sitarek & Bednarek 2010). Hence, the emission region must have been compact (Giannios et al. 2009; Marscher 2014) and probably did not fill the width of the jet. All of this is contrary to the standard expectation that emission regions producing such energetic flares should be located within ∼1 pc of the black hole and fill the width of the jet.

Acknowledgments

The authors wish to thank the anonymous referee for a constructive report, which helped to improve the manuscript. The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of H.E.S.S. is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG), the Helmholtz Association, the Alexander von Humboldt Foundation, the French Ministry of Higher Education, Research and Innovation, the Centre National de la Recherche Scientifique (CNRS/IN2P3 and CNRS/INSU), the Commissariat à l’énergie atomique et aux énergies alternatives (CEA), the U.K. Science and Technology Facilities Council (STFC), the Knut and Alice Wallenberg Foundation, the National Science Centre, Poland grant no. 2016/22/M/ST9/00382, the South African Department of Science and Technology and National Research Foundation, the University of Namibia, the National Commission on Research, Science & Technology of Namibia (NCRST), the Austrian Federal Ministry of Education, Science and Research and the Austrian Science Fund (FWF), the Australian Research Council (ARC), the Japan Society for the Promotion of Science and by the University of Amsterdam. We appreciate the excellent work of the technical support staff in Berlin, Zeuthen, Heidelberg, Palaiseau, Paris, Saclay, Tübingen and in Namibia in the construction and operation of the equipment. This work benefited from services provided by the H.E.S.S. Virtual Organisation, supported by the national resource providers of the EGI Federation. 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 (FPA2017-87859-P, FPA2017-85668-P, FPA2017-82729-C6-2-R, FPA2017-82729-C6-6-R, FPA2017-82729-C6-5-R, AYA2015-71042-P, AYA2016-76012-C3-1-P, ESP2017-87055-C2-2-P, FPA2017-90566-REDC); 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-153/28.08.2018 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 SEV-2015-0548, the Unidad de Excelencia “María de Maeztu” MDM-2014-0369 and the “la Caixa” Foundation (fellowship LCF/BQ/PI18/11630012), by the Croatian Science Foundation (HrZZ) Project IP-2016-06-9782 and the University of Rijeka Project 13.12.1.3.02, 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. The VLBA is an instrument of the National Radio Astronomy Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated by Associated Universities, Inc. The research at Boston University was supported by NASA Fermi Guest Investigator grant 80NSSC17K0649. This research has made use of data obtained with the Global Millimeter VLBI Array (GMVA), which consists of telescopes operated by the MPIfR, IRAM, Onsala, Metsahovi, Yebes, the Korean VLBI Network, the Greenland Telescope, the Green Bank Observatory and the VLBA. The GMVA VLBI data were correlated at the correlator of the MPIfR in Bonn, Germany.

References

Appendix A: Linear variability time scale

As an additional test to the exponential time scale of variability, a linear variability time scale between subsequent flux points (following Zhang et al. 1999) is used:

$$\begin{array}{c}\hfill t_{\mathrm{lin}}=\frac{F_i+F_{i+1}}{2}\frac{t_{i+1}-t_i}{|F_{i+1}-F_i|}·\end{array}$$(A.1)

The time scales defined by Eq. (A.1) corroborate the results obtained with the exponential time scale, namely Eq. (3). The step with the fastest variability according to the exponential time scale (t_exp = (16 ± 5_stat) min deviating 2.4σ from the harmonic mean) exhibits a linear time scale of t_lin = (17 ± 5_stat) min deviating 3σ from the corresponding harmonic mean.

In the optical domain, Eq. (A.1) results in the same values as the exponential time scale. In turn, the same conclusion can be drawn about the deviating time steps from the harmonic mean.

Noteworthy, the cadence of observations with H.E.S.S., MAGIC, and ATOM is comparable within a factor of 2. Therefore, the presented light curves in Fig. 10 probe the same frequencies of variability at VHE γ-ray and optical ranges.

All Tables

All Figures

	Fig. 1. Measured VHE γ-ray spectra from H.E.S.S. with spectral type and time frames as indicated. Points have a significance of at least 3σ, while upper limits are at a 99% confidence level. The ‘full data set’ spectrum is the average spectrum of all observations from JD 2457536 – 2457546.
In the text

Fig. 2. Light curves of PKS 1510–089 during the observation period. a: nightly averaged VHE γ-ray light curve from H.E.S.S. (red points) and MAGIC (green open squares). The grey bars mark the systematic uncertainties. b: 24 h-average HE γ-ray light curve from Fermi-LAT. The dashed line marks the 11 years average with F̅ = 3.92 × 10−7 ph cm−2 s−1. Arrows mark 95% CL upper limits for bins with test statistics (TS) < 9.0. c: 24 h-binned HE γ-ray spectral index assuming a power-law spectrum. The dashed line marks the 11-year average with the index = 2.402. d: nightly averaged optical R-band light curve (de-reddened) from ATOM. The horizontal bars mark the observation time.

In the text

Fig. 3. Light curves of PKS 1510–089 of the flare night, JD 2457539. a: VHE γ-ray light curve from H.E.S.S. (red points) and MAGIC (green open squares). The binning is 28 min and 20 min for the H.E.S.S. and MAGIC light curves, respectively. The grey bars mark the systematic uncertainties. b: detected energies from Fermi-LAT above 1 GeV. The clustering of points in two bunches is a result of the orbital motion of the Fermi satellite. c: optical R-band light curve from ATOM showing individual exposures of about 8 min duration each.

In the text

	Fig. 4. Measured VHE γ-ray spectrum from MAGIC for JD 2457539 in green open squares, and the VHE γ-ray low state spectrum from MAGIC Collaboration (2018a) in grey filled squares.
In the text

Fig. 5. Sequence of total (contours) and polarised (colour-scale) intensity images of PKS 1510–089 at 43 GHz, convolved with a beam of full width at half maximum (FWHM) dimensions 0.38 × 0.15 mas2 along PA = −10°. The global total intensity peak is 2677 mJy beam−1 and the global polarised intensity peak is 147 mJy beam−1. Black line segments within each image show the direction of the polarisation electric vector; the length of each segment is proportional to the polarised intensity value; the black vertical line indicates the position of the core, A0; circles within the images mark positions of superluminal knots K16 (red) as well as K15 (blue), and quasi-stationary features A1 (yellow) and A2 (gray), according to the modelling.

In the text

Fig. 6. Sequence of total (contours) and polarised (colour-scale) intensity images of PKS 1510–089 at 86 GHz, with a beam of FWHM dimensions 0.10 × 0.07 mas2 along PA = −15°. The global total intensity peak is 1137 mJy beam−1 and the global polarised intensity peak is 74 mJy beam−1; the contours are 0.5%, 1%, 2%, ...64%, and 95% of the global intensity peak. Black line segments within each image show the direction of the polarisation electric vector; the length of the segment is proportional to the polarised intensity values; the black vertical line indicates the position of the core, A0; the blue, red, yellow, and grey circles mark positions of knots K15, K16, A2, and A1, respectively, according to the modelling.

In the text

Fig. 7. Left: separation of knots K16 (red/brown circles at 7 mm/3 mm, respectively) from the core A0 (black dotted line) according to the 43 GHz and 86 GHz maps; the red lines approximate the motion of K16; dotted yellow and blue lines show the average positions of A1 (yellow/gray triangles at 7 mm/3 mm) and A2 (cyan/black triangles at 7 mm/3 mm); the red line segment at the position of A0 indicates the 1σ uncertainty of the ejection time of K16. Right: light curves of the core A0 (black), A1 (yellow), A2 (cyan), K15 (blue), K16 (red), the sum of all components (magenta) at 43 GHz, and from the entire source at 37 GHz (gray). The red vertical line indicates the ejection time (and its 1σ uncertainty) of K16 from A0. In both panels, the black vertical line marks the VHE γ-ray flare.

In the text

	Fig. 8. Top: degree of linear polarisation of the core, A0, computed from the VLBA data at 43 GHz. Bottom: position angle of polarisation of the core. The red vertical line marks the time of ejection of K16 with its 1σ uncertainty. The black vertical line marks the VHE γ-ray flare.
In the text

Fig. 9. Left: total (contours) and polarised (colour-scale) intensity images of PKS 1510–089 at 43 GHz. The total intensity peak is 2357 mJy beam−1 and the polarised intensity peak is 69 mJy beam−1; the contours are 0.5%, 1%, 2%...64%, and 95% of the total intensity peak. The red circle indicates the position of K16 according to modelling; the blue line across K16 and perpendicular to the jet direction shows the profile line for calculating polarisation parameters. Right: profiles of degree (top) and position angle (bottom) of polarisation at different epochs; the magenta horizontal line indicates the average jet direction.

In the text

Fig. 10. Light curves and time scales of PKS 1510–089 of the flare night, JD 2457539. a: VHE γ-ray light curve from H.E.S.S. (red points) and MAGIC (green open squares). The binning is 28 min and 20 min for the H.E.S.S. and MAGIC light curves, respectively. The grey bars mark the systematic uncertainties. b: inverse exponential time scales between subsequent points of the VHE γ-ray light curve. Systematic uncertainties have been considered for the open symbol, as it marks the step from H.E.S.S. data points to MAGIC data points. The grey dash-dotted line marks the harmonic mean with ${\overline{t}}_{\exp }^{-1}\sim 0.56$ h−1. c: inverse exponential time scales between subsequent points of the R-band light curve. The grey dash-dotted line marks the harmonic mean with ${\overline{t}}_{\exp }^{-1}\sim 0.03$ h−1. d: optical R-band light curve from ATOM showing individual exposures of 8 min duration each.

In the text

Fig. 11. Various fits to the MAGIC light curve above 100 GeV on JD 2457539. Olive green line shows the exponential fit, magenta the broken exponential and cyan the linear fit. The exponential fit to the first seven points of the light curve is shown with red line, with the red shaded region showing the 68% C.L. uncertainty band and the green line shows its extrapolation to the end of the light curve.

In the text

Fig. 12. Intrinsic spectra as measured by MAGIC on the night of the flare. Left: evolution of the spectral index with the line indicating the average value. Right: MAGIC spectra of the source obtained before (red points and magenta shaded region) and after (blue points and cyan shaded region) the break in the light curve. The dashed line shows the Crab nebula spectrum for a reference (Aleksić et al. 2015).

In the text

Fig. 13. Observed HE and VHE γ-ray spectra of PKS 1510–089 during the flare night, JD 2457539. The Fermi-LAT butterflies are integrated over the precise H.E.S.S. (red) and MAGIC (green) observation windows and are plotted until the highest detected photon energy. The VHE γ-ray spectra are corrected for EBL absorption using the model of Franceschini et al. (2008), with the H.E.S.S. spectrum in red and the MAGIC spectrum in green. The grey data points are the low-state spectra from MAGIC Collaboration (2018a). For this spectrum, HE points above 1 GeV are omitted due to potential bias (for details, see MAGIC Collaboration 2018a).

In the text

	Fig. 14. Two-dimensional likelihood surface of a power-law fit to the Fermi-LAT data contemporaneous to the H.E.S.S. observations (left) and MAGIC observations (right). The best-fit values are marked in red for H.E.S.S. and green for MAGIC. The likelihood surfaces are used in a simultaneous fit with IACT data to search for BLR absorption features. See main text for further details.
In the text

Fig. 15. HE and VHE γ-ray spectra of PKS 1510–089 during the flare night, JD 2457539. The Fermi-LAT butterflies are integrated over the precise H.E.S.S. (red, left panel) and MAGIC (green, right panel) observation windows and are plotted up to the highest detected photon energy. The VHE γ-ray spectra are corrected for the EBL absorption with the model of Franceschini et al. (2008), with the H.E.S.S. spectrum in red and the MAGIC spectrum in green. The lines mark the best-fit spectra for combined Fermi-LAT and H.E.S.S. and MAGIC fits for the BLR ring geometry and different intrinsic spectra: a power law (PL), power law with sub-exponential cut-off (EPL), and a broken power law (BPL). In the PL case, curvature is only provided by the BLR absorption, whereas the EPL and BPL models include intrinsic spectral curvature. The best-fit values of the distance r from the black hole are provided in the legend.

In the text

Fig. 16. Likelihood and Δχ2 profiles for the distance r in the BLR ring geometry for the combined Fermi-LAT and H.E.S.S. (left) and Fermi-LAT and MAGIC fits (right). The spread of the profiles shows the effect of the IACT systematic uncertainties on the measured fluxes. For the H.E.S.S. and Fermi-LAT fits, the TS value is defined as twice the difference between the best-fit log-likelihood and the log-likelihood profiled over r. The dashed lines show where the TS or Δχ2 values are equal to 2.71, which corresponds to lower limits at 95% confidence.

In the text

Fig. 17. Top: HE γ-ray light curve of PKS 1510–089 in 24 h bins. Grey arrows are upper limits (95% confidence level for TS < 9.0), and the grey dashed line indicates the 11-year average. Bottom: separation of knot K16 (red) from the core A0 (black line) according to the 43 GHz maps. The red line approximates the motion of K16 with the 1σ uncertainty given by the shaded region. The same is shown for A1 (blue) and A2 (magenta). In both panels the green dashed line indicates the time of the VHE γ-ray flare.

In the text