© ESO 2020

1. Introduction

Blazars (Urry & Padovani 1995) are a sub-class of active galactic nuclei (AGNs) that exhibit relativistic jets closely aligned to the line of sight of an observer on Earth. Due to strong relativistic beaming effects in this geometry, the intensity from these sources appears greatly boosted in the observer frame and is dominated by the non-thermal continuum produced within the jet. Blazars are characterised by rapid variability across the entire non-thermal waveband that spans a wide energy range from radio to very-high-energy (VHE, E >  100 GeV) γ rays. BL Lacs are a special class of blazar showing extremely weak or no emission lines in their optical and ultraviolet (UV) spectra.

Multi-wavelength (MWL) observations show that the non-thermal emission spectral energy distribution (SED) of blazars usually exhibits a double-peaked structure (e.g. Pian et al. 1998). The first SED peak is commonly attributed to synchrotron radiation of relativistic electrons located inside an emitting region within the jet and moving relativistically towards the observer with bulk Lorentz factor Γ_bulk. The origin of the high-energy peak in the SED is debatable. Within the framework of leptonic models, the origin of this component is often ascribed to inverse Compton (IC) up-scattering of low-energy photons by high-energy electrons. For the BL Lacs the synchrotron photons present within the jet are commonly believed to serve as seeds for IC up-scattering, known as the synchrotron self-Compton (SSC) scenario (Konigl 1981; Maraschi et al. 1992; Dermer & Schlickeiser 1993; Ghisellini et al. 1998). Alternatively, in hadronic scenarios the second SED peak is attributed to relativistic protons accelerated within the jet, either via synchrotron radiation (Aharonian 2000) or via secondary emission from electron-positron pairs generated in inelastic collisions between a high-energy proton and ambient matter (pp interactions; Ackermann et al. 2013) or (internal and/or external) low-energy photon fields (p-γ interactions, Mannheim 1993; Mücke et al. 2003; Böttcher 2007). The pp interaction channel is generally neglected for blazars because the particle density inside the emitting region is considered to be too low. According to the location of the first SED peak, blazars are further classified (Padovani & Giommi 1995) into high-frequency-peaked BL Lac objects (HBLs, peaking at X-ray frequencies), intermediate- and low-frequency-peaked BL Lac objects (IBLs and LBLs, peaking at optical-infrared frequencies). Usually the TeV-loud blazars have the first peak of the SED at UV–X-ray energies and belong to the class of HBLs.

1ES 1959+650 is a well-known HBL object located nearby with a redshift z = 0.048 (Perlman et al. 1996). It was first detected in the radio band by the NRAO Green Bank Telescope (Gregory & Condon 1991) and in the X-ray band by the Einstein IPC Slew Survey (Elvis et al. 1992). Its first detection at TeV energies was by the Utah Seven Telescope Array experiment (Nishiyama 1999). This source has also been detected in high-energy (HE; 100 MeV < E< 10 GeV) γ rays with the Fermi Large Area Telescope (Fermi-LAT; Acero et al. 2015). During May–July 2002 the source exhibited strong flaring activities and flux variations in the VHE band (Aharonian et al. 2003; Holder et al. 2002; Daniel et al. 2005). Krawczynski et al. (2004) performed a MWL campaign during this period including TeV γ-ray, X-ray, optical, and radio observations and reported the detection of an “orphan flare” on 4 June 2002, a strong outburst in the VHE γ-ray band without a simultaneous X-ray counterpart. The authors reported a correlation between the X-ray and γ-ray fluxes in general except during the orphan flare. Correlated variability between these two energy bands can usually be explained by standard leptonic models whereas the lack of such a correlation challenges the SSC interpretation of the VHE flux. Hence, investigations of correlation between these two energy bands are particularly interesting for bright TeV HBLs such as 1ES 1959+650. The origin of the TeV orphan flare detected in 2002 was explained by Böttcher (2005) using a hadronic synchrotron mirror model where the flare is produced due to the interaction of relativistic protons inside the jet with external photons supplied by the reflected electron-synchrotron emission from nearby gas clouds. The source was later detected in a low VHE state by the MAGIC telescopes in 2004 and during the 2006 MWL campaign. The integral flux above 180 GeV is (4.7 ± 0.5 ± 1.6)×10⁻¹¹ cm² ⋅ s⁻¹, which is equal to ∼20% of the Crab Nebula flux¹ above the same energy threshold (Albert et al. 2006) during the former observation. The integral flux above 300 GeV is ∼10% of the Crab Nebula flux during the latter observation (Tagliaferri et al. 2008). Aliu et al. (2013) report the source detection by VERITAS with a significance of 16.4σ in a low VHE flux state and higher X-ray variability compared to other energy bands. Another strong VHE outburst during 2012 was reported by Aliu et al. (2014) where an increased VHE flux was observed without a simultaneous activity in the X-ray and UV fluxes that could be explained by the reflected emission model similar to Böttcher (2005). Regarding the long-term behaviour in X-rays, the photon spectral index scatters below and above 2 when the spectra are fitted by a power-law (PL) function (Krawczynski et al. 2004; Patel et al. 2018). This behaviour of the PL index suggests that the low-energy peak of the SED moves around the X-ray band because the photon index of a PL is ∼2 at the peak energy of the SED. Significant flux variabilities are also seen in HE γ rays with the Fermi-LAT (Ciprini & Fermi Large Area Telescope Collaboration 2015; Patel et al. 2018), although the determination of the spectral shape requires an observation longer than one day. The γ-ray spectrum in 100 MeV–100 GeV is expressed by a PL function according to the LAT 4-year point source catalogue (3FGL; Acero et al. 2015), and the photon power-law index is 1.88 ± 0.02. In the 8-year point source catalogue (4FGL²), the power-law index is reported to be 1.82 ± 0.01 in 50 MeV–1 TeV although a log-parabola (LogP) is preferred to a PL as the spectral type. The index in 10 GeV–2 TeV is 1.94 ± 0.06 for 7-year data (the 3FHL catalogue; Ajello et al. 2017). Therefore, the high-energy peak of the SED is anticipated to be located above 10 GeV

The first potential association between a HE neutrino event and the flaring blazar TXS 0506+056 (IceCube Collaboration 2018) has inaugurated a new era in multi-messenger astronomy and triggered many studies related to the neutrino-blazar coincidence (Ansoldi et al. 2018; Cerruti et al. 2018; Keivani et al. 2018; Gao et al. 2019). 1ES 1959+650 is also an interesting candidate for neutrino point-source search by IceCube. In 2005 the AMANDA neutrino telescope reported the detection of neutrinos with a hint of spatial correlation with the source direction (Halzen & Hooper 2005), although the observations were not statistically significant. The most recent IceCube analysis, spanning 8 years of data however, results in a local p-value at the position of 1ES 1959+650 of p ∼ 0.25 (IceCube Collaboration 2019), which is consistent with the background-only hypothesis.

Since 2015 the source entered into an active state across several energy bands, most notably in optical (Baliyan et al. 2016), X-rays (Kapanadze 2015; Kapanadze et al. 2016) and also γ rays, as reported by the preliminary data analysis from the MAGIC, Fermi-LAT, FACT, and VERITAS collaborations (e.g. Buson et al. 2016; Biland et al. 2016; Biland & Mirzoyan 2016). In this paper we report the results of a MWL campaign led by the MAGIC collaboration during April–November 2016 when the source was in an active state. The MAGIC telescopes observed three major VHE γ-ray flares from this source on 13 and 14 June, and 1 July 2016. The flaring events of 1ES 1950+650 are particularly interesting because the close proximity and brightness of the source allow us to perform detailed spectral measurements up to TeV energies, study the flux variability patterns, test emission models related to the origin of the VHE γ rays, and investigate their connection to cosmic ray and neutrino production. Apart from 1ES 1959+650, such a detailed analysis with short-timescale variations is only possible for very bright and nearby sources such as Mrk 421 and Mrk 501 (e.g. Tavecchio et al. 2001; Ahnen et al. 2018). The main focus of this paper is devoted to the extreme flaring events of 1ES 1959+650 from 2016, their MWL spectral and temporal properties, and the investigation of their broadband characteristics to understand the underlying physical conditions inside the source during the flares.

The paper is organised as follows. Section 2 describes the details of data analysis methods across all wavebands. Section 3 reports the results from the spectral and temporal analysis in the VHE band and from instruments observing at lower frequencies along with an investigation of the intra-night variability behaviour. Section 4 discusses the dimension of the emission region and broadband SED modelling of the flaring events. Section 5 ends with the summary and conclusions.

2. Observations and data analysis

For the present study, we performed two kinds of multi-wavelength data analysis. One is a long-term analysis of the source flux in four wavebands and the other is a quasi-simultaneous spectral analysis at the two flux peaks of 13 and 14 June 2016. A Swift observation was performed during the MAGIC observation on each of these days, and hence we have quasi-simultaneous data from UV to VHE only for them amongst the three high VHE-flux days. In the following we report a brief explanation of the instrumentation involved in these campaigns; we describe the observations performed and the data analysis techniques adopted.

