© ESO, 2017

1. Introduction

PKS 1510−089 is a bright flat spectrum radio quasar (FSRQ) located at a redshift of z = 0.36 (Tanner et al. 1996). The source is one of only six objects firmly classified as a FSRQ from which gamma-ray emission has been detected in the very-high-energy (VHE, > 100 GeV) range (Abramowski et al. 2013). Moreover, one of the highest recorded apparent speeds of superluminal motion, up to ~ 46 c, has been seen in the ultrarelativistic jet of PKS 1510−089 (Jorstad et al. 2005). As it is in many other FSRQs, the GeV gamma-ray emission of PKS 1510−089 is strongly variable (Abdo et al. 2010; Saito et al. 2013; Aleksić et al. 2014). The doubling time of the PKS 1510−089 flares observed in the GeV range go down to 1 h (Saito et al. 2013).

Most of the FSRQs have been detected in the VHE gamma-ray range during (usually short) flares (see, e.g., MAGIC Collaboration et al. 2008; Aleksić et al. 2011; Ahnen et al. 2015). Since 2013, the Major Atmospheric Gamma Imaging Cherenkov (MAGIC) telescope system has performed regular monitoring of PKS 1510−089. Interestingly, until 2015, no variability was seen in PKS 1510−089 in VHE gamma rays in H.E.S.S. (Abramowski et al. 2013) or in MAGIC (Aleksić et al. 2014) observations. It should be noted, however, that both VHE gamma-ray detections happened during long periods of enhanced optical and GeV gamma-ray activity. Hence no low-state VHE gamma-ray emission has been established so far from PKS 1510−089.

In May 2015, a strong flare of PKS 1510−089 was observed in GeV gamma rays by the Large Area Telescope (LAT) on board the Fermi satellite. At this time the source also showed high activity in the optical (Jankowsky et al. 2015; Mirzoyan 2015) and IR bands (Sameer et al. 2015; Carrasco et al. 2015). The high state triggered further MAGIC observations, which led to the detection of an enhanced VHE gamma-ray activity from the source (Mirzoyan 2015). The VHE gamma-ray emission has been also observed by the H.E.S.S. telescope (Zacharias et al. 2016). In May 2016 another flare occurred (de Naurois 2016; Mirzoyan 2016) with an even stronger VHE gamma-ray flux than in May 2015. The May 2016 flare will be discussed in a separate paper.

In this paper we report on the observations of PKS 1510−089 during the May 2015 flare. In Sect. 2 we briefly introduce the instruments that provided multiwavelength data and describe the data reduction procedures. In Sect. 3 we present the multiwavelength behavior of the source. Section 4 is devoted to the interpretation of the data in the framework of an external Compton model. The most important results are summarized in Sect. 5.

2. Instruments, observations and data analysis

During the May 2015 outburst PKS 1510−089 was observed by various instruments in a broad range of frequencies (from radio up to VHE gamma rays). In this section we introduce the different instruments and data sets and explain the data analysis procedure.

2.1. MAGIC

MAGIC is a system of two imaging atmospheric Cherenkov telescopes with a mirror dish diameter of 17 m each. They are located in the Canary Islands, on La Palma (28.7° N, 17.9° W), at a height of 2200 m above sea level (Aleksić et al. 2016a). As PKS 1510−089 is a southern source, only observable at zenith angle > 38°, the corresponding trigger threshold is ≳ 90 GeV (Aleksić et al. 2016b), about 1.7 times larger than for the low zenith observations.

The MAGIC telescopes observed PKS 1510−089 for 5.4 h between 18 and 24 of May, 2015 (MJD 57 160−57 166). The data were analyzed using MARS, the standard analysis package of MAGIC (Zanin et al. 2013; Aleksić et al. 2016b). As part of the data set was affected by Calima¹ we applied a correction for the atmosphere transmission based on LIDAR information (Fruck & Gaug 2015).

2.2. Fermi-LAT

Fermi-LAT monitors the gamma-ray sky every 3 h in the energy range from 20 MeV to beyond 300 GeV (Atwood et al. 2009). An analysis of the publicly available Pass 8 SOURCE class events was performed for a region of interest (ROI) of 10° radius centered at the position of PKS 1510−089. In order to reduce contamination from the Earth’s limb, a zenith angle cut of < 90° was applied. The analysis was performed with the ScienceTools software package version v10r0p5 using the P8R2_SOURCE_V6² instrument response function and the gll_iem_v06 and iso_P8R2_SOURCE_V6_v06 models³ for the Galactic and isotropic diffuse emission (Acero et al. 2016), respectively.

An unbinned likelihood analysis was applied using gtlike, including in the model all 3FGL sources (Acero et al. 2015) within 20° from PKS 1510−089. The spectral indices and fluxes were left free for sources within 10°, while sources from 10° to 20° have their parameters fixed to their catalog value. A first unbinned likelihood fit was performed for the events collected within almost four months of data from 22 March 2015 to 19 July 2015 (MJD 57 103−57 223) in the energy range between 100 MeV and 800 GeV. The sources with a test statistic (TS; Mattox et al. 1996) below 5 were removed from the model. Next, the optimized output model was used to produce the light curves and spectra of PKS 1510−089 in different time bins (from 1 day to 3 h) and energy ranges (E> 100 MeV, E> 1 GeV). For the calculation of the light curves, all sources were fixed in the model except PKS 1510−089 for which both the flux normalization and the spectral index were left free and modeled as a power law. For the calculation of the spectral points, the spectral index of PKS 1510−089 was also fixed to its best-fit value during the considered time period for which the spectral points are estimated. The normalization of the Galactic and isotropic diffuse emission models was left to vary freely during the calculation of both the light curves and the spectra.

2.3. Swift-XRT and UVOT