2.1. MAGIC Telescopes

MAGIC is an array of two Imaging Atmospheric Cherenkov Telescopes (IACTs) located at La Palma in the Canary Islands, Spain (Aleksić et al. 2016). The location coordinates are 28°.7N, 17°.9W and the altitude is 2 200 m above sea level. The diameter of each telescope dish is 17 m. The standard trigger energy threshold for low zenith angle observations is ∼50 GeV (Aleksić et al. 2016).

Our dataset was collected by applying two kinds of observation strategy. One is dedicated monitoring of 1ES 1959+650 in 2016. Intensive observations were also triggered by high flux states of the source in optical, X-ray, HE γ ray, and VHE γ ray. The effective observation time reached 72 h over 67 nights between 29 April (MJD 57507) and 29 November (MJD 57721) 2016, including high and low flux states of the source. The majority of data from this source was taken with a zenith angle ranging from 35° to 50°. For such a zenith angle the energy threshold is higher than the value mentioned in the previous paragraph, for instance ∼100 GeV for a zenith angle of 40° (Aleksić et al. 2016). Some of the data were taken under moonlight. For these data the level of background noise detected by every pixel in the MAGIC cameras increases, which affects the overall shape and parameters of the Cherenkov shower image. This mainly imposes non-triviality to discriminate the γ rays from the background, and to reconstruct their energy and arrival direction. Consequently, the analysis energy threshold increases. The main peculiarities and details of such data analysis are described in Ahnen et al. (2017).

Observation of Cherenkov light is affected by the atmospheric transparency. The MAGIC collaboration has a self-made LIDAR system for monitoring the atmospheric transmission (Fruck & Gaug 2015). For some of the observation nights, the LIDAR information was not available due to technical problems. The transmission condition can also be estimated with thermal radiation measured with a pyrometer (Will 2017) and with the number of stars detected by a CCD camera, which is installed mainly for monitoring the telescope mis-pointing (Riegel et al. 2005; Bretz et al. 2009).

We analysed the data using the MAGIC Standard Analysis Software (MARS; Moralejo et al. 2009; Zanin 2013). Data with aerosol transmission levels from a height of 9 km above the ground level of MAGIC lower than 75%, with too high background light rate due to the moon (above 4.5 times brighter than dark conditions), and with zenith angles larger than 45° were discarded to keep a low analysis energy threshold for the spectral analyses, ∼100 GeV for dark conditions. After these quality cuts, ∼62 h of data over 61 nights from 29 April to 21 November remained for further analysis. For data selection based on the atmospheric transparency, the pyrometer data and the number of stars were used in addition to the LIDAR data to validate the data quality of the nights without LIDAR. For all the nights with the LIDAR data, atmospheric transmission correction based on the information was applied (Fruck & Gaug 2015).

For the long-term analysis, we derived the night-wise γ-ray flux with energies above 300 GeV. The energy threshold was set to 300 GeV for variability studies of the integral flux to reduce its relative error. Next we fitted the observed spectrum at every night in the range 150 GeV–1 TeV with a LogP function (compatible with most of the data) to study the relation between the spectral hardness and the flux amplitude. The fitting function was folded by the energy dispersion and a model of the γ-ray absorption by extragalactic background light (EBL). Additional details about the analysis are given in Appendix A.

Three major flares were observed on 13 June, 14 June, and 1 July 2016. We performed a detailed analysis of the data during these highest VHE flux states. In order to reconstruct their intrinsic source spectra, the observed spectra were unfolded by the energy dispersion and the EBL absorption. In addition, we tested whether each of the following four functions can describe the intrinsic spectra: a simple PL, a PL with an exponential cutoff, a LogP, and a LogP with an exponential cutoff. The explicit formulae of these functions are defined in Appendix B. In order to determine the peak energy of the second SED component that is constrained by the MAGIC observations, we defined two additional functions having the forms given in Eqs. (B.5) and (B.6). For these functions, the local spectral index at E_peak was set to −2 in the expressions for $\frac{\mathrm{d}F}{\mathrm{d}E}$, thus enabling us to measure the peak location, where $E^2\times \frac{\mathrm{d}F}{\mathrm{d}E}$ tends to become flat. The fitting energy range is from 100 GeV to 9 TeV. For the spectral analyses on the three nights, we restricted the time window in order to avoid relatively strong moonlight and to get precise spectral measurements at the low-energy threshold. In addition to the spectral analyses, the intra-night variability was inspected. We produced light curves above 300 GeV with a fixed time-binning of 10 min for these nights and evaluated a characteristic timescale of the flux variabilities.

2.2. Fermi large area telescope

The LAT is one of two instruments on board the Fermi Gamma-ray Space Telescope (Atwood et al. 2009; Ackermann et al. 2012). The LAT covers an energy range from a few tens of MeV to more than 300 GeV. It covered the whole sky every three hours during its standard survey mode. Its point spread function (PSF) is about 0$\stackrel{°}{.}$8 with 68% containment at 1 GeV³. In order to suppress contamination of gamma rays from Galactic diffuse emission and nearby sources to 1ES 1959+650, we set the analysis range from 300 MeV to 300 GeV and used only three-fourths of the P8R2_SOURCE_V6 data with better PSF, namely the event types PSF1, PSF2, and PSF3 (Atwood et al. 2013).

We performed the standard binned likelihood analysis of the data from 28 April (MJD 57506.0) to 24 November (MJD 57716.0) 2016, binning by 3 days. In addition, we focused on the data for 1.5 days from 12 June 21:00 to 14 June 9:00 in 2016, which are quasi-simultaneous with the MAGIC data on 13 and 14 June 2016 and the central time roughly coincides with that of the MAGIC analysis period composed of those two days. In addition to the version 11-05-03 of the Fermi standard ScienceTools⁴, we used the version 0.15.1 of the Fermipy python package (Wood et al. 2017). The region of interest (RoI) for the likelihood fitting is taken to be a square of width 20°. The likelihood model includes the sources within a square region whose width is 40°. The included sources were taken from the LAT 4-year Point Source Catalogue (3FGL; Acero et al. 2015). In addition, an FSRQ CGRaBS J1933+6540 (also known as TXS 1933+655), which was detected by the LAT after the release of the 3FGL catalogue (Cheung 2018), the Galactic diffuse component and the isotropic background were also included in the model. The models for the Galactic diffuse and the isotropic component were given by the files of gll_iem_v06.fits and iso_P8R2_SOURCE_V6_PSF{1-3}_v06.txt, respectively⁵. We modelled CGRaBS J1933+6540 using a power-law spectrum and point-like spatial distribution at the position catalogued by Beasley et al. (2002). The background components with the detection significance lower than 1σ in each analysis period were removed from the model. The normalisation of the components with the significance > 3σ and all spectral parameters of the components with the significance > 4σ were kept free to vary for the fitting. The other parameters were fixed to the values in the 3FGL catalogue. The spectral index of 1ES 1959+650 was fixed to 1.88 ± 0.02 in a case when the significance was lower than 4σ. For the binned likelihood fit, the energy dispersion correction was enabled.

2.3. Neil Gehrels Swift Observatory

We used the publicly available data of two instruments onboard Neil Gehrels Swift Observatory (hereafter Swift), namely the X-ray telescope (XRT) and the UV/Optical Telescope (UVOT) in the LAT analysis period from 28 April to 22 November 2016, covering our MAGIC analysis period.

2.3.1. X-ray Telescope

The XRT is sensitive in the energy range of 0.2–10 keV (Burrows et al. 2004). The observations were performed in Window Timing mode and the exposure time distributes between 250 s and 1 981 s. 1ES 1959+650 was observed with the XRT 80 times during the period. The data were reduced with the version 0.13.3 of the standard software xrtpipeline. For the calibration files, the version 20160609 of the XRT CALDB were used. We extracted the counts within 47.1462 arcsec from the source. We used the version 12.9.1 of XSPEC for high-level analysis. We re-binned the data to have at least 25 counts per bin and fitted the spectra only above 0.5 keV with a PL and a LogP function. For the fitting, the equivalent hydrogen column density is fixed at 10²¹ cm⁻². For each observation of the XRT, we performed a fit of the spectrum with XSPEC. As the long term analysis, for each observation we produced the energy-flux light curve in the energy range 0.5 to 5 keV. The results of the fit with the LogP were also used to trace the relation between the spectral hardness and the flux. The details of the analysis are described in Appendix A. For the spectral analysis at the major flares, we found observations from 02:44 to 04:00, 13 June 2016 (ID: 35025243) and from 02:16 to 03:20, 14 June 2016 (ID: 35025245). These X-ray times are subsets of the MAGIC observation time. The exposure time is 972 s and 865 s, respectively. We derived the differential energy flux from 0.6 keV to 7.5 keV.

2.3.2. UV/Optical telescope

The UVOT has six filters that have a narrow effective waveband ranging from 170 nm to 600 nm (Roming et al. 2005; Poole et al. 2007). We measured the energy flux in an aperture with a radius of 5″ for one of the filters, fixing the Galactic extinction as E(B − V) = 0.1478 (Schlafly & Finkbeiner 2011). As the background region, we took an annulus from 27.5″ to 35″. In order to produce a long-term light curve, we used the data of the filter W1 centred at 260 nm; we did so because the source was observed most frequently with W1 among the filters of the UVOT, 81 times from 28 April to 21 November 2016. For the simultaneous multi-waveband analysis on 13 and 14 June 2016, the data of the filters [W1, W2: centred at 192.8 nm] and [W1, M1: centred at 224.6 nm, W2] were available, respectively.

3. Results

3.1. Long-term light curves in the very-high-energy γ-ray and other wavebands

The long-term flux light curve of γ rays with energies above 300 GeV from 29 April to 21 November 2016 is displayed on the top panel of Fig. 1. It exhibits a large flux variability of more than one decade. Such an erratic trend in the light curves is a common feature of HBLs. On 13 June, 14 June, and 1 July 2016, the flux above 300 GeV reached ∼3 Crab Units (C.U.)⁶ (13 and 14 June are treated as separate flares in our work mainly due to difference in their finer-scale temporal variability; see Sect. 3.5). This is the highest flux level observed from this source since 2002. On the other hand, the lowest flux level is ∼0.2 C.U., which is comparable to the value in the past low states, as mentioned in Sect. 1. There have been several other smaller flares with varying rise and fall times (e.g. two flares with flux level ∼2 C.U. around MJD 57550 and 57570). The light curves obtained with LAT, XRT and UVOT are plotted in the second, third, and last panel of Fig. 1, respectively. The flux in the UV and HE bands changes less than the one observed in the X-ray and VHE bands. The relation among the flux variability in VHE, HE, and X-rays looks complicated with varying rising and falling trends amongst different wavebands. The flux in the UV band exhibits an overall increase throughout the observation period.

Fig. 1. Long-term light curves of 1ES 1959+650 in 2016 with four instruments. From top to bottom: VHE gamma-ray flux (>300 GeV) from MAGIC, HE gamma-ray flux (0.3–300 GeV) measured with Fermi-LAT, X-ray energy flux (0.5–5.0 keV) from Swift-XRT, and UV energy flux (W1 filter, ∼260 nm) from Swift-UVOT. The flux in Crab Units is indicated with an additional y-axis in the top panel.

3.2. X-ray to very-high-energy γ-ray long-term correlation

X-rays and VHE γ rays are the most variable energy bands of 1ES 1959+650, as shown in Fig. 1. It is anticipated that the low-energy and high-energy peaks of the SED are located in the X-ray and VHE γ-ray bands, respectively, and the transition between different flux states of the source mostly affect the spectra in these two bands. In order to study the correlation between these two energy bands, the method of discrete correlation function (DCF; Edelson & Krolik 1988) was used. Figure 2 shows a plot of the correlation coefficients as a function of time lag in the range [−100, 100] days between X rays and VHE γ rays. A 2.5-day time-binning lag was used to keep sufficient statistics. An overall good correlation was found in the long-term with correlation coefficient of 0.76 ± 0.1 and zero time lag, as can be seen from Fig. 2.

Fig. 2. Discrete correlation function as a function of time lag for the VHE γ-ray and X-ray light curves of the 1ES 1959+650 2016 multi-wavelength monitoring campaign. In the long term, the VHE flux shows a correlation (DCF ∼ 0.76 ± 0.1) with the X-ray flux with zero time lag.

It should be noted that the long-term X-ray to VHE γ-ray correlation does not necessarily apply to the extreme flaring episodes of the source (e.g. during the 2002 MWL campaign, a general correlation was found between these two bands, except during the orphan flaring behaviour as reported in Krawczynski et al. 2004). For the mid-June 2016 high states of 1ES 1959+650, it is difficult to estimate the degree of correlation between X rays and VHE γ rays due to the limited MWL statistics for the short duration flares. Hence the correlation information is unconstrained for the modelling of the flares. Our calculations can infer that the derived correlation at least holds for the long-term behaviour of the source during 2016, barring its exceptional activities.

In general, the X-ray to γ-ray cross-correlation can be quite complex with different trends between the higher and lower energy bands (Ahnen et al. 2018) or between different observation epochs. A precise cross-correlation study requires a dense long-term multi-band dataset that can provide sufficient statistics for binning the data into multiple observation periods under varying source conditions and finer energy intervals. This will be followed up in a future paper with a denser and longer multi-wavelength campaign.

3.3. Spectral index vs. flux correlation in the very-high-energy γ rays and X-rays

In order to further investigate the long-term source behaviour in the X-ray and VHE γ-ray bands, we studied variations of the spectral indices as a function of the source fluxes in the above two energy bands. Figure 3 shows the spectral index as a function of flux in the VHE and X-ray bands for individual daily measurements with MAGIC and single observations with XRT, respectively. The value of α in the LogP function given in Eq. (B.3) was used as a measure of the spectral index. It corresponds to the value of the energy-dependent PL index at the normalisation energy E₀, which is fixed at 1 keV and 300 GeV for the XRT and the MAGIC bands respectively. The integrated photon flux above E₀ = 300 GeV and the differential energy flux at E₀ = 1 keV were used for the VHE and X-ray observations, respectively.

Fig. 3. Correlation between α of LogP fitting to the MAGIC spectrum for each night (left panel) and the XRT spectrum of each observation (right). The normalisation energy is fixed at 300 GeV for MAGIC and at 1 keV for XRT. The MAGIC spectra are deabsorbed with the model of Franceschini et al. (2008). More details can be found in the text.

In order to quantify the correlation between the spectral index and the flux, we used the weighted Pearson correlation coefficient (see Appendix C). The X-ray index variation is clearly not compatible with a constant function and shows the typical harder-when-brighter behaviour throughout 2016, confirming the past trends (e.g. Krawczynski et al. 2004). The weighted Pearson coefficient between the index and the flux in the X-ray band is $r=0.{75}_{-0.05}^{+0.05}$. In the VHE γ-ray band, the correlation coefficient $r=0.{64}_{-0.08}^{+0.09}$, suggests a harder-when-brighter behaviour in the VHE band during 2016. A richer dataset will be necessary to reinforce our claim regarding the strength of correlation.

3.4. Spectra during the highest flux nights

3.4.1. Very-high-energy γ ray

A spectrum at the VHE γ rays is especially important for this study because the high-energy peak of the SED is located in this energy range as explained in the following text, and therefore it is used to constrain the spectrum of the emitting particles. The SEDs in the VHE band on 13 June, 14 June, and 1 July 2016 are plotted in Fig. 4. They are unfolded with the instrument response function of MAGIC. The spectra are quite flat and extend beyond a few TeV. We fit the SEDs during the three nights with four functions (given by Eqs. (B.1)–(B.4); see also Eqs. (B.5) and (B.6)) and the results are documented in Table 1. In all cases, the EBL-corrected VHE spectra are more compatible with a curvature in the spectra rather than a simple PL (Eq. (B.1)), but no clear preference among Eqs. (B.2)–(B.4) were found on 13 and 14 June 2016. For 1 July, LogP with cutoff-type of spectrum is preferred over a pure LogP spectrum with a significance of ∼3σ. However, the curvature parameter β for the LogP with cutoff spectrum is consistent with zero and the best-fit function is essentially the same as a PL with cutoff-type of spectrum. In the case of either of these curved functions, the power-law index is harder than 2 around 300 GeV. The high flux and our intensive observations enabled us to determine the cutoff energy of the PL with a cutoff-type function and the peak energy of the LogP function to ∼10% and ∼20% statistical uncertainty, respectively. The SEDs peak at ∼0.4–0.7 TeV, and they have a cutoff above a few TeV when fitted with Eq. (B.2). The cutoff energy and the peak energy on 13 June 2016 are higher than those on the other two nights. These peak energies are similar to the peak values of Mrk 501 SED observed in April 1997 and in June and July 2005 (Djannati-Atai et al. 1999; Albert et al. 2007). Compared to the peak energy of Mrk 501 SEDs observed on some nights in 2012 (Ahnen et al. 2018), those peaks are a few times lower.

Fig. 4. VHE SEDs during the nights of highest flux: 13 June, 14 June, and 1 July 2016 (from top to bottom). The black circle and red square markers represent the observed and the EBL-deabsorbed spectra, respectively. These spectra have been unfolded with the instrument response function of MAGIC. The absorption by the EBL has been corrected via the model of Franceschini et al. (2008).

3.4.2. Other wavebands

Here we report the results of the SED analysis for the other bands from 13 and 14 June 2016, and 1 July 2016 excluded due to the lack of simultaneous X-ray and optical data. The spectral points are reported in the broadband SEDs plotted in Figs. 6–8.