The multi-epochs (16 individual pointings) event list obtained by the X-ray Telescope (XRT) (Burrows et al. 2004) on board the Swift satellite during the period from 11 May to 25 May 2015 (MJD 57 153−57 167) with the total exposure time of 26.6 ks were downloaded from the publicly available database table SWIFTXRLOG (Swift-XRT Instrument Log). The individual exposures ranged from 0.6 to 4 ks. They were processed using the HEASOFT package version 6.18. All the observations from this period were performed in photon counting (PC) mode. The source region was defined as a circle of 20 pixel (~47′′) radius at the center of the source, while the background region was defined by an annulus centered at the source with inner and outer radii of 40 (~94′′) and 80 pixels (~188′′), respectively. The source and background spectra were extracted using XSELECT task (v2.4c). The source spectrum count rate did not exceed 0.5 counts/s in any of the observation epochs. Therefore, no pile-up correction was needed. For the light curve analysis we combined three pairs of epochs (MJD 57 157, 57 161, 57 163) separated by ~2 h.

The xrtexpomap task (v0.2.7) was used to correct the flux loss caused because some of the CCD pixels were not used during the data collection. The xrtmkarf task (v0.6.3) took into account vignetting and bad pixels. The grppha task was used to group source spectra in such a way that each bin contains 20 counts. XSPEC task (v12.9.0i) was used to calculate the flux and power-law model spectral parameters using a fixed equivalent Galactic hydrogen column density of n_H = 6.89 × 10²⁰ [cm^-2] (Kalberla et al. 2005).

The second instrument on board the Swift satellite, the Ultraviolet/Optical Telescope (UVOT, Poole et al. 2008) was used to monitor the flux of the source in the 180−600 nm wavelength range. Following Raiteri et al. (2010) we used an iterative procedure for the data calibration, where the effective wavelength, counts-to-flux conversion factor, and Galactic extinction for each filter were calculated by taking into account the filter’s effective area and the source’s spectral shape. Of the 16 pointings in the investigated time period, 8 were taken with a full set of filters (v, b, u, w1, m2, w2). For the SED modeling we used pointings on MJD 57 160 (with u and w2 filters) contemporaneous with the MAGIC flaring state, and on MJD 57 165 (all filters available) during the post-flare MAGIC observations.

2.4. Optical photometry and polarization

PKS 1510−089 is regularly monitored as part of the Tuorla blazar monitoring program⁴ in the R band using a 35 cm Celestron telescope attached to the KVA (Kunglinga Vetenskapsakademi) telescope located at La Palma. The data analysis was performed with the semi-automatic pipeline using the standard analysis procedures (Nilsson et al., in prep.). The differential photometry was performed using the comparison star magnitudes from Villata et al. (1997). The magnitudes were corrected for the Galactic extinction using values from Schlafly & Finkbeiner (2011).

The optical polarization observations were performed with a number of instruments: Nordic Optical Telescope (NOT), Steward Observatory, Perkins Telescopes, RINGO3, AZT-8, and LX-200. The NOT polarimetric observations were done with ALFOSC in the R band using the standard setup for linear polarization observations (lambda/2 retarder followed by a calcite). The data were analyzed using the standard procedures with semi-automatic software as in Hovatta et al. (2016). The Steward polarimetric observations were obtained as part of an ongoing monitoring program of gamma-ray-bright blazars in support of the Fermi mission⁵. The observations were performed in the 5000−7000 Å  band. The data analysis pipeline is described in Smith et al. (2009). Polarimetric R-band observations were also provided by the 1.8 m Perkins telescope at Lowell Observatory equipped with PRISM (Perkins Reimaging System). The data analysis was done following the standard procedures as in Chatterjee et al. (2008). Polarization observations were also taken with the RINGO3 polarimeter (Arnold et al. 2012) on the fully robotic and autonomous Liverpool Telescope on La Palma, Canary Islands (Steele et al. 2004), as part of the Liverpool blazar monitoring campaign (see Jermak et al. 2016), in collaboration with the Monitoring AGN Polarimetry at the LIverpool Telescope (MAPLIT) program. Simultaneous observations in the blue, 350−640 nm; green, 650−760 nm; and red, 770−1000 nm, passbands were taken using a rapidly rotating (once per 4 s) polaroid which modulates the incoming beam of light in eight rotor positions. For this work we use only the 650−760 nm measurements, which are the closest to the R band used in the above-mentioned polarization instruments. The beam is simultaneously split by two dichroic mirrors into three electron multiplying CCD cameras (EMCCDs). The combination of the flux from the eight rotor positions using equations from Clarke & Neumayer (2002) can be used to find the linear Stokes parameters, which were used to calculate the degree and angle of polarization. Finally, additional optical polarimetric data are reported here from the 70 cm AZT-8 telescope (Crimea) and the 40 cm LX-200 telescope (St. Petersburg), both equipped with nearly identical imaging photometers-polarimeters (Larionov et al. 2008). Polarimetric observations were performed using two Savart plates rotated by 45° relative to each another. By swapping the plates, the observer can obtain the relative Stokes q and u parameters from the two split images of each source in the field. Instrumental polarization was found via stars located near the object under the assumption that their radiation is unpolarized. The electric vector position angle (EVPA) was corrected for the n × 180° ambiguity by minimizing the difference to the closest data point unless there is a gap of over 14 days.

2.5. Infrared

PKS 1510−089 is monitored by a number of IR instruments. We used the publicly available data in the B, V, R, J, and K bands from the Small and Moderate Aperture Research Telescope System (SMARTS) instrument located at Cerro Tololo Interamerican Observatory (CTIO) in Chile. The data reduction and calibration is described in Bonning et al. (2012). We converted the magnitudes into flux units using Bessell et al. (1998) and corrected for the interstellar dust absorption following Schlafly & Finkbeiner (2011).

The observations at Teide Observatory (Canary Islands) were obtained with the 1.52 m Carlos Sanchez Telescope (TCS), using the near-infrared camera CAIN during the nights of MJD 57 162−57 174. This camera is equipped with a 256 × 256 pixel NICMOS-3 detector providing a scale of 1′′/pixel. Data were acquired in the three filters J, H, and Ks. Observations were performed using a 5-point dither pattern (repeated twice) in order to facilitate a proper sky background subtraction. At each point, the exposure time was about 1 min, split into individual exposures of 10 s in the J filter and 6 s in the H and Ks filters to avoid saturation by sky brightness. Image reduction was performed with the caindr package under the IRAF environment⁶. Data reduction includes flat-fielding, sky subtraction, and the shift and combination of all frames taken in the same dither cycle. Photometric calibration was made based on field stars from the 2MASS catalog (Cutri et al. 2003). The photometric zero point was determined for each frame by averaging the offset between the instrumental and the 2MASS magnitudes of the catalog. Deviant stars were excluded and typical errors remained below 5%.