The fitting result of the LAT spectrum for the 1.5 days is documented in Table 2. The power-law index is 1.56 ± 0.20. The index is marginally harder than the values reported in the 3FGL and the 4FGL catalogue (1.88 ± 0.02 and 1.82 ± 0.01, respectively) by less than 2σ. It should be noted that the analysis energy ranges are not identical to that of our analysis. We cannot find a significant curvature or break in the spectrum because of the small photon statistics. The parameters of the XRT datasets, which were simultaneous with the MAGIC observations, are listed in Table 2. The spectra are fitted by PL and a LogP function does not improve the goodness of the fit. The power-law index is 1.81 ± 0.01 on 13 June and 1.82 ± 0.01 on 14 June 2016. It is clearly harder than when the source was in a low state in 2006 as Tagliaferri et al. (2008) reported a PL index 2.197 ± 0.001 with the data of Suzaku-XIS for 0.7–10 keV. Our results can be compared to the SED of Tagliaferri et al. (2008) with the XIS and Suzaku-HXD/PIN also in Figs. 6–8.

3.5. Intra-night variability

The investigation of intra-night variability in the VHE band is essential not only to constrain the size of the emission region, but it also plays an important role in replicating the physical conditions inside the source leading to the origin of the second SED peak. The observed VHE γ-ray flux exhibited fast variations for some nights in 2016, particularly for the nights with the highest VHE flux levels. We analysed the light curves with a fixed time-binning of 10 min and found that the light curves above 300 GeV for 13 June and 1 July 2016 show significant intra-night variability over short timescales, as shown in Fig. 5. The flux level above 300 GeV of our standard candle, the Crab Nebula is also shown for comparison purposes with a red dotted line. No significant intra-night variability was observed on 14 June 2016. A common method to quantify the mean variability of the source is given by the fractional variability amplitude (Vaughan et al. 2003). For a set of N flux points x_i with corresponding errors σ_i, err, having mean flux x_mean and mean squared error ${\sigma }_{\text{mean},\text{err}}^2$, the fractional variability is defined as

$$\begin{array}{c}\hfill F_{\mathrm{var}}=\sqrt{\frac{S^2-{\sigma }_{\mathrm{mean},\mathrm{err}}^2}{x_{\mathrm{mean}}^2}},\end{array}$$(1)

Fig. 5. Fast intra-night variability observed in the VHE band on 13 June (left panel) and 1 July (right panel) 2016. The light curves above 300 GeV are constructed with a fixed time-binning of 10 min. Green dashed curve: fit with the function given in Eq. (3); solid red curve: fit with the function given in Eq. (4). See the figure legend for the rise and decay times inferred from the fit. The blue dashed line represents the steady flux of the Crab Nebula above 300 GeV, shown for comparison purposes.

where S² denotes the sample variance. The error in F_var is calculated following Eq. (B2) in Vaughan et al. (2003). The fractional variability amplitude for 13 June, 14 June, and 1 July 2016 are 0.20 ± 0.02, 0.06 ± 0.05, and 0.16 ± 0.02, respectively. Another approach to give a quantitative measure of the variability is to calculate the power of variability from power spectral density (PSD; Vaughan et al. 2003). The analysis of our data points shows that the power-law index obtained from a fit to the PSD has the hardest value for 13 June 2016, followed by 1 July; 14 June 2016 has the softest index of the three nights, which is a result similar to the one obtained from the fractional variability amplitude. However, due to a limited number of data points, determining the slope of the PSD is not very meaningful and the fractional variability amplitude gives a more reliable measure of the flux variations.

An estimate of the fastest variability timescale can be obtained from the doubling time which is defined using the formulae from Zhang et al. (1999),

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

where F_i, F_i + 1 and t_i, t_i + 1 denote the fluxes and corresponding observation times for two consecutive data points in the light curve. The errors of the doubling timescale are propagated through the errors in the flux measurement. For the night of 13 June 2016, the pair-wise shortest variability timescale was found between the 8th and the 9th data points (t_var = 36 ± 14 min). The same quantity calculated for the night of 1 July was found to have the minimum value between the 2nd and the 3rd data points with flux doubling time, t_var = 36 ± 15 min.

The rise and decay times of the individual substructures in the light curves can be obtained by fitting the peaks with an exponential or sum of two exponential functions represented by the formulae

$$\begin{array}{cc}\hfill F(t)& =A_0{e}^{-|t-t_0|/t_r},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill F(t)& =A_0/(e^{\frac{t_0-t}{t_r}}+e^{\frac{t-t_0}{t_f}}),\hfill \end{array}$$(4)

where A₀ is defined as the flux and two times the flux at t₀ for Eqs. (3) and (4), respectively; t_r, t_f are respectively the rise and decay times of the flare, which are left as free parameters; and the flux doubling time in this formalism is defined as t_rise/fall = t_r/f × ln(2). For 13 June and 1 July 2016, the results of the double-exponential fit (solid red curves in Fig. 5) and the single-exponential fit (green dashed curve in Fig. 5) are summarised in Table 3. The doubling timescales obtained from the fitting are comparable to the results from the Zhang et al. (1999) formulation, the former being slightly biased by the choice of the fitting range. For the theoretical discussions in the next section, the timescales obtained from the Zhang et al. (1999) formulation are used.

4. Discussions

4.1. Size of the emission region

The variability of blazars can act as a powerful probe to characterise the emission region and to investigate the undergoing physical processes inside the source. The emitting region is assumed to be a spherical blob of radius R in our broadband SED models (discussed below). The observed variability timescale t_var can be used to derive an upper limit (UL) to the size of the radiating blob in the co-moving frame of the jet as (Tavecchio et al. 2010a)

$$\begin{array}{c}\hfill R\le \frac{c{t}_{\mathrm{var}}\delta }{1+z},\end{array}$$(5)

where δ represents the Doppler factor and z represents the redshift of the source. The fastest variability timescale of the source in the VHE γ-ray band observed by MAGIC is used in our calculation for this purpose. However, we note that the spectra at the VHE γ rays used in the broadband SED modelling described in the next section represent an average emission state for the entire night, and thus are not a true representative of the finer scale variability observed in the light curve. Under the assumption of t_var ∼ 35 min as derived in the previous section, the upper limit to R with δ = 20−60 can be given in the range 10¹⁵ − 3 × 10¹⁵ cm.

An additional constraint to the value of R can be provided from the condition that the radius of the emission region should always be greater than the gyro-radius of the highest energy protons, valid only for the hadronic-dominated solutions discussed below. This condition is always respected in our modelling and is given by the formula (Böttcher et al. 2013)

$$\begin{array}{c}\hfill B\text{(G)}\ge 30\frac{E_{p,\max }}{{10}^{19}(\mathrm{eV})}\frac{{10}^{15}}{R(\mathrm{cm})},\end{array}$$(6)

where E_p, max represents the maximum energy of the protons and B represents the strength of the magnetic field inside the emission volume.

4.2. Multi-wavelength spectral characteristics

We assembled quasi-simultaneous MWL data during the VHE outbursts over optical–VHE γ-ray wavebands, as shown in Figs. 6–8, in order to investigate the broadband spectral behaviour of the source. For comparison, we also plot the MWL spectra from 2006 in Fig. 6 during a low VHE flux state measured by MAGIC (grey points; Tagliaferri et al. 2008). A clear shift of both SED peaks towards higher energies (in the X-ray and VHE γ-ray domain) was observed during the flares with respect to the historical data. In addition, the spectral indices are harder during the flares in both these bands. This behaviour is reported as being typical for HBLs in Tagliaferri et al. (2008) and has been mentioned by many other authors in the past (e.g. Tagliaferri et al. 2003; Pian et al. 2014). The lower energy SED peak during the flares, although not well constrained, lies above a few 10¹⁸ Hz. This indicates that the synchrotron peak frequency shifted towards the extreme-HBL (EHBL; Costamante et al. 2001) regime during the major flaring episodes of the source. A similar behaviour was also previously observed for Mrk 501 (Ahnen et al. 2018). The corresponding EBL-corrected γ-ray luminosity is above 10⁴⁵ erg s⁻¹ which is slightly higher than the usual expectation from the luminosity-peak-frequency anti-correlation behaviour predicted in the so-called blazar sequence (Ghisellini et al. 2017). However, we note that the blazar sequence is constructed using average SEDs that also contain data collected during quiescent states. Moreover, the ratio between the peak luminosities of the higher and lower energy SED components (defined by the Compton dominance parameter in the blazar sequence) changed by a factor of ∼4 between the historical data and the flaring γ-ray states in 2016 corresponding to ∼1 order of magnitude change in the γ-ray luminosity. In the blazar sequence the Compton dominance changes by ∼1 order of magnitude for a ∼4 order of magnitude change in the luminosity.

Fig. 6. One-zone SSC models applied to 13 June (left panel) and 14 June (right panel) 2016. The symbols corresponding to the data from different instruments are given in the legend. The historical data are taken from Tagliaferri et al. (2008). The black (solid), brown (dot-dashed), and blue (dashed) curves represent the summed emission component in increasing order of Doppler factor δ. We found a satisfactory explanation of the MWL data with high values of δ >  45 for 13 June 2016. The data from 14 June 2016 do not strictly require such high values and can be modelled with moderate values of δ >  30. For more details see the discussion in Sect. 4.3.1 and the parameters in Table 4.

4.3. Broadband emission modelling

We modelled the broadband SEDs during the flares (13 and 14 June 2016; 1 July discarded due to lack of simultaneous X-ray and optical data) using three different theoretical frameworks considering one-zone leptonic, hadronic, and lepto-hadronic models. For the modelling, a modified version of the code described in Ansoldi et al. (2018) was used. The emitting region is assumed to be a spherical blob of radius R filled with a tangled magnetic field of strength B, moving down the jet with bulk Lorentz factor Γ_bulk. The viewing angle of the radiated photons in the jet frame are at an angle θ with respect to an observer on Earth. The radiative output is calculated in the jet co-moving frame and then transformed to the observer frame via the Doppler factor δ = [Γ_bulk(1 − βcos(θ))]⁻¹.

4.3.1. Leptonic model

First, we investigated a one-zone SSC model for the SEDs by assuming a stationary population of primary electrons within the emitting region. The primary electron distribution is assumed to follow a broken power-law described by two slopes, n₁ and n₂; a break Lorentz factor γ_e, brk; and a minimum and maximum Lorentz factor γ_e, min and γ_e, max, respectively. The break energy is calculated by balancing the synchrotron cooling and the electron escape timescale using the following condition

$$\begin{array}{c}\hfill t_{e,\mathrm{sync}}=\frac{7.75\times {10}^8}{B^2\times {\gamma }_{e,\mathrm{brk}}}=\frac{R}{c},\end{array}$$(7)

where t_e, sync represents the synchrotron cooling timescale of the electrons. From our modelling, values of n₁ between 2.25–2.3 were found to provide a satisfactory description of the Swift-XRT and UVOT data. Under the assumption that the break in the electron spectrum is induced by radiative cooling, the second index is constrained as n₂ = n₁ + 1. Although the peak of the first SED is not well defined, due to lack of simultaneous hard X-ray data, it helps to constrain the value of the magnetic field strength, which is then used to derive γ_e, brk for a given value of R (from Eq. (7)). The deabsorbed VHE spectrum is quite flat and extends up to several TeV, especially for the 13 June flare. However, due to the fast radiative cooling of the electrons and the Klein-Nishina effect, the inverse Compton component is generally suppressed and cannot easily explain the flat photon spectrum observed at TeV energies. To overcome this effect our model requires large values of Doppler factor and low magnetic field strength to generate the broadband spectra up to VHE in the simple SSC solutions. The results from the modelling are shown in Fig. 6. The MWL SED of 13 June 2016 can be satisfactorily explained with δ≥45–50, whereas that of 14 June requires smaller values of δ≥30, which mainly arises from differences in spectral hardness and/or cutoff in the VHE data measured by MAGIC for 13 and 14 June 2016. Smaller values of the Doppler factor are ruled out for the range of magnetic field strength considered in this work. The complete list of parameters for these models is reported in Table 4.

Next, we compare these results to those of three previous flares of 1ES 1959+650 with one-zone SSC models. Krawczynski et al. (2004) applied an SSC model to the VHE high state of 1ES 1959+650 observed in 2002. The authors averaged the VHE spectra during six nights with the flux greater than 1 C.U. above 2 TeV and estimated the averaged X-ray spectrum corresponding to it. The VHE flux (here ν >  3 × 10²⁶ Hz) and the X-ray flux (here ν ∼ 10¹⁸ Hz) in the averaged spectra are comparable to those in our data. The X-ray and the VHE γ-ray flux values in 2002 were lower by roughly 20–40% than those on 13 June 2016. The authors concluded that the estimated MWL SED is reproduced by δ = 20, R = 5.8 × 10¹⁵ cm while other parameters are comparable to ours. Our dataset covers a much wider energy range than the one used by Krawczynski et al. (2004), and the VHE spectra extending up to TeV energies is flatter than their model. The flat SED requires large values of the Doppler factors because such a spectrum is only reproduced by the SSC emission radiated from electrons with energy below γ_e, brk. Since γ_e, brk is constrained in the previous paragraph, δ must be large so that sub-TeV γ-rays are dominantly radiated by those electrons.

In Fig. 6, the data points taken from Tagliaferri et al. (2008) exhibit a high state in X-ray and a low state in the VHE γ ray between 24 and 29 May 2006. The flux ratio of these two energy bands differs from that of our data by a factor of ∼4. Tagliaferri et al. (2008) described the SED by a one-zone SSC model with δ = 18, R = 7.3 × 10¹⁵ cm and B = 0.25 G. The difference in the luminosity ratio, which is determined by L_SSC/L_sync = U_sync/U_B, arises due to the difference in Compton dominance. The increased Compton dominance in our modelling is due to a combination of 7–10 times smaller R compensated by 2–3 times higher δ compared to the modelling of Tagliaferri et al. (2008).

Aliu et al. (2014) reported a VHE flare on 20 May 2012 (MJD 56067) without a simultaneously observed high X-ray state. The UV and X-ray spectra observed during the high and low VHE states are very similar, while at the same time being significantly different from those of the flares in June 2016. The authors applied a time-independent SSC model to the high state and an averaged low state. According to these models, the synchrotron peak is located at ∼10^16.5 Hz, which is more than two orders of magnitude lower than our models. Their model peak is produced by the minimum electron Lorentz factor γ_e, min of an order of 10⁶. This is much higher than that of our models, γ_e, min = 3–7 × 10². Such high γ_e, min values were also suggested by Patel et al. (2018) in the context of two-zone SSC models for several high-state periods and a low-state period in 2016.

The high-energy SED peak especially on 13 June 2016 lies close to the regime of the so-called hard-TeV BL Lac objects (EHBLs with high-energy peak above ∼2 TeV). Generally, their models require high electron spectral break energy γ_e, brk ∼ 10⁶ and large Doppler factors δ = 20–60 (Costamante et al. 2018), similar to ours. However, the magnetic field strength is extremely weak, at mG level or even lower (see e.g. Kaufmann et al. 2011). In addition to typical EHBLs, some HBLs exhibit serendipitous EHBL-like nature temporarily. The most deeply investigated amongst them, Mrk 501 had shown harder VHE spectra or a shift of the second SED peak up to ∼1 TeV during flaring activities (Albert et al. 2007; Aliu et al. 2016). The one-zone SED modelling of such states given in Albert et al. (2007) implies [B, δ] = [0.05 G, 50], [0.23 G, 25] and R ∼ 10¹⁵ cm, roughly compatible with our models. A temporary transition of Mrk 501 towards the EHBL regime was also observed in 2012 (Ahnen et al. 2018). For the strongest VHE flare on 9 June 2012, a two-zone SSC model applied by Ahnen et al. (2018) yielded B = 6.8 × 10⁻² G, R = 3.3 × 10¹⁵ cm, and γ_e, min = 2 × 10³. Consequently, our SSC model parameters are close to the range predicted for temporal or standard EHBLs. This might further corroborate the transition of 1ES 1959+650 towards an EHBL-like state during the mid-June 2016 VHE outbursts.

4.3.2. Hadronic model

We also investigated an alternative scenario, in which the high-energy component of the SED is associated with relativistic protons additionally injected into the emission region along with the primary leptons. The proton distribution is described with a power law with an exponential cutoff function having proton spectral index n_p and exponential cutoff Lorentz factor γ_p, max. We fixed the minimum proton Lorentz factor to γ_p, min = 1 in order to get a conservative estimate of the proton luminosity budget. γ_p, max in the co-moving frame is determined from the condition t_acc = minimum[t_esc,  t_psync,  t_p − γ], where t_acc = 10η_accE_p/eBc (Cerruti et al. 2015) denotes the acceleration timescale and t_esc, t_psync, and t_p − γ denote the particle escape, proton-synchrotron, and photo-meson cooling timescales, respectively (see Fig. 9). The particle escape is parametrised by an efficiency factor η_esc such that t_esc = η_escR/c (Aliu et al. 2014).