We also used IR photometry data obtained with the 1.2 m telescope of Mt Abu InfraRed Observatory (MIRO), India, mounted with the Near Infrared Camera and Spectrograph (NICS) equipped with 1024 × 1024 HgCdTe Hawaii array detector. The field of view is 8′ × 8′ with a pixel scale of 0.5′′/pixel. The observations on PKS 1510−089 were performed with a 4-position dither with offsets of 30 arcsec, keeping the comparison stars⁷ 1, 2, 3, 4, and 6 in the field of the source. The sky and dark contributions were removed using these dithered images and aperture photometry was performed using standard procedures under IRAF (see Banerjee & Ashok 2012, for details on data reduction and analysis). The source magnitudes in the J, H, and Ks bands were calibrated using correction factors obtained using the weighted average of the standard values of comparison stars mentioned above.

2.6. Radio

Radio monitoring observations were performed with Metsähovi Radio Telescope operating at 37 GHz frequency. The instrument and data reduction procedures are described in Teräesranta et al. (1998) and Aleksić et al. (2014).

The quasar PKS 1510−089 was observed within a sample of gamma-ray blazars that the Boston University (BU) blazar group monitors with the Very Long Baseline Array (VLBA) approximately monthly at 43 GHz (the VLBA-BU-BLAZAR project). The observations of PKS 1510−089 are usually performed via nine short scans (~ 5 min each) within a span of 7−8 h. The data were calibrated at the VLBA DiFX correlator and reduced using the Astronomical Image Process System (AIPS) and Difmap software packages, as described in Jorstad et al. (2005). The calibrated data are available online⁸. We analyzed the data obtained from February 2015 to April 2016 (12 epochs). We modeled the total intensity images by components with circular Gaussian brightness distributions. For each component we determined flux density, distance and position angle (PA) with respect to the VLBI core⁹, and size. The 43 GHz core is expected to be located at a distance of ~6.5 pc from the central engine of PKS 1510−089 (see Pushkarev et al. 2012; Aleksić et al. 2014). A map of the parsec scale jet of the quasar formed from 20 stacked images over six years of VLBA observations at 43 GHz is plotted in Fig. 1 (individual images can be found at the BU blazar group website http://www.bu.edu/blazars/VLBA_GLAST/1510.html).

Fig. 1 Stacked map of 43 GHz total (black contours) and polarized (color scale) intensity images of the inner pc-scale jet of PKS 1510−089 with the direction of electric field vector polarization denoted by black line segments (contour levels are indicated at the top, polarized flux levels are shown in the color bar to the right). All images have been convolved with the same Gaussian beam, shown in the lower right corner.

The image shows the VLBI core, which is the brightest compact feature located at the southeast end of the jet. The core is used as a reference point in the stacking procedure since it is assumed to be stationary. The stacked image reveals the full opening angle of the jet, as well as the location of the jet axis. As can be inferred from Fig. 1, the jet axis is along PA ~ −30deg, while the projected opening angle is ~ 60deg. The core at 43 GHz is only partially optically thick, which is supported by synchronous optical and radio core polarization variability in a number of blazars (D’Arcangelo et al. 2007; D’arcangelo et al. 2009). Figure 1 reveals a shift between the total and polarized intensity peaks and the complex structure of the polarized emission in the core. This favors the hypothesis that the core of PKS 1510−089 is a recollimation shock, which was inferred previously from the polarization structure in other blazars (e.g., Cawthorne et al. 2013).

3. Results

In Fig. 2 we show the multiwavelength light curve of PKS 1510−089 during the May 2015 outburst in the investigated period of MJD 57 151−57 174.

Fig. 2 Multiwavelength light curve of PKS 1510−089 during the May 2015 flare. From top to bottom: nightly gamma-ray flux above 150 GeV from MAGIC (the dashed line shows the average emission in Feb.−Apr. 2012, Aleksić et al. 2014); Fermi-LAT flux above 0.1 GeV in 6 h binning, and the corresponding spectral index (the dashed line shows the average emission from the 3FGL catalog, Acero et al. 2015); X-ray spectral flux (filled circles) and spectral index (empty circles) measured by Swift-XRT; polarization percentage and polarization angle measured by NOT, Steward, Perkins, RINGO3, AZT-8, and LX-200 (see legend); optical emission in R band (KVA, SMARTS) and UV emission in w2-band (Swift-UVOT); IR emission in J band (SMARTS, MIRO-NICS, TCS); radio observations by Metsähovi at 37 GHz. Data from IR to UV are corrected for Galactic absorption. The red and blue shaded regions show Period A and Period B, respectively, for which the spectral modeling is performed.

Following the observations of MAGIC, we define two observation periods: Period A (MJD ~ 57 160−57 161) and Period B (MJD ~ 57 164−57 166). The multiwavelength SED of both periods is investigated.

3.1. MAGIC

The MAGIC light curve (top panel of Fig. 2) shows clear variability, with the highest flux observed during the two nights of Period A. The hypothesis of constant flux during all five observation nights of MAGIC can be clearly rejected, with a chance probability of 7.7 × 10^-8. Even allowing for a 20% variable systematic uncertainty on individual night fluxes (motivated by the variable systematic uncertainty estimate given in Aleksić et al. 2016b, rescaled to a softer source) we still obtain a small value of chance probability of 3.2 × 10^-4 that the flux is constant. The flux during Period A is ~ 5 times higher than that observed during the previous detection by MAGIC in 2012 (Aleksić et al. 2014). In the following observations during Period B, the VHE gamma-ray flux decreased to a level consistent with the detection in 2012.

In order to search for a possible short time variability, we binned the light curve during Period A into 20 min bins (see Fig. 3).

Fig. 3 Light curve above 150 GeV obtained with the MAGIC telescopes during the flare in Period A. The fluxes are computed in 20 min bins. The black line shows the constant flux fit, (4.0 ± 0.5) × 10-11 cm-2s-1. The gray band shows the 95% confidence level interval allowing for 20% variable systematic uncertainty (see text for details).

No variability is detected at such a time scale. Fitting the light curve with a constant flux hypothesis we obtain χ²/N_d.o.f. = 5.9/7. We estimate the maximum variability that can be hidden by the uncertainties of the measurement by computing for each 20 min light curve bin a 95% confidence level interval on the flux using the Rolke et al. (2005) prescription. We include a 20% variable systematic uncertainty in these calculations. By comparing the least constraining upper edge of the 95% confidence level interval with the average flux from those two nights we obtained that the flux did not increase by more than a factor of 3.5 on time scales of 20 min.

For spectral analysis we combined the two nights of Period A. The obtained VHE gamma-ray spectrum of PKS 1510−089 is shown in Fig. 4.

Fig. 4 Spectral energy distribution of PKS 1510−089 constructed from the MAGIC data gathered on MJD 57 160 and 57 161 (Period A). The observed spectrum is shown by red empty squares with inclined hashing, and the EBL-deabsorbed spectrum according to the Domínguez et al. (2011) model by magenta filled circles with vertical hashing. The hashing shows the uncertainty of the forward folding with a power law. The SED constructed from the data taken on MJD between 57 164 and 57 166 (Period B) is shown as empty (observed) and filled (EBL-deabsorbed) blue circles. For comparison, MAGIC measurements performed in Feb.−Mar. 2012 and H.E.S.S. measurements from Mar. 2009 (Abramowski et al. 2013) are shown as gray diamonds and stars, respectively.

It can be described by a power law, dN/ dE = f × (E/ 200 GeV)^{− α} [cm^-2 s^-1 TeV^-1], with f = (1.96 ± 0.55_stat ± 0.36_syst) × 10^-10 and α = 4.59 ± 0.75_stat. The spectral parameters are obtained using a forward folding method (Albert et al. 2007). The statistical uncertainty on the spectral index is much larger than the typical systematic uncertainty of ± 0.15 of the observations performed with the MAGIC telescopes. The systematic uncertainty on the flux normalization does not include the uncertainty of the energy scale of MAGIC, which for this data set we estimate as ≲ 19%, slightly larger than the ≲ 15% given in Aleksić et al. (2016b), due to the need of LIDAR correction of Calima-affected data. Correcting for the absorption of TeV gamma rays due to the interaction with the extragalactic background light according to the Domínguez et al. (2011) model, an intrinsic spectrum with normalization of f_deabs = (4.2 ± 1.0_stat ± 0.76_syst) × 10^-10 and index α_deabs = 3.17 ± 0.80_stat are obtained. The spectral shape is marginally consistent (but with large uncertainties) with the previous measurements by H.E.S.S. (observed slope 5.4 ± 0.7_stat ± 0.3_syst, Abramowski et al. 2013) and MAGIC (intrinsic slope 2.5 ± 0.6_stat, Aleksić et al. 2014).

For comparison, we have also reconstructed the average flux from the MAGIC measurements performed during Period B. The combination of weaker emission and observations performed during higher atmospheric transmission results in the energy range of the reconstructed spectrum shifted to lower energies. The observed flux in Period B is at a similar level as the one detected by MAGIC in 2012. We obtained the observed and intrinsic VHE spectral slopes of PKS 1510−089 during Period B of 4.75 ± 0.62_stat and 4.33 ± 0.75_stat, respectively.

3.2. Fermi-LAT

The GeV gamma-ray flux of PKS 1510−089 is highly variable in the investigated period. A few individual flares are visible, with time scales of a few days. Owing to the short-term gamma-ray variability, in order to get a spectrum comparable to the MAGIC observations, for SED analysis we selected the events observed by the LAT within 6 h centered on each of the MAGIC observations. We calculated two GeV spectra, which correspond to the two different states of the source contemporaneously to the MAGIC observations in Period A and Period B, respectively. The best description of the GeV spectrum measured by the LAT in Period A is a power law with spectral index 2.20 ± 0.07 and a flux above 100 MeV of (6.8 ± 0.5) × 10^-6 ph cm^-2 s^-1 (TS = 842). During Period B, the Fermi-LAT measured spectrum is best described by a power law with a similar spectral index of 2.17 ± 0.08, but a significantly lower flux, (3.7 ± 0.3) × 10^-6 ph cm^-2 s^-1 (TS = 564). For comparison, the 3FGL flux above 100 MeV is (0.94 ± 0.01) × 10^-6 ph cm^-2 s^-1 (Acero et al. 2015). Therefore, PKS 1510−089 reached a factor of ~7 to 4 times its average flux over the first four years of Fermi-LAT observations during the two epochs for which the spectral analyses was performed. The GeV flux was, however, still smaller by a factor of about 2 than the daily peak flux observed in 2011 (Saito et al. 2015). The Fermi-LAT spectrum during both Period A and B was slightly harder than in the neighboring days (see also Sect. 3.5). In the case of Period B there is a weak hint of a spectral curvature. The likelihood ratio test gives 2.8σ preference (ΔTS = 8.0 for one additional degree of freedom) of the log-parabola shape over the power-law spectral model. The corresponding photon index can be described as (2.05 ± 0.10) + (0.21 ± 0.09) × ln(E/ 0.45 GeV).

To characterize the variability in the Fermi-LAT light curve we calculated the power density spectrum (PDS) for the 2015 flare epoch (MJD 57 143.9375−57 182.9375) with a 3 h binning and for a mission lifetime at the time of the analysis (MJD 54 682.655−57 484.655) with 1 day binning. In both cases the flux calculation was performed in the 0.1−800 GeV energy range.

The estimated white noise level, based on the data error values was subtracted from each PDS. The PDS is rebinned in logarithmic frequency intervals, and normalized to variance per frequency unit divided by the square of the mean flux.

Fig. 5 Power density spectrum normalized to variance per frequency unit for the mission-long Fermi-LAT light curve (black) and for the 2015 flare epoch (magenta).

The PDS level, computed with this normalization, is similar for the 2015 flare and for the mission lifetime light curve (see Fig. 5). This suggests that the fractional variability is the same and presumably driven by the same variability process. The shape of the overall PDS can be described by a power law for frequencies above 0.01 day^-1. A power law fitted to the 0.007−0.5 day^-1 frequency range PDS gives an index of 1.14 ± 0.07 for full data set and 0.97 ± 0.30 for the 2015 flare. The uncertainty value of the fit is based on the scatter of the measurements with respect to the fitted line only and therefore may be underestimated. The power-law index is similar to that of other FSRQs (Ackermann et al. 2011). The stochastic nature of the variability together with the limited observation length lead to a large uncertainty in the PDS at the lowest frequencies.

3.3. Swift-XRT

The X-ray flux, measured by Swift-XRT, shows a gradual decrease in the observed flux during the studied period. Except for the first point at MJD 57 153.6, which happened before the two large Fermi-LAT flares, the X-ray flux can be much better described by an exponential decline (χ²/N_d.o.f. = 22.6/10) than a constant value (χ²/N_dof = 47.3/11). Similarly, the corresponding X-ray spectral indices are better described by a linear softening of the spectrum (χ²/N_dof = 41.0/10) than a constant (χ²/N_d.o.f. = 49.9/11).

The flux is marginally higher during MJD 57 156−57 162, when a broad Fermi-LAT flare occurred. XRT observations on MJD 57 166−57 168, during the following broad Fermi flare did not show a clear increase in the X-ray flux. For spectral analysis we combined the four pointings taken during Period A, each of them within 6−8 h of the MAGIC observations. The X-ray spectrum in Period A can be marginally well described (χ²/N_d.o.f. = 71/53) with a power law with an index of 1.48 ± 0.05. Swift-XRT observations performed during Period B resulted in a spectrum that can be well fitted with a power law with a significantly softer index of 1.70 ± 0.04 with χ²/N_d.o.f. = 110/89, or by a curved spectrum with an index of 2.23 ± 0.14−(0.366 ± 0.093) × log (E/ 0.3 keV) with χ²/N_dof = 96/88.

3.4. IR, optical, and UVOT

The optical-UV SED of PKS 1510−089 consists of an unpolarized, quasi-stable accretion disk component and a non-thermal jet emission. The variability and polarization of the jet component might be diluted by the accretion disk component. The effect is strongest in the UV bands, and weak in the R band where the polarization is measured. It does not affect the timing of the observed polarization variability.

The optical emission of PKS 1510−089 during the investigated period is clearly variable (with a factor of 2 difference between the lowest and highest flux). The optical variability does not strictly follow the gamma-ray variability. In the optical R-band, and to a lesser extent also in the UV w2-band, the flux slowly increased throughout May 2015. Similarly to the optical range, the IR flux doubled during the ~ 8 days before Period A. At the end of the observation period it returned to the pre-flare state.

Throughout the investigated period, a smooth rotation of optical EVPA by ~ 100° occurred. The rotations of optical polarization angle also accompanied the 2009 and 2012 gamma-ray flaring states (Marscher et al. 2010; Aleksić et al. 2014). However, the rotations of the EVPA seem very common in PKS 1510−089, e.g., recent work by Jermak et al. (2016) also identified a rotation in the 2011 data; therefore, further data are needed to firmly associate them with the emission of VHE gamma rays. The low percentage of polarization seems to be typical for this source (Jermak et al. 2016). Indeed, during Period A the polarization percentage is low (~ 5%). It triples between period A and B. In the few days before Period A (i.e., during the raising part of the Fermi-LAT flare that culminated with detection of VHE gamma-ray emission by MAGIC), a higher polarization was also observed. The polarization behavior during the 2015 flaring period agrees with what is expected from a knot following a spiral path through a mainly toroidal magnetic field (Marscher et al. 2010). An alternative explanation might be the light travel time effects within an axisymmetric emission region pervaded by a predominately helical magnetic field (Zhang et al. 2015). The evolution of the polarization is probably further complicated by the superposition of individual flares seen in Fermi-LAT.

3.5. Radio

PKS 1510−089 shows moderate variability in the Metsähovi observations performed at 37 GHz in May 2015. The sampling is, however, rather sparse, and especially the local peaks of the GeV flux are mostly not covered by the observations. The observations on MJD 57 161 during Period A are burdened with a larger uncertainty, due to adverse atmospheric conditions.

Figure 6 shows the light curve of the core from February 2015 to April 2016.

Fig. 6 VLBA light curves of the core (circles) and knot K15 (stars) at 43 GHz. The vertical dashed line and gray shaded region show the zero separation epoch of K15 and its uncertainty.

The light curve reveals a significant gradual increase in flux in the second half of 2015 that is most likely connected with a disturbance (knot K15) detected in the VLBA images starting in December 2015 (Fig. 7).

Fig. 7 Total intensity images of the PKS1510-089 core region at 43 GHz, with a global peak intensity of Ipeak = 3.566 Jy/beam and 0.15 mas FWHM circular Gaussian restoring beam (bottom right circle). The solid and dashed lines designate positions of the VLBI core and K15, respectively, across the epochs.

Knot K15 is bright and relatively slow, with an apparent speed β_app = (5.3 ± 1.4) c. According to currently available data, the knot was ejected on MJD 57 230 ± 52 (see Fig. 8).

Fig. 8 Separation of knot K15 from the core.

A similar behavior has also been observed during a high gamma-ray state in Feb.−Apr. 2012, when the emergence of a new radio knot (K12) from the core was associated with a VHE outburst (Aleksić et al. 2014). The large uncertainty in the ejection time of K15 does not allow us to associate it firmly with a particular GeV flare as the source showed activity in Fermi-LAT close to the time of K15 separation (e.g., the hard flare on MJD 57 245 in Fig. 9).

Fig. 9 Fermi-LAT light curve above 0.1 GeV (top panel) binned in 3-day intervals and the corresponding evolution of spectral index (middle panel) and the optical R-band flux observed by KVA (bottom panel). The red and blue hashed bands show Period A and B, respectively, while the gray band shows the 1σ uncertainty on the K15 knot separation time.

Interestingly, the position angle of K15 is ~+50°, while the typical projected direction of the pc-scale jet is ~−30° (see Fig. 1), where seven knots have been observed to date (see, e.g., Jorstad et al. 2005; Aleksić et al. 2014). On the other hand, knot K11, ejected in October 2011 (before the 2012 VHE outburst), had a similar PA to the most recently observed one for K15 close to the core, but turned toward the usual pc-scale jet direction a few months later. We note that K15 is a factor of a few times slower than K11, and so may eventually follow a similar trajectory. The slower apparent speed and brighter flux of K15 suggest that the velocity vector of this disturbance is closer to the line of sight than K11, causing it to have a higher Doppler factor.

4. Modeling the spectral energy distribution

The gamma-ray emission of FSRQs is usually explained as the effect of the inverse Compton scattering of electrons on a radiation field external to the jet (see, e.g., Sikora et al. 1994; Ghisellini et al. 2010), the so-called external Compton (EC) scenario. The radiation field can originate from the accretion disk, broad-line region (BLR) or the dust torus (DT). This scenario has been applied to explain the emission of PKS 1510−089 in its previous flaring episodes (Abdo et al. 2010). The origin of the radiation field is closely connected to the location of the emission region. Moreover, the observed VHE gamma rays escaping from the emission region suggest that the emission region is located outside the BLR in order to escape the absorption by e⁺e⁻ pair production process on BLR photons (Abramowski et al. 2013; Aleksić et al. 2014). Dotson et al. (2015) investigated the 2009 GeV flares of PKS 1510−089 using the energy dependence of the flare decay time as the diagnostic for the emission zone location. They claimed that two of the flares happened around the distance of the DT, while the other two came from the vicinity of VLBI core. On the other hand, the modeling of GeV flares seen from PKS 1510−089 in 2011 placed the emission region at the distance of 0.3−3 pc from the black hole, with the EC process happening on a mixture of BLR and DT radiation fields (Saito et al. 2015). A similar location (~ 1 pc from the black hole) with EC mainly on DT was invoked to explain the high optical and gamma-ray state observed at the beginning of 2012 (Aleksić et al. 2014). The broadband emission could be also explained by a much more distant region (~ 20 pc) for a two-zone model where the jet consists of the inner spine and outer sheath layer (Aleksić et al. 2014; MacDonald et al. 2015).

Fig. 10 Multiwavelength spectral energy distribution of PKS 1510−089 in Periods A (red symbols) and B (blue symbols). The red and the blue curves show the result of the emission model for the two periods. The black dashed and long-dashed lines show the adopted emission for the accretion disk and the dusty torus, respectively. For comparison, the dashed gray line shows the model derived for the SED in 2012 (from Aleksić et al. 2014). Historical measurements (ASDC, see http://www.asdc.asi.it/) are shown as gray points.

In Fig. 10 we present the two spectral energy distributions of PKS 1510−089 constructed from the data taken in Periods A and B, corresponding to high and low gamma-ray flux, respectively. As can be seen, most of the flux variation (≈a factor of 2−3) occurs in the Fermi-LAT and MAGIC bands. The low-energy flux (optical, X-rays) is almost constant between the two periods. Remarkably, the high-energy peak during the period B has a very similar level to the 2012 high state, even though the IR−UV emission is a factor of ~ 3 higher.

We model these SEDs of PKS 1510−089 in the framework of the same one-zone model, as discussed in Aleksić et al. (2014). To explain the sub-TeV emission observed by MAGIC we assume that the emission region is located beyond the BLR radius where the external photon field is dominated by the thermal IR radiation of the DT.

For the setup of the model we assume the scaling laws and the prescriptions given in Ghisellini & Tavecchio (2009). Specifically, the radius of the BLR is given by cm and that of the torus by cm, where is the accretion disk luminosity in units of 10⁴⁵ erg s^-1. We calculate the IR radiation energy density assuming that a fraction f_IR = 0.6 of the disk radiation is intercepted and reprocessed by the torus heated to 1000 K. Similarly, the BLR intercepts f_BLR = 0.1 of the disk radiation. With these prescriptions, we note that the energy densities of the BLR and torus radiation fields do not depend on the disk luminosity since they depend on the constant ratio . Assuming the same disk luminosity L_disk = 6.7 × 10⁴⁵ erg s^-1 as in Aleksić et al. (2014) the scaling law of Ghisellini & Tavecchio (2009) (based on reverberation mapping measurements of BLR size; see, e.g., Bentz et al. 2009) allows us to infer a BLR radius of R_BLR = 2.6 × 10¹⁷ cm.

We fix the distance of the emission region from the base of the jet to r = 6 × 10¹⁷ cm. The emission region is filling the whole cross section of the jet, which for an assumed jet semi-aperture angle θ_j = 0.047 rad, results in the emission region radius R = 2.8 × 10¹⁶ cm. This size for the emission region is in line, even for moderate values of the Doppler factor, with the constraints set by the time scale variability of a few days observed by MAGIC.

The apparent superluminal motion of radio component K15 puts limits on Γ ≳ 5. The large uncertainty in the separation time from the radio core does not allow us to firmly associate such limits with the speed of the emission region responsible for the emission in investigated Periods A and B. Moreover, the apparent speed is measured over a much longer time scale, when the emission region might have decelerated. On the other hand, the beaming of the emission is constrained by the observed luminosity of the dominating high-energy peak. It can be estimated as , where is the energy density of the external radiation field measured in the frame of reference of the blob, n_e is the number density of the electrons, V = 4/3πR³ is the volume of the emission region measured in its own frame of reference, and δ is the Doppler factor of the blob. The average squared Lorentz factor of the electrons ⟨ γ² ⟩ can be approximated as γ_b if the distribution starts from γ_min ≈ 1 and follows an index of 2 up to the break of γ_b. Assuming Γ ≈ δ, the total kinetic power of a jet composed of cold protons with number density of n_p is P_jet = πR²n_pm_pc³Γ². Combining the two formulas, for a jet with n_e ≈ n_p we obtain .

The observed luminosity of ~ 4 × 10⁴⁷ erg s^-1 of the EC peak requires the jet Lorentz factor to be at least 10 even for the case of P_jet = 10L_disk. For the modeling we apply the same values as used in Aleksić et al. (2014) for the jet bulk Lorentz factor Γ = 20 and δ = 25. For such a bulk Lorentz factor the assumed jet semi-aperture angle of 0.047 rad is broader than 0.01(Γ/20)^-1 suggested by radio observations (see, e.g., Jorstad et al. 2005). The radio observations are sensitive to the jet opening angle at a distance equal to or greater than the location of the radio core, i.e., a few pc from the base of the jet. On the other hand, the emission region assumed in the modeling is located closer to the base of the jet, where the jet opening angle can be larger, as observed, for example in the case of the M87 radio galaxy (Asada & Nakamura 2012).

Having fixed these values, the free parameters of the model are only the intensity of the magnetic field B and those describing the electron energy distribution. Hence, the observed variability according to this scenario is caused by changes in the conditions of the plasma flowing through the shock region. Since we assume that the emission occurs outside the intense radiation field of the BLR, the IC emission and the radiative losses of the emitting electrons are dominated by the scattering of the IR radiation field of the torus. As the energy density of this radiation field is relatively low, the cooling of the electrons is not very effective. A simple calculation shows that the Lorentz factor at which the cooling time equals the dynamical time R/c is on the order of . If we assume that the electrons are injected starting from γ_inj,min = 1, with a power-law distribution with slope n_inn, in equilibrium we would expect a break at γ_cool, above which the slope of the distribution would steepen to n_inj + 1. However, such a break could not properly describe the SED, since the required break (estimated using the X-ray and the MAGIC slopes) is larger (even taking into account Klein-Nishina effects on the spectrum in the MAGIC energy range). To reproduce the SED we therefore assume a scenario in which the electron energy distribution is a double broken power law; there is a cooling break at γ_c, and a second break (most probably connected with the acceleration process) at γ_b. The particles are injected into the emission region with a broken power-law energy distribution with slopes and , and break Lorentz factor γ_b. In equilibrium conditions, the electron energy distribution displays three power laws with slopes from γ_min to the cooling electron Lorentz factor γ_c, from γ_c to γ_b, and above γ_b (see, e.g., Ghisellini et al. 2002). We note that similar hard spectra have been also postulated when modeling blue blazars (Ghisellini et al. 2012).

The values of these parameters required to reproduce the SEDs are listed in Table 1. In fact, the difference in the emission between Period A and B can be explained with a relatively small change in the fit parameters, namely a slightly stronger magnetic field and lower maximum and break energies of the electrons

The model discussed here has some caveats. As typically happens in blazar modeling, the radio points overshoot the model line, which has a strong low-energy cutoff due to synchrotron self-absorption. This emission is normally thought to occur from much larger regions farther along the jet. Moreover, the variations in the optical emission occur on longer time scales than the flares observed in GeV. Hence, additional optical emission might be produced by the high-energy electrons swept with the flow farther along the jet (up to a few pc from the base of the jet). In those regions the external radiation field density would be too low to efficiently produce high-energy photons via the inverse Compton process, turning synchrotron emission into the dominant radiation process. Even though electrons can suffer adiabatic energy losses, the total observed synchrotron radiation might still slowly increase, due to aggregation of electrons from multiple individual flares. In fact, if the electrons reach the radio core (and beyond) located at the distance of d_core they could be responsible for the emission of a new radio knot. The separation time of such a knot from the core could be estimated to occur (1 + z)d_core/ (cδΓ) = 21(d_core/ 6.5pc)(δ/ 25)^-1(Γ/20)^-1 days after the gamma-ray flare. The extent of the core might shorten the time delay before the knot–core interaction starts. We encourage further trials of modeling of the observed high state with scenarios employing emission from a greater length of the jet than a single active region.

5. Discussion and conclusions

Using the MAGIC telescopes data, we detected enhanced VHE gamma-ray emission from the direction of PKS 1510−089 during the high optical and GeV state of the source in May 2015. It was the first time that VHE gamma-ray variability was detected for this source. The spectral shape is, however, consistent within the statistical uncertainties with the previous measurements of the source. During May 2015 the IR through UV data showed a gradual increase in flux, while the flux in the X-ray range was slowly decreasing.

The May 2015 data revealed, similarly to the 2012 data, that the enhanced VHE gamma-ray emission occurred during the rotation of the optical polarization angle. Also, similarly to other gamma-ray flares, an ejection of a new radio component was observed; however, it had an unusual position angle. This suggests that the association of VHE gamma-ray emission with the rotation of EVPA and ejection of a new radio component might be a common feature of PKS 1510−089.

The source was modeled with the external Compton scenario. The evolution of the state of the source from the VHE gamma-ray flare to a lower emission (at the 2012 high state level reported in Aleksić et al. 2014) can be explained by relatively small changes in the conditions of the plasma flowing through the emission region. The presented scenario is, however, only one possible solution. As discussed in Aleksić et al. (2014), if we assume that the VHE flaring is indeed connected to the ejection of the new component (in this case K15) from the VLBA core and the rotation of the optical polarization angle, it would be natural to assume a single emission region located far outside the dusty torus. In this case the seed photons could be provided by the sheath of the jet and this scenario has been shown to provide a feasible description of the previous flaring epochs of PKS 1510−089 (Aleksić et al. 2014; MacDonald et al. 2015). The VHE gamma-ray variability with time scale τ seen during the 2015 outburst puts constraints on the size, and therefore also on the location of the emission region. Assuming that the spine of the jet fills a significant fraction of the jet (as in Aleksić et al. 2014), the location of the emission region cannot be farther than d = τδc/⁽(1 + z)θ_j⁾ = 2.7(t/ 3 days)(δ/ 25)(θ_j/ 1°)^-1 pc. Therefore, a high Doppler factor and a narrow jet are needed in order to place the emission region at the radio core. In fact, such low values of the jet extension, (0.2 ± 0.2)° (Jorstad et al. 2005) and 0.9° (Pushkarev et al. 2009) at the radio core are reported by the radio observations. Alternatively, the inner spine can be much narrower than the whole jet, as suggested by MacDonald et al. (2015).

To further study the connection of VHE emission with events at lower frequencies, long-term monitoring data are needed and this question will be addressed in a future publication. With the detection of this flare from PKS 1510−089, VHE gamma-ray variability (on time scales varying from tens of minutes to days) has been observed in all FSRQs known in VHE gamma rays. Fast-varying VHE gamma-ray emission is common among the brightest gamma-ray FSRQs. As it seems that most of the gamma-ray FSRQs can only be detected during these flares, it is not surprising at all that only a handful have been detected in VHE gamma-rays.

Acknowledgments

The MAGIC collaboration 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 (FPA2015-69818-P, FPA2012-36668, FPA2015-68278-P, FPA2015-69210-C6-2-R, FPA2015-69210-C6-4-R, FPA2015-69210-C6-6-R, AYA2013-47447-C3-1-P, AYA2015-71042-P, ESP2015-71662-C2-2-P, CSD2009-00064), and the Japanese JSPS and MEXT is gratefully acknowledged. This work was also supported by the Spanish Centro de Excelencia “Severo Ochoa” SEV-2012-0234 and SEV-2015-0548, and Unidad de Excelencia “María de Maeztu” MDM-2014-0369, by grant 268740 of the Academy of Finland, by the Croatian Science Foundation (HrZZ) Project 09/176 and the University of Rijeka Project 13.12.1.3.02, by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3, and by the Polish MNiSzW grant 745/N-HESS-MAGIC/2010/0. The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT, as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States; the Commissariat à l’Énergie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France; the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy; the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK), and Japan Aerospace Exploration Agency (JAXA) in Japan; and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden. The Metsähovi team acknowledges the support from the Academy of Finland to our observing projects (numbers 212656, 210338, 121148, and others). This paper has made use of up-to-date SMARTS optical/near-infrared light curves that are available at www.astro.yale.edu/smarts/glast/home.php. MIRO is operated by Physical Research Laboratory, Ahmedabad with support from Dept. of Space, government of India. This article is based on observations made with the 1.5 m telescope Carlos Sánchez operated by the Instituto de Astrofisica de Canarias in the Teide Observatory. The data presented here were obtained in part with ALFOSC, which is provided by the Instituto de Astrofisica de Andalucia (IAA) under a joint agreement with the University of Copenhagen and NOTSA. The Liverpool Telescope is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council. St. Petersburg University team acknowledges support from Russian RFBR grant 15-02-00949 and St.Petersburg University research grant 6.38.335.2015. The BU group acknowledges support by NASA under Fermi Guest Investigator grants NNX11AQ03G and NNX14AQ58G. 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 under cooperative agreement by Associated Universities, Inc. We would like to thank the anonymous referee for the comments on the manuscript.

References

All Tables

All Figures

Fig. 1 Stacked map of 43 GHz total (black contours) and polarized (color scale) intensity images of the inner pc-scale jet of PKS 1510−089 with the direction of electric field vector polarization denoted by black line segments (contour levels are indicated at the top, polarized flux levels are shown in the color bar to the right). All images have been convolved with the same Gaussian beam, shown in the lower right corner.

In the text

Fig. 2 Multiwavelength light curve of PKS 1510−089 during the May 2015 flare. From top to bottom: nightly gamma-ray flux above 150 GeV from MAGIC (the dashed line shows the average emission in Feb.−Apr. 2012, Aleksić et al. 2014); Fermi-LAT flux above 0.1 GeV in 6 h binning, and the corresponding spectral index (the dashed line shows the average emission from the 3FGL catalog, Acero et al. 2015); X-ray spectral flux (filled circles) and spectral index (empty circles) measured by Swift-XRT; polarization percentage and polarization angle measured by NOT, Steward, Perkins, RINGO3, AZT-8, and LX-200 (see legend); optical emission in R band (KVA, SMARTS) and UV emission in w2-band (Swift-UVOT); IR emission in J band (SMARTS, MIRO-NICS, TCS); radio observations by Metsähovi at 37 GHz. Data from IR to UV are corrected for Galactic absorption. The red and blue shaded regions show Period A and Period B, respectively, for which the spectral modeling is performed.

In the text

	Fig. 3 Light curve above 150 GeV obtained with the MAGIC telescopes during the flare in Period A. The fluxes are computed in 20 min bins. The black line shows the constant flux fit, (4.0 ± 0.5) × 10-11 cm-2s-1. The gray band shows the 95% confidence level interval allowing for 20% variable systematic uncertainty (see text for details).
In the text

Fig. 4 Spectral energy distribution of PKS 1510−089 constructed from the MAGIC data gathered on MJD 57 160 and 57 161 (Period A). The observed spectrum is shown by red empty squares with inclined hashing, and the EBL-deabsorbed spectrum according to the Domínguez et al. (2011) model by magenta filled circles with vertical hashing. The hashing shows the uncertainty of the forward folding with a power law. The SED constructed from the data taken on MJD between 57 164 and 57 166 (Period B) is shown as empty (observed) and filled (EBL-deabsorbed) blue circles. For comparison, MAGIC measurements performed in Feb.−Mar. 2012 and H.E.S.S. measurements from Mar. 2009 (Abramowski et al. 2013) are shown as gray diamonds and stars, respectively.

In the text

	Fig. 5 Power density spectrum normalized to variance per frequency unit for the mission-long Fermi-LAT light curve (black) and for the 2015 flare epoch (magenta).
In the text

	Fig. 6 VLBA light curves of the core (circles) and knot K15 (stars) at 43 GHz. The vertical dashed line and gray shaded region show the zero separation epoch of K15 and its uncertainty.
In the text

	Fig. 7 Total intensity images of the PKS1510-089 core region at 43 GHz, with a global peak intensity of Ipeak = 3.566 Jy/beam and 0.15 mas FWHM circular Gaussian restoring beam (bottom right circle). The solid and dashed lines designate positions of the VLBI core and K15, respectively, across the epochs.
In the text

	Fig. 8 Separation of knot K15 from the core.
In the text

	Fig. 9 Fermi-LAT light curve above 0.1 GeV (top panel) binned in 3-day intervals and the corresponding evolution of spectral index (middle panel) and the optical R-band flux observed by KVA (bottom panel). The red and blue hashed bands show Period A and B, respectively, while the gray band shows the 1σ uncertainty on the K15 knot separation time.
In the text

Fig. 10 Multiwavelength spectral energy distribution of PKS 1510−089 in Periods A (red symbols) and B (blue symbols). The red and the blue curves show the result of the emission model for the two periods. The black dashed and long-dashed lines show the adopted emission for the accretion disk and the dusty torus, respectively. For comparison, the dashed gray line shows the model derived for the SED in 2012 (from Aleksić et al. 2014). Historical measurements (ASDC, see http://www.asdc.asi.it/) are shown as gray points.

In the text