In the hadronic scenario that we investigated, direct synchrotron radiation by the highest energy relativistic protons (a few EeV in the co-moving frame) can satisfactorily reproduce the second SED peak located at few hundreds of GeV (i.e. in the VHE regime constrained by the MAGIC observations). The lower energy peak is still associated with synchrotron radiation by the primary leptons. The hadronic solutions are shown in Fig. 7. The photo-meson cascade component arises due to emission by secondary leptons that are generated when a high-energy proton interacts with the low-energy synchrotron photon field. It gives a sub-dominant contribution to the overall SED in the chosen parameter space. In the proton-synchrotron solutions, the protons have to be accelerated up to few EeV energies, which can be achieved if the source possesses very high acceleration efficiency (η_acc = 1) under magnetically dense environments (see the timescale plots in Fig. 9 (left panel). Large magnetic fields of the order of 100 G are adopted in our purely hadronic solutions in order to overcome the slow cooling timescale of protons, which is generally insufficient to explain the variability timescales of less than an hour observed from this source during 2016. Under these conditions, the protons can cool down with timescale t_psync ∼ 2.5 × 10⁴ s, shorter than the co-moving frame variability timescales (Δt_jet = δΔt_var) exhibited by the source in the VHE band. The requirement of somewhat larger values of the magnetic field is typical for proton-synchrotron models. Moreover, for our choice of R, a few times 10¹⁴ cm and assuming a jet-opening angle close to 1 degree, the distance from the central core d becomes a few times 10¹⁶ cm, where B ∼100 G can be expected (e.g. Barkov et al. 2012). No spectral break due to cooling is assumed in the propagated spectrum of the protons in our simple formalism since they remain un-cooled before escaping (i.e. t_psync and t_esc are competing processes having almost equal values for the highest energy protons). In this high-B domain, the electrons are in the fast-cooling regime and parametrised by a simple power-law distribution. The complete list of model parameters for 13 and 14 June 2016 can be found in Table 4. The Doppler factor required for the hadronic solutions is considerably smaller (δ ∼ 25) compared to that required for the purely leptonic models, especially for 13 June 2016 and represents fairly typical values. In the domain of such typical values of δ, magnetic field strengths of less than 100 G (which also implies lower values of the maximum proton energy) are insufficient to explain the flat TeV spectra of 1ES 1959+650 in purely hadronic solutions. The difference between the VHE spectra from 13 and 14 June 2016 can be mainly attributed to small differences in the values of γ_p, max in our hadronic solutions (14 June requires slightly smaller values of γ_p, max than 13 June: ∼5 × 10⁹ and ∼7 × 10⁹, respectively). The total jet power is evaluated as

$$\begin{array}{c}\hfill L_j=\pi R^{2}c{\mathrm{\Gamma }}_{\mathrm{bulk}}^2(u_p+u_B+u_e),\end{array}$$(8)

Fig. 7. One-zone hadronic models applied to 13 June (left panel) and 14 June (right panel) 2016. The symbols corresponding to the data from different instruments are given in the legend. Solid black line: summed; dashed blue line: electron-synchrotron; dotted green line: SSC; dot-dashed sea green line: proton-synchrotron; dot-dot-dashed orange line: p − γ cascade. The higher energy peak in the SED is dominated by synchrotron radiation by relativistic protons, which can be achieved with B ∼ 100 G and Ep, max >  1018 eV and jet power Lj ∼ 1046 erg s−1 (∼LEdd). For more details see the discussion in Sect. 4.3.2 and the parameters in Table 4.

where u_p, u_B, and u_e represent the energy densities carried by the protons, magnetic field, and electrons, respectively. In the proton-synchrotron-dominated solutions the required jet power amounts to L_j ∼ 10⁴⁶ erg s⁻¹, comparable to the Eddington luminosity (L_Edd ∼ 10⁴⁶ erg s⁻¹) of the source (assuming M_BH = 10⁸ M_⊙; Falomo et al. 2002).

Aharonian (2000) applied a proton-synchrotron scenario to a TeV spectrum of Mrk 501 during extraordinary flares measured in 1997. The EBL-corrected SED peaks at 1–2 TeV with an exponential-like cutoff around 6 TeV is similar to the spectrum of 1ES 1959+650 on 13 June 2016. The author argued that the characteristics of the TeV flares were explained by synchrotron emission from ultra-relativistic protons with γ_p.max ≥ 10¹⁰, strong magnetic field B = 30–100 G, and the Doppler factor δ = 10–30. These values are comparable to the parameters in our hadronic solution. Aharonian (2000) also noted that the spectral shape is stable regardless of any possible changes in R and B, provided δ and η_acc remain unchanged. This agrees with our findings that the spectral shape in the VHE γ-ray band of the two sources is similar to each other during the flares.

4.3.3. Lepto-hadronic model and implications for neutrino emission

In general, the proton-synchrotron models predict neutrino fluxes below the sensitivity of the current generation of neutrino telescopes. In order to further investigate the potential of neutrino emission we also studied a lepto-hadronic model for both 13 June and 14 June 2016. The high-energy SED peak is different between the two nights mainly at the VHE γ-ray band (13 June 2016 has a slightly harder spectrum and higher VHE flux). The photo-meson cascade component (the main hadronic component in our lepto-hadronic model) can take into account the differences between the VHE spectra of 13 June 2016 and 14 June 2016 with a slight difference in the particle energetics between the two nights. Our conclusions about the level of neutrino emission from the source remain the same for both nights taking into account such small differences in the γ-ray spectra. Hence, we only take the night of 13 June 2016 as a reference in our paper.

We assume an additional proton population with a power-law spectrum characterised by the same functional form as described in the hadronic modelling subsection, along with the relativistic electrons inside the emission region. In the lepto-hadronic solutions, the second SED peak is comprised of contributions from both the SSC component and the p − γ cascade component as shown in Fig. 8 (left panel). In this case, the required maximum proton energy is governed by the particle escape timescale, as can be seen from the timescale plot in Fig. 9 (right panel) (γ_p, max ∼ 6 × 10⁷). The values of the other model parameters are given in Table 4.

Fig. 8. One-zone lepto-hadronic models (left panel) and the predicted neutrino flux (right panel) for 13 June 2016. The definition of symbols and lines in the SED model in the left panel is the same as Fig. 7. The higher energy peak in the SED in this case, is a combination of the SSC and photo-meson cascade component which can be achieved with B ∼ 1 G, Ep, max >  1016 eV at the cost of high jet power Lj >  1048 erg s−1 (≫LEdd). The meaning of the different curves in the neutrino spectra (right panel) is mentioned in the legend. The IceCube sensitivity curve is taken from IceCube Collaboration (2019) corresponding to declination 60°. Neutrino spectra predicted in the proton-synchrotron solutions of Fig. 7 peak at very high energies and provide low neutrino flux in the range 0.1–100 PeV. The neutrino peak is shifted to lower energies in the lepto-hadronic solutions providing slightly higher flux at the cost of very high values of the jet luminosity. For more details see the discussion in Sect. 4.3.3 and the parameters in Table 4.

Fig. 9. Comparison between acceleration timescale (solid black line with acceleration efficiency ηacc  =  1) and different cooling timescales (tesc, tpsync, tpγ) for the hadronic (left panel) and lepto-hadronic (right panel) scenarios of 13 June 2016. Also shown is the energy at which the proton gyro-radius becomes equal to the radius of the emission region (dot-dashed red line).

The peak of the neutrino spectra is mainly governed by the maximum proton energy. In our proton-synchrotron solutions, due to the requirement of high values of the maximum proton energy to explain the electromagnetic SED, the inferred neutrino spectrum peaks at energies above a few EeV in the observer frame and the flux at 0.1–100 PeV is quite low. In the case of the lepto-hadronic solution due to the requirement of much lower values of γ_p, max the neutrino spectra peak about two orders of magnitude lower in energy. The inferred individual and summed components of the neutrino spectra predicted from the lepto-hadronic solution are shown in Fig. 8 (right panel); for comparison, the neutrino spectrum from the proton-synchrotron model is also shown (brown dashed line). Additionally, the IceCube sensitivity curve, calculated for 8 years of operation and at the declination of 1ES 1959+650 (IceCube Collaboration 2019) is overlayed in the figure. We note that a direct comparison of our model-derived neutrino spectra to the IceCube sensitivity is difficult due to the variable nature of the 1ES 1959+650 electro-magnetic emission. The neutrino spectra, calculated for a short-lasting high emission state, hardly reaches the limit of the IceCube sensitivity. Therefore, we can expect that on average the neutrino emission from this source will be much lower. From the lepto-hadronic solution, the integrated neutrino flux in the range 600 GeV–100 TeV (i.e. the central 90% neutrino energy range for the declination of the source ∼65°, calculated from Fig. 1, bottom panel, in IceCube Collaboration 2019) is ∼5.5 × 10⁻¹³  TeV cm⁻² s⁻¹, which is comparable to the upper limit flux for 1ES 1959+650 obtained by IceCube (9.86 × 10⁻¹³  TeV cm⁻² s⁻¹ at 90% C.L.⁷). Moreover, the lepto-hadronic solutions require very high values of jet power (L_j >  10⁴⁸ erg s⁻¹) exceeding L_Edd by more than 2 orders of magnitude, and hence are energetically less favourable (but see Barkov et al. 2012). Although relaxing the condition on the minimum proton Lorentz factor γ_p, min = 1 can reduce the luminosity to some extent, it is still insufficient to achieve sub-Eddington values. Based on the conclusions from our one-zone electromagnetic emission modelling we infer that it is difficult to produce detectable neutrino emission during the 2016 flares of 1ES 1959+650.

5. Summary and conclusions

In this work, we have reported the spectral and temporal properties of the 1ES 1959+650 MWL emission in 2016 with a special emphasis on the major VHE γ-ray flares observed by MAGIC during this period. During the 2016 long-term MWL monitoring campaign, the X-ray flux was found to be correlated in general with the VHE γ-ray flux having a discrete correlation coefficient of 0.76 ± 0.1 with no lag. For the individual extreme flaring episodes of 13 June and 14 June 2016, a correlation behaviour could not be quantified due to a lack of sufficient number of data points for the short duration flares. Hence the correlation information was not used in the broadband SED modelling. In the long-term, the X-ray spectral index hardens with increasing flux level and a hint of similar behaviour is also visible in the VHE band. In our follow-up paper with more quasi-simultaneous multi-wavelength data, we will provide a detailed cross-correlation study between the X-ray and the VHE bands. The absence of long-term correlation between these two bands might imply two different emission regions (e.g. Patel et al. 2018) or two completely independent particle populations giving rise to the emission components in the different wavebands. The blazar has shown extreme flaring episodes in the VHE band especially on 13 June, 14 June, and 1 July 2016 with the highest flux > 300 GeV reaching 2.5–3 times the Crab Nebula flux above the same energy. This is the highest flux observed from this source since the orphan flares in 2002 (the flares presented in this paper are not orphan). The VHE spectra during the flares are quite flat extending to several TeVs. MAGIC plays a crucial role in constraining the location of the second peak in the broadband SED. The first SED peak is not well constrained due to the lack of simultaneous hard X-ray data and is treated as a somewhat free parameter in our emission modelling scenarios.

The nights of 13 June and 1 July 2016 showed fast intra-night variability in the VHE band on sub-hour timescales, which indicates the presence of small compact emission regions inside the jet. The VHE flux on the night of 14 June 2016 was fairly constant without any signs of variation over short timescales. Owing to the different temporal characteristics of the two flares on 13 June and 14 June 2016, we consider them independent in the modelling assuming that they arise from two different independent emission zones within the jet. If the same emission region were to be the cause for both the flares, then in time T ∼ 2 days, the blob would be travelling a distance z ∼ cTΓ²∼1–2 parsec (assuming Γ∼20–40). After travelling such a large distance, the blob would expand and lose energy adiabatically, and its magnetic field strength would decrease (Tagliaferri et al. 2008). As a consequence this would weaken the produced flux; however, this is not observed in the VHE band (13 June and 14 June have comparable flux levels).

We discussed the broadband spectral characteristics of the source in the framework of leptonic, hadronic and lepto-hadronic emission scenarios. In all cases we used the formulation of isotropic target photon fields for inverse Compton or photo-meson interactions that is provided by the synchrotron radiation of the primary electrons responsible for the first SED peak. To explain the broadband spectra up to TeV energies, the SSC models require large values of Doppler factor in general (δ >  30), which indicates highly relativistic small regions in the jet responsible for the γ-ray emission. The requirement of larger δ on 13 June 2016 compared to 14 June (see Table 4) can be mainly attributed to the difference in the spectral hardness and/or cutoff in the VHE band. We note, however, that our assumption of cooling break (ΔN = n₂ − n₁ = 1) comes from the simplest expectations of a break due to radiative cooling and in reality the acceleration and cooling processes can be more complex leading to a different spectral behaviour (see e.g. Tavecchio et al. 2010b). In the context of this model, the unusually high flux in VHE γ rays is considered to be produced by high Compton dominance related to the small emission region and the strong relativistic Doppler boosting compared with those of studies of different periods. In addition, some obtained parameters in the SSC model (such as B, δ, γ_e, min) are similar to those predicted for Mrk 501 during an EHBL-like state. This might imply the transition of 1ES 1959+650 to such a state during extreme flaring periods.

We have also investigated alternative solutions where the jet is composed of relativistic protons in addition to the accelerated electrons. In the first scenario, the so-called proton-synchrotron model requires high values of the magnetic field strength (of the order 100 G) and acceleration efficiency close to the theoretical maximum (η_acc ∼ 1) to explain the γ-ray observations. In this parameter regime, the electrons cool down very rapidly. A hard injection spectrum of the electrons (<2) is thus required to explain the X-ray observations. Such a spectrum can be generated by acceleration mechanisms such as stochastic acceleration (Virtanen & Vainio 2005) or magnetic reconnection (Cerutti et al. 2012; Sironi & Spitkovsky 2014).

The position of the second SED peak strongly depends on the maximum energy of the protons, which in our model is determined by a balance of the acceleration and cooling timescales (proton-synchrotron, escape, photo-meson). Compared to the SSC models, the proton-synchrotron solutions require smaller values of the Doppler factor (δ ∼ 25). We have also investigated mixed lepto-hadronic models where the high-energy SED peak is a combination of the SSC and proton-induced cascade emission. The required jet power for the proton-synchrotron solutions is comparable to the Eddington luminosity of the source (∼10⁴⁶ erg s⁻¹) and that for the lepto-hadronic solutions exceeds L_Edd by about 2 orders of magnitude. However, super-Eddington values of jet power in blazars have been predicted by various other authors (e.g. Barkov et al. 2012; Basumallick & Gupta 2017 and references therein). We also note that the jet power can be significantly reduced by assuming external photon fields inside the emission region as in the structured jet scenario discussed in Tavecchio et al. (2014; see also Righi et al. 2017).

The neutrino spectra predicted from the proton-synchrotron model peak at very high neutrino energies (i.e. > 10¹⁸ eV in the observer frame, which is a consequence of the requirement of the high maximum proton energy in such solutions). It provides low neutrino flux in the range 0.1–100 PeV. The neutrino flux in this range can be boosted by choosing a lower value of the maximum proton energy as shown in the lepto-hadronic solutions. However, such a scenario is also energetically less favoured due to the requirement of high values of jet power as discussed above. Our predicted neutrino spectra during the brightest 2016 flare do not significantly exceed the IceCube sensitivity limit (calculated using 8 years of IceCube active observing time) in all cases. The model-predicted integrated neutrino flux in the range 600 GeV–100 TeV (90% energy confidence interval) is comparable to the flux upper limit in the location of 1ES 1959+650 derived from 8 years of IceCube data. Our conclusions are in agreement with the non-detection of significant neutrino excess in the IceCube data analysis following the 2016 γ-ray flares (Kintscher et al. 2018).

In this work, a comparative study was done for different classes of SED models to demonstrate the multiple possibilities, which naturally leads to some degeneracy in the parameter space. Future multi-messenger and multi-wavelength observations can play a very crucial role to distinguish between the hadronic and leptonic scenarios and constrain the model parameters. For example, a multi-year multi-waveband monitoring campaign can help to follow the transition between high and low emission states. Such a long-term dataset is of paramount importance in order to understand the nature of the emitting particles, follow the evolution of the model parameters, and characterise the undergoing physical conditions which might change rapidly with the changing state of the source. These studies will be followed up in our future publication with a long-term monitoring campaign.

Acknowledgments

We would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF and MPG, the Italian INFN and INAF, the Swiss National Fund SNF, the ERDF under the Spanish MINECO (FPA2015-69818-P, FPA2012-36668, FPA2015-68378-P, FPA2015-69210-C6-2-R, FPA2015-69210-C6-4-R, FPA2015-69210-C6-6-R, AYA2015-71042-P, AYA2016-76012-C3-1-P, ESP2015-71662-C2-2-P, FPA2017-90566-REDC), the Indian Department of Atomic Energy, the Japanese JSPS and MEXT and the Bulgarian Ministry of Education and Science, National RI Roadmap Project DO1-153/28.08.2018 is gratefully acknowledged. This work was also supported by the Spanish Centro de Excelencia “Severo Ochoa” SEV-2016-0588 and SEV-2015-0548, and Unidad de Excelencia “María de Maeztu” MDM-2014-0369, 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. We also appreciate a helpful framework which we used for the Fermi-LAT data analysis, Fermipy (Wood et al. 2017).

References

Appendix A: Analysis details of the spectral index vs. flux correlation

We describe the details of the analysis reported in Sect. 3.3 here. The spectral fit is done with a LogP function given in Eq. (B.3), which is simple and compatible with most of these spectra. To fit the MAGIC spectra, the LogP function was folded by two functions of the photon energy. One is the dispersion of the reconstructed energy from the true value, and the other is a correction factor for photon absorption due to the EBL (using the model of Franceschini et al. 2008). On the other hand, to reconstruct the intrinsic source spectra shown in Fig. 4 we unfolded the observed spectra by the energy dispersion and the EBL absorption.

The following procedure is common for the MAGIC and XRT analyses. The energy-dependent photon index is defined as Γ(E) = α + 2β log₁₀(E/E₀) (see Eq. (4) in Massaro et al. 2004). The fitting range of the VHE spectra is restricted to 150 GeV–1 TeV in order to avoid a possible high-energy cutoff. For each fit, the value of χ² is calculated. If the LogP function deviates from the spectrum so that the fit probability is smaller than 5%, the night (for MAGIC) or observation (for XRT) is removed from the sample. The 5% cut corresponds to about 2σ and leaves most of the data points compatible with the LogP spectral shape. For the nights and observations that satisfy the above criteria, we adopted the value of α as a measure of the local spectral index at the normalization energy E₀.

Appendix B: Spectral fitting functions

The functions used for the spectral fitting are defined as follows: a simple power law (PL)

$$\begin{array}{c}\hfill \frac{\mathrm{d}F}{\mathrm{d}E}=F_0{\left(\frac{E}{{\mathrm{E}}_0}\right)}^{-\mathrm{\Gamma }},\end{array}$$(B.1)

a PL with an exponential cutoff

$$\begin{array}{c}\hfill \frac{\mathrm{d}F}{\mathrm{d}E}=F_0{\left(\frac{E}{{\mathrm{E}}_0}\right)}^{-\mathrm{\Gamma }}exp(-\frac{E}{E_{\mathrm{cut}}}),\end{array}$$(B.2)

a log-parabola (LogP)

$$\begin{array}{c}\hfill \frac{\mathrm{d}F}{\mathrm{d}E}=F_0{\left(\frac{E}{{\mathrm{E}}_0}\right)}^{-\alpha -\beta [{log}_{10}(E/{\mathrm{E}}_0)]},\end{array}$$(B.3)

and a LogP with an exponential cutoff

$$\begin{array}{c}\hfill \frac{\mathrm{d}F}{\mathrm{d}E}=F_0{\left(\frac{E}{{\mathrm{E}}_0}\right)}^{-\alpha -\beta [{log}_{10}(E/{\mathrm{E}}_0)]}exp(-\frac{E}{E_{\mathrm{cut}}}),\end{array}$$(B.4)

where dF/dE is the differential γ-ray flux as a function of the energy E. The value of the normalization energy E₀ is fixed at 300 GeV. The expressions of the LogP and the LogP with a cutoff with E_peak are

$$\begin{array}{c}\hfill \frac{\mathrm{d}F}{\mathrm{d}E}=F_0{\left(\frac{E}{E_{\mathrm{peak}}}\right)}^{-2}{10}^{-\beta {[log(E/E_{\mathrm{peak}})]}^2}\end{array}$$(B.5)

and

$$\begin{array}{c}\hfill \frac{\mathrm{d}F}{\mathrm{d}E}=F_0{\left(\frac{E}{E_{\mathrm{peak}}}\right)}^{-2}{10}^{-\beta {[{log}_{10}(E/E_{\mathrm{peak}})]}^2}exp(-\frac{E}{E_{\mathrm{cut}}}),\end{array}$$(B.6)

respectively. The parameter E_peak corresponds to the peak energy of a SED.

Appendix C: Weighted Pearson correlation coefficient calculation

We calculate the weighted mean and weighted covariance for two quantities x and y (in our case spectral index and flux) with errors σ_x, σ_y respectively using the following formulae

$$\begin{array}{c}\hfill \mathrm{mean}(x;w)=\frac{{\mathrm{\Sigma }}_i{w}_i{x}_i}{{\mathrm{\Sigma }}_i{w}_i}(similarexpressionfory),\end{array}$$(C.1)

where $w_i=1/{\sigma }_{x\mathbf{i}}^2$ or $w_i=1/{\sigma }_{y\mathbf{i}}^2$;

$$\begin{array}{c}\hfill \mathrm{cov}(x,y;w)=\frac{{\mathrm{\Sigma }}_i{w}_i(x_i-\mathrm{mean}(x;w)).(y_i-\mathrm{mean}(y;w))}{{\mathrm{\Sigma }}_i{w}_i},\end{array}$$(C.2)

where w_i = 1/(σ_xiσ_yi).

Using the above definitions, the weighted Pearson correlation coefficient can be calculated as

$$\begin{array}{c}\hfill \mathrm{corr}(x,y;w)=\frac{\mathrm{cov}(x,y;w)}{\sqrt{\mathrm{cov}(x,x;w).\mathrm{cov}(y,y;w)}}\end{array}$$(C.3)

The errors of the correlation coefficient have been calculated using the z-transformed Discrete Correlation Function (Alexander 2013; see also Edelson & Krolik 1988 and Anderhub et al. 2009).

All Tables

All Figures

	Fig. 1. Long-term light curves of 1ES 1959+650 in 2016 with four instruments. From top to bottom: VHE gamma-ray flux (>300 GeV) from MAGIC, HE gamma-ray flux (0.3–300 GeV) measured with Fermi-LAT, X-ray energy flux (0.5–5.0 keV) from Swift-XRT, and UV energy flux (W1 filter, ∼260 nm) from Swift-UVOT. The flux in Crab Units is indicated with an additional y-axis in the top panel.
In the text

	Fig. 2. Discrete correlation function as a function of time lag for the VHE γ-ray and X-ray light curves of the 1ES 1959+650 2016 multi-wavelength monitoring campaign. In the long term, the VHE flux shows a correlation (DCF ∼ 0.76 ± 0.1) with the X-ray flux with zero time lag.
In the text

	Fig. 3. Correlation between α of LogP fitting to the MAGIC spectrum for each night (left panel) and the XRT spectrum of each observation (right). The normalisation energy is fixed at 300 GeV for MAGIC and at 1 keV for XRT. The MAGIC spectra are deabsorbed with the model of Franceschini et al. (2008). More details can be found in the text.
In the text

	Fig. 4. VHE SEDs during the nights of highest flux: 13 June, 14 June, and 1 July 2016 (from top to bottom). The black circle and red square markers represent the observed and the EBL-deabsorbed spectra, respectively. These spectra have been unfolded with the instrument response function of MAGIC. The absorption by the EBL has been corrected via the model of Franceschini et al. (2008).
In the text

Fig. 5. Fast intra-night variability observed in the VHE band on 13 June (left panel) and 1 July (right panel) 2016. The light curves above 300 GeV are constructed with a fixed time-binning of 10 min. Green dashed curve: fit with the function given in Eq. (3); solid red curve: fit with the function given in Eq. (4). See the figure legend for the rise and decay times inferred from the fit. The blue dashed line represents the steady flux of the Crab Nebula above 300 GeV, shown for comparison purposes.

In the text

Fig. 6. One-zone SSC models applied to 13 June (left panel) and 14 June (right panel) 2016. The symbols corresponding to the data from different instruments are given in the legend. The historical data are taken from Tagliaferri et al. (2008). The black (solid), brown (dot-dashed), and blue (dashed) curves represent the summed emission component in increasing order of Doppler factor δ. We found a satisfactory explanation of the MWL data with high values of δ >  45 for 13 June 2016. The data from 14 June 2016 do not strictly require such high values and can be modelled with moderate values of δ >  30. For more details see the discussion in Sect. 4.3.1 and the parameters in Table 4.

In the text

Fig. 7. One-zone hadronic models applied to 13 June (left panel) and 14 June (right panel) 2016. The symbols corresponding to the data from different instruments are given in the legend. Solid black line: summed; dashed blue line: electron-synchrotron; dotted green line: SSC; dot-dashed sea green line: proton-synchrotron; dot-dot-dashed orange line: p − γ cascade. The higher energy peak in the SED is dominated by synchrotron radiation by relativistic protons, which can be achieved with B ∼ 100 G and Ep, max >  1018 eV and jet power Lj ∼ 1046 erg s−1 (∼LEdd). For more details see the discussion in Sect. 4.3.2 and the parameters in Table 4.

In the text

Fig. 8. One-zone lepto-hadronic models (left panel) and the predicted neutrino flux (right panel) for 13 June 2016. The definition of symbols and lines in the SED model in the left panel is the same as Fig. 7. The higher energy peak in the SED in this case, is a combination of the SSC and photo-meson cascade component which can be achieved with B ∼ 1 G, Ep, max >  1016 eV at the cost of high jet power Lj >  1048 erg s−1 (≫LEdd). The meaning of the different curves in the neutrino spectra (right panel) is mentioned in the legend. The IceCube sensitivity curve is taken from IceCube Collaboration (2019) corresponding to declination 60°. Neutrino spectra predicted in the proton-synchrotron solutions of Fig. 7 peak at very high energies and provide low neutrino flux in the range 0.1–100 PeV. The neutrino peak is shifted to lower energies in the lepto-hadronic solutions providing slightly higher flux at the cost of very high values of the jet luminosity. For more details see the discussion in Sect. 4.3.3 and the parameters in Table 4.

In the text

	Fig. 9. Comparison between acceleration timescale (solid black line with acceleration efficiency ηacc  =  1) and different cooling timescales (tesc, tpsync, tpγ) for the hadronic (left panel) and lepto-hadronic (right panel) scenarios of 13 June 2016. Also shown is the energy at which the proton gyro-radius becomes equal to the radius of the emission region (dot-dashed red line).
In the text